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;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
19: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
22: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
23: N = m * n;
25: /* Generate matrix */
26: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
27: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
28: PetscCall(MatSetFromOptions(C));
29: PetscCall(MatSetUp(C));
30: PetscCall(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: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
38: }
39: if (i < m - 1) {
40: J = Ii + n;
41: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
42: }
43: if (j > 0) {
44: J = Ii - 1;
45: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
46: }
47: if (j < n - 1) {
48: J = Ii + 1;
49: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, ADD_VALUES));
50: }
51: v = 4.0;
52: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
53: }
54: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
55: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
57: /* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */
58: /* PetscCall(MatShift(C,alpha)); */
59: /* PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); */
61: /* Setup and solve for system */
62: /* Create vectors. */
63: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
64: PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
65: PetscCall(VecSetFromOptions(x));
66: PetscCall(VecDuplicate(x, &b));
67: PetscCall(VecDuplicate(x, &u));
68: PetscCall(VecDuplicate(x, &u_tmp));
69: /* Set exact solution u; then compute right-hand-side vector b. */
70: PetscCall(VecSet(u, 1.0));
71: PetscCall(MatMult(C, u, b));
73: for (k = 0; k < 3; k++) {
74: if (k == 0) { /* CG */
75: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
76: PetscCall(KSPSetOperators(ksp, C, C));
77: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n CG: \n"));
78: PetscCall(KSPSetType(ksp, KSPCG));
79: } else if (k == 1) { /* MINRES */
80: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
81: PetscCall(KSPSetOperators(ksp, C, C));
82: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n MINRES: \n"));
83: PetscCall(KSPSetType(ksp, KSPMINRES));
84: } else { /* SYMMLQ */
85: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
86: PetscCall(KSPSetOperators(ksp, C, C));
87: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n SYMMLQ: \n"));
88: PetscCall(KSPSetType(ksp, KSPSYMMLQ));
89: }
90: PetscCall(KSPGetPC(ksp, &pc));
91: /* PetscCall(PCSetType(pc,PCICC)); */
92: PetscCall(PCSetType(pc, PCJACOBI));
93: PetscCall(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: PetscCall(KSPSetFromOptions(ksp));
104: /* Solve linear system; */
105: PetscCall(KSPSetUp(ksp));
106: PetscCall(KSPSolve(ksp, b, x));
108: PetscCall(KSPGetIterationNumber(ksp, &its));
109: /* Check error */
110: PetscCall(VecCopy(u, u_tmp));
111: PetscCall(VecAXPY(u_tmp, none, x));
112: PetscCall(VecNorm(u_tmp, NORM_2, &err_norm));
113: PetscCall(MatMult(C, x, u_tmp));
114: PetscCall(VecAXPY(u_tmp, none, b));
115: PetscCall(VecNorm(u_tmp, NORM_2, &res_norm));
117: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
118: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g;", (double)res_norm));
119: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm %g.\n", (double)err_norm));
120: PetscCall(KSPDestroy(&ksp));
121: }
123: /*
124: Free work space. All PETSc objects should be destroyed when they
125: are no longer needed.
126: */
127: PetscCall(VecDestroy(&b));
128: PetscCall(VecDestroy(&u));
129: PetscCall(VecDestroy(&x));
130: PetscCall(VecDestroy(&u_tmp));
131: PetscCall(MatDestroy(&C));
133: PetscCall(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*/