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