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