1 static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 PetscMPIInt size;
8 Vec x, y, b, s1, s2;
9 Mat A; /* linear system matrix */
10 Mat sA, sB, sFactor, B, C; /* symmetric matrices */
11 PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 1, i, j, k1, k2, col[3], lf, block, row, Ii, J, n1, inc;
12 PetscReal norm1, norm2, rnorm, tol = 10 * PETSC_SMALL;
13 PetscScalar neg_one = -1.0, four = 4.0, value[3];
14 IS perm, iscol;
15 PetscRandom rdm;
16 PetscBool doIcc = PETSC_TRUE, equal;
17 MatInfo minfo1, minfo2;
18 MatFactorInfo factinfo;
19 MatType type;
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscInitialize(&argc, &args, NULL, help));
23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
24 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
27
28 n = mbs * bs;
29 PetscCall(MatCreate(PETSC_COMM_SELF, &A));
30 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
31 PetscCall(MatSetType(A, MATSEQBAIJ));
32 PetscCall(MatSetFromOptions(A));
33 PetscCall(MatSeqBAIJSetPreallocation(A, bs, nz, NULL));
34
35 PetscCall(MatCreate(PETSC_COMM_SELF, &sA));
36 PetscCall(MatSetSizes(sA, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
37 PetscCall(MatSetType(sA, MATSEQSBAIJ));
38 PetscCall(MatSetFromOptions(sA));
39 PetscCall(MatGetType(sA, &type));
40 PetscCall(PetscObjectTypeCompare((PetscObject)sA, MATSEQSBAIJ, &doIcc));
41 PetscCall(MatSeqSBAIJSetPreallocation(sA, bs, nz, NULL));
42 PetscCall(MatSetOption(sA, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
43
44 /* Test MatGetOwnershipRange() */
45 PetscCall(MatGetOwnershipRange(A, &Ii, &J));
46 PetscCall(MatGetOwnershipRange(sA, &i, &j));
47 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n"));
48
49 /* Assemble matrix */
50 if (bs == 1) {
51 PetscCall(PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL));
52 if (prob == 1) { /* tridiagonal matrix */
53 value[0] = -1.0;
54 value[1] = 2.0;
55 value[2] = -1.0;
56 for (i = 1; i < n - 1; i++) {
57 col[0] = i - 1;
58 col[1] = i;
59 col[2] = i + 1;
60 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
61 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
62 }
63 i = n - 1;
64 col[0] = 0;
65 col[1] = n - 2;
66 col[2] = n - 1;
67
68 value[0] = 0.1;
69 value[1] = -1;
70 value[2] = 2;
71
72 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
73 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
74
75 i = 0;
76 col[0] = n - 1;
77 col[1] = 1;
78 col[2] = 0;
79 value[0] = 0.1;
80 value[1] = -1.0;
81 value[2] = 2;
82
83 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
84 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
85
86 } else if (prob == 2) { /* matrix for the five point stencil */
87 n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
88 PetscCheck(n1 * n1 == n, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "sqrt(n) must be a positive integer!");
89 for (i = 0; i < n1; i++) {
90 for (j = 0; j < n1; j++) {
91 Ii = j + n1 * i;
92 if (i > 0) {
93 J = Ii - n1;
94 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
95 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
96 }
97 if (i < n1 - 1) {
98 J = Ii + n1;
99 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
100 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
101 }
102 if (j > 0) {
103 J = Ii - 1;
104 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
105 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
106 }
107 if (j < n1 - 1) {
108 J = Ii + 1;
109 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
110 PetscCall(MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
111 }
112 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
113 PetscCall(MatSetValues(sA, 1, &Ii, 1, &Ii, &four, INSERT_VALUES));
114 }
115 }
116 }
117
118 } else { /* bs > 1 */
119 for (block = 0; block < n / bs; block++) {
120 /* diagonal blocks */
121 value[0] = -1.0;
122 value[1] = 4.0;
123 value[2] = -1.0;
124 for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
125 col[0] = i - 1;
126 col[1] = i;
127 col[2] = i + 1;
128 PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
129 PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
130 }
131 i = bs - 1 + block * bs;
132 col[0] = bs - 2 + block * bs;
133 col[1] = bs - 1 + block * bs;
134
135 value[0] = -1.0;
136 value[1] = 4.0;
137
138 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
139 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
140
141 i = 0 + block * bs;
142 col[0] = 0 + block * bs;
143 col[1] = 1 + block * bs;
144
145 value[0] = 4.0;
146 value[1] = -1.0;
147
148 PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
149 PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
150 }
151 /* off-diagonal blocks */
152 value[0] = -1.0;
153 for (i = 0; i < (n / bs - 1) * bs; i++) {
154 col[0] = i + bs;
155
156 PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
157 PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES));
158
159 col[0] = i;
160 row = i + bs;
161
162 PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
163 PetscCall(MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES));
164 }
165 }
166 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
167 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
168
169 PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY));
170 PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY));
171
172 /* Test MatGetInfo() of A and sA */
173 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1));
174 PetscCall(MatGetInfo(sA, MAT_LOCAL, &minfo2));
175 i = (int)(minfo1.nz_used - minfo2.nz_used);
176 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
177 k1 = (int)(minfo1.nz_allocated - minfo1.nz_used);
178 k2 = (int)(minfo2.nz_allocated - minfo2.nz_used);
179 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error (compare A and sA): MatGetInfo()\n"));
180
181 /* Test MatDuplicate() */
182 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1));
183 PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB));
184 PetscCall(MatEqual(sA, sB, &equal));
185 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDuplicate()");
186
187 /* Test MatNorm() */
188 PetscCall(MatNorm(A, NORM_FROBENIUS, &norm1));
189 PetscCall(MatNorm(sB, NORM_FROBENIUS, &norm2));
190 rnorm = PetscAbsReal(norm1 - norm2) / norm2;
191 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2));
192 PetscCall(MatNorm(A, NORM_INFINITY, &norm1));
193 PetscCall(MatNorm(sB, NORM_INFINITY, &norm2));
194 rnorm = PetscAbsReal(norm1 - norm2) / norm2;
195 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2));
196 PetscCall(MatNorm(A, NORM_1, &norm1));
197 PetscCall(MatNorm(sB, NORM_1, &norm2));
198 rnorm = PetscAbsReal(norm1 - norm2) / norm2;
199 if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n", (double)norm1, (double)norm2));
200
201 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
202 PetscCall(MatGetInfo(A, MAT_LOCAL, &minfo1));
203 PetscCall(MatGetInfo(sB, MAT_LOCAL, &minfo2));
204 i = (int)(minfo1.nz_used - minfo2.nz_used);
205 j = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
206 k1 = (int)(minfo1.nz_allocated - minfo1.nz_used);
207 k2 = (int)(minfo2.nz_allocated - minfo2.nz_used);
208 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error(compare A and sB): MatGetInfo()\n"));
209
210 PetscCall(MatGetSize(A, &Ii, &J));
211 PetscCall(MatGetSize(sB, &i, &j));
212 if (i - Ii || j - J) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n"));
213
214 PetscCall(MatGetBlockSize(A, &Ii));
215 PetscCall(MatGetBlockSize(sB, &i));
216 if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n"));
217
218 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
219 PetscCall(PetscRandomSetFromOptions(rdm));
220 PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
221 PetscCall(VecDuplicate(x, &s1));
222 PetscCall(VecDuplicate(x, &s2));
223 PetscCall(VecDuplicate(x, &y));
224 PetscCall(VecDuplicate(x, &b));
225 PetscCall(VecSetRandom(x, rdm));
226
227 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
228 #if !defined(PETSC_USE_COMPLEX)
229 /* Scaling matrix with complex numbers results non-spd matrix,
230 causing crash of MatForwardSolve() and MatBackwardSolve() */
231 PetscCall(MatDiagonalScale(A, x, x));
232 PetscCall(MatDiagonalScale(sB, x, x));
233 PetscCall(MatMultEqual(A, sB, 10, &equal));
234 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale");
235
236 PetscCall(MatGetDiagonal(A, s1));
237 PetscCall(MatGetDiagonal(sB, s2));
238 PetscCall(VecAXPY(s2, neg_one, s1));
239 PetscCall(VecNorm(s2, NORM_1, &norm1));
240 if (norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal(), ||s1-s2||=%g\n", (double)norm1));
241
242 {
243 PetscScalar alpha = 0.1;
244 PetscCall(MatScale(A, alpha));
245 PetscCall(MatScale(sB, alpha));
246 }
247 #endif
248
249 /* Test MatGetRowMaxAbs() */
250 PetscCall(MatGetRowMaxAbs(A, s1, NULL));
251 PetscCall(MatGetRowMaxAbs(sB, s2, NULL));
252 PetscCall(VecNorm(s1, NORM_1, &norm1));
253 PetscCall(VecNorm(s2, NORM_1, &norm2));
254 norm1 -= norm2;
255 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatGetRowMaxAbs() \n"));
256
257 /* Test MatMult() */
258 for (i = 0; i < 40; i++) {
259 PetscCall(VecSetRandom(x, rdm));
260 PetscCall(MatMult(A, x, s1));
261 PetscCall(MatMult(sB, x, s2));
262 PetscCall(VecNorm(s1, NORM_1, &norm1));
263 PetscCall(VecNorm(s2, NORM_1, &norm2));
264 norm1 -= norm2;
265 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), norm1-norm2: %g\n", (double)norm1));
266 }
267
268 /* MatMultAdd() */
269 for (i = 0; i < 40; i++) {
270 PetscCall(VecSetRandom(x, rdm));
271 PetscCall(VecSetRandom(y, rdm));
272 PetscCall(MatMultAdd(A, x, y, s1));
273 PetscCall(MatMultAdd(sB, x, y, s2));
274 PetscCall(VecNorm(s1, NORM_1, &norm1));
275 PetscCall(VecNorm(s2, NORM_1, &norm2));
276 norm1 -= norm2;
277 if (norm1 < -tol || norm1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), norm1-norm2: %g\n", (double)norm1));
278 }
279
280 /* Test MatMatMult() for sbaij and dense matrices */
281 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, 5 * n, NULL, &B));
282 PetscCall(MatSetRandom(B, rdm));
283 PetscCall(MatMatMult(sA, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
284 PetscCall(MatMatMultEqual(sA, B, C, 5 * n, &equal));
285 PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error: MatMatMult()");
286 PetscCall(MatDestroy(&C));
287 PetscCall(MatDestroy(&B));
288
289 /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
290 PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iscol));
291 PetscCall(ISDestroy(&iscol));
292 norm1 = tol;
293 inc = bs;
294
295 /* initialize factinfo */
296 PetscCall(PetscMemzero(&factinfo, sizeof(MatFactorInfo)));
297
298 for (lf = -1; lf < 10; lf += inc) {
299 if (lf == -1) { /* Cholesky factor of sB (duplicate sA) */
300 factinfo.fill = 5.0;
301
302 PetscCall(MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sFactor));
303 PetscCall(MatCholeskyFactorSymbolic(sFactor, sB, perm, &factinfo));
304 } else if (!doIcc) break;
305 else { /* incomplete Cholesky factor */ factinfo.fill = 5.0;
306 factinfo.levels = lf;
307
308 PetscCall(MatGetFactor(sB, MATSOLVERPETSC, MAT_FACTOR_ICC, &sFactor));
309 PetscCall(MatICCFactorSymbolic(sFactor, sB, perm, &factinfo));
310 }
311 PetscCall(MatCholeskyFactorNumeric(sFactor, sB, &factinfo));
312 /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */
313
314 /* test MatGetDiagonal on numeric factor */
315 /*
316 if (lf == -1) {
317 PetscCall(MatGetDiagonal(sFactor,s1));
318 printf(" in ex74.c, diag: \n");
319 PetscCall(VecView(s1,PETSC_VIEWER_STDOUT_SELF));
320 }
321 */
322
323 PetscCall(MatMult(sB, x, b));
324
325 /* test MatForwardSolve() and MatBackwardSolve() */
326 if (lf == -1) {
327 PetscCall(MatForwardSolve(sFactor, b, s1));
328 PetscCall(MatBackwardSolve(sFactor, s1, s2));
329 PetscCall(VecAXPY(s2, neg_one, x));
330 PetscCall(VecNorm(s2, NORM_2, &norm2));
331 if (10 * norm1 < norm2) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n", (double)norm2, bs));
332 }
333
334 /* test MatSolve() */
335 PetscCall(MatSolve(sFactor, b, y));
336 PetscCall(MatDestroy(&sFactor));
337 /* Check the error */
338 PetscCall(VecAXPY(y, neg_one, x));
339 PetscCall(VecNorm(y, NORM_2, &norm2));
340 if (10 * norm1 < norm2 && lf - inc != -1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n", lf - inc, lf, (double)norm1, (double)norm2));
341 norm1 = norm2;
342 if (norm2 < tol && lf != -1) break;
343 }
344
345 #if defined(PETSC_HAVE_MUMPS)
346
347 #if defined(PETSC_USE_REAL___FLOAT128)
348 tol = 1e-10; // since MUMPS is run in double
349 #endif
350
351 PetscCall(MatGetFactor(sA, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &sFactor));
352 PetscCall(MatCholeskyFactorSymbolic(sFactor, sA, NULL, NULL));
353 PetscCall(MatCholeskyFactorNumeric(sFactor, sA, NULL));
354 for (i = 0; i < 10; i++) {
355 PetscCall(VecSetRandom(b, rdm));
356 PetscCall(MatSolve(sFactor, b, y));
357 /* Check the error */
358 PetscCall(MatMult(sA, y, x));
359 PetscCall(VecAXPY(x, neg_one, b));
360 PetscCall(VecNorm(x, NORM_2, &norm2));
361 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve(), norm2: %g\n", (double)norm2));
362 }
363 PetscCall(MatDestroy(&sFactor));
364 #endif
365
366 PetscCall(ISDestroy(&perm));
367
368 PetscCall(MatDestroy(&A));
369 PetscCall(MatDestroy(&sB));
370 PetscCall(MatDestroy(&sA));
371 PetscCall(VecDestroy(&x));
372 PetscCall(VecDestroy(&y));
373 PetscCall(VecDestroy(&s1));
374 PetscCall(VecDestroy(&s2));
375 PetscCall(VecDestroy(&b));
376 PetscCall(PetscRandomDestroy(&rdm));
377
378 PetscCall(PetscFinalize());
379 return 0;
380 }
381
382 /*TEST
383
384 test:
385 args: -bs {{1 2 3 4 5 6 7 8}}
386 output_file: output/empty.out
387
388 TEST*/
389