Actual source code: ex24.c
2: static char help[] = "Tests CG, MINRES and SYMMLQ on symmetric matrices with SBAIJ format. The preconditioner ICC only works on sequential SBAIJ format. \n\n";
4: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Mat C;
9: PetscScalar v, none = -1.0;
10: PetscInt i, j, Ii, J, Istart, Iend, N, m = 4, n = 4, its, k;
11: PetscMPIInt size, rank;
12: PetscReal err_norm, res_norm;
13: Vec x, b, u, u_tmp;
14: PC pc;
15: KSP ksp;
18: PetscInitialize(&argc, &args, (char *)0, help);
19: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
20: MPI_Comm_size(PETSC_COMM_WORLD, &size);
21: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
22: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
23: N = m * n;
25: /* Generate matrix */
26: MatCreate(PETSC_COMM_WORLD, &C);
27: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
28: MatSetFromOptions(C);
29: MatSetUp(C);
30: MatGetOwnershipRange(C, &Istart, &Iend);
31: for (Ii = Istart; Ii < Iend; Ii++) {
32: v = -1.0;
33: i = Ii / n;
34: j = Ii - i * n;
35: if (i > 0) {
36: J = Ii - n;
37: MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES);
38: }
39: if (i < m - 1) {
40: J = Ii + n;
41: MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES);
42: }
43: if (j > 0) {
44: J = Ii - 1;
45: MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES);
46: }
47: if (j < n - 1) {
48: J = Ii + 1;
49: MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES);
50: }
51: v = 4.0;
52: MatSetValues(C, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
53: }
54: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
57: /* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */
58: /* MatShift(C,alpha); */
59: /* MatView(C,PETSC_VIEWER_STDOUT_WORLD); */
61: /* Setup and solve for system */
62: /* Create vectors. */
63: VecCreate(PETSC_COMM_WORLD, &x);
64: VecSetSizes(x, PETSC_DECIDE, N);
65: VecSetFromOptions(x);
66: VecDuplicate(x, &b);
67: VecDuplicate(x, &u);
68: VecDuplicate(x, &u_tmp);
69: /* Set exact solution u; then compute right-hand-side vector b. */
70: VecSet(u, 1.0);
71: MatMult(C, u, b);
73: for (k = 0; k < 3; k++) {
74: if (k == 0) { /* CG */
75: KSPCreate(PETSC_COMM_WORLD, &ksp);
76: KSPSetOperators(ksp, C, C);
77: PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n");
78: KSPSetType(ksp, KSPCG);
79: } else if (k == 1) { /* MINRES */
80: KSPCreate(PETSC_COMM_WORLD, &ksp);
81: KSPSetOperators(ksp, C, C);
82: PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n");
83: KSPSetType(ksp, KSPMINRES);
84: } else { /* SYMMLQ */
85: KSPCreate(PETSC_COMM_WORLD, &ksp);
86: KSPSetOperators(ksp, C, C);
87: PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n");
88: KSPSetType(ksp, KSPSYMMLQ);
89: }
90: KSPGetPC(ksp, &pc);
91: /* PCSetType(pc,PCICC); */
92: PCSetType(pc, PCJACOBI);
93: KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
95: /*
96: Set runtime options, e.g.,
97: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
98: These options will override those specified above as long as
99: KSPSetFromOptions() is called _after_ any other customization
100: routines.
101: */
102: KSPSetFromOptions(ksp);
104: /* Solve linear system; */
105: KSPSetUp(ksp);
106: KSPSolve(ksp, b, x);
108: KSPGetIterationNumber(ksp, &its);
109: /* Check error */
110: VecCopy(u, u_tmp);
111: VecAXPY(u_tmp, none, x);
112: VecNorm(u_tmp, NORM_2, &err_norm);
113: MatMult(C, x, u_tmp);
114: VecAXPY(u_tmp, none, b);
115: VecNorm(u_tmp, NORM_2, &res_norm);
117: PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its);
118: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g;", (double)res_norm);
119: PetscPrintf(PETSC_COMM_WORLD, " Error norm %g.\n", (double)err_norm);
120: KSPDestroy(&ksp);
121: }
123: /*
124: Free work space. All PETSc objects should be destroyed when they
125: are no longer needed.
126: */
127: VecDestroy(&b);
128: VecDestroy(&u);
129: VecDestroy(&x);
130: VecDestroy(&u_tmp);
131: MatDestroy(&C);
133: PetscFinalize();
134: return 0;
135: }
137: /*TEST
139: test:
140: args: -pc_type icc -mat_type seqsbaij -mat_ignore_lower_triangular
142: test:
143: suffix: 2
144: args: -pc_type icc -pc_factor_levels 2 -mat_type seqsbaij -mat_ignore_lower_triangular
146: test:
147: suffix: 3
148: nsize: 2
149: args: -pc_type bjacobi -sub_pc_type icc -mat_type mpisbaij -mat_ignore_lower_triangular -ksp_max_it 8
151: test:
152: suffix: 4
153: nsize: 2
154: args: -pc_type bjacobi -sub_pc_type icc -sub_pc_factor_levels 1 -mat_type mpisbaij -mat_ignore_lower_triangular
156: TEST*/