Actual source code: ex75.c
2: static char help[] = "Tests various routines in MatMPISBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Vec x, y, u, s1, s2;
9: Mat A, sA, sB;
10: PetscRandom rctx;
11: PetscReal r1, r2, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON;
12: PetscScalar one = 1.0, neg_one = -1.0, value[3], four = 4.0, alpha = 0.1;
13: PetscInt n, col[3], n1, block, row, i, j, i2, j2, Ii, J, rstart, rend, bs = 1, mbs = 16, d_nz = 3, o_nz = 3, prob = 1;
14: PetscMPIInt size, rank;
15: PetscBool flg;
18: PetscInitialize(&argc, &args, (char *)0, help);
19: PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL);
20: PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
21: PetscOptionsGetInt(NULL, NULL, "-prob", &prob, NULL);
23: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
24: MPI_Comm_size(PETSC_COMM_WORLD, &size);
26: /* Create a BAIJ matrix A */
27: n = mbs * bs;
28: MatCreate(PETSC_COMM_WORLD, &A);
29: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
30: MatSetType(A, MATBAIJ);
31: MatSetFromOptions(A);
32: MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL);
33: MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL);
34: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
36: if (bs == 1) {
37: if (prob == 1) { /* tridiagonal matrix */
38: value[0] = -1.0;
39: value[1] = 0.0;
40: value[2] = -1.0;
41: for (i = 1; i < n - 1; i++) {
42: col[0] = i - 1;
43: col[1] = i;
44: col[2] = i + 1;
45: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
46: }
47: i = n - 1;
48: col[0] = 0;
49: col[1] = n - 2;
50: col[2] = n - 1;
51: value[0] = 0.1;
52: value[1] = -1.0;
53: value[2] = 0.0;
54: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
56: i = 0;
57: col[0] = 0;
58: col[1] = 1;
59: col[2] = n - 1;
60: value[0] = 0.0;
61: value[1] = -1.0;
62: value[2] = 0.1;
63: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
64: } else if (prob == 2) { /* matrix for the five point stencil */
65: n1 = (int)PetscSqrtReal((PetscReal)n);
66: for (i = 0; i < n1; i++) {
67: for (j = 0; j < n1; j++) {
68: Ii = j + n1 * i;
69: if (i > 0) {
70: J = Ii - n1;
71: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
72: }
73: if (i < n1 - 1) {
74: J = Ii + n1;
75: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
76: }
77: if (j > 0) {
78: J = Ii - 1;
79: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
80: }
81: if (j < n1 - 1) {
82: J = Ii + 1;
83: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
84: }
85: MatSetValues(A, 1, &Ii, 1, &Ii, &four, INSERT_VALUES);
86: }
87: }
88: }
89: /* end of if (bs == 1) */
90: } else { /* bs > 1 */
91: for (block = 0; block < n / bs; block++) {
92: /* diagonal blocks */
93: value[0] = -1.0;
94: value[1] = 4.0;
95: value[2] = -1.0;
96: for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
97: col[0] = i - 1;
98: col[1] = i;
99: col[2] = i + 1;
100: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
101: }
102: i = bs - 1 + block * bs;
103: col[0] = bs - 2 + block * bs;
104: col[1] = bs - 1 + block * bs;
105: value[0] = -1.0;
106: value[1] = 4.0;
107: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
109: i = 0 + block * bs;
110: col[0] = 0 + block * bs;
111: col[1] = 1 + block * bs;
112: value[0] = 4.0;
113: value[1] = -1.0;
114: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
115: }
116: /* off-diagonal blocks */
117: value[0] = -1.0;
118: for (i = 0; i < (n / bs - 1) * bs; i++) {
119: col[0] = i + bs;
120: MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
121: col[0] = i;
122: row = i + bs;
123: MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
124: }
125: }
126: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
128: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
130: /* Get SBAIJ matrix sA from A */
131: MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA);
133: /* Test MatGetSize(), MatGetLocalSize() */
134: MatGetSize(sA, &i, &j);
135: MatGetSize(A, &i2, &j2);
136: i -= i2;
137: j -= j2;
138: if (i || j) {
139: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetSize()\n", rank);
140: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
141: }
143: MatGetLocalSize(sA, &i, &j);
144: MatGetLocalSize(A, &i2, &j2);
145: i2 -= i;
146: j2 -= j;
147: if (i2 || j2) {
148: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetLocalSize()\n", rank);
149: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
150: }
152: /* vectors */
153: /*--------------------*/
154: /* i is obtained from MatGetLocalSize() */
155: VecCreate(PETSC_COMM_WORLD, &x);
156: VecSetSizes(x, i, PETSC_DECIDE);
157: VecSetFromOptions(x);
158: VecDuplicate(x, &y);
159: VecDuplicate(x, &u);
160: VecDuplicate(x, &s1);
161: VecDuplicate(x, &s2);
163: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
164: PetscRandomSetFromOptions(rctx);
165: VecSetRandom(x, rctx);
166: VecSet(u, one);
168: /* Test MatNorm() */
169: MatNorm(A, NORM_FROBENIUS, &r1);
170: MatNorm(sA, NORM_FROBENIUS, &r2);
171: rnorm = PetscAbsReal(r1 - r2) / r2;
172: if (rnorm > tol && rank == 0) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs);
173: MatNorm(A, NORM_INFINITY, &r1);
174: MatNorm(sA, NORM_INFINITY, &r2);
175: rnorm = PetscAbsReal(r1 - r2) / r2;
176: if (rnorm > tol && rank == 0) PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs);
177: MatNorm(A, NORM_1, &r1);
178: MatNorm(sA, NORM_1, &r2);
179: rnorm = PetscAbsReal(r1 - r2) / r2;
180: if (rnorm > tol && rank == 0) PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs);
182: /* Test MatGetOwnershipRange() */
183: MatGetOwnershipRange(sA, &rstart, &rend);
184: MatGetOwnershipRange(A, &i2, &j2);
185: i2 -= rstart;
186: j2 -= rend;
187: if (i2 || j2) {
188: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank);
189: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
190: }
192: /* Test MatDiagonalScale() */
193: MatDiagonalScale(A, x, x);
194: MatDiagonalScale(sA, x, x);
195: MatMultEqual(A, sA, 10, &flg);
198: /* Test MatGetDiagonal(), MatScale() */
199: MatGetDiagonal(A, s1);
200: MatGetDiagonal(sA, s2);
201: VecNorm(s1, NORM_1, &r1);
202: VecNorm(s2, NORM_1, &r2);
203: r1 -= r2;
204: if (r1 < -tol || r1 > tol) {
205: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n", rank, (double)r1);
206: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
207: }
209: MatScale(A, alpha);
210: MatScale(sA, alpha);
212: /* Test MatGetRowMaxAbs() */
213: MatGetRowMaxAbs(A, s1, NULL);
214: MatGetRowMaxAbs(sA, s2, NULL);
216: VecNorm(s1, NORM_1, &r1);
217: VecNorm(s2, NORM_1, &r2);
218: r1 -= r2;
219: if (r1 < -tol || r1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetRowMaxAbs() \n");
221: /* Test MatMult(), MatMultAdd() */
222: MatMultEqual(A, sA, 10, &flg);
223: if (!flg) {
224: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale()\n", rank);
225: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
226: }
228: MatMultAddEqual(A, sA, 10, &flg);
229: if (!flg) {
230: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd()\n", rank);
231: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
232: }
234: /* Test MatMultTranspose(), MatMultTransposeAdd() */
235: for (i = 0; i < 10; i++) {
236: VecSetRandom(x, rctx);
237: MatMultTranspose(A, x, s1);
238: MatMultTranspose(sA, x, s2);
239: VecNorm(s1, NORM_1, &r1);
240: VecNorm(s2, NORM_1, &r2);
241: r1 -= r2;
242: if (r1 < -tol || r1 > tol) {
243: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale(), err=%g\n", rank, (double)r1);
244: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
245: }
246: }
247: for (i = 0; i < 10; i++) {
248: VecSetRandom(x, rctx);
249: VecSetRandom(y, rctx);
250: MatMultTransposeAdd(A, x, y, s1);
251: MatMultTransposeAdd(sA, x, y, s2);
252: VecNorm(s1, NORM_1, &r1);
253: VecNorm(s2, NORM_1, &r2);
254: r1 -= r2;
255: if (r1 < -tol || r1 > tol) {
256: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd(), err=%g \n", rank, (double)r1);
257: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
258: }
259: }
261: /* Test MatDuplicate() */
262: MatDuplicate(sA, MAT_COPY_VALUES, &sB);
263: MatEqual(sA, sB, &flg);
264: if (!flg) PetscPrintf(PETSC_COMM_WORLD, " Error in MatDuplicate(), sA != sB \n");
265: MatMultEqual(sA, sB, 5, &flg);
266: if (!flg) {
267: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMult()\n", rank);
268: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
269: }
270: MatMultAddEqual(sA, sB, 5, &flg);
271: if (!flg) {
272: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMultAdd(()\n", rank);
273: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
274: }
275: MatDestroy(&sB);
276: VecDestroy(&u);
277: VecDestroy(&x);
278: VecDestroy(&y);
279: VecDestroy(&s1);
280: VecDestroy(&s2);
281: MatDestroy(&sA);
282: MatDestroy(&A);
283: PetscRandomDestroy(&rctx);
284: PetscFinalize();
285: return 0;
286: }
288: /*TEST
290: test:
291: nsize: {{1 3}}
292: args: -bs {{1 2 3 5 7 8}} -mat_ignore_lower_triangular -prob {{1 2}}
294: TEST*/