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