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";
2
3 #include <petscksp.h>
4
main(int argc,char ** args)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;
15
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;
23
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));
55
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)); */
59
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));
71
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));
93
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));
102
103 /* Solve linear system; */
104 PetscCall(KSPSetUp(ksp));
105 PetscCall(KSPSolve(ksp, b, x));
106
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));
115
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 }
121
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));
131
132 PetscCall(PetscFinalize());
133 return 0;
134 }
135
136 /*TEST
137
138 test:
139 args: -pc_type icc -mat_type seqsbaij -mat_ignore_lower_triangular
140
141 test:
142 suffix: 2
143 args: -pc_type icc -pc_factor_levels 2 -mat_type seqsbaij -mat_ignore_lower_triangular
144
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
149
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
154
155 TEST*/
156