1 static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Vec x, y, b;
8 Mat A; /* linear system matrix */
9 Mat sA, sC; /* symmetric part of the matrices */
10 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, col[3], block, row, Ii, J, n1, lvl;
11 PetscMPIInt size;
12 PetscReal norm2;
13 PetscScalar neg_one = -1.0, four = 4.0, value[3];
14 IS perm, cperm;
15 PetscRandom rdm;
16 PetscBool reorder = PETSC_FALSE, displ = PETSC_FALSE;
17 MatFactorInfo factinfo;
18 PetscBool equal;
19 PetscBool TestAIJ = PETSC_FALSE, TestBAIJ = PETSC_TRUE;
20 PetscInt TestShift = 0;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscInitialize(&argc, &args, NULL, help));
24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
25 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
28 PetscCall(PetscOptionsGetBool(NULL, NULL, "-reorder", &reorder, NULL));
29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-testaij", &TestAIJ, NULL));
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-testShift", &TestShift, NULL));
31 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));
32
33 n = mbs * bs;
34 if (TestAIJ) { /* A is in aij format */
35 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, n, n, nz, NULL, &A));
36 TestBAIJ = PETSC_FALSE;
37 } else { /* A is in baij format */
38 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A));
39 TestAIJ = PETSC_FALSE;
40 }
41
42 /* Assemble matrix */
43 if (bs == 1) {
44 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL));
45 if (prob == 1) { /* tridiagonal matrix */
46 value[0] = -1.0;
47 value[1] = 2.0;
48 value[2] = -1.0;
49 for (i = 1; i < n - 1; i++) {
50 col[0] = i - 1;
51 col[1] = i;
52 col[2] = i + 1;
53 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
54 }
55 i = n - 1;
56 col[0] = 0;
57 col[1] = n - 2;
58 col[2] = n - 1;
59
60 value[0] = 0.1;
61 value[1] = -1;
62 value[2] = 2;
63 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
64
65 i = 0;
66 col[0] = 0;
67 col[1] = 1;
68 col[2] = n - 1;
69
70 value[0] = 2.0;
71 value[1] = -1.0;
72 value[2] = 0.1;
73 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
74 } else if (prob == 2) { /* matrix for the five point stencil */
75 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
76 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
77 for (i = 0; i < n1; i++) {
78 for (j = 0; j < n1; j++) {
79 Ii = j + n1 * i;
80 if (i > 0) {
81 J = Ii - n1;
82 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
83 }
84 if (i < n1 - 1) {
85 J = Ii + n1;
86 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
87 }
88 if (j > 0) {
89 J = Ii - 1;
90 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
91 }
92 if (j < n1 - 1) {
93 J = Ii + 1;
94 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
95 }
96 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
97 }
98 }
99 }
100 } else { /* bs > 1 */
101 for (block = 0; block < n / bs; block++) {
102 /* diagonal blocks */
103 value[0] = -1.0;
104 value[1] = 4.0;
105 value[2] = -1.0;
106 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
107 col[0] = i - 1;
108 col[1] = i;
109 col[2] = i + 1;
110 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
111 }
112 i = bs - 1 + block * bs;
113 col[0] = bs - 2 + block * bs;
114 col[1] = bs - 1 + block * bs;
115
116 value[0] = -1.0;
117 value[1] = 4.0;
118 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
119
120 i = 0 + block * bs;
121 col[0] = 0 + block * bs;
122 col[1] = 1 + block * bs;
123
124 value[0] = 4.0;
125 value[1] = -1.0;
126 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
127 }
128 /* off-diagonal blocks */
129 value[0] = -1.0;
130 for (i = 0; i < (n / bs - 1) * bs; i++) {
131 col[0] = i + bs;
132 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
133 col[0] = i;
134 row = i + bs;
135 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
136 }
137 }
138
139 if (TestShift) {
140 /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
141 for (i = 0; i < bs; i++) {
142 row = i;
143 col[0] = i;
144 value[0] = 0.0;
145 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
146 }
147 }
148
149 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
150 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
151
152 /* Test MatConvert */
153 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
154 PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));
155 PetscCall(MatMultEqual(A, sA, 20, &equal));
156 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_USER, "A != sA");
157
158 /* Test MatGetOwnershipRange() */
159 PetscCall(MatGetOwnershipRange(A, &Ii, &J));
160 PetscCall(MatGetOwnershipRange(sA, &i, &j));
161 PetscCheck(i == Ii && j == J, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetOwnershipRange() in MatSBAIJ format");
162
163 /* Vectors */
164 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
165 PetscCall(PetscRandomSetFromOptions(rdm));
166 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
167 PetscCall(VecDuplicate(x, &b));
168 PetscCall(VecDuplicate(x, &y));
169 PetscCall(VecSetRandom(x, rdm));
170
171 /* Test MatReordering() - not work on sbaij matrix */
172 if (reorder) {
173 PetscCall(MatGetOrdering(A, MATORDERINGRCM, &perm, &cperm));
174 } else {
175 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &cperm));
176 }
177 PetscCall(ISDestroy(&cperm));
178
179 /* initialize factinfo */
180 PetscCall(MatFactorInfoInitialize(&factinfo));
181 if (TestShift == 1) {
182 factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO;
183 factinfo.shiftamount = 0.1;
184 } else if (TestShift == 2) {
185 factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
186 }
187
188 /* Test MatCholeskyFactor(), MatICCFactor() */
189 /*------------------------------------------*/
190 /* Test aij matrix A */
191 if (TestAIJ) {
192 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "AIJ: \n"));
193 i = 0;
194 for (lvl = -1; lvl < 10; lvl++) {
195 if (lvl == -1) { /* Cholesky factor */
196 factinfo.fill = 5.0;
197
198 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
199 PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
200 } else { /* incomplete Cholesky factor */
201 factinfo.fill = 5.0;
202 factinfo.levels = lvl;
203
204 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
205 PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
206 }
207 PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));
208
209 PetscCall(MatMult(A, x, b));
210 PetscCall(MatSolve(sC, b, y));
211 PetscCall(MatDestroy(&sC));
212
213 /* Check the residual */
214 PetscCall(VecAXPY(y, neg_one, x));
215 PetscCall(VecNorm(y, NORM_2, &norm2));
216
217 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
218 }
219 }
220
221 /* Test baij matrix A */
222 if (TestBAIJ) {
223 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "BAIJ: \n"));
224 i = 0;
225 for (lvl = -1; lvl < 10; lvl++) {
226 if (lvl == -1) { /* Cholesky factor */
227 factinfo.fill = 5.0;
228
229 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
230 PetscCall(MatCholeskyFactorSymbolic(sC, A, perm, &factinfo));
231 } else { /* incomplete Cholesky factor */
232 factinfo.fill = 5.0;
233 factinfo.levels = lvl;
234
235 PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
236 PetscCall(MatICCFactorSymbolic(sC, A, perm, &factinfo));
237 }
238 PetscCall(MatCholeskyFactorNumeric(sC, A, &factinfo));
239
240 PetscCall(MatMult(A, x, b));
241 PetscCall(MatSolve(sC, b, y));
242 PetscCall(MatDestroy(&sC));
243
244 /* Check the residual */
245 PetscCall(VecAXPY(y, neg_one, x));
246 PetscCall(VecNorm(y, NORM_2, &norm2));
247 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
248 }
249 }
250
251 /* Test sbaij matrix sA */
252 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SBAIJ: \n"));
253 i = 0;
254 for (lvl = -1; lvl < 10; lvl++) {
255 if (lvl == -1) { /* Cholesky factor */
256 factinfo.fill = 5.0;
257
258 PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sC));
259 PetscCall(MatCholeskyFactorSymbolic(sC, sA, perm, &factinfo));
260 } else { /* incomplete Cholesky factor */
261 factinfo.fill = 5.0;
262 factinfo.levels = lvl;
263
264 PetscCall(MatGetFactor(sA, MATSOLVERPETSC, MAT_FACTOR_ICC, &sC));
265 PetscCall(MatICCFactorSymbolic(sC, sA, perm, &factinfo));
266 }
267 PetscCall(MatCholeskyFactorNumeric(sC, sA, &factinfo));
268
269 if (lvl == 0 && bs == 1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
270 /*
271 Mat B;
272 PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&B));
273 PetscCall(MatICCFactor(B,perm,&factinfo));
274 PetscCall(MatEqual(sC,B,&equal));
275 PetscCheck(equal, PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
276 PetscCall(MatDestroy(&B));
277 */
278 }
279
280 PetscCall(MatMult(sA, x, b));
281 PetscCall(MatSolve(sC, b, y));
282
283 /* Test MatSolves() */
284 if (bs == 1) {
285 Vecs xx, bb;
286 PetscCall(VecsCreateSeq(PETSC_COMM_SELF, n, 4, &xx));
287 PetscCall(VecsDuplicate(xx, &bb));
288 PetscCall(MatSolves(sC, bb, xx));
289 PetscCall(VecsDestroy(xx));
290 PetscCall(VecsDestroy(bb));
291 }
292 PetscCall(MatDestroy(&sC));
293
294 /* Check the residual */
295 PetscCall(VecAXPY(y, neg_one, x));
296 PetscCall(VecNorm(y, NORM_2, &norm2));
297 if (displ) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " lvl: %" PetscInt_FMT ", residual: %g\n", lvl, (double)norm2));
298 }
299
300 PetscCall(ISDestroy(&perm));
301 PetscCall(MatDestroy(&A));
302 PetscCall(MatDestroy(&sA));
303 PetscCall(VecDestroy(&x));
304 PetscCall(VecDestroy(&y));
305 PetscCall(VecDestroy(&b));
306 PetscCall(PetscRandomDestroy(&rdm));
307
308 PetscCall(PetscFinalize());
309 return 0;
310 }
311
312 /*TEST
313
314 test:
315 args: -bs {{1 2 3 4 5 6 7 8}}
316 output_file: output/empty.out
317
318 test:
319 suffix: 3
320 args: -testaij
321 output_file: output/empty.out
322
323 TEST*/
324