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, PetscBool);
6: PetscErrorCode CheckMat(Mat, Mat, PetscBool, const char *);
7: PetscErrorCode ISL2GMapNoNeg(ISLocalToGlobalMapping, IS, IS *);
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, lis, is2, reven, rodd, ceven, codd;
18: IS *rows, *cols;
19: IS irow[2], icol[2];
20: PetscLayout rlayout, clayout;
21: const PetscInt *rrange, *crange, *idxs1, *idxs2;
22: MatType lmtype;
23: PetscScalar diag = 2., *vals;
24: PetscInt n, m, i, lm, ln;
25: PetscInt rst, ren, cst, cen, nr, nc, rbs = 1, cbs = 1;
26: PetscMPIInt rank, size, lrank, rrank;
27: PetscBool testT, squaretest, isaij;
28: PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE, allow_repeated = PETSC_TRUE;
29: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric, test_matlab = PETSC_FALSE, test_setvalues = PETSC_TRUE;
31: PetscFunctionBeginUser;
32: PetscCall(PetscInitialize(&argc, &args, NULL, 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: PetscCall(PetscOptionsGetBool(NULL, NULL, "-allow_repeated", &allow_repeated, NULL));
44: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matlab", &test_matlab, NULL));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_setvalues", &test_setvalues, NULL));
46: PetscCall(PetscOptionsGetInt(NULL, NULL, "-rbs", &rbs, NULL));
47: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cbs", &cbs, NULL));
48: PetscCheck(size == 1 || m >= 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 4 for parallel runs");
49: PetscCheck(size != 1 || m >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 2 for uniprocessor runs");
50: PetscCheck(n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of cols should be larger or equal 2");
52: /* create a MATIS matrix */
53: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
54: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
55: PetscCall(MatSetType(A, MATIS));
56: PetscCall(MatSetFromOptions(A));
57: if (!negmap && !repmap) {
58: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
59: Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
60: Equivalent to passing NULL for the mapping */
61: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is));
62: } else if (negmap && !repmap) { /* non repeated but with negative indices */
63: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is));
64: } else if (!negmap && repmap) { /* non negative but repeated indices */
65: IS isl[2];
67: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]));
68: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]));
69: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
70: PetscCall(ISDestroy(&isl[0]));
71: PetscCall(ISDestroy(&isl[1]));
72: } else { /* negative and repeated indices */
73: IS isl[2];
75: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]));
76: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]));
77: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
78: PetscCall(ISDestroy(&isl[0]));
79: PetscCall(ISDestroy(&isl[1]));
80: }
81: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmap));
82: PetscCall(ISDestroy(&is));
84: if (m != n || diffmap) {
85: PetscCall(ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is));
86: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmap));
87: PetscCall(ISDestroy(&is));
88: } else {
89: PetscCall(PetscObjectReference((PetscObject)cmap));
90: rmap = cmap;
91: }
93: PetscCall(MatISSetAllowRepeated(A, allow_repeated));
94: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
95: PetscCall(MatSetBlockSizes(A, rbs, cbs));
96: PetscCall(MatISStoreL2L(A, PETSC_FALSE));
97: PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));
98: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
99: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
100: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
101: for (i = 0; i < lm; i++) {
102: PetscScalar v[3];
103: PetscInt cols[3];
105: cols[0] = (i - 1 + n) % n;
106: cols[1] = i % n;
107: cols[2] = (i + 1) % n;
108: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
109: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
110: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
111: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
112: PetscCall(MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES));
113: }
114: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
115: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
117: /* activate tests for square matrices with same maps only */
118: PetscCall(MatHasCongruentLayouts(A, &squaretest));
119: if (squaretest && rmap != cmap) {
120: PetscInt nr, nc;
122: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nr));
123: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nc));
124: if (nr != nc) squaretest = PETSC_FALSE;
125: else {
126: PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
127: PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
128: PetscCall(PetscArraycmp(idxs1, idxs2, nr, &squaretest));
129: PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
130: PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
131: }
132: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
133: }
134: if (negmap && repmap) squaretest = PETSC_FALSE;
136: /* test MatISGetLocalMat */
137: PetscCall(MatISGetLocalMat(A, &B));
138: PetscCall(MatGetType(B, &lmtype));
140: /* test MatGetInfo */
141: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n"));
142: PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
143: PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
144: 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,
145: (PetscInt)info.nz_unneeded, (PetscInt)info.assemblies, (PetscInt)info.mallocs));
146: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
147: PetscCall(MatGetInfo(A, MAT_GLOBAL_MAX, &info));
148: 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,
149: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
150: PetscCall(MatGetInfo(A, MAT_GLOBAL_SUM, &info));
151: 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,
152: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
154: /* test MatIsSymmetric */
155: PetscCall(MatIsSymmetric(A, 0.0, &issymmetric));
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric));
158: /* Create a MPIAIJ matrix, same as A */
159: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
160: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
161: PetscCall(MatSetBlockSizes(B, rbs, cbs));
162: PetscCall(MatSetType(B, MATAIJ));
163: PetscCall(MatSetFromOptions(B));
164: PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
165: PetscCall(MatMPIAIJSetPreallocation(B, 3, NULL, 3, NULL));
166: PetscCall(MatMPIBAIJSetPreallocation(B, 1, 3, NULL, 3, NULL));
167: #if defined(PETSC_HAVE_HYPRE)
168: PetscCall(MatHYPRESetPreallocation(B, 3, NULL, 3, NULL));
169: #endif
170: PetscCall(MatISSetPreallocation(B, 3, NULL, 3, NULL));
171: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
172: for (i = 0; i < lm; i++) {
173: PetscScalar v[3];
174: PetscInt cols[3];
176: cols[0] = (i - 1 + n) % n;
177: cols[1] = i % n;
178: cols[2] = (i + 1) % n;
179: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
180: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
181: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
182: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
183: PetscCall(MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES));
184: }
185: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
188: /* test MatView */
189: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n"));
190: PetscCall(MatView(A, NULL));
191: PetscCall(MatView(B, NULL));
193: /* test MATLAB ASCII view */
194: if (test_matlab) { /* output is different when using real or complex numbers */
195: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView ASCII MATLAB\n"));
196: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB));
197: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
198: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
199: }
201: /* test CheckMat */
202: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n"));
203: PetscCall(CheckMat(A, B, PETSC_FALSE, "CheckMat"));
205: /* test binary MatView/MatLoad */
206: {
207: PetscMPIInt color = rank % 2;
208: MPI_Comm comm;
209: char name[PETSC_MAX_PATH_LEN];
210: PetscViewer wview, cview, sview, view;
211: Mat A2;
213: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, rank, &comm));
215: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "world_is"));
216: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_WRITE, &wview));
217: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "seq_is_%d", rank));
218: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_WRITE, &sview));
219: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "color_is_%d", color));
220: PetscCall(PetscViewerBinaryOpen(comm, name, FILE_MODE_WRITE, &cview));
221: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary world\n"));
222: PetscCall(MatView(A, wview));
223: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary self\n"));
224: PetscCall(MatView(A, sview));
225: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary subcomm\n"));
226: PetscCall(MatView(A, cview));
227: PetscCall(PetscViewerDestroy(&wview));
228: PetscCall(PetscViewerDestroy(&cview));
229: PetscCall(PetscViewerDestroy(&sview));
231: /* Load a world matrix */
232: PetscCall(MatCreate(PETSC_COMM_WORLD, &A2));
233: PetscCall(MatSetType(A2, MATIS));
234: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
236: /* Read back the same matrix and check */
237: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from world\n"));
238: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "world_is", FILE_MODE_READ, &view));
239: PetscCall(MatLoad(A2, view));
240: PetscCall(CheckMat(A, A2, PETSC_TRUE, "Load"));
241: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
242: PetscCall(PetscViewerDestroy(&view));
244: /* Read the matrix from rank 0 only */
245: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from self\n"));
246: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "seq_is_0", FILE_MODE_READ, &view));
247: PetscCall(MatLoad(A2, view));
248: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
249: PetscCall(PetscViewerDestroy(&view));
251: /* Read the matrix from subcomm */
252: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from subcomm\n"));
253: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "color_is_0", FILE_MODE_READ, &view));
254: PetscCall(MatLoad(A2, view));
255: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
256: PetscCall(PetscViewerDestroy(&view));
258: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
259: PetscCall(MatDestroy(&A2));
261: /* now load the original matrix from color 0 only processes */
262: if (!color) {
263: PetscCall(PetscPrintf(comm, "Test subcomm MatLoad from world\n"));
264: PetscCall(MatCreate(comm, &A2));
265: PetscCall(MatSetType(A2, MATIS));
266: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm), PETSC_VIEWER_ASCII_INFO_DETAIL));
267: PetscCall(PetscViewerBinaryOpen(comm, "world_is", FILE_MODE_READ, &view));
268: PetscCall(MatLoad(A2, view));
269: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_(comm)));
270: PetscCall(PetscViewerDestroy(&view));
271: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm)));
272: PetscCall(MatDestroy(&A2));
273: }
275: PetscCallMPI(MPI_Comm_free(&comm));
276: }
278: /* test MatDuplicate and MatAXPY */
279: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n"));
280: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
281: PetscCall(CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY"));
283: /* test MatConvert */
284: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n"));
285: PetscCall(MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2));
286: PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX"));
287: PetscCall(MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2));
288: PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX"));
289: PetscCall(MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2));
290: PetscCall(CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX"));
291: PetscCall(MatDestroy(&A2));
292: PetscCall(MatDestroy(&B2));
293: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n"));
294: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
295: PetscCall(MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2));
296: PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX"));
297: PetscCall(MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2));
298: PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX"));
299: PetscCall(MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2));
300: PetscCall(CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX"));
301: PetscCall(MatDestroy(&A2));
302: PetscCall(MatDestroy(&B2));
303: PetscCall(PetscStrcmp(lmtype, MATSEQAIJ, &isaij));
304: if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
305: 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};
306: ISLocalToGlobalMapping tcmap, trmap;
308: for (ri = 0; ri < 2; ri++) {
309: PetscInt *r;
311: r = (PetscInt *)(ri == 0 ? rr : rk);
312: for (ci = 0; ci < 2; ci++) {
313: PetscInt *c, rb, cb;
315: c = (PetscInt *)(ci == 0 ? cr : ck);
316: for (rb = 1; rb < 4; rb++) {
317: PetscCall(ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is));
318: PetscCall(ISLocalToGlobalMappingCreateIS(is, &trmap));
319: PetscCall(ISDestroy(&is));
320: for (cb = 1; cb < 4; cb++) {
321: Mat T, lT, T2;
322: char testname[256];
324: PetscCall(PetscSNPrintf(testname, sizeof(testname), "MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")", ri, ci, rb, cb));
325: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test %s\n", testname));
327: PetscCall(ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is));
328: PetscCall(ISLocalToGlobalMappingCreateIS(is, &tcmap));
329: PetscCall(ISDestroy(&is));
331: PetscCall(MatCreate(PETSC_COMM_SELF, &T));
332: PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4));
333: PetscCall(MatSetType(T, MATIS));
334: PetscCall(MatSetLocalToGlobalMapping(T, trmap, tcmap));
335: PetscCall(ISLocalToGlobalMappingDestroy(&tcmap));
336: PetscCall(MatISGetLocalMat(T, &lT));
337: PetscCall(MatSetType(lT, MATSEQAIJ));
338: PetscCall(MatSeqAIJSetPreallocation(lT, cb * 4, NULL));
339: PetscCall(MatSetRandom(lT, NULL));
340: PetscCall(MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT));
341: PetscCall(MatISRestoreLocalMat(T, &lT));
342: PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
343: PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
345: PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2));
346: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX"));
347: PetscCall(MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2));
348: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX"));
349: PetscCall(MatDestroy(&T2));
350: PetscCall(MatDuplicate(T, MAT_COPY_VALUES, &T2));
351: PetscCall(MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2));
352: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX"));
353: PetscCall(MatDestroy(&T));
354: PetscCall(MatDestroy(&T2));
355: }
356: PetscCall(ISLocalToGlobalMappingDestroy(&trmap));
357: }
358: }
359: }
360: }
362: /* test MatDiagonalScale */
363: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n"));
364: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
365: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
366: PetscCall(MatCreateVecs(A, &x, &y));
367: PetscCall(VecSetRandom(x, NULL));
368: if (issymmetric) {
369: PetscCall(VecCopy(x, y));
370: } else {
371: PetscCall(VecSetRandom(y, NULL));
372: PetscCall(VecScale(y, 8.));
373: }
374: PetscCall(MatDiagonalScale(A2, y, x));
375: PetscCall(MatDiagonalScale(B2, y, x));
376: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale"));
377: PetscCall(MatDestroy(&A2));
378: PetscCall(MatDestroy(&B2));
379: PetscCall(VecDestroy(&x));
380: PetscCall(VecDestroy(&y));
382: /* test MatPtAP (A IS and B AIJ) */
383: if (isaij && m == n) {
384: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n"));
385: /* There's a bug in MatCreateSubMatrices_MPIAIJ I cannot figure out */
386: if (!allow_repeated || !repmap || size == 1) {
387: PetscCall(MatISStoreL2L(A, PETSC_TRUE));
388: PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A2));
389: PetscCall(MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B2));
390: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX"));
391: PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &A2));
392: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX"));
393: PetscCall(MatDestroy(&A2));
394: PetscCall(MatDestroy(&B2));
395: }
396: }
398: /* test MatGetLocalSubMatrix */
399: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n"));
400: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2));
401: PetscCall(ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven));
402: PetscCall(ISComplement(reven, 0, lm, &rodd));
403: PetscCall(ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven));
404: PetscCall(ISComplement(ceven, 0, ln, &codd));
405: PetscCall(MatGetLocalSubMatrix(A2, reven, ceven, &Aee));
406: PetscCall(MatGetLocalSubMatrix(A2, reven, codd, &Aeo));
407: PetscCall(MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe));
408: PetscCall(MatGetLocalSubMatrix(A2, rodd, codd, &Aoo));
409: for (i = 0; i < lm; i++) {
410: PetscInt j, je, jo, colse[3], colso[3];
411: PetscScalar ve[3], vo[3];
412: PetscScalar v[3];
413: PetscInt cols[3];
414: PetscInt row = i / 2;
416: cols[0] = (i - 1 + n) % n;
417: cols[1] = i % n;
418: cols[2] = (i + 1) % n;
419: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
420: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
421: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
422: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
423: for (j = 0, je = 0, jo = 0; j < 3; j++) {
424: if (cols[j] % 2) {
425: vo[jo] = v[j];
426: colso[jo++] = cols[j] / 2;
427: } else {
428: ve[je] = v[j];
429: colse[je++] = cols[j] / 2;
430: }
431: }
432: if (i % 2) {
433: PetscCall(MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES));
434: PetscCall(MatSetValuesBlockedLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES));
435: } else {
436: PetscCall(MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES));
437: PetscCall(MatSetValuesBlockedLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES));
438: }
439: }
440: PetscCall(MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee));
441: PetscCall(MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo));
442: PetscCall(MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe));
443: PetscCall(MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo));
444: PetscCall(ISDestroy(&reven));
445: PetscCall(ISDestroy(&ceven));
446: PetscCall(ISDestroy(&rodd));
447: PetscCall(ISDestroy(&codd));
448: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
449: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
450: PetscCall(MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN));
451: PetscCall(CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix"));
452: PetscCall(MatDestroy(&A2));
454: /* test MatConvert_Nest_IS */
455: testT = PETSC_FALSE;
456: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL));
458: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n"));
459: nr = 2;
460: nc = 2;
461: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL));
462: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
463: if (testT) {
464: PetscCall(MatGetOwnershipRange(A, &cst, &cen));
465: PetscCall(MatGetOwnershipRangeColumn(A, &rst, &ren));
466: } else {
467: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
468: PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
469: }
470: PetscCall(PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats));
471: for (i = 0; i < nr * nc; i++) {
472: if (testT) {
473: PetscCall(MatCreateTranspose(A, &mats[i]));
474: PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]));
475: } else {
476: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &mats[i]));
477: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]));
478: }
479: }
480: for (i = 0; i < nr; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]));
481: for (i = 0; i < nc; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]));
482: PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2));
483: PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2));
484: for (i = 0; i < nr; i++) PetscCall(ISDestroy(&rows[i]));
485: for (i = 0; i < nc; i++) PetscCall(ISDestroy(&cols[i]));
486: for (i = 0; i < 2 * nr * nc; i++) PetscCall(MatDestroy(&mats[i]));
487: PetscCall(PetscFree3(rows, cols, mats));
488: PetscCall(MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T));
489: PetscCall(MatDestroy(&B2));
490: PetscCall(MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2));
491: PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX"));
492: PetscCall(MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2));
493: PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX"));
494: PetscCall(MatDestroy(&B2));
495: PetscCall(MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2));
496: PetscCall(CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX"));
497: PetscCall(MatDestroy(&T));
498: PetscCall(MatDestroy(&A2));
500: /* test MatCreateSubMatrix */
501: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n"));
502: if (rank == 0) {
503: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is));
504: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2));
505: } else if (rank == 1) {
506: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
507: if (n > 3) {
508: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2));
509: } else {
510: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
511: }
512: } else if (rank == 2 && n > 4) {
513: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
514: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2));
515: } else {
516: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
517: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
518: }
519: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2));
520: PetscCall(MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2));
521: PetscCall(CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix"));
523: PetscCall(MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2));
524: PetscCall(MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2));
525: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix"));
526: PetscCall(MatDestroy(&A2));
527: PetscCall(MatDestroy(&B2));
529: if (!issymmetric) {
530: PetscCall(MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2));
531: PetscCall(MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2));
532: PetscCall(MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2));
533: PetscCall(MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2));
534: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix"));
535: }
537: PetscCall(MatDestroy(&A2));
538: PetscCall(MatDestroy(&B2));
539: PetscCall(ISDestroy(&is));
540: PetscCall(ISDestroy(&is2));
542: /* test MatCreateSubMatrices */
543: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n"));
544: PetscCall(MatGetLayouts(A, &rlayout, &clayout));
545: PetscCall(PetscLayoutGetRanges(rlayout, &rrange));
546: PetscCall(PetscLayoutGetRanges(clayout, &crange));
547: lrank = (size + rank - 1) % size;
548: rrank = (rank + 1) % size;
549: PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[lrank + 1] - rrange[lrank], rrange[lrank], 1, &irow[0]));
550: PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[rrank + 1] - crange[rrank], crange[rrank], 1, &icol[0]));
551: PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[rrank + 1] - rrange[rrank], rrange[rrank], 1, &irow[1]));
552: PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[lrank + 1] - crange[lrank], crange[lrank], 1, &icol[1]));
553: PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub));
554: PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub));
555: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
556: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
557: PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub));
558: PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub));
559: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
560: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
561: PetscCall(MatDestroySubMatrices(2, &Asub));
562: PetscCall(MatDestroySubMatrices(2, &Bsub));
563: PetscCall(ISDestroy(&irow[0]));
564: PetscCall(ISDestroy(&irow[1]));
565: PetscCall(ISDestroy(&icol[0]));
566: PetscCall(ISDestroy(&icol[1]));
568: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
569: if (size > 1) {
570: if (rank == 0) {
571: PetscInt st, len;
573: st = (m + 1) / 2;
574: len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
575: PetscCall(ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is));
576: } else {
577: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
578: }
579: } else {
580: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
581: }
582: /* local IS for local zero operations */
583: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
584: PetscCall(ISCreateStride(PETSC_COMM_WORLD, lm ? 1 : 0, 0, 1, &lis));
586: if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
587: PetscInt *idx0, *idx1, n0, n1;
588: IS Ais[2], Bis[2];
590: /* test MatDiagonalSet */
591: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n"));
592: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
593: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
594: PetscCall(MatCreateVecs(A, NULL, &x));
595: PetscCall(VecSetRandom(x, NULL));
596: PetscCall(MatDiagonalSet(A2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
597: PetscCall(MatDiagonalSet(B2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
598: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet"));
599: PetscCall(VecDestroy(&x));
600: PetscCall(MatDestroy(&A2));
601: PetscCall(MatDestroy(&B2));
603: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
604: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n"));
605: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
606: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
607: PetscCall(MatShift(A2, 2.0));
608: PetscCall(MatShift(B2, 2.0));
609: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatShift"));
610: PetscCall(MatDestroy(&A2));
611: PetscCall(MatDestroy(&B2));
613: /* nonzero diag value is supported for square matrices only */
614: PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, is, diag, PETSC_FALSE));
615: PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, lis, diag, PETSC_TRUE));
617: /* test MatIncreaseOverlap */
618: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n"));
619: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
620: n0 = (ren - rst) / 2;
621: n1 = (ren - rst) / 3;
622: PetscCall(PetscMalloc1(n0, &idx0));
623: PetscCall(PetscMalloc1(n1, &idx1));
624: for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
625: for (i = 0; i < n1; i++) idx1[i] = rst + i;
626: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]));
627: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]));
628: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]));
629: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]));
630: PetscCall(MatIncreaseOverlap(A, 2, Ais, 3));
631: PetscCall(MatIncreaseOverlap(B, 2, Bis, 3));
632: /* Non deterministic output! */
633: PetscCall(ISSort(Ais[0]));
634: PetscCall(ISSort(Ais[1]));
635: PetscCall(ISSort(Bis[0]));
636: PetscCall(ISSort(Bis[1]));
637: PetscCall(ISView(Ais[0], NULL));
638: PetscCall(ISView(Bis[0], NULL));
639: PetscCall(ISView(Ais[1], NULL));
640: PetscCall(ISView(Bis[1], NULL));
641: PetscCall(MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub));
642: PetscCall(MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub));
643: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]"));
644: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]"));
645: PetscCall(MatDestroySubMatrices(2, &Asub));
646: PetscCall(MatDestroySubMatrices(2, &Bsub));
647: PetscCall(ISDestroy(&Ais[0]));
648: PetscCall(ISDestroy(&Ais[1]));
649: PetscCall(ISDestroy(&Bis[0]));
650: PetscCall(ISDestroy(&Bis[1]));
651: }
652: PetscCall(TestMatZeroRows(A, B, squaretest, is, 0.0, PETSC_FALSE));
653: PetscCall(TestMatZeroRows(A, B, squaretest, lis, 0.0, PETSC_TRUE));
654: PetscCall(ISDestroy(&is));
655: PetscCall(ISDestroy(&lis));
657: /* test MatTranspose */
658: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n"));
659: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
660: PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &B2));
661: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose"));
663: PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &A2));
664: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose"));
665: PetscCall(MatDestroy(&A2));
667: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
668: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
669: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose"));
670: PetscCall(MatDestroy(&A2));
672: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
673: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose"));
674: PetscCall(MatDestroy(&A2));
675: PetscCall(MatDestroy(&B2));
677: /* test MatISFixLocalEmpty */
678: if (isaij) {
679: PetscInt r[2];
681: r[0] = 0;
682: r[1] = PetscMin(m, n) - 1;
683: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n"));
684: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
686: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
687: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
688: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
689: PetscCall(CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)"));
691: PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
692: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
693: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
694: PetscCall(MatZeroRows(B2, 2, r, 0.0, NULL, NULL));
695: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
696: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
697: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
698: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
699: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)"));
700: PetscCall(MatDestroy(&A2));
702: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
703: PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
704: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
705: PetscCall(MatTranspose(B2, MAT_INPLACE_MATRIX, &B2));
706: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
707: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
708: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
709: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
710: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
711: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)"));
713: PetscCall(MatDestroy(&A2));
714: PetscCall(MatDestroy(&B2));
716: if (squaretest) {
717: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
718: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
719: PetscCall(MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL));
720: PetscCall(MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL));
721: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
722: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
723: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
724: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
725: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
726: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)"));
727: PetscCall(MatDestroy(&A2));
728: PetscCall(MatDestroy(&B2));
729: }
730: }
732: /* test MatInvertBlockDiagonal
733: special cases for block-diagonal matrices */
734: if (m == n) {
735: ISLocalToGlobalMapping map;
736: Mat Abd, Bbd;
737: IS is, bis;
738: const PetscScalar *isbd, *aijbd;
739: const PetscInt *sts, *idxs;
740: PetscInt *idxs2, diff, perm, nl, bs, st, en, in;
741: PetscBool ok;
743: for (diff = 0; diff < 3; diff++) {
744: for (perm = 0; perm < 3; perm++) {
745: for (bs = 1; bs < 4; bs++) {
746: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs));
747: PetscCall(PetscMalloc1(bs * bs, &vals));
748: PetscCall(MatGetOwnershipRanges(A, &sts));
749: switch (diff) {
750: case 1: /* inverted layout by processes */
751: in = 1;
752: st = sts[size - rank - 1];
753: en = sts[size - rank];
754: nl = en - st;
755: break;
756: case 2: /* round-robin layout */
757: in = size;
758: st = rank;
759: nl = n / size;
760: if (rank < n % size) nl++;
761: break;
762: default: /* same layout */
763: in = 1;
764: st = sts[rank];
765: en = sts[rank + 1];
766: nl = en - st;
767: break;
768: }
769: PetscCall(ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is));
770: PetscCall(ISGetLocalSize(is, &nl));
771: PetscCall(ISGetIndices(is, &idxs));
772: PetscCall(PetscMalloc1(nl, &idxs2));
773: for (i = 0; i < nl; i++) {
774: switch (perm) { /* invert some of the indices */
775: case 2:
776: idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
777: break;
778: case 1:
779: idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
780: break;
781: default:
782: idxs2[i] = idxs[i];
783: break;
784: }
785: }
786: PetscCall(ISRestoreIndices(is, &idxs));
787: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis));
788: PetscCall(ISLocalToGlobalMappingCreateIS(bis, &map));
789: PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd));
790: PetscCall(ISLocalToGlobalMappingDestroy(&map));
791: PetscCall(MatISSetPreallocation(Abd, bs, NULL, 0, NULL));
792: for (i = 0; i < nl; i++) {
793: PetscInt b1, b2;
795: for (b1 = 0; b1 < bs; b1++) {
796: for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
797: }
798: PetscCall(MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES));
799: }
800: PetscCall(MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY));
801: PetscCall(MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY));
802: PetscCall(MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd));
803: PetscCall(MatInvertBlockDiagonal(Abd, &isbd));
804: PetscCall(MatInvertBlockDiagonal(Bbd, &aijbd));
805: PetscCall(MatGetLocalSize(Bbd, &nl, NULL));
806: ok = PETSC_TRUE;
807: for (i = 0; i < nl / bs; i++) {
808: PetscInt b1, b2;
810: for (b1 = 0; b1 < bs; b1++) {
811: for (b2 = 0; b2 < bs; b2++) {
812: if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
813: if (!ok) {
814: 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])));
815: break;
816: }
817: }
818: if (!ok) break;
819: }
820: if (!ok) break;
821: }
822: PetscCall(MatDestroy(&Abd));
823: PetscCall(MatDestroy(&Bbd));
824: PetscCall(PetscFree(vals));
825: PetscCall(ISDestroy(&is));
826: PetscCall(ISDestroy(&bis));
827: }
828: }
829: }
830: }
832: /* test MatGetDiagonalBlock */
833: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n"));
834: PetscCall(MatGetDiagonalBlock(A, &A2));
835: PetscCall(MatGetDiagonalBlock(B, &B2));
836: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
837: PetscCall(MatScale(A, 2.0));
838: PetscCall(MatScale(B, 2.0));
839: PetscCall(MatGetDiagonalBlock(A, &A2));
840: PetscCall(MatGetDiagonalBlock(B, &B2));
841: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
843: /* test MatISSetAllowRepeated on a MATIS */
844: PetscCall(MatISSetAllowRepeated(A, allow_repeated));
845: if (allow_repeated) { /* original MATIS maybe with repeated entries, test assembling of local matrices */
846: Mat lA, lA2;
848: for (PetscInt i = 0; i < 1; i++) { /* TODO: make MatScatter inherit from MATSHELL and support MatProducts */
849: PetscBool usemult = PETSC_FALSE;
851: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
852: if (i) {
853: Mat tA;
855: usemult = PETSC_TRUE;
856: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries and local shell matrices\n"));
857: PetscCall(MatISGetLocalMat(A2, &lA2));
858: PetscCall(MatConvert(lA2, MATSHELL, MAT_INITIAL_MATRIX, &tA));
859: PetscCall(MatISRestoreLocalMat(A2, &lA2));
860: PetscCall(MatISSetLocalMat(A2, tA));
861: PetscCall(MatDestroy(&tA));
862: } else {
863: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries\n"));
864: }
865: PetscCall(MatISSetAllowRepeated(A2, PETSC_FALSE));
866: PetscCall(MatISGetLocalMat(A, &lA));
867: PetscCall(MatISGetLocalMat(A2, &lA2));
868: if (!repmap) PetscCall(CheckMat(lA, lA2, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
869: PetscCall(MatISRestoreLocalMat(A, &lA));
870: PetscCall(MatISRestoreLocalMat(A2, &lA2));
871: if (repmap) PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with repeated entries"));
872: else PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
873: PetscCall(MatDestroy(&A2));
874: }
875: } else { /* original matis with non-repeated entries, this should only recreate the local matrices */
876: Mat lA;
877: PetscBool flg;
879: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(true) with non repeated entries\n"));
880: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
881: PetscCall(MatISSetAllowRepeated(A2, PETSC_TRUE));
882: PetscCall(MatISGetLocalMat(A2, &lA));
883: PetscCall(MatAssembled(lA, &flg));
884: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local mat should be unassembled");
885: PetscCall(MatISRestoreLocalMat(A2, &lA));
886: PetscCall(MatDestroy(&A2));
887: }
889: /* Test MatZeroEntries */
890: PetscCall(MatZeroEntries(A));
891: PetscCall(MatZeroEntries(B));
892: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatZeroEntries"));
894: /* Test MatSetValues and MatSetValuesBlocked */
895: if (test_setvalues) {
896: PetscCall(PetscMalloc1(lm * ln, &vals));
897: for (i = 0; i < lm * ln; i++) vals[i] = i + 1.0;
898: PetscCall(MatGetLocalSize(A, NULL, &ln));
899: PetscCall(MatISSetPreallocation(A, ln, NULL, n - ln, NULL));
900: PetscCall(MatSeqAIJSetPreallocation(B, ln, NULL));
901: PetscCall(MatMPIAIJSetPreallocation(B, ln, NULL, n - ln, NULL));
902: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
903: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
905: PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
906: PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
907: PetscCall(MatSetValues(A, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
908: PetscCall(MatSetValues(B, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
909: PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
910: PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
911: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
912: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
913: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
914: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
915: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValues"));
917: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rmap, &idxs1));
918: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cmap, &idxs2));
919: PetscCall(MatSetValuesBlocked(A, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
920: PetscCall(MatSetValuesBlocked(B, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
921: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rmap, &idxs1));
922: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cmap, &idxs2));
923: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
924: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
925: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
926: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
927: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValuesBlocked"));
928: PetscCall(PetscFree(vals));
929: }
931: /* free testing matrices */
932: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
933: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
934: PetscCall(MatDestroy(&A));
935: PetscCall(MatDestroy(&B));
936: PetscCall(PetscFinalize());
937: return 0;
938: }
940: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
941: {
942: Mat Bcheck;
943: PetscReal error;
945: PetscFunctionBeginUser;
946: if (!usemult && B) {
947: PetscBool hasnorm;
949: PetscCall(MatHasOperation(B, MATOP_NORM, &hasnorm));
950: if (!hasnorm) usemult = PETSC_TRUE;
951: }
952: if (!usemult) {
953: if (B) {
954: MatType Btype;
956: PetscCall(MatGetType(B, &Btype));
957: PetscCall(MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck));
958: } else {
959: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
960: }
961: if (B) { /* if B is present, subtract it */
962: PetscCall(MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN));
963: }
964: PetscCall(MatNorm(Bcheck, NORM_INFINITY, &error));
965: if (error > PETSC_SQRT_MACHINE_EPSILON) {
966: ISLocalToGlobalMapping rl2g, cl2g;
968: PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Bcheck"));
969: PetscCall(MatView(Bcheck, NULL));
970: if (B) {
971: PetscCall(PetscObjectSetName((PetscObject)B, "B"));
972: PetscCall(MatView(B, NULL));
973: PetscCall(MatDestroy(&Bcheck));
974: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
975: PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Assembled A"));
976: PetscCall(MatView(Bcheck, NULL));
977: }
978: PetscCall(MatDestroy(&Bcheck));
979: PetscCall(PetscObjectSetName((PetscObject)A, "A"));
980: PetscCall(MatView(A, NULL));
981: PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
982: if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g, NULL));
983: if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g, NULL));
984: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
985: }
986: PetscCall(MatDestroy(&Bcheck));
987: } else {
988: PetscBool ok, okt;
990: PetscCall(MatMultEqual(A, B, 3, &ok));
991: PetscCall(MatMultTransposeEqual(A, B, 3, &okt));
992: PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ? %d, multtranspose ok ? %d", func, ok, okt);
993: }
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag, PetscBool local)
998: {
999: Mat B, Bcheck, B2 = NULL, lB;
1000: Vec x = NULL, b = NULL, b2 = NULL;
1001: ISLocalToGlobalMapping l2gr, l2gc;
1002: PetscReal error;
1003: char diagstr[16];
1004: const PetscInt *idxs;
1005: PetscInt rst, ren, i, n, N, d;
1006: PetscMPIInt rank;
1007: PetscBool miss, haszerorows;
1008: IS gis;
1010: PetscFunctionBeginUser;
1011: if (diag == 0.) {
1012: PetscCall(PetscStrncpy(diagstr, "zero", sizeof(diagstr)));
1013: } else {
1014: PetscCall(PetscStrncpy(diagstr, "nonzero", sizeof(diagstr)));
1015: }
1016: PetscCall(ISView(is, NULL));
1017: PetscCall(MatGetLocalToGlobalMapping(A, &l2gr, &l2gc));
1018: /* tests MatDuplicate and MatCopy */
1019: if (diag == 0.) {
1020: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
1021: } else {
1022: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
1023: PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
1024: }
1025: PetscCall(MatISGetLocalMat(B, &lB));
1026: PetscCall(MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows));
1027: if (squaretest && haszerorows) {
1028: PetscCall(MatCreateVecs(B, &x, &b));
1029: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
1030: PetscCall(VecSetLocalToGlobalMapping(b, l2gr));
1031: PetscCall(VecSetLocalToGlobalMapping(x, l2gc));
1032: PetscCall(VecSetRandom(x, NULL));
1033: PetscCall(VecSetRandom(b, NULL));
1034: /* mimic b[is] = x[is] */
1035: PetscCall(VecDuplicate(b, &b2));
1036: PetscCall(VecSetLocalToGlobalMapping(b2, l2gr));
1037: PetscCall(VecCopy(b, b2));
1038: if (local) {
1039: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1040: PetscCall(ISGetLocalSize(gis, &n));
1041: PetscCall(ISGetIndices(gis, &idxs));
1042: } else {
1043: PetscCall(ISGetLocalSize(is, &n));
1044: PetscCall(ISGetIndices(is, &idxs));
1045: }
1046: PetscCall(VecGetSize(x, &N));
1047: for (i = 0; i < n; i++) {
1048: if (0 <= idxs[i] && idxs[i] < N) {
1049: PetscCall(VecSetValue(b2, idxs[i], diag, INSERT_VALUES));
1050: PetscCall(VecSetValue(x, idxs[i], 1., INSERT_VALUES));
1051: }
1052: }
1053: if (local) {
1054: PetscCall(ISRestoreIndices(gis, &idxs));
1055: PetscCall(ISDestroy(&gis));
1056: } else {
1057: PetscCall(ISRestoreIndices(is, &idxs));
1058: }
1059: PetscCall(VecAssemblyBegin(b2));
1060: PetscCall(VecAssemblyEnd(b2));
1061: PetscCall(VecAssemblyBegin(x));
1062: PetscCall(VecAssemblyEnd(x));
1063: /* test ZeroRows on MATIS */
1064: if (local) {
1065: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsLocal (diag %s)\n", diagstr));
1066: PetscCall(MatZeroRowsLocalIS(B, is, diag, x, b));
1067: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumnsLocal (diag %s)\n", diagstr));
1068: PetscCall(MatZeroRowsColumnsLocalIS(B2, is, diag, NULL, NULL));
1069: } else {
1070: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1071: PetscCall(MatZeroRowsIS(B, is, diag, x, b));
1072: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr));
1073: PetscCall(MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL));
1074: }
1075: } else if (haszerorows) {
1076: /* test ZeroRows on MATIS */
1077: if (local) {
1078: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsLocal (diag %s)\n", diagstr));
1079: PetscCall(MatZeroRowsLocalIS(B, is, diag, NULL, NULL));
1080: } else {
1081: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1082: PetscCall(MatZeroRowsIS(B, is, diag, NULL, NULL));
1083: }
1084: b = b2 = x = NULL;
1085: } else {
1086: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr));
1087: b = b2 = x = NULL;
1088: }
1090: if (squaretest && haszerorows) {
1091: PetscCall(VecAXPY(b2, -1., b));
1092: PetscCall(VecNorm(b2, NORM_INFINITY, &error));
1093: PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr);
1094: }
1096: /* test MatMissingDiagonal */
1097: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n"));
1098: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1099: PetscCall(MatMissingDiagonal(B, &miss, &d));
1100: PetscCall(MatGetOwnershipRange(B, &rst, &ren));
1101: PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
1102: 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));
1103: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1104: PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));
1106: PetscCall(VecDestroy(&x));
1107: PetscCall(VecDestroy(&b));
1108: PetscCall(VecDestroy(&b2));
1110: /* check the result of ZeroRows with that from MPIAIJ routines
1111: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
1112: if (haszerorows) {
1113: PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck));
1114: PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1115: if (local) {
1116: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1117: PetscCall(MatZeroRowsIS(Bcheck, gis, diag, NULL, NULL));
1118: PetscCall(ISDestroy(&gis));
1119: } else {
1120: PetscCall(MatZeroRowsIS(Bcheck, is, diag, NULL, NULL));
1121: }
1122: PetscCall(CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows"));
1123: PetscCall(MatDestroy(&Bcheck));
1124: }
1125: PetscCall(MatDestroy(&B));
1127: if (B2) { /* test MatZeroRowsColumns */
1128: PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &B));
1129: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1130: if (local) {
1131: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1132: PetscCall(MatZeroRowsColumnsIS(B, gis, diag, NULL, NULL));
1133: PetscCall(ISDestroy(&gis));
1134: } else {
1135: PetscCall(MatZeroRowsColumnsIS(B, is, diag, NULL, NULL));
1136: }
1137: PetscCall(CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns"));
1138: PetscCall(MatDestroy(&B));
1139: PetscCall(MatDestroy(&B2));
1140: }
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: PetscErrorCode ISL2GMapNoNeg(ISLocalToGlobalMapping mapping, IS is, IS *newis)
1145: {
1146: PetscInt n, *idxout, nn = 0;
1147: const PetscInt *idxin;
1149: PetscFunctionBegin;
1150: PetscCall(ISGetLocalSize(is, &n));
1151: PetscCall(ISGetIndices(is, &idxin));
1152: PetscCall(PetscMalloc1(n, &idxout));
1153: PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
1154: PetscCall(ISRestoreIndices(is, &idxin));
1155: for (PetscInt i = 0; i < n; i++)
1156: if (idxout[i] > -1) idxout[nn++] = idxout[i];
1157: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), nn, idxout, PETSC_OWN_POINTER, newis));
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: /*TEST
1163: test:
1164: requires: !complex
1165: args: -test_matlab -test_trans
1167: test:
1168: suffix: 2
1169: nsize: 4
1170: args: -mat_is_convert_local_nest -nr 3 -nc 4
1172: test:
1173: suffix: 3
1174: nsize: 5
1175: args: -m 11 -n 10 -mat_is_convert_local_nest -nr 2 -nc 1 -cbs 2
1177: test:
1178: suffix: 4
1179: nsize: 6
1180: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
1182: test:
1183: suffix: 5
1184: nsize: 6
1185: args: -m 12 -n 12 -test_trans -nr 3 -nc 1 -rbs 2
1187: test:
1188: suffix: 6
1189: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -rbs 6 -cbs 3
1191: test:
1192: suffix: 7
1193: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
1195: test:
1196: suffix: 8
1197: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
1199: test:
1200: suffix: 9
1201: nsize: 5
1202: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
1204: test:
1205: suffix: 10
1206: nsize: 5
1207: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
1209: test:
1210: suffix: vscat_default
1211: nsize: 5
1212: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
1213: output_file: output/ex23_11.out
1215: test:
1216: suffix: 12
1217: nsize: 3
1218: args: -m 12 -n 12 -symmetric -mat_is_localmat_type sbaij -test_trans -nr 2 -nc 3 -test_setvalues 0
1220: testset:
1221: output_file: output/ex23_13.out
1222: nsize: 3
1223: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
1224: filter: grep -v "type:" | grep -v "block size is 1" | grep -v "not using I-node routines"
1225: test:
1226: suffix: baij
1227: args: -mat_is_localmat_type baij
1228: test:
1229: requires: viennacl
1230: suffix: viennacl
1231: args: -mat_is_localmat_type aijviennacl
1232: test:
1233: requires: cuda
1234: suffix: cusparse
1235: args: -mat_is_localmat_type aijcusparse
1236: test:
1237: requires: kokkos_kernels
1238: suffix: kokkos
1239: args: -mat_is_localmat_type aijkokkos
1241: test:
1242: suffix: negrep
1243: nsize: {{1 3}separate output}
1244: 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} -allow_repeated 0
1246: test:
1247: suffix: negrep_allowrep
1248: nsize: {{1 3}separate output}
1249: 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} -allow_repeated
1251: TEST*/