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;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
20: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-prob", &prob, NULL));
23: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
24: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
26: /* Create a BAIJ matrix A */
27: n = mbs * bs;
28: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
29: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
30: PetscCall(MatSetType(A, MATBAIJ));
31: PetscCall(MatSetFromOptions(A));
32: PetscCall(MatMPIBAIJSetPreallocation(A, bs, d_nz, NULL, o_nz, NULL));
33: PetscCall(MatSeqBAIJSetPreallocation(A, bs, d_nz, NULL));
34: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
72: }
73: if (i < n1 - 1) {
74: J = Ii + n1;
75: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
76: }
77: if (j > 0) {
78: J = Ii - 1;
79: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
80: }
81: if (j < n1 - 1) {
82: J = Ii + 1;
83: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES));
84: }
85: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES));
121: col[0] = i;
122: row = i + bs;
123: PetscCall(MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES));
124: }
125: }
126: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
127: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
128: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
130: /* Get SBAIJ matrix sA from A */
131: PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
133: /* Test MatGetSize(), MatGetLocalSize() */
134: PetscCall(MatGetSize(sA, &i, &j));
135: PetscCall(MatGetSize(A, &i2, &j2));
136: i -= i2;
137: j -= j2;
138: if (i || j) {
139: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetSize()\n", rank));
140: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
141: }
143: PetscCall(MatGetLocalSize(sA, &i, &j));
144: PetscCall(MatGetLocalSize(A, &i2, &j2));
145: i2 -= i;
146: j2 -= j;
147: if (i2 || j2) {
148: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatGetLocalSize()\n", rank));
149: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
150: }
152: /* vectors */
153: /*--------------------*/
154: /* i is obtained from MatGetLocalSize() */
155: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
156: PetscCall(VecSetSizes(x, i, PETSC_DECIDE));
157: PetscCall(VecSetFromOptions(x));
158: PetscCall(VecDuplicate(x, &y));
159: PetscCall(VecDuplicate(x, &u));
160: PetscCall(VecDuplicate(x, &s1));
161: PetscCall(VecDuplicate(x, &s2));
163: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
164: PetscCall(PetscRandomSetFromOptions(rctx));
165: PetscCall(VecSetRandom(x, rctx));
166: PetscCall(VecSet(u, one));
168: /* Test MatNorm() */
169: PetscCall(MatNorm(A, NORM_FROBENIUS, &r1));
170: PetscCall(MatNorm(sA, NORM_FROBENIUS, &r2));
171: rnorm = PetscAbsReal(r1 - r2) / r2;
172: if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
173: PetscCall(MatNorm(A, NORM_INFINITY, &r1));
174: PetscCall(MatNorm(sA, NORM_INFINITY, &r2));
175: rnorm = PetscAbsReal(r1 - r2) / r2;
176: if (rnorm > tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%" PetscInt_FMT "\n", (double)r1, (double)r2, bs));
177: PetscCall(MatNorm(A, NORM_1, &r1));
178: PetscCall(MatNorm(sA, NORM_1, &r2));
179: rnorm = PetscAbsReal(r1 - r2) / r2;
180: if (rnorm > tol && rank == 0) PetscCall(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: PetscCall(MatGetOwnershipRange(sA, &rstart, &rend));
184: PetscCall(MatGetOwnershipRange(A, &i2, &j2));
185: i2 -= rstart;
186: j2 -= rend;
187: if (i2 || j2) {
188: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MaGetOwnershipRange()\n", rank));
189: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
190: }
192: /* Test MatDiagonalScale() */
193: PetscCall(MatDiagonalScale(A, x, x));
194: PetscCall(MatDiagonalScale(sA, x, x));
195: PetscCall(MatMultEqual(A, sA, 10, &flg));
196: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Error in MatDiagonalScale");
198: /* Test MatGetDiagonal(), MatScale() */
199: PetscCall(MatGetDiagonal(A, s1));
200: PetscCall(MatGetDiagonal(sA, s2));
201: PetscCall(VecNorm(s1, NORM_1, &r1));
202: PetscCall(VecNorm(s2, NORM_1, &r2));
203: r1 -= r2;
204: if (r1 < -tol || r1 > tol) {
205: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n", rank, (double)r1));
206: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
207: }
209: PetscCall(MatScale(A, alpha));
210: PetscCall(MatScale(sA, alpha));
212: /* Test MatGetRowMaxAbs() */
213: PetscCall(MatGetRowMaxAbs(A, s1, NULL));
214: PetscCall(MatGetRowMaxAbs(sA, s2, NULL));
216: PetscCall(VecNorm(s1, NORM_1, &r1));
217: PetscCall(VecNorm(s2, NORM_1, &r2));
218: r1 -= r2;
219: if (r1 < -tol || r1 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetRowMaxAbs() \n"));
221: /* Test MatMult(), MatMultAdd() */
222: PetscCall(MatMultEqual(A, sA, 10, &flg));
223: if (!flg) {
224: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale()\n", rank));
225: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
226: }
228: PetscCall(MatMultAddEqual(A, sA, 10, &flg));
229: if (!flg) {
230: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd()\n", rank));
231: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
232: }
234: /* Test MatMultTranspose(), MatMultTransposeAdd() */
235: for (i = 0; i < 10; i++) {
236: PetscCall(VecSetRandom(x, rctx));
237: PetscCall(MatMultTranspose(A, x, s1));
238: PetscCall(MatMultTranspose(sA, x, s2));
239: PetscCall(VecNorm(s1, NORM_1, &r1));
240: PetscCall(VecNorm(s2, NORM_1, &r2));
241: r1 -= r2;
242: if (r1 < -tol || r1 > tol) {
243: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMult() or MatScale(), err=%g\n", rank, (double)r1));
244: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
245: }
246: }
247: for (i = 0; i < 10; i++) {
248: PetscCall(VecSetRandom(x, rctx));
249: PetscCall(VecSetRandom(y, rctx));
250: PetscCall(MatMultTransposeAdd(A, x, y, s1));
251: PetscCall(MatMultTransposeAdd(sA, x, y, s2));
252: PetscCall(VecNorm(s1, NORM_1, &r1));
253: PetscCall(VecNorm(s2, NORM_1, &r2));
254: r1 -= r2;
255: if (r1 < -tol || r1 > tol) {
256: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatMultAdd(), err=%g \n", rank, (double)r1));
257: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
258: }
259: }
261: /* Test MatDuplicate() */
262: PetscCall(MatDuplicate(sA, MAT_COPY_VALUES, &sB));
263: PetscCall(MatEqual(sA, sB, &flg));
264: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error in MatDuplicate(), sA != sB \n"));
265: PetscCall(MatMultEqual(sA, sB, 5, &flg));
266: if (!flg) {
267: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMult()\n", rank));
268: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
269: }
270: PetscCall(MatMultAddEqual(sA, sB, 5, &flg));
271: if (!flg) {
272: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d], Error: MatDuplicate() or MatMultAdd(()\n", rank));
273: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
274: }
275: PetscCall(MatDestroy(&sB));
276: PetscCall(VecDestroy(&u));
277: PetscCall(VecDestroy(&x));
278: PetscCall(VecDestroy(&y));
279: PetscCall(VecDestroy(&s1));
280: PetscCall(VecDestroy(&s2));
281: PetscCall(MatDestroy(&sA));
282: PetscCall(MatDestroy(&A));
283: PetscCall(PetscRandomDestroy(&rctx));
284: PetscCall(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*/