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 TestMatZeroRowsShell(Mat, Mat, IS, PetscScalar);
7: PetscErrorCode TestMatZeroRowsMissingDiagonal(MPI_Comm);
8: PetscErrorCode CheckMat(Mat, Mat, PetscBool, const char *);
9: PetscErrorCode CheckVariableBlockSizes(Mat);
10: PetscErrorCode CheckRepeatedMapFiltering(MPI_Comm, PetscInt, PetscInt, InsertMode);
11: PetscErrorCode ISL2GMapNoNeg(ISLocalToGlobalMapping, IS, IS *);
13: int main(int argc, char **args)
14: {
15: Mat A, B, A2, B2, T;
16: Mat Aee, Aeo, Aoe, Aoo;
17: Mat *mats, *Asub, *Bsub;
18: Vec x, y;
19: MatInfo info;
20: ISLocalToGlobalMapping cmap, rmap;
21: IS is, lis, is2, reven, rodd, ceven, codd;
22: IS *rows, *cols;
23: IS irow[2], icol[2];
24: PetscLayout rlayout, clayout;
25: const PetscInt *rrange, *crange, *idxs1, *idxs2;
26: MatType lmtype;
27: PetscScalar diag = 2., *vals;
28: PetscInt n, m, i, lm, ln;
29: PetscInt rst, ren, cst, cen, nr, nc, rbs = 1, cbs = 1;
30: PetscMPIInt rank, size, lrank, rrank;
31: PetscBool testT, squaretest, isaij;
32: PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE, allow_repeated = PETSC_TRUE;
33: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric, test_matlab = PETSC_FALSE, test_setvalues = PETSC_TRUE;
34: PetscBool test_variableblocksizes = PETSC_FALSE;
36: PetscFunctionBeginUser;
37: PetscCall(PetscInitialize(&argc, &args, NULL, help));
38: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
39: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
40: m = n = 2 * size;
41: PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL));
42: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
43: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
44: PetscCall(PetscOptionsGetBool(NULL, NULL, "-negmap", &negmap, NULL));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-repmap", &repmap, NULL));
46: PetscCall(PetscOptionsGetBool(NULL, NULL, "-permmap", &permute, NULL));
47: PetscCall(PetscOptionsGetBool(NULL, NULL, "-diffmap", &diffmap, NULL));
48: PetscCall(PetscOptionsGetBool(NULL, NULL, "-allow_repeated", &allow_repeated, NULL));
49: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matlab", &test_matlab, NULL));
50: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_setvalues", &test_setvalues, NULL));
51: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_variableblocksizes", &test_variableblocksizes, NULL));
52: PetscCall(PetscOptionsGetInt(NULL, NULL, "-rbs", &rbs, NULL));
53: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cbs", &cbs, NULL));
54: PetscCheck(size == 1 || m >= 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 4 for parallel runs");
55: PetscCheck(size != 1 || m >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 2 for uniprocessor runs");
56: PetscCheck(n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of cols should be larger or equal 2");
58: /* create a MATIS matrix */
59: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
60: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
61: PetscCall(MatSetType(A, MATIS));
62: PetscCall(MatSetFromOptions(A));
63: if (!negmap && !repmap) {
64: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
65: Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
66: Equivalent to passing NULL for the mapping */
67: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is));
68: } else if (negmap && !repmap) { /* non repeated but with negative indices */
69: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is));
70: } else if (!negmap && repmap) { /* non negative but repeated indices */
71: IS isl[2];
73: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]));
74: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]));
75: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
76: PetscCall(ISDestroy(&isl[0]));
77: PetscCall(ISDestroy(&isl[1]));
78: } else { /* negative and repeated indices */
79: IS isl[2];
81: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]));
82: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]));
83: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
84: PetscCall(ISDestroy(&isl[0]));
85: PetscCall(ISDestroy(&isl[1]));
86: }
87: PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmap));
88: PetscCall(ISDestroy(&is));
90: if (m != n || diffmap) {
91: PetscCall(ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is));
92: PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmap));
93: PetscCall(ISDestroy(&is));
94: } else {
95: PetscCall(PetscObjectReference((PetscObject)cmap));
96: rmap = cmap;
97: }
99: PetscCall(MatISSetAllowRepeated(A, allow_repeated));
100: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
101: PetscCall(MatSetBlockSizes(A, rbs, cbs));
102: PetscCall(MatISStoreL2L(A, PETSC_FALSE));
103: PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));
104: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
105: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
106: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
107: for (i = 0; i < lm; i++) {
108: PetscScalar v[3];
109: PetscInt cols[3];
111: cols[0] = (i - 1 + n) % n;
112: cols[1] = i % n;
113: cols[2] = (i + 1) % n;
114: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
115: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
116: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
117: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
118: PetscCall(MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES));
119: }
120: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
121: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
123: /* activate tests for square matrices with same maps only */
124: PetscCall(MatHasCongruentLayouts(A, &squaretest));
125: if (squaretest && rmap != cmap) {
126: PetscInt nr, nc;
128: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nr));
129: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nc));
130: if (nr != nc) squaretest = PETSC_FALSE;
131: else {
132: PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
133: PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
134: PetscCall(PetscArraycmp(idxs1, idxs2, nr, &squaretest));
135: PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
136: PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
137: }
138: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
139: }
140: if (negmap && repmap) squaretest = PETSC_FALSE;
141: PetscCheck(squaretest || !test_variableblocksizes, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Variable block size test only for square matrices");
143: /* test MatISGetLocalMat */
144: PetscCall(MatISGetLocalMat(A, &B));
145: PetscCall(MatGetType(B, &lmtype));
146: if (test_variableblocksizes) {
147: PetscInt bsizes[2] = {rbs, lm - rbs};
149: PetscCall(MatSetVariableBlockSizes(B, PETSC_STATIC_ARRAY_LENGTH(bsizes), bsizes));
150: }
152: /* test MatGetInfo */
153: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n"));
154: PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
155: PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
156: 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,
157: (PetscInt)info.nz_unneeded, (PetscInt)info.assemblies, (PetscInt)info.mallocs));
158: PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
159: PetscCall(MatGetInfo(A, MAT_GLOBAL_MAX, &info));
160: 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,
161: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
162: PetscCall(MatGetInfo(A, MAT_GLOBAL_SUM, &info));
163: 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,
164: (PetscInt)info.assemblies, (PetscInt)info.mallocs));
166: /* test MatIsSymmetric */
167: PetscCall(MatIsSymmetric(A, 0.0, &issymmetric));
168: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric));
170: /* Create a MPIAIJ matrix, same as A */
171: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
172: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
173: PetscCall(MatSetBlockSizes(B, rbs, cbs));
174: PetscCall(MatSetType(B, MATAIJ));
175: PetscCall(MatSetFromOptions(B));
176: PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
177: PetscCall(MatMPIAIJSetPreallocation(B, 3, NULL, 3, NULL));
178: PetscCall(MatMPIBAIJSetPreallocation(B, 1, 3, NULL, 3, NULL));
179: #if PetscDefined(HAVE_HYPRE)
180: PetscCall(MatHYPRESetPreallocation(B, 3, NULL, 3, NULL));
181: #endif
182: PetscCall(MatISSetPreallocation(B, 3, NULL, 3, NULL));
183: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
184: for (i = 0; i < lm; i++) {
185: PetscScalar v[3];
186: PetscInt cols[3];
188: cols[0] = (i - 1 + n) % n;
189: cols[1] = i % n;
190: cols[2] = (i + 1) % n;
191: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
192: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
193: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
194: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
195: PetscCall(MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES));
196: }
197: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
198: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
200: /* test MatView */
201: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n"));
202: PetscCall(MatView(A, NULL));
203: PetscCall(MatView(B, NULL));
205: /* test MATLAB ASCII view */
206: if (test_matlab) { /* output is different when using real or complex numbers */
207: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView ASCII MATLAB\n"));
208: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB));
209: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
210: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
211: }
213: /* test CheckMat */
214: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n"));
215: PetscCall(CheckMat(A, B, PETSC_FALSE, "CheckMat"));
217: /* test binary MatView/MatLoad */
218: {
219: PetscMPIInt color = rank % 2;
220: MPI_Comm comm;
221: char name[PETSC_MAX_PATH_LEN];
222: PetscViewer wview, cview, sview, view;
223: Mat A2;
225: PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, rank, &comm));
227: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "world_is"));
228: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_WRITE, &wview));
229: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "seq_is_%d", rank));
230: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_WRITE, &sview));
231: PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "color_is_%d", color));
232: PetscCall(PetscViewerBinaryOpen(comm, name, FILE_MODE_WRITE, &cview));
233: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary world\n"));
234: PetscCall(MatView(A, wview));
235: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary self\n"));
236: PetscCall(MatView(A, sview));
237: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary subcomm\n"));
238: PetscCall(MatView(A, cview));
239: PetscCall(PetscViewerDestroy(&wview));
240: PetscCall(PetscViewerDestroy(&cview));
241: PetscCall(PetscViewerDestroy(&sview));
243: /* Load a world matrix */
244: PetscCall(MatCreate(PETSC_COMM_WORLD, &A2));
245: PetscCall(MatSetType(A2, MATIS));
246: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
248: /* Read back the same matrix and check */
249: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from world\n"));
250: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "world_is", FILE_MODE_READ, &view));
251: PetscCall(MatLoad(A2, view));
252: if (test_variableblocksizes) PetscCall(CheckVariableBlockSizes(A2));
253: PetscCall(CheckMat(A, A2, PETSC_TRUE, "Load"));
254: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
255: PetscCall(PetscViewerDestroy(&view));
257: /* Read the matrix from rank 0 only */
258: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from self\n"));
259: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "seq_is_0", FILE_MODE_READ, &view));
260: PetscCall(MatLoad(A2, view));
261: if (test_variableblocksizes) PetscCall(CheckVariableBlockSizes(A2));
262: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
263: PetscCall(PetscViewerDestroy(&view));
265: /* Read the matrix from subcomm */
266: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from subcomm\n"));
267: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "color_is_0", FILE_MODE_READ, &view));
268: PetscCall(MatLoad(A2, view));
269: if (test_variableblocksizes) PetscCall(CheckVariableBlockSizes(A2));
270: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
271: PetscCall(PetscViewerDestroy(&view));
273: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
274: PetscCall(MatDestroy(&A2));
276: /* now load the original matrix from color 0 only processes */
277: if (!color) {
278: PetscCall(PetscPrintf(comm, "Test subcomm MatLoad from world\n"));
279: PetscCall(MatCreate(comm, &A2));
280: PetscCall(MatSetType(A2, MATIS));
281: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm), PETSC_VIEWER_ASCII_INFO_DETAIL));
282: PetscCall(PetscViewerBinaryOpen(comm, "world_is", FILE_MODE_READ, &view));
283: PetscCall(MatLoad(A2, view));
284: if (test_variableblocksizes) PetscCall(CheckVariableBlockSizes(A2));
285: PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_(comm)));
286: PetscCall(PetscViewerDestroy(&view));
287: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm)));
288: PetscCall(MatDestroy(&A2));
289: }
291: PetscCallMPI(MPI_Comm_free(&comm));
292: }
294: /* test MatDuplicate and MatAXPY */
295: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n"));
296: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
297: PetscCall(CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY"));
299: /* test MatConvert */
300: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n"));
301: PetscCall(MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2));
302: PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX"));
303: PetscCall(MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2));
304: PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX"));
305: PetscCall(MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2));
306: PetscCall(CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX"));
307: PetscCall(MatDestroy(&A2));
308: PetscCall(MatDestroy(&B2));
309: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n"));
310: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
311: PetscCall(MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2));
312: PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX"));
313: PetscCall(MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2));
314: PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX"));
315: PetscCall(MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2));
316: PetscCall(CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX"));
317: PetscCall(MatDestroy(&A2));
318: PetscCall(MatDestroy(&B2));
319: PetscCall(PetscStrcmp(lmtype, MATSEQAIJ, &isaij));
320: if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
321: 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};
322: ISLocalToGlobalMapping tcmap, trmap;
324: for (ri = 0; ri < 2; ri++) {
325: PetscInt *r;
327: r = (PetscInt *)(ri == 0 ? rr : rk);
328: for (ci = 0; ci < 2; ci++) {
329: PetscInt *c, rb, cb;
331: c = (PetscInt *)(ci == 0 ? cr : ck);
332: for (rb = 1; rb < 4; rb++) {
333: PetscCall(ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is));
334: PetscCall(ISLocalToGlobalMappingCreateIS(is, &trmap));
335: PetscCall(ISDestroy(&is));
336: for (cb = 1; cb < 4; cb++) {
337: Mat T, lT, T2;
338: char testname[256];
340: PetscCall(PetscSNPrintf(testname, sizeof(testname), "MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")", ri, ci, rb, cb));
341: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test %s\n", testname));
343: PetscCall(ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is));
344: PetscCall(ISLocalToGlobalMappingCreateIS(is, &tcmap));
345: PetscCall(ISDestroy(&is));
347: PetscCall(MatCreate(PETSC_COMM_SELF, &T));
348: PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4));
349: PetscCall(MatSetType(T, MATIS));
350: PetscCall(MatSetLocalToGlobalMapping(T, trmap, tcmap));
351: PetscCall(ISLocalToGlobalMappingDestroy(&tcmap));
352: PetscCall(MatISGetLocalMat(T, &lT));
353: PetscCall(MatSetType(lT, MATSEQAIJ));
354: PetscCall(MatSeqAIJSetPreallocation(lT, cb * 4, NULL));
355: PetscCall(MatSetRandom(lT, NULL));
356: PetscCall(MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT));
357: PetscCall(MatISRestoreLocalMat(T, &lT));
358: PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
359: PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
361: PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2));
362: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX"));
363: PetscCall(MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2));
364: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX"));
365: PetscCall(MatDestroy(&T2));
366: PetscCall(MatDuplicate(T, MAT_COPY_VALUES, &T2));
367: PetscCall(MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2));
368: PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX"));
369: PetscCall(MatDestroy(&T));
370: PetscCall(MatDestroy(&T2));
371: }
372: PetscCall(ISLocalToGlobalMappingDestroy(&trmap));
373: }
374: }
375: }
376: }
378: /* test MatDiagonalScale */
379: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n"));
380: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
381: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
382: PetscCall(MatCreateVecs(A, &x, &y));
383: PetscCall(VecSetRandom(x, NULL));
384: if (issymmetric) {
385: PetscCall(VecCopy(x, y));
386: } else {
387: PetscCall(VecSetRandom(y, NULL));
388: PetscCall(VecScale(y, 8.));
389: }
390: PetscCall(MatDiagonalScale(A2, y, x));
391: PetscCall(MatDiagonalScale(B2, y, x));
392: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale"));
393: PetscCall(MatDestroy(&A2));
394: PetscCall(MatDestroy(&B2));
395: PetscCall(VecDestroy(&x));
396: PetscCall(VecDestroy(&y));
398: /* test MatPtAP (A IS and B AIJ) */
399: if (isaij && m == n) {
400: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n"));
401: /* There's a bug in MatCreateSubMatrices_MPIAIJ I cannot figure out */
402: if (!allow_repeated || !repmap || size == 1) {
403: PetscCall(MatISStoreL2L(A, PETSC_TRUE));
404: PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A2));
405: PetscCall(MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B2));
406: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX"));
407: PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &A2));
408: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX"));
409: PetscCall(MatDestroy(&A2));
410: PetscCall(MatDestroy(&B2));
411: }
412: }
414: /* test MatGetLocalSubMatrix */
415: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n"));
416: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2));
417: PetscCall(ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven));
418: PetscCall(ISComplement(reven, 0, lm, &rodd));
419: PetscCall(ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven));
420: PetscCall(ISComplement(ceven, 0, ln, &codd));
421: PetscCall(MatGetLocalSubMatrix(A2, reven, ceven, &Aee));
422: PetscCall(MatGetLocalSubMatrix(A2, reven, codd, &Aeo));
423: PetscCall(MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe));
424: PetscCall(MatGetLocalSubMatrix(A2, rodd, codd, &Aoo));
425: for (i = 0; i < lm; i++) {
426: PetscInt j, je, jo, colse[3], colso[3];
427: PetscScalar ve[3], vo[3];
428: PetscScalar v[3];
429: PetscInt cols[3];
430: PetscInt row = i / 2;
432: cols[0] = (i - 1 + n) % n;
433: cols[1] = i % n;
434: cols[2] = (i + 1) % n;
435: v[0] = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
436: v[1] = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
437: v[2] = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
438: PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
439: for (j = 0, je = 0, jo = 0; j < 3; j++) {
440: if (cols[j] % 2) {
441: vo[jo] = v[j];
442: colso[jo++] = cols[j] / 2;
443: } else {
444: ve[je] = v[j];
445: colse[je++] = cols[j] / 2;
446: }
447: }
448: if (i % 2) {
449: PetscCall(MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES));
450: PetscCall(MatSetValuesBlockedLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES));
451: } else {
452: PetscCall(MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES));
453: PetscCall(MatSetValuesBlockedLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES));
454: }
455: }
456: PetscCall(MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee));
457: PetscCall(MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo));
458: PetscCall(MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe));
459: PetscCall(MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo));
460: PetscCall(ISDestroy(&reven));
461: PetscCall(ISDestroy(&ceven));
462: PetscCall(ISDestroy(&rodd));
463: PetscCall(ISDestroy(&codd));
464: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
465: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
466: PetscCall(MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN));
467: PetscCall(CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix"));
468: PetscCall(MatDestroy(&A2));
470: /* test MatConvert_Nest_IS */
471: testT = PETSC_FALSE;
472: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL));
474: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n"));
475: nr = 2;
476: nc = 2;
477: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL));
478: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
479: if (testT) {
480: PetscCall(MatGetOwnershipRange(A, &cst, &cen));
481: PetscCall(MatGetOwnershipRangeColumn(A, &rst, &ren));
482: } else {
483: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
484: PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
485: }
486: PetscCall(PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats));
487: for (i = 0; i < nr * nc; i++) {
488: if (testT) {
489: PetscCall(MatCreateTranspose(A, &mats[i]));
490: PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]));
491: } else {
492: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &mats[i]));
493: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]));
494: }
495: }
496: for (i = 0; i < nr; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]));
497: for (i = 0; i < nc; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]));
498: PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2));
499: PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2));
500: for (i = 0; i < nr; i++) PetscCall(ISDestroy(&rows[i]));
501: for (i = 0; i < nc; i++) PetscCall(ISDestroy(&cols[i]));
502: for (i = 0; i < 2 * nr * nc; i++) PetscCall(MatDestroy(&mats[i]));
503: PetscCall(PetscFree3(rows, cols, mats));
504: PetscCall(MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T));
505: PetscCall(MatDestroy(&B2));
506: PetscCall(MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2));
507: PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX"));
508: PetscCall(MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2));
509: PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX"));
510: PetscCall(MatDestroy(&B2));
511: PetscCall(MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2));
512: PetscCall(CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX"));
513: PetscCall(MatDestroy(&T));
514: PetscCall(MatDestroy(&A2));
516: /* test MatCreateSubMatrix */
517: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n"));
518: if (rank == 0) {
519: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is));
520: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2));
521: } else if (rank == 1) {
522: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
523: if (n > 3) {
524: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2));
525: } else {
526: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
527: }
528: } else if (rank == 2 && n > 4) {
529: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
530: PetscCall(ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2));
531: } else {
532: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
533: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
534: }
535: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2));
536: PetscCall(MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2));
537: PetscCall(CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix"));
539: PetscCall(MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2));
540: PetscCall(MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2));
541: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix"));
542: PetscCall(MatDestroy(&A2));
543: PetscCall(MatDestroy(&B2));
545: if (!issymmetric) {
546: PetscCall(MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2));
547: PetscCall(MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2));
548: PetscCall(MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2));
549: PetscCall(MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2));
550: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix"));
551: }
553: PetscCall(MatDestroy(&A2));
554: PetscCall(MatDestroy(&B2));
555: PetscCall(ISDestroy(&is));
556: PetscCall(ISDestroy(&is2));
558: /* test MatCreateSubMatrices */
559: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n"));
560: PetscCall(MatGetLayouts(A, &rlayout, &clayout));
561: PetscCall(PetscLayoutGetRanges(rlayout, &rrange));
562: PetscCall(PetscLayoutGetRanges(clayout, &crange));
563: lrank = (size + rank - 1) % size;
564: rrank = (rank + 1) % size;
565: PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[lrank + 1] - rrange[lrank], rrange[lrank], 1, &irow[0]));
566: PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[rrank + 1] - crange[rrank], crange[rrank], 1, &icol[0]));
567: PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[rrank + 1] - rrange[rrank], rrange[rrank], 1, &irow[1]));
568: PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[lrank + 1] - crange[lrank], crange[lrank], 1, &icol[1]));
569: PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub));
570: PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub));
571: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
572: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
573: PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub));
574: PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub));
575: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
576: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
577: PetscCall(MatDestroySubMatrices(2, &Asub));
578: PetscCall(MatDestroySubMatrices(2, &Bsub));
579: PetscCall(ISDestroy(&irow[0]));
580: PetscCall(ISDestroy(&irow[1]));
581: PetscCall(ISDestroy(&icol[0]));
582: PetscCall(ISDestroy(&icol[1]));
584: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
585: if (size > 1) {
586: if (rank == 0) {
587: PetscInt st, len;
589: st = (m + 1) / 2;
590: len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
591: PetscCall(ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is));
592: } else {
593: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
594: }
595: } else {
596: PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
597: }
598: /* local IS for local zero operations */
599: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
600: PetscCall(ISCreateStride(PETSC_COMM_WORLD, lm ? 1 : 0, 0, 1, &lis));
602: if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
603: PetscInt *idx0, *idx1, n0, n1;
604: IS Ais[2], Bis[2];
606: /* test MatDiagonalSet */
607: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n"));
608: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
609: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
610: PetscCall(MatCreateVecs(A, NULL, &x));
611: PetscCall(VecSetRandom(x, NULL));
612: PetscCall(MatDiagonalSet(A2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
613: PetscCall(MatDiagonalSet(B2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
614: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet"));
615: PetscCall(VecDestroy(&x));
616: PetscCall(MatDestroy(&A2));
617: PetscCall(MatDestroy(&B2));
619: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
620: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n"));
621: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
622: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
623: PetscCall(MatShift(A2, 2.0));
624: PetscCall(MatShift(B2, 2.0));
625: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatShift"));
626: PetscCall(MatDestroy(&A2));
627: PetscCall(MatDestroy(&B2));
629: /* nonzero diag value is supported for square matrices only */
630: PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, is, diag, PETSC_FALSE));
631: PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, lis, diag, PETSC_TRUE));
632: PetscCall(TestMatZeroRowsShell(A, B, is, diag));
633: PetscCall(TestMatZeroRowsMissingDiagonal(PETSC_COMM_WORLD));
635: /* test MatIncreaseOverlap */
636: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n"));
637: PetscCall(MatGetOwnershipRange(A, &rst, &ren));
638: n0 = (ren - rst) / 2;
639: n1 = (ren - rst) / 3;
640: PetscCall(PetscMalloc1(n0, &idx0));
641: PetscCall(PetscMalloc1(n1, &idx1));
642: for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
643: for (i = 0; i < n1; i++) idx1[i] = rst + i;
644: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]));
645: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]));
646: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]));
647: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]));
648: PetscCall(MatIncreaseOverlap(A, 2, Ais, 3));
649: PetscCall(MatIncreaseOverlap(B, 2, Bis, 3));
650: /* Non deterministic output! */
651: PetscCall(ISSort(Ais[0]));
652: PetscCall(ISSort(Ais[1]));
653: PetscCall(ISSort(Bis[0]));
654: PetscCall(ISSort(Bis[1]));
655: PetscCall(ISView(Ais[0], NULL));
656: PetscCall(ISView(Bis[0], NULL));
657: PetscCall(ISView(Ais[1], NULL));
658: PetscCall(ISView(Bis[1], NULL));
659: PetscCall(MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub));
660: PetscCall(MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub));
661: PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]"));
662: PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]"));
663: PetscCall(MatDestroySubMatrices(2, &Asub));
664: PetscCall(MatDestroySubMatrices(2, &Bsub));
665: PetscCall(ISDestroy(&Ais[0]));
666: PetscCall(ISDestroy(&Ais[1]));
667: PetscCall(ISDestroy(&Bis[0]));
668: PetscCall(ISDestroy(&Bis[1]));
669: }
670: PetscCall(TestMatZeroRows(A, B, squaretest, is, 0.0, PETSC_FALSE));
671: PetscCall(TestMatZeroRows(A, B, squaretest, lis, 0.0, PETSC_TRUE));
672: PetscCall(ISDestroy(&is));
673: PetscCall(ISDestroy(&lis));
675: /* test MatTranspose */
676: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n"));
677: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
678: PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &B2));
679: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose"));
681: PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &A2));
682: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose"));
683: PetscCall(MatDestroy(&A2));
685: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
686: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
687: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose"));
688: PetscCall(MatDestroy(&A2));
690: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
691: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose"));
692: PetscCall(MatDestroy(&A2));
693: PetscCall(MatDestroy(&B2));
695: /* test MatISFixLocalEmpty */
696: if (isaij) {
697: PetscInt r[2];
699: r[0] = 0;
700: r[1] = PetscMin(m, n) - 1;
701: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n"));
702: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
704: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
705: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
706: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
707: PetscCall(CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)"));
709: PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
710: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
711: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
712: PetscCall(MatZeroRows(B2, 2, r, 0.0, NULL, NULL));
713: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
714: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
715: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
716: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
717: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)"));
718: PetscCall(MatDestroy(&A2));
720: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
721: PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
722: PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
723: PetscCall(MatTranspose(B2, MAT_INPLACE_MATRIX, &B2));
724: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
725: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
726: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
727: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
728: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
729: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)"));
731: PetscCall(MatDestroy(&A2));
732: PetscCall(MatDestroy(&B2));
734: if (squaretest) {
735: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
736: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
737: PetscCall(MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL));
738: PetscCall(MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL));
739: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
740: PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
741: PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
742: PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
743: PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
744: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)"));
745: PetscCall(MatDestroy(&A2));
746: PetscCall(MatDestroy(&B2));
747: }
748: }
750: /* test MatInvertBlockDiagonal
751: special cases for block-diagonal matrices */
752: if (m == n) {
753: ISLocalToGlobalMapping map;
754: Mat Abd, Bbd;
755: IS is, bis;
756: const PetscScalar *isbd, *aijbd;
757: const PetscInt *sts, *idxs;
758: PetscInt *idxs2, diff, perm, nl, bs, st, en, in;
759: PetscBool ok;
761: for (diff = 0; diff < 3; diff++) {
762: for (perm = 0; perm < 3; perm++) {
763: for (bs = 1; bs < 4; bs++) {
764: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs));
765: PetscCall(PetscMalloc1(bs * bs, &vals));
766: PetscCall(MatGetOwnershipRanges(A, &sts));
767: switch (diff) {
768: case 1: /* inverted layout by processes */
769: in = 1;
770: st = sts[size - rank - 1];
771: en = sts[size - rank];
772: nl = en - st;
773: break;
774: case 2: /* round-robin layout */
775: in = size;
776: st = rank;
777: nl = n / size;
778: if (rank < n % size) nl++;
779: break;
780: default: /* same layout */
781: in = 1;
782: st = sts[rank];
783: en = sts[rank + 1];
784: nl = en - st;
785: break;
786: }
787: PetscCall(ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is));
788: PetscCall(ISGetLocalSize(is, &nl));
789: PetscCall(ISGetIndices(is, &idxs));
790: PetscCall(PetscMalloc1(nl, &idxs2));
791: for (i = 0; i < nl; i++) {
792: switch (perm) { /* invert some of the indices */
793: case 2:
794: idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
795: break;
796: case 1:
797: idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
798: break;
799: default:
800: idxs2[i] = idxs[i];
801: break;
802: }
803: }
804: PetscCall(ISRestoreIndices(is, &idxs));
805: PetscCall(ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis));
806: PetscCall(ISLocalToGlobalMappingCreateIS(bis, &map));
807: PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd));
808: PetscCall(ISLocalToGlobalMappingDestroy(&map));
809: PetscCall(MatISSetPreallocation(Abd, bs, NULL, 0, NULL));
810: for (i = 0; i < nl; i++) {
811: for (PetscInt b1 = 0; b1 < bs; b1++)
812: for (PetscInt b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
813: PetscCall(MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES));
814: }
815: PetscCall(MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY));
816: PetscCall(MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY));
817: PetscCall(MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd));
818: PetscCall(MatInvertBlockDiagonal(Abd, &isbd));
819: PetscCall(MatInvertBlockDiagonal(Bbd, &aijbd));
820: PetscCall(MatGetLocalSize(Bbd, &nl, NULL));
821: ok = PETSC_TRUE;
822: for (i = 0; i < nl / bs; i++) {
823: for (PetscInt b1 = 0; b1 < bs; b1++) {
824: for (PetscInt b2 = 0; b2 < bs; b2++) {
825: if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
826: if (!ok) {
827: 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])));
828: break;
829: }
830: }
831: if (!ok) break;
832: }
833: if (!ok) break;
834: }
835: PetscCall(MatDestroy(&Abd));
836: PetscCall(MatDestroy(&Bbd));
837: PetscCall(PetscFree(vals));
838: PetscCall(ISDestroy(&is));
839: PetscCall(ISDestroy(&bis));
840: }
841: }
842: }
843: }
845: /* test MatGetDiagonalBlock */
846: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n"));
847: PetscCall(MatGetDiagonalBlock(A, &A2));
848: PetscCall(MatGetDiagonalBlock(B, &B2));
849: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
850: PetscCall(MatScale(A, 2.0));
851: PetscCall(MatScale(B, 2.0));
852: PetscCall(MatGetDiagonalBlock(A, &A2));
853: PetscCall(MatGetDiagonalBlock(B, &B2));
854: PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
856: /* test MatISSetAllowRepeated on a MATIS */
857: PetscCall(MatISSetAllowRepeated(A, allow_repeated));
858: if (allow_repeated) { /* original MATIS maybe with repeated entries, test assembling of local matrices */
859: Mat lA, lA2;
861: for (PetscInt i = 0; i < 1; i++) { /* TODO: make MatScatter inherit from MATSHELL and support MatProducts */
862: PetscBool usemult = PETSC_FALSE;
864: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
865: if (i) {
866: Mat tA;
868: usemult = PETSC_TRUE;
869: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries and local shell matrices\n"));
870: PetscCall(MatISGetLocalMat(A2, &lA2));
871: PetscCall(MatConvert(lA2, MATSHELL, MAT_INITIAL_MATRIX, &tA));
872: PetscCall(MatISRestoreLocalMat(A2, &lA2));
873: PetscCall(MatISSetLocalMat(A2, tA));
874: PetscCall(MatDestroy(&tA));
875: } else {
876: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries\n"));
877: }
878: PetscCall(MatISSetAllowRepeated(A2, PETSC_FALSE));
879: PetscCall(MatISGetLocalMat(A, &lA));
880: PetscCall(MatISGetLocalMat(A2, &lA2));
881: if (!repmap) PetscCall(CheckMat(lA, lA2, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
882: PetscCall(MatISRestoreLocalMat(A, &lA));
883: PetscCall(MatISRestoreLocalMat(A2, &lA2));
884: if (repmap) PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with repeated entries"));
885: else PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
886: PetscCall(MatDestroy(&A2));
887: }
888: } else { /* original matis with non-repeated entries, this should only recreate the local matrices */
889: Mat lA;
890: PetscBool flg;
892: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(true) with non repeated entries\n"));
893: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
894: PetscCall(MatISSetAllowRepeated(A2, PETSC_TRUE));
895: PetscCall(MatISGetLocalMat(A2, &lA));
896: PetscCall(MatAssembled(lA, &flg));
897: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local mat should be unassembled");
898: PetscCall(MatISRestoreLocalMat(A2, &lA));
899: PetscCall(MatDestroy(&A2));
900: }
902: /* Test MatZeroEntries */
903: PetscCall(MatZeroEntries(A));
904: PetscCall(MatZeroEntries(B));
905: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatZeroEntries"));
907: /* Test MatSetValues and MatSetValuesBlocked */
908: if (test_setvalues) {
909: PetscCall(PetscMalloc1(lm * ln, &vals));
910: for (i = 0; i < lm * ln; i++) vals[i] = i + 1.0;
911: PetscCall(MatGetLocalSize(A, NULL, &ln));
912: PetscCall(MatISSetPreallocation(A, ln, NULL, n - ln, NULL));
913: PetscCall(MatSeqAIJSetPreallocation(B, ln, NULL));
914: PetscCall(MatMPIAIJSetPreallocation(B, ln, NULL, n - ln, NULL));
915: PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
916: PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
918: PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
919: PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
920: PetscCall(MatSetValues(A, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
921: PetscCall(MatSetValues(B, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
922: PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
923: PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
924: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
925: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
926: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
927: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
928: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValues"));
930: PetscCall(ISLocalToGlobalMappingGetBlockIndices(rmap, &idxs1));
931: PetscCall(ISLocalToGlobalMappingGetBlockIndices(cmap, &idxs2));
932: PetscCall(MatSetValuesBlocked(A, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
933: PetscCall(MatSetValuesBlocked(B, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
934: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rmap, &idxs1));
935: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cmap, &idxs2));
936: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
937: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
938: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
939: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
940: PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValuesBlocked"));
941: PetscCall(PetscFree(vals));
942: }
944: /* free testing matrices */
945: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
946: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
947: PetscCall(MatDestroy(&A));
948: PetscCall(MatDestroy(&B));
949: for (PetscInt bs = 1; bs <= 2; bs++) {
950: for (PetscInt mask = 0; mask < 5; mask++) {
951: PetscCall(CheckRepeatedMapFiltering(PETSC_COMM_WORLD, bs, mask, INSERT_VALUES));
952: PetscCall(CheckRepeatedMapFiltering(PETSC_COMM_WORLD, bs, mask, ADD_VALUES));
953: }
954: }
955: PetscCall(PetscFinalize());
956: return 0;
957: }
959: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
960: {
961: Mat Bcheck;
962: PetscReal error;
964: PetscFunctionBeginUser;
965: if (!usemult && B) {
966: PetscBool hasnorm;
968: PetscCall(MatHasOperation(B, MATOP_NORM, &hasnorm));
969: if (!hasnorm) usemult = PETSC_TRUE;
970: }
971: if (!usemult) {
972: if (B) {
973: MatType Btype;
975: PetscCall(MatGetType(B, &Btype));
976: PetscCall(MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck));
977: } else {
978: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
979: }
980: if (B) { /* if B is present, subtract it */
981: PetscCall(MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN));
982: }
983: PetscCall(MatNorm(Bcheck, NORM_INFINITY, &error));
984: if (error > PETSC_SQRT_MACHINE_EPSILON) {
985: ISLocalToGlobalMapping rl2g, cl2g;
987: PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Bcheck"));
988: PetscCall(MatView(Bcheck, NULL));
989: if (B) {
990: PetscCall(PetscObjectSetName((PetscObject)B, "B"));
991: PetscCall(MatView(B, NULL));
992: PetscCall(MatDestroy(&Bcheck));
993: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
994: PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Assembled A"));
995: PetscCall(MatView(Bcheck, NULL));
996: }
997: PetscCall(MatDestroy(&Bcheck));
998: PetscCall(PetscObjectSetName((PetscObject)A, "A"));
999: PetscCall(MatView(A, NULL));
1000: PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
1001: if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g, NULL));
1002: if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g, NULL));
1003: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
1004: }
1005: PetscCall(MatDestroy(&Bcheck));
1006: } else {
1007: PetscBool ok, okt;
1009: PetscCall(MatMultEqual(A, B, 3, &ok));
1010: PetscCall(MatMultTransposeEqual(A, B, 3, &okt));
1011: PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ? %d, multtranspose ok ? %d", func, ok, okt);
1012: }
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: PetscErrorCode CheckRepeatedMapFiltering(MPI_Comm comm, PetscInt bs, PetscInt mask, InsertMode mode)
1017: {
1018: Mat A, B, local, expected;
1019: ISLocalToGlobalMapping rmap, cmap;
1020: PetscInt rows[3], cols[3], nr = 0, nc = 0, i, j, r, c;
1021: PetscMPIInt rank, size;
1023: PetscFunctionBeginUser;
1024: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1025: PetscCallMPI(MPI_Comm_size(comm, &size));
1026: /* Only rank zero masks entries, so the other ranks retain every local occurrence. */
1027: for (i = 0; i < 3; i++) {
1028: rows[i] = cols[i] = rank;
1029: if (!rank && (mask == 4 || (i == 1 && (mask == 1 || mask == 3)))) rows[i] = -1;
1030: if (!rank && (mask == 4 || (i == 1 && (mask == 2 || mask == 3)))) cols[i] = -1;
1031: nr += rows[i] >= 0;
1032: nc += cols[i] >= 0;
1033: }
1034: PetscCall(ISLocalToGlobalMappingCreate(comm, bs, 3, rows, PETSC_COPY_VALUES, &rmap));
1035: if (mask >= 3) {
1036: PetscCall(PetscObjectReference((PetscObject)rmap));
1037: cmap = rmap;
1038: } else PetscCall(ISLocalToGlobalMappingCreate(comm, bs, 3, cols, PETSC_COPY_VALUES, &cmap));
1039: PetscCall(MatCreate(comm, &A));
1040: PetscCall(MatSetSizes(A, bs, bs, bs * size, bs * size));
1041: PetscCall(MatSetType(A, MATIS));
1042: PetscCall(MatISSetAllowRepeated(A, PETSC_TRUE));
1043: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
1044: PetscCall(MatISSetPreallocation(A, 3 * bs, NULL, 3 * bs, NULL));
1045: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, nr * bs, nc * bs, nc * bs, NULL, &expected));
1046: PetscCall(MatCreateAIJ(comm, bs, bs, bs * size, bs * size, bs, NULL, 0, NULL, &B));
1047: for (i = 0, r = 0; i < 3 * bs; i++) {
1048: for (j = 0, c = 0; j < 3 * bs; j++) {
1049: PetscScalar value = 1 + i * 3 * bs + j;
1051: PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &value, mode));
1052: if (rows[i / bs] >= 0 && cols[j / bs] >= 0) {
1053: PetscCall(MatSetValue(expected, r, c, value, INSERT_VALUES));
1054: PetscCall(MatSetValue(B, bs * rows[i / bs] + i % bs, bs * cols[j / bs] + j % bs, value, ADD_VALUES));
1055: }
1056: c += cols[j / bs] >= 0;
1057: }
1058: r += rows[i / bs] >= 0;
1059: }
1060: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1061: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1062: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1063: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1064: PetscCall(MatAssemblyBegin(expected, MAT_FINAL_ASSEMBLY));
1065: PetscCall(MatAssemblyEnd(expected, MAT_FINAL_ASSEMBLY));
1066: PetscCall(MatISGetLocalMat(A, &local));
1067: PetscCall(CheckMat(local, expected, PETSC_FALSE, "local matrix with repeated and masked indices"));
1068: PetscCall(MatISRestoreLocalMat(A, &local));
1069: PetscCall(CheckMat(A, B, PETSC_FALSE, "global matrix with repeated and masked indices"));
1070: PetscCall(MatDestroy(&expected));
1071: PetscCall(MatDestroy(&B));
1072: PetscCall(MatDestroy(&A));
1073: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
1074: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: PetscErrorCode CheckVariableBlockSizes(Mat A)
1079: {
1080: Mat lA;
1081: const PetscInt *bsizes;
1082: PetscInt m, nblocks, sum = 0;
1084: PetscFunctionBeginUser;
1085: PetscCall(MatISGetLocalMat(A, &lA));
1086: PetscCall(MatGetLocalSize(lA, &m, NULL));
1087: PetscCall(MatGetVariableBlockSizes(lA, &nblocks, &bsizes));
1088: for (PetscInt i = 0; i < nblocks; i++) sum += bsizes[i];
1089: PetscCheck(sum == m && (!m || nblocks), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid variable block sizes on loaded local matrix");
1090: PetscCall(MatISRestoreLocalMat(A, &lA));
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: static PetscErrorCode MatMult_Local(Mat A, Vec x, Vec y)
1095: {
1096: Mat local;
1098: PetscFunctionBeginUser;
1099: PetscCall(MatShellGetContext(A, &local));
1100: PetscCall(MatMult(local, x, y));
1101: PetscFunctionReturn(PETSC_SUCCESS);
1102: }
1104: static PetscErrorCode MatMultTranspose_Local(Mat A, Vec x, Vec y)
1105: {
1106: Mat local;
1108: PetscFunctionBeginUser;
1109: PetscCall(MatShellGetContext(A, &local));
1110: PetscCall(MatMultTranspose(local, x, y));
1111: PetscFunctionReturn(PETSC_SUCCESS);
1112: }
1114: PetscErrorCode TestMatZeroRowsShell(Mat A, Mat Afull, IS rows, PetscScalar diag)
1115: {
1116: Mat B, Bcheck, lB, local, shell;
1117: PetscInt m, n;
1119: PetscFunctionBeginUser;
1120: for (PetscInt columns = 0; columns < 2; columns++) {
1121: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
1122: PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck));
1123: PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1124: PetscCall(MatISGetLocalMat(B, &lB));
1125: local = lB;
1126: PetscCall(PetscObjectReference((PetscObject)local));
1127: PetscCall(MatGetLocalSize(local, &m, &n));
1128: PetscCall(MatISRestoreLocalMat(B, &lB));
1129: PetscCall(MatCreateShell(PETSC_COMM_SELF, m, n, m, n, local, &shell));
1130: PetscCall(MatShellSetOperation(shell, MATOP_MULT, (PetscErrorCodeFn *)MatMult_Local));
1131: PetscCall(MatShellSetOperation(shell, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_Local));
1132: PetscCall(MatISSetLocalMat(B, shell));
1133: PetscCall(MatDestroy(&shell));
1134: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1135: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1136: if (columns) {
1137: PetscCall(MatZeroRowsColumnsIS(B, rows, diag, NULL, NULL));
1138: PetscCall(MatZeroRowsColumnsIS(Bcheck, rows, diag, NULL, NULL));
1139: } else {
1140: PetscCall(MatZeroRowsIS(B, rows, diag, NULL, NULL));
1141: PetscCall(MatZeroRowsIS(Bcheck, rows, diag, NULL, NULL));
1142: }
1143: PetscCall(CheckMat(B, Bcheck, PETSC_TRUE, "Zero rows with a local shell matrix"));
1144: PetscCall(MatDestroy(&B));
1145: PetscCall(MatDestroy(&Bcheck));
1146: PetscCall(MatDestroy(&local));
1147: }
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: PetscErrorCode TestMatZeroRowsMissingDiagonal(MPI_Comm comm)
1152: {
1153: Mat A, B, local;
1154: ISLocalToGlobalMapping l2g;
1155: const PetscInt idxs[] = {0, 1};
1156: const PetscInt *cols;
1157: PetscInt ncols;
1158: const PetscScalar values[] = {2., 1.};
1159: PetscMPIInt rank;
1161: PetscFunctionBeginUser;
1162: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1163: PetscCall(ISLocalToGlobalMappingCreate(comm, 1, 2, idxs, PETSC_COPY_VALUES, &l2g));
1164: for (PetscInt columns = 0; columns < 2; columns++) {
1165: PetscCall(MatCreateIS(comm, 1, PETSC_DECIDE, PETSC_DECIDE, 2, 2, l2g, l2g, &A));
1166: PetscCall(MatISSetPreallocation(A, 2, NULL, 0, NULL));
1167: // Leave (1,1) structurally absent only on rank 0.
1168: PetscCall(MatSetValuesLocal(A, 1, idxs, 2, idxs, values, INSERT_VALUES));
1169: PetscCall(MatSetValueLocal(A, 1, 0, 1., INSERT_VALUES));
1170: if (rank) PetscCall(MatSetValueLocal(A, 1, 1, 0., INSERT_VALUES));
1171: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1172: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1173: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
1174: if (columns) {
1175: PetscCall(MatZeroRowsColumns(A, rank ? 0 : 1, idxs, 2., NULL, NULL));
1176: PetscCall(MatZeroRowsColumns(B, rank ? 0 : 1, idxs, 2., NULL, NULL));
1177: } else {
1178: PetscCall(MatZeroRows(A, rank ? 0 : 1, idxs, 2., NULL, NULL));
1179: PetscCall(MatZeroRows(B, rank ? 0 : 1, idxs, 2., NULL, NULL));
1180: }
1181: PetscCall(MatISGetLocalMat(A, &local));
1182: PetscCall(MatGetRow(local, 1, &ncols, &cols, NULL));
1183: PetscCheck(ncols == (rank ? 2 : 1) && cols[0] == 0 && (!rank || cols[1] == 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Zeroing row 0 changed the nonzero pattern of row 1");
1184: PetscCall(MatRestoreRow(local, 1, &ncols, &cols, NULL));
1185: PetscCall(MatISRestoreLocalMat(A, &local));
1186: PetscCall(CheckMat(A, B, PETSC_FALSE, "Zero rows with a missing local diagonal entry"));
1187: PetscCall(MatDestroy(&A));
1188: PetscCall(MatDestroy(&B));
1189: }
1190: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
1191: PetscFunctionReturn(PETSC_SUCCESS);
1192: }
1194: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag, PetscBool local)
1195: {
1196: Mat B, Bcheck, B2 = NULL, lB;
1197: Vec x = NULL, b = NULL, b2 = NULL;
1198: ISLocalToGlobalMapping l2gr, l2gc;
1199: PetscReal error;
1200: char diagstr[16];
1201: const PetscInt *idxs;
1202: PetscInt i, n, N;
1203: PetscBool haszerorows;
1204: IS gis;
1206: PetscFunctionBeginUser;
1207: if (diag == 0.) {
1208: PetscCall(PetscStrncpy(diagstr, "zero", sizeof(diagstr)));
1209: } else {
1210: PetscCall(PetscStrncpy(diagstr, "nonzero", sizeof(diagstr)));
1211: }
1212: PetscCall(ISView(is, NULL));
1213: PetscCall(MatGetLocalToGlobalMapping(A, &l2gr, &l2gc));
1214: /* tests MatDuplicate and MatCopy */
1215: if (diag == 0.) PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
1216: else {
1217: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
1218: PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
1219: }
1220: PetscCall(MatISGetLocalMat(B, &lB));
1221: PetscCall(MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows));
1222: if (squaretest && haszerorows) {
1223: PetscCall(MatCreateVecs(B, &x, &b));
1224: PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
1225: PetscCall(VecSetLocalToGlobalMapping(b, l2gr));
1226: PetscCall(VecSetLocalToGlobalMapping(x, l2gc));
1227: PetscCall(VecSetRandom(x, NULL));
1228: PetscCall(VecSetRandom(b, NULL));
1229: /* mimic b[is] = x[is] */
1230: PetscCall(VecDuplicate(b, &b2));
1231: PetscCall(VecSetLocalToGlobalMapping(b2, l2gr));
1232: PetscCall(VecCopy(b, b2));
1233: if (local) {
1234: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1235: PetscCall(ISGetLocalSize(gis, &n));
1236: PetscCall(ISGetIndices(gis, &idxs));
1237: } else {
1238: PetscCall(ISGetLocalSize(is, &n));
1239: PetscCall(ISGetIndices(is, &idxs));
1240: }
1241: PetscCall(VecGetSize(x, &N));
1242: for (i = 0; i < n; i++) {
1243: if (0 <= idxs[i] && idxs[i] < N) {
1244: PetscCall(VecSetValue(b2, idxs[i], diag, INSERT_VALUES));
1245: PetscCall(VecSetValue(x, idxs[i], 1., INSERT_VALUES));
1246: }
1247: }
1248: if (local) {
1249: PetscCall(ISRestoreIndices(gis, &idxs));
1250: PetscCall(ISDestroy(&gis));
1251: } else {
1252: PetscCall(ISRestoreIndices(is, &idxs));
1253: }
1254: PetscCall(VecAssemblyBegin(b2));
1255: PetscCall(VecAssemblyEnd(b2));
1256: PetscCall(VecAssemblyBegin(x));
1257: PetscCall(VecAssemblyEnd(x));
1258: /* test ZeroRows on MATIS */
1259: if (local) {
1260: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsLocal (diag %s)\n", diagstr));
1261: PetscCall(MatZeroRowsLocalIS(B, is, diag, x, b));
1262: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumnsLocal (diag %s)\n", diagstr));
1263: PetscCall(MatZeroRowsColumnsLocalIS(B2, is, diag, NULL, NULL));
1264: } else {
1265: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1266: PetscCall(MatZeroRowsIS(B, is, diag, x, b));
1267: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr));
1268: PetscCall(MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL));
1269: }
1270: } else if (haszerorows) {
1271: /* test ZeroRows on MATIS */
1272: if (local) {
1273: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsLocal (diag %s)\n", diagstr));
1274: PetscCall(MatZeroRowsLocalIS(B, is, diag, NULL, NULL));
1275: } else {
1276: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1277: PetscCall(MatZeroRowsIS(B, is, diag, NULL, NULL));
1278: }
1279: b = b2 = x = NULL;
1280: } else {
1281: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr));
1282: b = b2 = x = NULL;
1283: }
1285: if (squaretest && haszerorows) {
1286: PetscCall(VecAXPY(b2, -1., b));
1287: PetscCall(VecNorm(b2, NORM_INFINITY, &error));
1288: PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr);
1289: }
1290: PetscCall(VecDestroy(&x));
1291: PetscCall(VecDestroy(&b));
1292: PetscCall(VecDestroy(&b2));
1294: /* check the result of ZeroRows with that from MPIAIJ routines
1295: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
1296: if (haszerorows) {
1297: PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck));
1298: PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1299: if (local) {
1300: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1301: PetscCall(MatZeroRowsIS(Bcheck, gis, diag, NULL, NULL));
1302: PetscCall(ISDestroy(&gis));
1303: } else {
1304: PetscCall(MatZeroRowsIS(Bcheck, is, diag, NULL, NULL));
1305: }
1306: PetscCall(CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows"));
1307: PetscCall(MatDestroy(&Bcheck));
1308: }
1309: PetscCall(MatDestroy(&B));
1311: if (B2) { /* test MatZeroRowsColumns */
1312: PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &B));
1313: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1314: if (local) {
1315: PetscCall(ISL2GMapNoNeg(l2gr, is, &gis));
1316: PetscCall(MatZeroRowsColumnsIS(B, gis, diag, NULL, NULL));
1317: PetscCall(ISDestroy(&gis));
1318: } else {
1319: PetscCall(MatZeroRowsColumnsIS(B, is, diag, NULL, NULL));
1320: }
1321: PetscCall(CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns"));
1322: PetscCall(MatDestroy(&B));
1323: PetscCall(MatDestroy(&B2));
1324: }
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }
1328: PetscErrorCode ISL2GMapNoNeg(ISLocalToGlobalMapping mapping, IS is, IS *newis)
1329: {
1330: PetscInt n, *idxout, nn = 0;
1331: const PetscInt *idxin;
1333: PetscFunctionBegin;
1334: PetscCall(ISGetLocalSize(is, &n));
1335: PetscCall(ISGetIndices(is, &idxin));
1336: PetscCall(PetscMalloc1(n, &idxout));
1337: PetscCall(ISLocalToGlobalMappingApply(mapping, n, idxin, idxout));
1338: PetscCall(ISRestoreIndices(is, &idxin));
1339: for (PetscInt i = 0; i < n; i++)
1340: if (idxout[i] > -1) idxout[nn++] = idxout[i];
1341: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)is), nn, idxout, PETSC_OWN_POINTER, newis));
1342: PetscFunctionReturn(PETSC_SUCCESS);
1343: }
1345: /*TEST
1347: test:
1348: requires: !complex
1349: args: -test_matlab -test_trans
1351: test:
1352: suffix: 2
1353: nsize: 4
1354: args: -mat_is_convert_local_nest -nr 3 -nc 4
1356: test:
1357: suffix: 3
1358: nsize: 5
1359: args: -m 11 -n 10 -mat_is_convert_local_nest -nr 2 -nc 1 -cbs 2
1361: test:
1362: suffix: 4
1363: nsize: 6
1364: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
1366: test:
1367: suffix: 5
1368: nsize: 6
1369: args: -m 12 -n 12 -test_trans -nr 3 -nc 1 -rbs 2 -test_variableblocksizes -mat_is_view_variableblocksizes
1371: test:
1372: suffix: 6
1373: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -rbs 6 -cbs 3
1375: test:
1376: suffix: 7
1377: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
1379: test:
1380: suffix: 8
1381: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
1383: test:
1384: suffix: 9
1385: nsize: 5
1386: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
1388: test:
1389: suffix: 10
1390: nsize: 5
1391: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
1393: test:
1394: suffix: vscat_default
1395: nsize: 5
1396: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
1397: output_file: output/ex23_11.out
1399: test:
1400: suffix: 12
1401: nsize: 3
1402: args: -m 12 -n 12 -symmetric -mat_is_localmat_type sbaij -test_trans -nr 2 -nc 3 -test_setvalues 0
1404: testset:
1405: output_file: output/ex23_13.out
1406: nsize: 3
1407: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
1408: filter: grep -v "type:" | grep -v "not using I-node routines"
1409: test:
1410: suffix: baij
1411: args: -mat_is_localmat_type baij
1412: test:
1413: requires: viennacl
1414: suffix: viennacl
1415: args: -mat_is_localmat_type aijviennacl
1416: test:
1417: requires: cuda
1418: suffix: cusparse
1419: args: -mat_is_localmat_type aijcusparse
1420: test:
1421: requires: kokkos_kernels
1422: suffix: kokkos
1423: args: -mat_is_localmat_type aijkokkos
1425: test:
1426: suffix: negrep
1427: nsize: {{1 3}separate output}
1428: 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
1430: test:
1431: suffix: negrep_allowrep
1432: nsize: {{1 3}separate output}
1433: 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
1435: TEST*/