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