Actual source code: ex74.c
1: static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n";
3: #include <petscmat.h>
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;
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &args, (char *)0, 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));
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));
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));
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"));
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;
68: value[0] = 0.1;
69: value[1] = -1;
70: value[2] = 2;
72: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
73: PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
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;
83: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
84: PetscCall(MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES));
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: }
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;
135: value[0] = -1.0;
136: value[1] = 4.0;
138: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
139: PetscCall(MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES));
141: i = 0 + block * bs;
142: col[0] = 0 + block * bs;
143: col[1] = 1 + block * bs;
145: value[0] = 4.0;
146: value[1] = -1.0;
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;
156: PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
157: PetscCall(MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES));
159: col[0] = i;
160: row = i + bs;
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));
169: PetscCall(MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY));
170: PetscCall(MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY));
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"));
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()");
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));
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"));
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"));
214: PetscCall(MatGetBlockSize(A, &Ii));
215: PetscCall(MatGetBlockSize(sB, &i));
216: if (i - Ii) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n"));
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));
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");
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));
242: {
243: PetscScalar alpha = 0.1;
244: PetscCall(MatScale(A, alpha));
245: PetscCall(MatScale(sB, alpha));
246: }
247: #endif
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"));
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: }
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: }
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_DEFAULT, &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));
289: /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
290: PetscCall(MatGetOrdering(A, MATORDERINGNATURAL, &perm, &iscol));
291: PetscCall(ISDestroy(&iscol));
292: norm1 = tol;
293: inc = bs;
295: /* initialize factinfo */
296: PetscCall(PetscMemzero(&factinfo, sizeof(MatFactorInfo)));
298: for (lf = -1; lf < 10; lf += inc) {
299: if (lf == -1) { /* Cholesky factor of sB (duplicate sA) */
300: factinfo.fill = 5.0;
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;
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); */
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: */
323: PetscCall(MatMult(sB, x, b));
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: }
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: }
345: #if defined(PETSC_HAVE_MUMPS)
346: PetscCall(MatGetFactor(sA, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &sFactor));
347: PetscCall(MatCholeskyFactorSymbolic(sFactor, sA, NULL, NULL));
348: PetscCall(MatCholeskyFactorNumeric(sFactor, sA, NULL));
349: for (i = 0; i < 10; i++) {
350: PetscCall(VecSetRandom(b, rdm));
351: PetscCall(MatSolve(sFactor, b, y));
352: /* Check the error */
353: PetscCall(MatMult(sA, y, x));
354: PetscCall(VecAXPY(x, neg_one, b));
355: PetscCall(VecNorm(x, NORM_2, &norm2));
356: if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve(), norm2: %g\n", (double)norm2));
357: }
358: PetscCall(MatDestroy(&sFactor));
359: #endif
361: PetscCall(ISDestroy(&perm));
363: PetscCall(MatDestroy(&A));
364: PetscCall(MatDestroy(&sB));
365: PetscCall(MatDestroy(&sA));
366: PetscCall(VecDestroy(&x));
367: PetscCall(VecDestroy(&y));
368: PetscCall(VecDestroy(&s1));
369: PetscCall(VecDestroy(&s2));
370: PetscCall(VecDestroy(&b));
371: PetscCall(PetscRandomDestroy(&rdm));
373: PetscCall(PetscFinalize());
374: return 0;
375: }
377: /*TEST
379: test:
380: args: -bs {{1 2 3 4 5 6 7 8}}
382: TEST*/