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