Actual source code: ex23.c

  1: static char help[] = "Tests the use of interface functions for MATIS matrices and conversion routines.\n";

  3: #include <petscmat.h>

  5: PetscErrorCode TestMatZeroRows(Mat, Mat, PetscBool, IS, PetscScalar);
  6: PetscErrorCode CheckMat(Mat, Mat, PetscBool, const char *);

  8: int main(int argc, char **args)
  9: {
 10:   Mat                    A, B, A2, B2, T;
 11:   Mat                    Aee, Aeo, Aoe, Aoo;
 12:   Mat                   *mats, *Asub, *Bsub;
 13:   Vec                    x, y;
 14:   MatInfo                info;
 15:   ISLocalToGlobalMapping cmap, rmap;
 16:   IS                     is, is2, reven, rodd, ceven, codd;
 17:   IS                    *rows, *cols;
 18:   IS                     irow[2], icol[2];
 19:   PetscLayout            rlayout, clayout;
 20:   const PetscInt        *rrange, *crange, *idxs1, *idxs2;
 21:   MatType                lmtype;
 22:   PetscScalar            diag = 2., *vals;
 23:   PetscInt               n, m, i, lm, ln;
 24:   PetscInt               rst, ren, cst, cen, nr, nc, rbs = 1, cbs = 1;
 25:   PetscMPIInt            rank, size, lrank, rrank;
 26:   PetscBool              testT, squaretest, isaij;
 27:   PetscBool              permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE, allow_repeated = PETSC_TRUE;
 28:   PetscBool              diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric, test_matlab = PETSC_FALSE, test_setvalues = PETSC_TRUE;

 30:   PetscFunctionBeginUser;
 31:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 32:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 33:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 34:   m = n = 2 * size;
 35:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL));
 36:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 37:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 38:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-negmap", &negmap, NULL));
 39:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-repmap", &repmap, NULL));
 40:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-permmap", &permute, NULL));
 41:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-diffmap", &diffmap, NULL));
 42:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-allow_repeated", &allow_repeated, NULL));
 43:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matlab", &test_matlab, NULL));
 44:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_setvalues", &test_setvalues, NULL));
 45:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-rbs", &rbs, NULL));
 46:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-cbs", &cbs, NULL));
 47:   PetscCheck(size == 1 || m >= 4, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 4 for parallel runs");
 48:   PetscCheck(size != 1 || m >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of rows should be larger or equal 2 for uniprocessor runs");
 49:   PetscCheck(n >= 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Number of cols should be larger or equal 2");

 51:   /* create a MATIS matrix */
 52:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 53:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
 54:   PetscCall(MatSetType(A, MATIS));
 55:   PetscCall(MatSetFromOptions(A));
 56:   if (!negmap && !repmap) {
 57:     /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
 58:        Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
 59:        Equivalent to passing NULL for the mapping */
 60:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is));
 61:   } else if (negmap && !repmap) { /* non repeated but with negative indices */
 62:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is));
 63:   } else if (!negmap && repmap) { /* non negative but repeated indices */
 64:     IS isl[2];

 66:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]));
 67:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]));
 68:     PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
 69:     PetscCall(ISDestroy(&isl[0]));
 70:     PetscCall(ISDestroy(&isl[1]));
 71:   } else { /* negative and repeated indices */
 72:     IS isl[2];

 74:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]));
 75:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]));
 76:     PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is));
 77:     PetscCall(ISDestroy(&isl[0]));
 78:     PetscCall(ISDestroy(&isl[1]));
 79:   }
 80:   PetscCall(ISLocalToGlobalMappingCreateIS(is, &cmap));
 81:   PetscCall(ISDestroy(&is));

 83:   if (m != n || diffmap) {
 84:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is));
 85:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rmap));
 86:     PetscCall(ISDestroy(&is));
 87:   } else {
 88:     PetscCall(PetscObjectReference((PetscObject)cmap));
 89:     rmap = cmap;
 90:   }

 92:   PetscCall(MatISSetAllowRepeated(A, allow_repeated));
 93:   PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
 94:   PetscCall(MatSetBlockSizes(A, rbs, cbs));
 95:   PetscCall(MatISStoreL2L(A, PETSC_FALSE));
 96:   PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));
 97:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
 98:   PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
 99:   PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));
100:   for (i = 0; i < lm; i++) {
101:     PetscScalar v[3];
102:     PetscInt    cols[3];

104:     cols[0] = (i - 1 + n) % n;
105:     cols[1] = i % n;
106:     cols[2] = (i + 1) % n;
107:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
108:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
109:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
110:     PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
111:     PetscCall(MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES));
112:   }
113:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
114:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

116:   /* activate tests for square matrices with same maps only */
117:   PetscCall(MatHasCongruentLayouts(A, &squaretest));
118:   if (squaretest && rmap != cmap) {
119:     PetscInt nr, nc;

121:     PetscCall(ISLocalToGlobalMappingGetSize(rmap, &nr));
122:     PetscCall(ISLocalToGlobalMappingGetSize(cmap, &nc));
123:     if (nr != nc) squaretest = PETSC_FALSE;
124:     else {
125:       PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
126:       PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
127:       PetscCall(PetscArraycmp(idxs1, idxs2, nr, &squaretest));
128:       PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
129:       PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
130:     }
131:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
132:   }
133:   if (negmap && repmap) squaretest = PETSC_FALSE;

135:   /* test MatISGetLocalMat */
136:   PetscCall(MatISGetLocalMat(A, &B));
137:   PetscCall(MatGetType(B, &lmtype));

139:   /* test MatGetInfo */
140:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n"));
141:   PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
142:   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
143:   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,
144:                                                (PetscInt)info.nz_unneeded, (PetscInt)info.assemblies, (PetscInt)info.mallocs));
145:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
146:   PetscCall(MatGetInfo(A, MAT_GLOBAL_MAX, &info));
147:   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,
148:                                    (PetscInt)info.assemblies, (PetscInt)info.mallocs));
149:   PetscCall(MatGetInfo(A, MAT_GLOBAL_SUM, &info));
150:   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,
151:                                    (PetscInt)info.assemblies, (PetscInt)info.mallocs));

153:   /* test MatIsSymmetric */
154:   PetscCall(MatIsSymmetric(A, 0.0, &issymmetric));
155:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric));

157:   /* Create a MPIAIJ matrix, same as A */
158:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
159:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
160:   PetscCall(MatSetBlockSizes(B, rbs, cbs));
161:   PetscCall(MatSetType(B, MATAIJ));
162:   PetscCall(MatSetFromOptions(B));
163:   PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
164:   PetscCall(MatMPIAIJSetPreallocation(B, 3, NULL, 3, NULL));
165:   PetscCall(MatMPIBAIJSetPreallocation(B, 1, 3, NULL, 3, NULL));
166: #if defined(PETSC_HAVE_HYPRE)
167:   PetscCall(MatHYPRESetPreallocation(B, 3, NULL, 3, NULL));
168: #endif
169:   PetscCall(MatISSetPreallocation(B, 3, NULL, 3, NULL));
170:   PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool)!(repmap || negmap))); /* I do not want to precompute the pattern */
171:   for (i = 0; i < lm; i++) {
172:     PetscScalar v[3];
173:     PetscInt    cols[3];

175:     cols[0] = (i - 1 + n) % n;
176:     cols[1] = i % n;
177:     cols[2] = (i + 1) % n;
178:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
179:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
180:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
181:     PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
182:     PetscCall(MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES));
183:   }
184:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
185:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

187:   /* test MatView */
188:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n"));
189:   PetscCall(MatView(A, NULL));
190:   PetscCall(MatView(B, NULL));

192:   /* test MATLAB ASCII view */
193:   if (test_matlab) { /* output is different when using real or complex numbers */
194:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView ASCII MATLAB\n"));
195:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB));
196:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
197:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
198:   }

200:   /* test CheckMat */
201:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n"));
202:   PetscCall(CheckMat(A, B, PETSC_FALSE, "CheckMat"));

204:   /* test binary MatView/MatLoad */
205:   {
206:     PetscMPIInt color = rank % 2;
207:     MPI_Comm    comm;
208:     char        name[PETSC_MAX_PATH_LEN];
209:     PetscViewer wview, cview, sview, view;
210:     Mat         A2;

212:     PetscCallMPI(MPI_Comm_split(PETSC_COMM_WORLD, color, rank, &comm));

214:     PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "world_is"));
215:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_WRITE, &wview));
216:     PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "seq_is_%d", rank));
217:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_WRITE, &sview));
218:     PetscCall(PetscSNPrintf(name, PETSC_STATIC_ARRAY_LENGTH(name), "color_is_%d", color));
219:     PetscCall(PetscViewerBinaryOpen(comm, name, FILE_MODE_WRITE, &cview));
220:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary world\n"));
221:     PetscCall(MatView(A, wview));
222:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary self\n"));
223:     PetscCall(MatView(A, sview));
224:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatView on binary subcomm\n"));
225:     PetscCall(MatView(A, cview));
226:     PetscCall(PetscViewerDestroy(&wview));
227:     PetscCall(PetscViewerDestroy(&cview));
228:     PetscCall(PetscViewerDestroy(&sview));

230:     /* Load a world matrix */
231:     PetscCall(MatCreate(PETSC_COMM_WORLD, &A2));
232:     PetscCall(MatSetType(A2, MATIS));
233:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));

235:     /* Read back the same matrix and check */
236:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from world\n"));
237:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "world_is", FILE_MODE_READ, &view));
238:     PetscCall(MatLoad(A2, view));
239:     PetscCall(CheckMat(A, A2, PETSC_TRUE, "Load"));
240:     PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
241:     PetscCall(PetscViewerDestroy(&view));

243:     /* Read the matrix from rank 0 only */
244:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from self\n"));
245:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "seq_is_0", FILE_MODE_READ, &view));
246:     PetscCall(MatLoad(A2, view));
247:     PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
248:     PetscCall(PetscViewerDestroy(&view));

250:     /* Read the matrix from subcomm */
251:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatLoad from subcomm\n"));
252:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "color_is_0", FILE_MODE_READ, &view));
253:     PetscCall(MatLoad(A2, view));
254:     PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_WORLD));
255:     PetscCall(PetscViewerDestroy(&view));

257:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
258:     PetscCall(MatDestroy(&A2));

260:     /* now load the original matrix from color 0 only processes */
261:     if (!color) {
262:       PetscCall(PetscPrintf(comm, "Test subcomm MatLoad from world\n"));
263:       PetscCall(MatCreate(comm, &A2));
264:       PetscCall(MatSetType(A2, MATIS));
265:       PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(comm), PETSC_VIEWER_ASCII_INFO_DETAIL));
266:       PetscCall(PetscViewerBinaryOpen(comm, "world_is", FILE_MODE_READ, &view));
267:       PetscCall(MatLoad(A2, view));
268:       PetscCall(MatView(A2, PETSC_VIEWER_STDOUT_(comm)));
269:       PetscCall(PetscViewerDestroy(&view));
270:       PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(comm)));
271:       PetscCall(MatDestroy(&A2));
272:     }

274:     PetscCallMPI(MPI_Comm_free(&comm));
275:   }

277:   /* test MatDuplicate and MatAXPY */
278:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n"));
279:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
280:   PetscCall(CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY"));

282:   /* test MatConvert */
283:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n"));
284:   PetscCall(MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2));
285:   PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX"));
286:   PetscCall(MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2));
287:   PetscCall(CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX"));
288:   PetscCall(MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2));
289:   PetscCall(CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX"));
290:   PetscCall(MatDestroy(&A2));
291:   PetscCall(MatDestroy(&B2));
292:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n"));
293:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
294:   PetscCall(MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2));
295:   PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX"));
296:   PetscCall(MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2));
297:   PetscCall(CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX"));
298:   PetscCall(MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2));
299:   PetscCall(CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX"));
300:   PetscCall(MatDestroy(&A2));
301:   PetscCall(MatDestroy(&B2));
302:   PetscCall(PetscStrcmp(lmtype, MATSEQAIJ, &isaij));
303:   if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
304:     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};
305:     ISLocalToGlobalMapping tcmap, trmap;

307:     for (ri = 0; ri < 2; ri++) {
308:       PetscInt *r;

310:       r = (PetscInt *)(ri == 0 ? rr : rk);
311:       for (ci = 0; ci < 2; ci++) {
312:         PetscInt *c, rb, cb;

314:         c = (PetscInt *)(ci == 0 ? cr : ck);
315:         for (rb = 1; rb < 4; rb++) {
316:           PetscCall(ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is));
317:           PetscCall(ISLocalToGlobalMappingCreateIS(is, &trmap));
318:           PetscCall(ISDestroy(&is));
319:           for (cb = 1; cb < 4; cb++) {
320:             Mat  T, lT, T2;
321:             char testname[256];

323:             PetscCall(PetscSNPrintf(testname, sizeof(testname), "MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")", ri, ci, rb, cb));
324:             PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test %s\n", testname));

326:             PetscCall(ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is));
327:             PetscCall(ISLocalToGlobalMappingCreateIS(is, &tcmap));
328:             PetscCall(ISDestroy(&is));

330:             PetscCall(MatCreate(PETSC_COMM_SELF, &T));
331:             PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4));
332:             PetscCall(MatSetType(T, MATIS));
333:             PetscCall(MatSetLocalToGlobalMapping(T, trmap, tcmap));
334:             PetscCall(ISLocalToGlobalMappingDestroy(&tcmap));
335:             PetscCall(MatISGetLocalMat(T, &lT));
336:             PetscCall(MatSetType(lT, MATSEQAIJ));
337:             PetscCall(MatSeqAIJSetPreallocation(lT, cb * 4, NULL));
338:             PetscCall(MatSetRandom(lT, NULL));
339:             PetscCall(MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT));
340:             PetscCall(MatISRestoreLocalMat(T, &lT));
341:             PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
342:             PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));

344:             PetscCall(MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2));
345:             PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX"));
346:             PetscCall(MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2));
347:             PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX"));
348:             PetscCall(MatDestroy(&T2));
349:             PetscCall(MatDuplicate(T, MAT_COPY_VALUES, &T2));
350:             PetscCall(MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2));
351:             PetscCall(CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX"));
352:             PetscCall(MatDestroy(&T));
353:             PetscCall(MatDestroy(&T2));
354:           }
355:           PetscCall(ISLocalToGlobalMappingDestroy(&trmap));
356:         }
357:       }
358:     }
359:   }

361:   /* test MatDiagonalScale */
362:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n"));
363:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
364:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
365:   PetscCall(MatCreateVecs(A, &x, &y));
366:   PetscCall(VecSetRandom(x, NULL));
367:   if (issymmetric) {
368:     PetscCall(VecCopy(x, y));
369:   } else {
370:     PetscCall(VecSetRandom(y, NULL));
371:     PetscCall(VecScale(y, 8.));
372:   }
373:   PetscCall(MatDiagonalScale(A2, y, x));
374:   PetscCall(MatDiagonalScale(B2, y, x));
375:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale"));
376:   PetscCall(MatDestroy(&A2));
377:   PetscCall(MatDestroy(&B2));
378:   PetscCall(VecDestroy(&x));
379:   PetscCall(VecDestroy(&y));

381:   /* test MatPtAP (A IS and B AIJ) */
382:   if (isaij && m == n) {
383:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n"));
384:     /* There's a bug in MatCreateSubMatrices_MPIAIJ I cannot figure out */
385:     if (!allow_repeated || !repmap || size == 1) {
386:       PetscCall(MatISStoreL2L(A, PETSC_TRUE));
387:       PetscCall(MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &A2));
388:       PetscCall(MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B2));
389:       PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX"));
390:       PetscCall(MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &A2));
391:       PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX"));
392:       PetscCall(MatDestroy(&A2));
393:       PetscCall(MatDestroy(&B2));
394:     }
395:   }

397:   /* test MatGetLocalSubMatrix */
398:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n"));
399:   PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2));
400:   PetscCall(ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven));
401:   PetscCall(ISComplement(reven, 0, lm, &rodd));
402:   PetscCall(ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven));
403:   PetscCall(ISComplement(ceven, 0, ln, &codd));
404:   PetscCall(MatGetLocalSubMatrix(A2, reven, ceven, &Aee));
405:   PetscCall(MatGetLocalSubMatrix(A2, reven, codd, &Aeo));
406:   PetscCall(MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe));
407:   PetscCall(MatGetLocalSubMatrix(A2, rodd, codd, &Aoo));
408:   for (i = 0; i < lm; i++) {
409:     PetscInt    j, je, jo, colse[3], colso[3];
410:     PetscScalar ve[3], vo[3];
411:     PetscScalar v[3];
412:     PetscInt    cols[3];
413:     PetscInt    row = i / 2;

415:     cols[0] = (i - 1 + n) % n;
416:     cols[1] = i % n;
417:     cols[2] = (i + 1) % n;
418:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
419:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
420:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
421:     PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
422:     for (j = 0, je = 0, jo = 0; j < 3; j++) {
423:       if (cols[j] % 2) {
424:         vo[jo]      = v[j];
425:         colso[jo++] = cols[j] / 2;
426:       } else {
427:         ve[je]      = v[j];
428:         colse[je++] = cols[j] / 2;
429:       }
430:     }
431:     if (i % 2) {
432:       PetscCall(MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES));
433:       PetscCall(MatSetValuesLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES));
434:     } else {
435:       PetscCall(MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES));
436:       PetscCall(MatSetValuesLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES));
437:     }
438:   }
439:   PetscCall(MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee));
440:   PetscCall(MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo));
441:   PetscCall(MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe));
442:   PetscCall(MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo));
443:   PetscCall(ISDestroy(&reven));
444:   PetscCall(ISDestroy(&ceven));
445:   PetscCall(ISDestroy(&rodd));
446:   PetscCall(ISDestroy(&codd));
447:   PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
448:   PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
449:   PetscCall(MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN));
450:   PetscCall(CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix"));
451:   PetscCall(MatDestroy(&A2));

453:   /* test MatConvert_Nest_IS */
454:   testT = PETSC_FALSE;
455:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL));

457:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n"));
458:   nr = 2;
459:   nc = 2;
460:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL));
461:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
462:   if (testT) {
463:     PetscCall(MatGetOwnershipRange(A, &cst, &cen));
464:     PetscCall(MatGetOwnershipRangeColumn(A, &rst, &ren));
465:   } else {
466:     PetscCall(MatGetOwnershipRange(A, &rst, &ren));
467:     PetscCall(MatGetOwnershipRangeColumn(A, &cst, &cen));
468:   }
469:   PetscCall(PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats));
470:   for (i = 0; i < nr * nc; i++) {
471:     if (testT) {
472:       PetscCall(MatCreateTranspose(A, &mats[i]));
473:       PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]));
474:     } else {
475:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &mats[i]));
476:       PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]));
477:     }
478:   }
479:   for (i = 0; i < nr; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]));
480:   for (i = 0; i < nc; i++) PetscCall(ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]));
481:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2));
482:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2));
483:   for (i = 0; i < nr; i++) PetscCall(ISDestroy(&rows[i]));
484:   for (i = 0; i < nc; i++) PetscCall(ISDestroy(&cols[i]));
485:   for (i = 0; i < 2 * nr * nc; i++) PetscCall(MatDestroy(&mats[i]));
486:   PetscCall(PetscFree3(rows, cols, mats));
487:   PetscCall(MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T));
488:   PetscCall(MatDestroy(&B2));
489:   PetscCall(MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2));
490:   PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX"));
491:   PetscCall(MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2));
492:   PetscCall(CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX"));
493:   PetscCall(MatDestroy(&B2));
494:   PetscCall(MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2));
495:   PetscCall(CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX"));
496:   PetscCall(MatDestroy(&T));
497:   PetscCall(MatDestroy(&A2));

499:   /* test MatCreateSubMatrix */
500:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n"));
501:   if (rank == 0) {
502:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is));
503:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2));
504:   } else if (rank == 1) {
505:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
506:     if (n > 3) {
507:       PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2));
508:     } else {
509:       PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
510:     }
511:   } else if (rank == 2 && n > 4) {
512:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
513:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2));
514:   } else {
515:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
516:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2));
517:   }
518:   PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2));
519:   PetscCall(MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2));
520:   PetscCall(CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix"));

522:   PetscCall(MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2));
523:   PetscCall(MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2));
524:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix"));
525:   PetscCall(MatDestroy(&A2));
526:   PetscCall(MatDestroy(&B2));

528:   if (!issymmetric) {
529:     PetscCall(MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2));
530:     PetscCall(MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2));
531:     PetscCall(MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2));
532:     PetscCall(MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2));
533:     PetscCall(CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix"));
534:   }

536:   PetscCall(MatDestroy(&A2));
537:   PetscCall(MatDestroy(&B2));
538:   PetscCall(ISDestroy(&is));
539:   PetscCall(ISDestroy(&is2));

541:   /* test MatCreateSubMatrices */
542:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n"));
543:   PetscCall(MatGetLayouts(A, &rlayout, &clayout));
544:   PetscCall(PetscLayoutGetRanges(rlayout, &rrange));
545:   PetscCall(PetscLayoutGetRanges(clayout, &crange));
546:   lrank = (size + rank - 1) % size;
547:   rrank = (rank + 1) % size;
548:   PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[lrank + 1] - rrange[lrank], rrange[lrank], 1, &irow[0]));
549:   PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[rrank + 1] - crange[rrank], crange[rrank], 1, &icol[0]));
550:   PetscCall(ISCreateStride(PETSC_COMM_SELF, rrange[rrank + 1] - rrange[rrank], rrange[rrank], 1, &irow[1]));
551:   PetscCall(ISCreateStride(PETSC_COMM_SELF, crange[lrank + 1] - crange[lrank], crange[lrank], 1, &icol[1]));
552:   PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub));
553:   PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub));
554:   PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
555:   PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
556:   PetscCall(MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub));
557:   PetscCall(MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub));
558:   PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]"));
559:   PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]"));
560:   PetscCall(MatDestroySubMatrices(2, &Asub));
561:   PetscCall(MatDestroySubMatrices(2, &Bsub));
562:   PetscCall(ISDestroy(&irow[0]));
563:   PetscCall(ISDestroy(&irow[1]));
564:   PetscCall(ISDestroy(&icol[0]));
565:   PetscCall(ISDestroy(&icol[1]));

567:   /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
568:   if (size > 1) {
569:     if (rank == 0) {
570:       PetscInt st, len;

572:       st  = (m + 1) / 2;
573:       len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
574:       PetscCall(ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is));
575:     } else {
576:       PetscCall(ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is));
577:     }
578:   } else {
579:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is));
580:   }

582:   if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
583:     PetscInt *idx0, *idx1, n0, n1;
584:     IS        Ais[2], Bis[2];

586:     /* test MatDiagonalSet */
587:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n"));
588:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
589:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
590:     PetscCall(MatCreateVecs(A, NULL, &x));
591:     PetscCall(VecSetRandom(x, NULL));
592:     PetscCall(MatDiagonalSet(A2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
593:     PetscCall(MatDiagonalSet(B2, x, allow_repeated ? ADD_VALUES : INSERT_VALUES));
594:     PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet"));
595:     PetscCall(VecDestroy(&x));
596:     PetscCall(MatDestroy(&A2));
597:     PetscCall(MatDestroy(&B2));

599:     /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
600:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n"));
601:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
602:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
603:     PetscCall(MatShift(A2, 2.0));
604:     PetscCall(MatShift(B2, 2.0));
605:     PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatShift"));
606:     PetscCall(MatDestroy(&A2));
607:     PetscCall(MatDestroy(&B2));

609:     /* nonzero diag value is supported for square matrices only */
610:     PetscCall(TestMatZeroRows(A, B, PETSC_TRUE, is, diag));

612:     /* test MatIncreaseOverlap */
613:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n"));
614:     PetscCall(MatGetOwnershipRange(A, &rst, &ren));
615:     n0 = (ren - rst) / 2;
616:     n1 = (ren - rst) / 3;
617:     PetscCall(PetscMalloc1(n0, &idx0));
618:     PetscCall(PetscMalloc1(n1, &idx1));
619:     for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
620:     for (i = 0; i < n1; i++) idx1[i] = rst + i;
621:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]));
622:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]));
623:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]));
624:     PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]));
625:     PetscCall(MatIncreaseOverlap(A, 2, Ais, 3));
626:     PetscCall(MatIncreaseOverlap(B, 2, Bis, 3));
627:     /* Non deterministic output! */
628:     PetscCall(ISSort(Ais[0]));
629:     PetscCall(ISSort(Ais[1]));
630:     PetscCall(ISSort(Bis[0]));
631:     PetscCall(ISSort(Bis[1]));
632:     PetscCall(ISView(Ais[0], NULL));
633:     PetscCall(ISView(Bis[0], NULL));
634:     PetscCall(ISView(Ais[1], NULL));
635:     PetscCall(ISView(Bis[1], NULL));
636:     PetscCall(MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub));
637:     PetscCall(MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub));
638:     PetscCall(CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]"));
639:     PetscCall(CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]"));
640:     PetscCall(MatDestroySubMatrices(2, &Asub));
641:     PetscCall(MatDestroySubMatrices(2, &Bsub));
642:     PetscCall(ISDestroy(&Ais[0]));
643:     PetscCall(ISDestroy(&Ais[1]));
644:     PetscCall(ISDestroy(&Bis[0]));
645:     PetscCall(ISDestroy(&Bis[1]));
646:   }
647:   PetscCall(TestMatZeroRows(A, B, squaretest, is, 0.0));
648:   PetscCall(ISDestroy(&is));

650:   /* test MatTranspose */
651:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n"));
652:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
653:   PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &B2));
654:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose"));

656:   PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &A2));
657:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose"));
658:   PetscCall(MatDestroy(&A2));

660:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
661:   PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
662:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose"));
663:   PetscCall(MatDestroy(&A2));

665:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &A2));
666:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose"));
667:   PetscCall(MatDestroy(&A2));
668:   PetscCall(MatDestroy(&B2));

670:   /* test MatISFixLocalEmpty */
671:   if (isaij) {
672:     PetscInt r[2];

674:     r[0] = 0;
675:     r[1] = PetscMin(m, n) - 1;
676:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n"));
677:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));

679:     PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
680:     PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
681:     PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
682:     PetscCall(CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)"));

684:     PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
685:     PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
686:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
687:     PetscCall(MatZeroRows(B2, 2, r, 0.0, NULL, NULL));
688:     PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
689:     PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
690:     PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
691:     PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
692:     PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)"));
693:     PetscCall(MatDestroy(&A2));

695:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
696:     PetscCall(MatZeroRows(A2, 2, r, 0.0, NULL, NULL));
697:     PetscCall(MatTranspose(A2, MAT_INPLACE_MATRIX, &A2));
698:     PetscCall(MatTranspose(B2, MAT_INPLACE_MATRIX, &B2));
699:     PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
700:     PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
701:     PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
702:     PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
703:     PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
704:     PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)"));

706:     PetscCall(MatDestroy(&A2));
707:     PetscCall(MatDestroy(&B2));

709:     if (squaretest) {
710:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
711:       PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
712:       PetscCall(MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL));
713:       PetscCall(MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL));
714:       PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
715:       PetscCall(MatISFixLocalEmpty(A2, PETSC_TRUE));
716:       PetscCall(MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY));
717:       PetscCall(MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY));
718:       PetscCall(MatViewFromOptions(A2, NULL, "-fixempty_view"));
719:       PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)"));
720:       PetscCall(MatDestroy(&A2));
721:       PetscCall(MatDestroy(&B2));
722:     }
723:   }

725:   /* test MatInvertBlockDiagonal
726:        special cases for block-diagonal matrices */
727:   if (m == n) {
728:     ISLocalToGlobalMapping map;
729:     Mat                    Abd, Bbd;
730:     IS                     is, bis;
731:     const PetscScalar     *isbd, *aijbd;
732:     const PetscInt        *sts, *idxs;
733:     PetscInt              *idxs2, diff, perm, nl, bs, st, en, in;
734:     PetscBool              ok;

736:     for (diff = 0; diff < 3; diff++) {
737:       for (perm = 0; perm < 3; perm++) {
738:         for (bs = 1; bs < 4; bs++) {
739:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs));
740:           PetscCall(PetscMalloc1(bs * bs, &vals));
741:           PetscCall(MatGetOwnershipRanges(A, &sts));
742:           switch (diff) {
743:           case 1: /* inverted layout by processes */
744:             in = 1;
745:             st = sts[size - rank - 1];
746:             en = sts[size - rank];
747:             nl = en - st;
748:             break;
749:           case 2: /* round-robin layout */
750:             in = size;
751:             st = rank;
752:             nl = n / size;
753:             if (rank < n % size) nl++;
754:             break;
755:           default: /* same layout */
756:             in = 1;
757:             st = sts[rank];
758:             en = sts[rank + 1];
759:             nl = en - st;
760:             break;
761:           }
762:           PetscCall(ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is));
763:           PetscCall(ISGetLocalSize(is, &nl));
764:           PetscCall(ISGetIndices(is, &idxs));
765:           PetscCall(PetscMalloc1(nl, &idxs2));
766:           for (i = 0; i < nl; i++) {
767:             switch (perm) { /* invert some of the indices */
768:             case 2:
769:               idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
770:               break;
771:             case 1:
772:               idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
773:               break;
774:             default:
775:               idxs2[i] = idxs[i];
776:               break;
777:             }
778:           }
779:           PetscCall(ISRestoreIndices(is, &idxs));
780:           PetscCall(ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis));
781:           PetscCall(ISLocalToGlobalMappingCreateIS(bis, &map));
782:           PetscCall(MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd));
783:           PetscCall(ISLocalToGlobalMappingDestroy(&map));
784:           PetscCall(MatISSetPreallocation(Abd, bs, NULL, 0, NULL));
785:           for (i = 0; i < nl; i++) {
786:             PetscInt b1, b2;

788:             for (b1 = 0; b1 < bs; b1++) {
789:               for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
790:             }
791:             PetscCall(MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES));
792:           }
793:           PetscCall(MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY));
794:           PetscCall(MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY));
795:           PetscCall(MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd));
796:           PetscCall(MatInvertBlockDiagonal(Abd, &isbd));
797:           PetscCall(MatInvertBlockDiagonal(Bbd, &aijbd));
798:           PetscCall(MatGetLocalSize(Bbd, &nl, NULL));
799:           ok = PETSC_TRUE;
800:           for (i = 0; i < nl / bs; i++) {
801:             PetscInt b1, b2;

803:             for (b1 = 0; b1 < bs; b1++) {
804:               for (b2 = 0; b2 < bs; b2++) {
805:                 if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
806:                 if (!ok) {
807:                   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])));
808:                   break;
809:                 }
810:               }
811:               if (!ok) break;
812:             }
813:             if (!ok) break;
814:           }
815:           PetscCall(MatDestroy(&Abd));
816:           PetscCall(MatDestroy(&Bbd));
817:           PetscCall(PetscFree(vals));
818:           PetscCall(ISDestroy(&is));
819:           PetscCall(ISDestroy(&bis));
820:         }
821:       }
822:     }
823:   }

825:   /* test MatGetDiagonalBlock */
826:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n"));
827:   PetscCall(MatGetDiagonalBlock(A, &A2));
828:   PetscCall(MatGetDiagonalBlock(B, &B2));
829:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));
830:   PetscCall(MatScale(A, 2.0));
831:   PetscCall(MatScale(B, 2.0));
832:   PetscCall(MatGetDiagonalBlock(A, &A2));
833:   PetscCall(MatGetDiagonalBlock(B, &B2));
834:   PetscCall(CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock"));

836:   /* test MatISSetAllowRepeated on a MATIS */
837:   PetscCall(MatISSetAllowRepeated(A, allow_repeated));
838:   if (allow_repeated) { /* original MATIS maybe with repeated entries, test assembling of local matrices */
839:     Mat lA, lA2;

841:     for (PetscInt i = 0; i < 1; i++) { /* TODO: make MatScatter inherit from MATSHELL and support MatProducts */
842:       PetscBool usemult = PETSC_FALSE;

844:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
845:       if (i) {
846:         Mat tA;

848:         usemult = PETSC_TRUE;
849:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries and local shell matrices\n"));
850:         PetscCall(MatISGetLocalMat(A2, &lA2));
851:         PetscCall(MatConvert(lA2, MATSHELL, MAT_INITIAL_MATRIX, &tA));
852:         PetscCall(MatISRestoreLocalMat(A2, &lA2));
853:         PetscCall(MatISSetLocalMat(A2, tA));
854:         PetscCall(MatDestroy(&tA));
855:       } else {
856:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(false) with possibly repeated entries\n"));
857:       }
858:       PetscCall(MatISSetAllowRepeated(A2, PETSC_FALSE));
859:       PetscCall(MatISGetLocalMat(A, &lA));
860:       PetscCall(MatISGetLocalMat(A2, &lA2));
861:       if (!repmap) PetscCall(CheckMat(lA, lA2, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
862:       PetscCall(MatISRestoreLocalMat(A, &lA));
863:       PetscCall(MatISRestoreLocalMat(A2, &lA2));
864:       if (repmap) PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with repeated entries"));
865:       else PetscCall(CheckMat(A2, B, usemult, "MatISSetAllowRepeated(false) with non-repeated entries"));
866:       PetscCall(MatDestroy(&A2));
867:     }
868:   } else { /* original matis with non-repeated entries, this should only recreate the local matrices */
869:     Mat       lA;
870:     PetscBool flg;

872:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatISSetAllowRepeated(true) with non repeated entries\n"));
873:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A2));
874:     PetscCall(MatISSetAllowRepeated(A2, PETSC_TRUE));
875:     PetscCall(MatISGetLocalMat(A2, &lA));
876:     PetscCall(MatAssembled(lA, &flg));
877:     PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local mat should be unassembled");
878:     PetscCall(MatISRestoreLocalMat(A2, &lA));
879:     PetscCall(MatDestroy(&A2));
880:   }

882:   /* Test MatZeroEntries */
883:   PetscCall(MatZeroEntries(A));
884:   PetscCall(MatZeroEntries(B));
885:   PetscCall(CheckMat(A, B, PETSC_FALSE, "MatZeroEntries"));

887:   /* Test MatSetValues and MatSetValuesBlocked */
888:   if (test_setvalues) {
889:     PetscCall(PetscMalloc1(lm * ln, &vals));
890:     for (i = 0; i < lm * ln; i++) vals[i] = i + 1.0;
891:     PetscCall(MatGetLocalSize(A, NULL, &ln));
892:     PetscCall(MatISSetPreallocation(A, ln, NULL, n - ln, NULL));
893:     PetscCall(MatSeqAIJSetPreallocation(B, ln, NULL));
894:     PetscCall(MatMPIAIJSetPreallocation(B, ln, NULL, n - ln, NULL));
895:     PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
896:     PetscCall(ISLocalToGlobalMappingGetSize(cmap, &ln));

898:     PetscCall(ISLocalToGlobalMappingGetIndices(rmap, &idxs1));
899:     PetscCall(ISLocalToGlobalMappingGetIndices(cmap, &idxs2));
900:     PetscCall(MatSetValues(A, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
901:     PetscCall(MatSetValues(B, lm, idxs1, ln, idxs2, vals, ADD_VALUES));
902:     PetscCall(ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1));
903:     PetscCall(ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2));
904:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
905:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
906:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
907:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
908:     PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValues"));

910:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(rmap, &idxs1));
911:     PetscCall(ISLocalToGlobalMappingGetBlockIndices(cmap, &idxs2));
912:     PetscCall(MatSetValuesBlocked(A, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
913:     PetscCall(MatSetValuesBlocked(B, lm / rbs, idxs1, ln / cbs, idxs2, vals, ADD_VALUES));
914:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rmap, &idxs1));
915:     PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cmap, &idxs2));
916:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
917:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
918:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
919:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
920:     PetscCall(CheckMat(A, B, PETSC_FALSE, "MatSetValuesBlocked"));
921:     PetscCall(PetscFree(vals));
922:   }

924:   /* free testing matrices */
925:   PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
926:   PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
927:   PetscCall(MatDestroy(&A));
928:   PetscCall(MatDestroy(&B));
929:   PetscCall(PetscFinalize());
930:   return 0;
931: }

933: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
934: {
935:   Mat       Bcheck;
936:   PetscReal error;

938:   PetscFunctionBeginUser;
939:   if (!usemult && B) {
940:     PetscBool hasnorm;

942:     PetscCall(MatHasOperation(B, MATOP_NORM, &hasnorm));
943:     if (!hasnorm) usemult = PETSC_TRUE;
944:   }
945:   if (!usemult) {
946:     if (B) {
947:       MatType Btype;

949:       PetscCall(MatGetType(B, &Btype));
950:       PetscCall(MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck));
951:     } else {
952:       PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
953:     }
954:     if (B) { /* if B is present, subtract it */
955:       PetscCall(MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN));
956:     }
957:     PetscCall(MatNorm(Bcheck, NORM_INFINITY, &error));
958:     if (error > PETSC_SQRT_MACHINE_EPSILON) {
959:       ISLocalToGlobalMapping rl2g, cl2g;

961:       PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Bcheck"));
962:       PetscCall(MatView(Bcheck, NULL));
963:       if (B) {
964:         PetscCall(PetscObjectSetName((PetscObject)B, "B"));
965:         PetscCall(MatView(B, NULL));
966:         PetscCall(MatDestroy(&Bcheck));
967:         PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck));
968:         PetscCall(PetscObjectSetName((PetscObject)Bcheck, "Assembled A"));
969:         PetscCall(MatView(Bcheck, NULL));
970:       }
971:       PetscCall(MatDestroy(&Bcheck));
972:       PetscCall(PetscObjectSetName((PetscObject)A, "A"));
973:       PetscCall(MatView(A, NULL));
974:       PetscCall(MatGetLocalToGlobalMapping(A, &rl2g, &cl2g));
975:       if (rl2g) PetscCall(ISLocalToGlobalMappingView(rl2g, NULL));
976:       if (cl2g) PetscCall(ISLocalToGlobalMappingView(cl2g, NULL));
977:       SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
978:     }
979:     PetscCall(MatDestroy(&Bcheck));
980:   } else {
981:     PetscBool ok, okt;

983:     PetscCall(MatMultEqual(A, B, 3, &ok));
984:     PetscCall(MatMultTransposeEqual(A, B, 3, &okt));
985:     PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ?  %d, multtranspose ok ? %d", func, ok, okt);
986:   }
987:   PetscFunctionReturn(PETSC_SUCCESS);
988: }

990: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
991: {
992:   Mat                    B, Bcheck, B2 = NULL, lB;
993:   Vec                    x = NULL, b = NULL, b2 = NULL;
994:   ISLocalToGlobalMapping l2gr, l2gc;
995:   PetscReal              error;
996:   char                   diagstr[16];
997:   const PetscInt        *idxs;
998:   PetscInt               rst, ren, i, n, N, d;
999:   PetscMPIInt            rank;
1000:   PetscBool              miss, haszerorows;

1002:   PetscFunctionBeginUser;
1003:   if (diag == 0.) {
1004:     PetscCall(PetscStrncpy(diagstr, "zero", sizeof(diagstr)));
1005:   } else {
1006:     PetscCall(PetscStrncpy(diagstr, "nonzero", sizeof(diagstr)));
1007:   }
1008:   PetscCall(ISView(is, NULL));
1009:   PetscCall(MatGetLocalToGlobalMapping(A, &l2gr, &l2gc));
1010:   /* tests MatDuplicate and MatCopy */
1011:   if (diag == 0.) {
1012:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
1013:   } else {
1014:     PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
1015:     PetscCall(MatCopy(A, B, SAME_NONZERO_PATTERN));
1016:   }
1017:   PetscCall(MatISGetLocalMat(B, &lB));
1018:   PetscCall(MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows));
1019:   if (squaretest && haszerorows) {
1020:     PetscCall(MatCreateVecs(B, &x, &b));
1021:     PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B2));
1022:     PetscCall(VecSetLocalToGlobalMapping(b, l2gr));
1023:     PetscCall(VecSetLocalToGlobalMapping(x, l2gc));
1024:     PetscCall(VecSetRandom(x, NULL));
1025:     PetscCall(VecSetRandom(b, NULL));
1026:     /* mimic b[is] = x[is] */
1027:     PetscCall(VecDuplicate(b, &b2));
1028:     PetscCall(VecSetLocalToGlobalMapping(b2, l2gr));
1029:     PetscCall(VecCopy(b, b2));
1030:     PetscCall(ISGetLocalSize(is, &n));
1031:     PetscCall(ISGetIndices(is, &idxs));
1032:     PetscCall(VecGetSize(x, &N));
1033:     for (i = 0; i < n; i++) {
1034:       if (0 <= idxs[i] && idxs[i] < N) {
1035:         PetscCall(VecSetValue(b2, idxs[i], diag, INSERT_VALUES));
1036:         PetscCall(VecSetValue(x, idxs[i], 1., INSERT_VALUES));
1037:       }
1038:     }
1039:     PetscCall(VecAssemblyBegin(b2));
1040:     PetscCall(VecAssemblyEnd(b2));
1041:     PetscCall(VecAssemblyBegin(x));
1042:     PetscCall(VecAssemblyEnd(x));
1043:     PetscCall(ISRestoreIndices(is, &idxs));
1044:     /*  test ZeroRows on MATIS */
1045:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1046:     PetscCall(MatZeroRowsIS(B, is, diag, x, b));
1047:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr));
1048:     PetscCall(MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL));
1049:   } else if (haszerorows) {
1050:     /*  test ZeroRows on MATIS */
1051:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr));
1052:     PetscCall(MatZeroRowsIS(B, is, diag, NULL, NULL));
1053:     b = b2 = x = NULL;
1054:   } else {
1055:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr));
1056:     b = b2 = x = NULL;
1057:   }

1059:   if (squaretest && haszerorows) {
1060:     PetscCall(VecAXPY(b2, -1., b));
1061:     PetscCall(VecNorm(b2, NORM_INFINITY, &error));
1062:     PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr);
1063:   }

1065:   /* test MatMissingDiagonal */
1066:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n"));
1067:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1068:   PetscCall(MatMissingDiagonal(B, &miss, &d));
1069:   PetscCall(MatGetOwnershipRange(B, &rst, &ren));
1070:   PetscCall(PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD));
1071:   PetscCall(PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n", rank, rst, ren, (int)miss, d, diagstr));
1072:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
1073:   PetscCall(PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD));

1075:   PetscCall(VecDestroy(&x));
1076:   PetscCall(VecDestroy(&b));
1077:   PetscCall(VecDestroy(&b2));

1079:   /* check the result of ZeroRows with that from MPIAIJ routines
1080:      assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
1081:   if (haszerorows) {
1082:     PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck));
1083:     PetscCall(MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1084:     PetscCall(MatZeroRowsIS(Bcheck, is, diag, NULL, NULL));
1085:     PetscCall(CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows"));
1086:     PetscCall(MatDestroy(&Bcheck));
1087:   }
1088:   PetscCall(MatDestroy(&B));

1090:   if (B2) { /* test MatZeroRowsColumns */
1091:     PetscCall(MatDuplicate(Afull, MAT_COPY_VALUES, &B));
1092:     PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
1093:     PetscCall(MatZeroRowsColumnsIS(B, is, diag, NULL, NULL));
1094:     PetscCall(CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns"));
1095:     PetscCall(MatDestroy(&B));
1096:     PetscCall(MatDestroy(&B2));
1097:   }
1098:   PetscFunctionReturn(PETSC_SUCCESS);
1099: }

1101: /*TEST

1103:    test:
1104:       requires: !complex
1105:       args: -test_matlab -test_trans

1107:    test:
1108:       suffix: 2
1109:       nsize: 4
1110:       args: -mat_is_convert_local_nest -nr 3 -nc 4

1112:    test:
1113:       suffix: 3
1114:       nsize: 5
1115:       args: -m 11 -n 10 -mat_is_convert_local_nest -nr 2 -nc 1 -cbs 2

1117:    test:
1118:       suffix: 4
1119:       nsize: 6
1120:       args: -m 9 -n 12 -test_trans -nr 2 -nc 7

1122:    test:
1123:       suffix: 5
1124:       nsize: 6
1125:       args: -m 12 -n 12 -test_trans -nr 3 -nc 1 -rbs 2

1127:    test:
1128:       suffix: 6
1129:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -rbs 6 -cbs 3

1131:    test:
1132:       suffix: 7
1133:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

1135:    test:
1136:       suffix: 8
1137:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap

1139:    test:
1140:       suffix: 9
1141:       nsize: 5
1142:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap

1144:    test:
1145:       suffix: 10
1146:       nsize: 5
1147:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

1149:    test:
1150:       suffix: vscat_default
1151:       nsize: 5
1152:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
1153:       output_file: output/ex23_11.out

1155:    test:
1156:       suffix: 12
1157:       nsize: 3
1158:       args: -m 12 -n 12 -symmetric -mat_is_localmat_type sbaij -test_trans -nr 2 -nc 3 -test_setvalues 0

1160:    testset:
1161:       output_file: output/ex23_13.out
1162:       nsize: 3
1163:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
1164:       filter: grep -v "type:" | grep -v "block size is 1" | grep -v "not using I-node routines"
1165:       test:
1166:         suffix: baij
1167:         args: -mat_is_localmat_type baij
1168:       test:
1169:         requires: viennacl
1170:         suffix: viennacl
1171:         args: -mat_is_localmat_type aijviennacl
1172:       test:
1173:         requires: cuda
1174:         suffix: cusparse
1175:         args: -mat_is_localmat_type aijcusparse
1176:       test:
1177:         requires: kokkos_kernels
1178:         suffix: kokkos
1179:         args: -mat_is_localmat_type aijkokkos

1181:    test:
1182:       suffix: negrep
1183:       nsize: {{1 3}separate output}
1184:       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

1186:    test:
1187:       suffix: negrep_allowrep
1188:       nsize: {{1 3}separate output}
1189:       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

1191: TEST*/