Actual source code: ex23.c
2: static char help[] = "Tests the use of interface functions for MATIS matrices and conversion routines.\n";
4: #include <petscmat.h>
6: PetscErrorCode TestMatZeroRows(Mat, Mat, PetscBool, IS, PetscScalar);
7: PetscErrorCode CheckMat(Mat, Mat, PetscBool, const char *);
9: int main(int argc, char **args)
10: {
11: Mat A, B, A2, B2, T;
12: Mat Aee, Aeo, Aoe, Aoo;
13: Mat *mats, *Asub, *Bsub;
14: Vec x, y;
15: MatInfo info;
16: ISLocalToGlobalMapping cmap, rmap;
17: IS is, is2, reven, rodd, ceven, codd;
18: IS *rows, *cols;
19: IS irow[2], icol[2];
20: PetscLayout rlayout, clayout;
21: const PetscInt *rrange, *crange;
22: MatType lmtype;
23: PetscScalar diag = 2.;
24: PetscInt n, m, i, lm, ln;
25: PetscInt rst, ren, cst, cen, nr, nc;
26: PetscMPIInt rank, size, lrank, rrank;
27: PetscBool testT, squaretest, isaij;
28: PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE;
29: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric;
32: PetscInitialize(&argc, &args, (char *)0, help);
33: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
34: MPI_Comm_size(PETSC_COMM_WORLD, &size);
35: m = n = 2 * size;
36: PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL);
37: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
38: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
39: PetscOptionsGetBool(NULL, NULL, "-negmap", &negmap, NULL);
40: PetscOptionsGetBool(NULL, NULL, "-repmap", &repmap, NULL);
41: PetscOptionsGetBool(NULL, NULL, "-permmap", &permute, NULL);
42: PetscOptionsGetBool(NULL, NULL, "-diffmap", &diffmap, NULL);
47: /* create a MATIS matrix */
48: MatCreate(PETSC_COMM_WORLD, &A);
49: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n);
50: MatSetType(A, MATIS);
51: MatSetFromOptions(A);
52: if (!negmap && !repmap) {
53: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
54: Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
55: Equivalent to passing NULL for the mapping */
56: ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is);
57: } else if (negmap && !repmap) { /* non repeated but with negative indices */
58: ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is);
59: } else if (!negmap && repmap) { /* non negative but repeated indices */
60: IS isl[2];
62: ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]);
63: ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]);
64: ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is);
65: ISDestroy(&isl[0]);
66: ISDestroy(&isl[1]);
67: } else { /* negative and repeated indices */
68: IS isl[2];
70: ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]);
71: ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]);
72: ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is);
73: ISDestroy(&isl[0]);
74: ISDestroy(&isl[1]);
75: }
76: ISLocalToGlobalMappingCreateIS(is, &cmap);
77: ISDestroy(&is);
79: if (m != n || diffmap) {
80: ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is);
81: ISLocalToGlobalMappingCreateIS(is, &rmap);
82: ISDestroy(&is);
83: } else {
84: PetscObjectReference((PetscObject)cmap);
85: rmap = cmap;
86: }
88: MatSetLocalToGlobalMapping(A, rmap, cmap);
89: MatISStoreL2L(A, PETSC_FALSE);
90: MatISSetPreallocation(A, 3, NULL, 3, NULL);
91: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap)); /* I do not want to precompute the pattern */
92: ISLocalToGlobalMappingGetSize(rmap, &lm);
93: ISLocalToGlobalMappingGetSize(cmap, &ln);
94: for (i = 0; i < lm; i++) {
95: PetscScalar v[3];
96: PetscInt cols[3];
98: cols[0] = (i - 1 + n) % n;
99: cols[1] = i % n;
100: cols[2] = (i + 1) % n;
101: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
102: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
103: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
104: ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
105: MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES);
106: }
107: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
110: /* activate tests for square matrices with same maps only */
111: MatHasCongruentLayouts(A, &squaretest);
112: if (squaretest && rmap != cmap) {
113: PetscInt nr, nc;
115: ISLocalToGlobalMappingGetSize(rmap, &nr);
116: ISLocalToGlobalMappingGetSize(cmap, &nc);
117: if (nr != nc) squaretest = PETSC_FALSE;
118: else {
119: const PetscInt *idxs1, *idxs2;
121: ISLocalToGlobalMappingGetIndices(rmap, &idxs1);
122: ISLocalToGlobalMappingGetIndices(cmap, &idxs2);
123: PetscArraycmp(idxs1, idxs2, nr, &squaretest);
124: ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1);
125: ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2);
126: }
127: MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
128: }
130: /* test MatISGetLocalMat */
131: MatISGetLocalMat(A, &B);
132: MatGetType(B, &lmtype);
134: /* test MatGetInfo */
135: PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n");
136: MatGetInfo(A, MAT_LOCAL, &info);
137: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
138: PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "Process %2d: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", PetscGlobalRank, (PetscInt)info.nz_used, (PetscInt)info.nz_allocated,
139: (PetscInt)info.nz_unneeded, (PetscInt)info.assemblies, (PetscInt)info.mallocs));
140: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
141: MatGetInfo(A, MAT_GLOBAL_MAX, &info);
142: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "GlobalMax : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded,
143: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
144: MatGetInfo(A, MAT_GLOBAL_SUM, &info);
145: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "GlobalSum : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", (PetscInt)info.nz_used, (PetscInt)info.nz_allocated, (PetscInt)info.nz_unneeded,
146: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
148: /* test MatIsSymmetric */
149: MatIsSymmetric(A, 0.0, &issymmetric);
150: PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric);
152: /* Create a MPIAIJ matrix, same as A */
153: MatCreate(PETSC_COMM_WORLD, &B);
154: MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n);
155: MatSetType(B, MATAIJ);
156: MatSetFromOptions(B);
157: MatSetUp(B);
158: MatSetLocalToGlobalMapping(B, rmap, cmap);
159: MatMPIAIJSetPreallocation(B, 3, NULL, 3, NULL);
160: MatMPIBAIJSetPreallocation(B, 1, 3, NULL, 3, NULL);
161: #if defined(PETSC_HAVE_HYPRE)
162: MatHYPRESetPreallocation(B, 3, NULL, 3, NULL);
163: #endif
164: MatISSetPreallocation(B, 3, NULL, 3, NULL);
165: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap)); /* I do not want to precompute the pattern */
166: for (i = 0; i < lm; i++) {
167: PetscScalar v[3];
168: PetscInt cols[3];
170: cols[0] = (i - 1 + n) % n;
171: cols[1] = i % n;
172: cols[2] = (i + 1) % n;
173: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
174: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
175: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
176: ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
177: MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES);
178: }
179: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
180: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
182: /* test MatView */
183: PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n");
184: MatView(A, NULL);
185: MatView(B, NULL);
187: /* test CheckMat */
188: PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n");
189: CheckMat(A, B, PETSC_FALSE, "CheckMat");
191: /* test MatDuplicate and MatAXPY */
192: PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n");
193: MatDuplicate(A, MAT_COPY_VALUES, &A2);
194: CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY");
196: /* test MatConvert */
197: PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n");
198: MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2);
199: CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX");
200: MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2);
201: CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX");
202: MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2);
203: CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX");
204: MatDestroy(&A2);
205: MatDestroy(&B2);
206: PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n");
207: MatDuplicate(B, MAT_COPY_VALUES, &B2);
208: MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2);
209: CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX");
210: MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2);
211: CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX");
212: MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2);
213: CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX");
214: MatDestroy(&A2);
215: MatDestroy(&B2);
216: PetscStrcmp(lmtype, MATSEQAIJ, &isaij);
217: if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
218: PetscInt ri, ci, rr[3] = {0, 1, 0}, cr[4] = {1, 2, 0, 1}, rk[3] = {0, 2, 1}, ck[4] = {1, 0, 3, 2};
219: ISLocalToGlobalMapping tcmap, trmap;
221: for (ri = 0; ri < 2; ri++) {
222: PetscInt *r;
224: r = (PetscInt *)(ri == 0 ? rr : rk);
225: for (ci = 0; ci < 2; ci++) {
226: PetscInt *c, rb, cb;
228: c = (PetscInt *)(ci == 0 ? cr : ck);
229: for (rb = 1; rb < 4; rb++) {
230: ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is);
231: ISLocalToGlobalMappingCreateIS(is, &trmap);
232: ISDestroy(&is);
233: for (cb = 1; cb < 4; cb++) {
234: Mat T, lT, T2;
235: char testname[256];
237: PetscSNPrintf(testname, sizeof(testname), "MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")", ri, ci, rb, cb);
238: PetscPrintf(PETSC_COMM_WORLD, "Test %s\n", testname);
240: ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is);
241: ISLocalToGlobalMappingCreateIS(is, &tcmap);
242: ISDestroy(&is);
244: MatCreate(PETSC_COMM_SELF, &T);
245: MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4);
246: MatSetType(T, MATIS);
247: MatSetLocalToGlobalMapping(T, trmap, tcmap);
248: ISLocalToGlobalMappingDestroy(&tcmap);
249: MatISGetLocalMat(T, &lT);
250: MatSetType(lT, MATSEQAIJ);
251: MatSeqAIJSetPreallocation(lT, cb * 4, NULL);
252: MatSetRandom(lT, NULL);
253: MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT);
254: MatISRestoreLocalMat(T, &lT);
255: MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
256: MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);
258: MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2);
259: CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX");
260: MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2);
261: CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX");
262: MatDestroy(&T2);
263: MatDuplicate(T, MAT_COPY_VALUES, &T2);
264: MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2);
265: CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX");
266: MatDestroy(&T);
267: MatDestroy(&T2);
268: }
269: ISLocalToGlobalMappingDestroy(&trmap);
270: }
271: }
272: }
273: }
275: /* test MatDiagonalScale */
276: PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n");
277: MatDuplicate(A, MAT_COPY_VALUES, &A2);
278: MatDuplicate(B, MAT_COPY_VALUES, &B2);
279: MatCreateVecs(A, &x, &y);
280: VecSetRandom(x, NULL);
281: if (issymmetric) {
282: VecCopy(x, y);
283: } else {
284: VecSetRandom(y, NULL);
285: VecScale(y, 8.);
286: }
287: MatDiagonalScale(A2, y, x);
288: MatDiagonalScale(B2, y, x);
289: CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale");
290: MatDestroy(&A2);
291: MatDestroy(&B2);
292: VecDestroy(&x);
293: VecDestroy(&y);
295: /* test MatPtAP (A IS and B AIJ) */
296: if (isaij && m == n) {
297: PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n");
298: MatISStoreL2L(A, PETSC_TRUE);
299: MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &A2);
300: MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B2);
301: CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX");
302: MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &A2);
303: CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX");
304: MatDestroy(&A2);
305: MatDestroy(&B2);
306: }
308: /* test MatGetLocalSubMatrix */
309: PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n");
310: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2);
311: ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven);
312: ISComplement(reven, 0, lm, &rodd);
313: ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven);
314: ISComplement(ceven, 0, ln, &codd);
315: MatGetLocalSubMatrix(A2, reven, ceven, &Aee);
316: MatGetLocalSubMatrix(A2, reven, codd, &Aeo);
317: MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe);
318: MatGetLocalSubMatrix(A2, rodd, codd, &Aoo);
319: for (i = 0; i < lm; i++) {
320: PetscInt j, je, jo, colse[3], colso[3];
321: PetscScalar ve[3], vo[3];
322: PetscScalar v[3];
323: PetscInt cols[3];
324: PetscInt row = i / 2;
326: cols[0] = (i - 1 + n) % n;
327: cols[1] = i % n;
328: cols[2] = (i + 1) % n;
329: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
330: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
331: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
332: ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
333: for (j = 0, je = 0, jo = 0; j < 3; j++) {
334: if (cols[j] % 2) {
335: vo[jo] = v[j];
336: colso[jo++] = cols[j] / 2;
337: } else {
338: ve[je] = v[j];
339: colse[je++] = cols[j] / 2;
340: }
341: }
342: if (i % 2) {
343: MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES);
344: MatSetValuesLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES);
345: } else {
346: MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES);
347: MatSetValuesLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES);
348: }
349: }
350: MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee);
351: MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo);
352: MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe);
353: MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo);
354: ISDestroy(&reven);
355: ISDestroy(&ceven);
356: ISDestroy(&rodd);
357: ISDestroy(&codd);
358: MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
359: MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
360: MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN);
361: CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix");
362: MatDestroy(&A2);
364: /* test MatConvert_Nest_IS */
365: testT = PETSC_FALSE;
366: PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL);
368: PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n");
369: nr = 2;
370: nc = 2;
371: PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL);
372: PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL);
373: if (testT) {
374: MatGetOwnershipRange(A, &cst, &cen);
375: MatGetOwnershipRangeColumn(A, &rst, &ren);
376: } else {
377: MatGetOwnershipRange(A, &rst, &ren);
378: MatGetOwnershipRangeColumn(A, &cst, &cen);
379: }
380: PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats);
381: for (i = 0; i < nr * nc; i++) {
382: if (testT) {
383: MatCreateTranspose(A, &mats[i]);
384: MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]);
385: } else {
386: MatDuplicate(A, MAT_COPY_VALUES, &mats[i]);
387: MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]);
388: }
389: }
390: for (i = 0; i < nr; i++) ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]);
391: for (i = 0; i < nc; i++) ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]);
392: MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2);
393: MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2);
394: for (i = 0; i < nr; i++) ISDestroy(&rows[i]);
395: for (i = 0; i < nc; i++) ISDestroy(&cols[i]);
396: for (i = 0; i < 2 * nr * nc; i++) MatDestroy(&mats[i]);
397: PetscFree3(rows, cols, mats);
398: MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T);
399: MatDestroy(&B2);
400: MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2);
401: CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX");
402: MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2);
403: CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX");
404: MatDestroy(&B2);
405: MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2);
406: CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX");
407: MatDestroy(&T);
408: MatDestroy(&A2);
410: /* test MatCreateSubMatrix */
411: PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n");
412: if (rank == 0) {
413: ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is);
414: ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2);
415: } else if (rank == 1) {
416: ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is);
417: if (n > 3) {
418: ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2);
419: } else {
420: ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2);
421: }
422: } else if (rank == 2 && n > 4) {
423: ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
424: ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2);
425: } else {
426: ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
427: ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2);
428: }
429: MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2);
430: MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2);
431: CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix");
433: MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2);
434: MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2);
435: CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix");
436: MatDestroy(&A2);
437: MatDestroy(&B2);
439: if (!issymmetric) {
440: MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2);
441: MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2);
442: MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2);
443: MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2);
444: CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix");
445: }
447: MatDestroy(&A2);
448: MatDestroy(&B2);
449: ISDestroy(&is);
450: ISDestroy(&is2);
452: /* test MatCreateSubMatrices */
453: PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n");
454: MatGetLayouts(A, &rlayout, &clayout);
455: PetscLayoutGetRanges(rlayout, &rrange);
456: PetscLayoutGetRanges(clayout, &crange);
457: lrank = (size + rank - 1) % size;
458: rrank = (rank + 1) % size;
459: ISCreateStride(PETSC_COMM_SELF, (rrange[lrank + 1] - rrange[lrank]), rrange[lrank], 1, &irow[0]);
460: ISCreateStride(PETSC_COMM_SELF, (crange[rrank + 1] - crange[rrank]), crange[rrank], 1, &icol[0]);
461: ISCreateStride(PETSC_COMM_SELF, (rrange[rrank + 1] - rrange[rrank]), rrange[rrank], 1, &irow[1]);
462: ISCreateStride(PETSC_COMM_SELF, (crange[lrank + 1] - crange[lrank]), crange[lrank], 1, &icol[1]);
463: MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub);
464: MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub);
465: CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]");
466: CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]");
467: MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub);
468: MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub);
469: CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]");
470: CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]");
471: MatDestroySubMatrices(2, &Asub);
472: MatDestroySubMatrices(2, &Bsub);
473: ISDestroy(&irow[0]);
474: ISDestroy(&irow[1]);
475: ISDestroy(&icol[0]);
476: ISDestroy(&icol[1]);
478: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
479: if (size > 1) {
480: if (rank == 0) {
481: PetscInt st, len;
483: st = (m + 1) / 2;
484: len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
485: ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is);
486: } else {
487: ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
488: }
489: } else {
490: ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is);
491: }
493: if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
494: PetscInt *idx0, *idx1, n0, n1;
495: IS Ais[2], Bis[2];
497: /* test MatDiagonalSet */
498: PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n");
499: MatDuplicate(A, MAT_COPY_VALUES, &A2);
500: MatDuplicate(B, MAT_COPY_VALUES, &B2);
501: MatCreateVecs(A, NULL, &x);
502: VecSetRandom(x, NULL);
503: MatDiagonalSet(A2, x, INSERT_VALUES);
504: MatDiagonalSet(B2, x, INSERT_VALUES);
505: CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet");
506: VecDestroy(&x);
507: MatDestroy(&A2);
508: MatDestroy(&B2);
510: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
511: PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n");
512: MatDuplicate(A, MAT_COPY_VALUES, &A2);
513: MatDuplicate(B, MAT_COPY_VALUES, &B2);
514: MatShift(A2, 2.0);
515: MatShift(B2, 2.0);
516: CheckMat(A2, B2, PETSC_FALSE, "MatShift");
517: MatDestroy(&A2);
518: MatDestroy(&B2);
520: /* nonzero diag value is supported for square matrices only */
521: TestMatZeroRows(A, B, PETSC_TRUE, is, diag);
523: /* test MatIncreaseOverlap */
524: PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n");
525: MatGetOwnershipRange(A, &rst, &ren);
526: n0 = (ren - rst) / 2;
527: n1 = (ren - rst) / 3;
528: PetscMalloc1(n0, &idx0);
529: PetscMalloc1(n1, &idx1);
530: for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
531: for (i = 0; i < n1; i++) idx1[i] = rst + i;
532: ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]);
533: ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]);
534: ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]);
535: ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]);
536: MatIncreaseOverlap(A, 2, Ais, 3);
537: MatIncreaseOverlap(B, 2, Bis, 3);
538: /* Non deterministic output! */
539: ISSort(Ais[0]);
540: ISSort(Ais[1]);
541: ISSort(Bis[0]);
542: ISSort(Bis[1]);
543: ISView(Ais[0], NULL);
544: ISView(Bis[0], NULL);
545: ISView(Ais[1], NULL);
546: ISView(Bis[1], NULL);
547: MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub);
548: MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub);
549: CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]");
550: CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]");
551: MatDestroySubMatrices(2, &Asub);
552: MatDestroySubMatrices(2, &Bsub);
553: ISDestroy(&Ais[0]);
554: ISDestroy(&Ais[1]);
555: ISDestroy(&Bis[0]);
556: ISDestroy(&Bis[1]);
557: }
558: TestMatZeroRows(A, B, squaretest, is, 0.0);
559: ISDestroy(&is);
561: /* test MatTranspose */
562: PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n");
563: MatTranspose(A, MAT_INITIAL_MATRIX, &A2);
564: MatTranspose(B, MAT_INITIAL_MATRIX, &B2);
565: CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose");
567: MatTranspose(A, MAT_REUSE_MATRIX, &A2);
568: CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose");
569: MatDestroy(&A2);
571: MatDuplicate(A, MAT_COPY_VALUES, &A2);
572: MatTranspose(A2, MAT_INPLACE_MATRIX, &A2);
573: CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose");
574: MatDestroy(&A2);
576: MatTranspose(A, MAT_INITIAL_MATRIX, &A2);
577: CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose");
578: MatDestroy(&A2);
579: MatDestroy(&B2);
581: /* test MatISFixLocalEmpty */
582: if (isaij) {
583: PetscInt r[2];
585: r[0] = 0;
586: r[1] = PetscMin(m, n) - 1;
587: PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n");
588: MatDuplicate(A, MAT_COPY_VALUES, &A2);
590: MatISFixLocalEmpty(A2, PETSC_TRUE);
591: MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
592: MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
593: CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)");
595: MatZeroRows(A2, 2, r, 0.0, NULL, NULL);
596: MatViewFromOptions(A2, NULL, "-fixempty_view");
597: MatDuplicate(B, MAT_COPY_VALUES, &B2);
598: MatZeroRows(B2, 2, r, 0.0, NULL, NULL);
599: MatISFixLocalEmpty(A2, PETSC_TRUE);
600: MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
601: MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
602: MatViewFromOptions(A2, NULL, "-fixempty_view");
603: CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)");
604: MatDestroy(&A2);
606: MatDuplicate(A, MAT_COPY_VALUES, &A2);
607: MatZeroRows(A2, 2, r, 0.0, NULL, NULL);
608: MatTranspose(A2, MAT_INPLACE_MATRIX, &A2);
609: MatTranspose(B2, MAT_INPLACE_MATRIX, &B2);
610: MatViewFromOptions(A2, NULL, "-fixempty_view");
611: MatISFixLocalEmpty(A2, PETSC_TRUE);
612: MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
613: MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
614: MatViewFromOptions(A2, NULL, "-fixempty_view");
615: CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)");
617: MatDestroy(&A2);
618: MatDestroy(&B2);
620: if (squaretest) {
621: MatDuplicate(A, MAT_COPY_VALUES, &A2);
622: MatDuplicate(B, MAT_COPY_VALUES, &B2);
623: MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL);
624: MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL);
625: MatViewFromOptions(A2, NULL, "-fixempty_view");
626: MatISFixLocalEmpty(A2, PETSC_TRUE);
627: MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
628: MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
629: MatViewFromOptions(A2, NULL, "-fixempty_view");
630: CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)");
631: MatDestroy(&A2);
632: MatDestroy(&B2);
633: }
634: }
636: /* test MatInvertBlockDiagonal
637: special cases for block-diagonal matrices */
638: if (m == n) {
639: ISLocalToGlobalMapping map;
640: Mat Abd, Bbd;
641: IS is, bis;
642: const PetscScalar *isbd, *aijbd;
643: PetscScalar *vals;
644: const PetscInt *sts, *idxs;
645: PetscInt *idxs2, diff, perm, nl, bs, st, en, in;
646: PetscBool ok;
648: for (diff = 0; diff < 3; diff++) {
649: for (perm = 0; perm < 3; perm++) {
650: for (bs = 1; bs < 4; bs++) {
651: PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs);
652: PetscMalloc1(bs * bs, &vals);
653: MatGetOwnershipRanges(A, &sts);
654: switch (diff) {
655: case 1: /* inverted layout by processes */
656: in = 1;
657: st = sts[size - rank - 1];
658: en = sts[size - rank];
659: nl = en - st;
660: break;
661: case 2: /* round-robin layout */
662: in = size;
663: st = rank;
664: nl = n / size;
665: if (rank < n % size) nl++;
666: break;
667: default: /* same layout */
668: in = 1;
669: st = sts[rank];
670: en = sts[rank + 1];
671: nl = en - st;
672: break;
673: }
674: ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is);
675: ISGetLocalSize(is, &nl);
676: ISGetIndices(is, &idxs);
677: PetscMalloc1(nl, &idxs2);
678: for (i = 0; i < nl; i++) {
679: switch (perm) { /* invert some of the indices */
680: case 2:
681: idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
682: break;
683: case 1:
684: idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
685: break;
686: default:
687: idxs2[i] = idxs[i];
688: break;
689: }
690: }
691: ISRestoreIndices(is, &idxs);
692: ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis);
693: ISLocalToGlobalMappingCreateIS(bis, &map);
694: MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd);
695: ISLocalToGlobalMappingDestroy(&map);
696: MatISSetPreallocation(Abd, bs, NULL, 0, NULL);
697: for (i = 0; i < nl; i++) {
698: PetscInt b1, b2;
700: for (b1 = 0; b1 < bs; b1++) {
701: for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
702: }
703: MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES);
704: }
705: MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY);
706: MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY);
707: MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd);
708: MatInvertBlockDiagonal(Abd, &isbd);
709: MatInvertBlockDiagonal(Bbd, &aijbd);
710: MatGetLocalSize(Bbd, &nl, NULL);
711: ok = PETSC_TRUE;
712: for (i = 0; i < nl / bs; i++) {
713: PetscInt b1, b2;
715: for (b1 = 0; b1 < bs; b1++) {
716: for (b2 = 0; b2 < bs; b2++) {
717: if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
718: if (!ok) {
719: PetscPrintf(PETSC_COMM_SELF, "[%d] ERROR block %" PetscInt_FMT ", entry %" PetscInt_FMT " %" PetscInt_FMT ": %g %g\n", rank, i, b1, b2, (double)PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2]), (double)PetscAbsScalar(aijbd[i * bs * bs + b1 * bs + b2]));
720: break;
721: }
722: }
723: if (!ok) break;
724: }
725: if (!ok) break;
726: }
727: MatDestroy(&Abd);
728: MatDestroy(&Bbd);
729: PetscFree(vals);
730: ISDestroy(&is);
731: ISDestroy(&bis);
732: }
733: }
734: }
735: }
737: /* test MatGetDiagonalBlock */
738: PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n");
739: MatGetDiagonalBlock(A, &A2);
740: MatGetDiagonalBlock(B, &B2);
741: CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock");
742: MatScale(A, 2.0);
743: MatScale(B, 2.0);
744: MatGetDiagonalBlock(A, &A2);
745: MatGetDiagonalBlock(B, &B2);
746: CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock");
748: /* free testing matrices */
749: ISLocalToGlobalMappingDestroy(&cmap);
750: ISLocalToGlobalMappingDestroy(&rmap);
751: MatDestroy(&A);
752: MatDestroy(&B);
753: PetscFinalize();
754: return 0;
755: }
757: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
758: {
759: Mat Bcheck;
760: PetscReal error;
763: if (!usemult && B) {
764: PetscBool hasnorm;
766: MatHasOperation(B, MATOP_NORM, &hasnorm);
767: if (!hasnorm) usemult = PETSC_TRUE;
768: }
769: if (!usemult) {
770: if (B) {
771: MatType Btype;
773: MatGetType(B, &Btype);
774: MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck);
775: } else {
776: MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck);
777: }
778: if (B) { /* if B is present, subtract it */
779: MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN);
780: }
781: MatNorm(Bcheck, NORM_INFINITY, &error);
782: if (error > PETSC_SQRT_MACHINE_EPSILON) {
783: ISLocalToGlobalMapping rl2g, cl2g;
785: PetscObjectSetName((PetscObject)Bcheck, "Bcheck");
786: MatView(Bcheck, NULL);
787: if (B) {
788: PetscObjectSetName((PetscObject)B, "B");
789: MatView(B, NULL);
790: MatDestroy(&Bcheck);
791: MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck);
792: PetscObjectSetName((PetscObject)Bcheck, "Assembled A");
793: MatView(Bcheck, NULL);
794: }
795: MatDestroy(&Bcheck);
796: PetscObjectSetName((PetscObject)A, "A");
797: MatView(A, NULL);
798: MatGetLocalToGlobalMapping(A, &rl2g, &cl2g);
799: if (rl2g) ISLocalToGlobalMappingView(rl2g, NULL);
800: if (cl2g) ISLocalToGlobalMappingView(cl2g, NULL);
801: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
802: }
803: MatDestroy(&Bcheck);
804: } else {
805: PetscBool ok, okt;
807: MatMultEqual(A, B, 3, &ok);
808: MatMultTransposeEqual(A, B, 3, &okt);
810: }
811: return 0;
812: }
814: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
815: {
816: Mat B, Bcheck, B2 = NULL, lB;
817: Vec x = NULL, b = NULL, b2 = NULL;
818: ISLocalToGlobalMapping l2gr, l2gc;
819: PetscReal error;
820: char diagstr[16];
821: const PetscInt *idxs;
822: PetscInt rst, ren, i, n, N, d;
823: PetscMPIInt rank;
824: PetscBool miss, haszerorows;
827: if (diag == 0.) {
828: PetscStrcpy(diagstr, "zero");
829: } else {
830: PetscStrcpy(diagstr, "nonzero");
831: }
832: ISView(is, NULL);
833: MatGetLocalToGlobalMapping(A, &l2gr, &l2gc);
834: /* tests MatDuplicate and MatCopy */
835: if (diag == 0.) {
836: MatDuplicate(A, MAT_COPY_VALUES, &B);
837: } else {
838: MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B);
839: MatCopy(A, B, SAME_NONZERO_PATTERN);
840: }
841: MatISGetLocalMat(B, &lB);
842: MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows);
843: if (squaretest && haszerorows) {
844: MatCreateVecs(B, &x, &b);
845: MatDuplicate(B, MAT_COPY_VALUES, &B2);
846: VecSetLocalToGlobalMapping(b, l2gr);
847: VecSetLocalToGlobalMapping(x, l2gc);
848: VecSetRandom(x, NULL);
849: VecSetRandom(b, NULL);
850: /* mimic b[is] = x[is] */
851: VecDuplicate(b, &b2);
852: VecSetLocalToGlobalMapping(b2, l2gr);
853: VecCopy(b, b2);
854: ISGetLocalSize(is, &n);
855: ISGetIndices(is, &idxs);
856: VecGetSize(x, &N);
857: for (i = 0; i < n; i++) {
858: if (0 <= idxs[i] && idxs[i] < N) {
859: VecSetValue(b2, idxs[i], diag, INSERT_VALUES);
860: VecSetValue(x, idxs[i], 1., INSERT_VALUES);
861: }
862: }
863: VecAssemblyBegin(b2);
864: VecAssemblyEnd(b2);
865: VecAssemblyBegin(x);
866: VecAssemblyEnd(x);
867: ISRestoreIndices(is, &idxs);
868: /* test ZeroRows on MATIS */
869: PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr);
870: MatZeroRowsIS(B, is, diag, x, b);
871: PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr);
872: MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL);
873: } else if (haszerorows) {
874: /* test ZeroRows on MATIS */
875: PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr);
876: MatZeroRowsIS(B, is, diag, NULL, NULL);
877: b = b2 = x = NULL;
878: } else {
879: PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr);
880: b = b2 = x = NULL;
881: }
883: if (squaretest && haszerorows) {
884: VecAXPY(b2, -1., b);
885: VecNorm(b2, NORM_INFINITY, &error);
887: }
889: /* test MatMissingDiagonal */
890: PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n");
891: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
892: MatMissingDiagonal(B, &miss, &d);
893: MatGetOwnershipRange(B, &rst, &ren);
894: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
895: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n", rank, rst, ren, (int)miss, d, diagstr);
896: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
897: PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
899: VecDestroy(&x);
900: VecDestroy(&b);
901: VecDestroy(&b2);
903: /* check the result of ZeroRows with that from MPIAIJ routines
904: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
905: if (haszerorows) {
906: MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck);
907: MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
908: MatZeroRowsIS(Bcheck, is, diag, NULL, NULL);
909: CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows");
910: MatDestroy(&Bcheck);
911: }
912: MatDestroy(&B);
914: if (B2) { /* test MatZeroRowsColumns */
915: MatDuplicate(Afull, MAT_COPY_VALUES, &B);
916: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
917: MatZeroRowsColumnsIS(B, is, diag, NULL, NULL);
918: CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns");
919: MatDestroy(&B);
920: MatDestroy(&B2);
921: }
922: return 0;
923: }
925: /*TEST
927: test:
928: args: -test_trans
930: test:
931: suffix: 2
932: nsize: 4
933: args: -matis_convert_local_nest -nr 3 -nc 4
935: test:
936: suffix: 3
937: nsize: 5
938: args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1
940: test:
941: suffix: 4
942: nsize: 6
943: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
945: test:
946: suffix: 5
947: nsize: 6
948: args: -m 12 -n 12 -test_trans -nr 3 -nc 1
950: test:
951: suffix: 6
952: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
954: test:
955: suffix: 7
956: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
958: test:
959: suffix: 8
960: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
962: test:
963: suffix: 9
964: nsize: 5
965: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
967: test:
968: suffix: 10
969: nsize: 5
970: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
972: test:
973: suffix: vscat_default
974: nsize: 5
975: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
976: output_file: output/ex23_11.out
978: test:
979: suffix: 12
980: nsize: 3
981: args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3
983: testset:
984: output_file: output/ex23_13.out
985: nsize: 3
986: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
987: filter: grep -v "type:"
988: test:
989: suffix: baij
990: args: -matis_localmat_type baij
991: test:
992: requires: viennacl
993: suffix: viennacl
994: args: -matis_localmat_type aijviennacl
995: test:
996: requires: cuda
997: suffix: cusparse
998: args: -matis_localmat_type aijcusparse
1000: test:
1001: suffix: negrep
1002: nsize: {{1 3}separate output}
1003: args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output}
1005: TEST*/