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