Actual source code: ex23.c


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

  4: #include <petscmat.h>

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

  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, is2, reven, rodd, ceven, codd;
 18:   IS                    *rows, *cols;
 19:   IS                     irow[2], icol[2];
 20:   PetscLayout            rlayout, clayout;
 21:   const PetscInt        *rrange, *crange;
 22:   MatType                lmtype;
 23:   PetscScalar            diag = 2.;
 24:   PetscInt               n, m, i, lm, ln;
 25:   PetscInt               rst, ren, cst, cen, nr, nc;
 26:   PetscMPIInt            rank, size, lrank, rrank;
 27:   PetscBool              testT, squaretest, isaij;
 28:   PetscBool              permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE;
 29:   PetscBool              diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric;

 32:   PetscInitialize(&argc, &args, (char *)0, help);
 33:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 34:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 35:   m = n = 2 * size;
 36:   PetscOptionsGetBool(NULL, NULL, "-symmetric", &symmetric, NULL);
 37:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 38:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
 39:   PetscOptionsGetBool(NULL, NULL, "-negmap", &negmap, NULL);
 40:   PetscOptionsGetBool(NULL, NULL, "-repmap", &repmap, NULL);
 41:   PetscOptionsGetBool(NULL, NULL, "-permmap", &permute, NULL);
 42:   PetscOptionsGetBool(NULL, NULL, "-diffmap", &diffmap, NULL);

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

 62:     ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &isl[0]);
 63:     ISCreateStride(PETSC_COMM_WORLD, n, n - 1, -1, &isl[1]);
 64:     ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is);
 65:     ISDestroy(&isl[0]);
 66:     ISDestroy(&isl[1]);
 67:   } else { /* negative and repeated indices */
 68:     IS isl[2];

 70:     ISCreateStride(PETSC_COMM_WORLD, n + 1, -1, 1, &isl[0]);
 71:     ISCreateStride(PETSC_COMM_WORLD, n + 1, n - 1, -1, &isl[1]);
 72:     ISConcatenate(PETSC_COMM_WORLD, 2, isl, &is);
 73:     ISDestroy(&isl[0]);
 74:     ISDestroy(&isl[1]);
 75:   }
 76:   ISLocalToGlobalMappingCreateIS(is, &cmap);
 77:   ISDestroy(&is);

 79:   if (m != n || diffmap) {
 80:     ISCreateStride(PETSC_COMM_WORLD, m, permute ? m - 1 : 0, permute ? -1 : 1, &is);
 81:     ISLocalToGlobalMappingCreateIS(is, &rmap);
 82:     ISDestroy(&is);
 83:   } else {
 84:     PetscObjectReference((PetscObject)cmap);
 85:     rmap = cmap;
 86:   }

 88:   MatSetLocalToGlobalMapping(A, rmap, cmap);
 89:   MatISStoreL2L(A, PETSC_FALSE);
 90:   MatISSetPreallocation(A, 3, NULL, 3, NULL);
 91:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap)); /* I do not want to precompute the pattern */
 92:   ISLocalToGlobalMappingGetSize(rmap, &lm);
 93:   ISLocalToGlobalMappingGetSize(cmap, &ln);
 94:   for (i = 0; i < lm; i++) {
 95:     PetscScalar v[3];
 96:     PetscInt    cols[3];

 98:     cols[0] = (i - 1 + n) % n;
 99:     cols[1] = i % n;
100:     cols[2] = (i + 1) % n;
101:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
102:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
103:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
104:     ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
105:     MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES);
106:   }
107:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
108:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

110:   /* activate tests for square matrices with same maps only */
111:   MatHasCongruentLayouts(A, &squaretest);
112:   if (squaretest && rmap != cmap) {
113:     PetscInt nr, nc;

115:     ISLocalToGlobalMappingGetSize(rmap, &nr);
116:     ISLocalToGlobalMappingGetSize(cmap, &nc);
117:     if (nr != nc) squaretest = PETSC_FALSE;
118:     else {
119:       const PetscInt *idxs1, *idxs2;

121:       ISLocalToGlobalMappingGetIndices(rmap, &idxs1);
122:       ISLocalToGlobalMappingGetIndices(cmap, &idxs2);
123:       PetscArraycmp(idxs1, idxs2, nr, &squaretest);
124:       ISLocalToGlobalMappingRestoreIndices(rmap, &idxs1);
125:       ISLocalToGlobalMappingRestoreIndices(cmap, &idxs2);
126:     }
127:     MPIU_Allreduce(MPI_IN_PLACE, &squaretest, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A));
128:   }

130:   /* test MatISGetLocalMat */
131:   MatISGetLocalMat(A, &B);
132:   MatGetType(B, &lmtype);

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

148:   /* test MatIsSymmetric */
149:   MatIsSymmetric(A, 0.0, &issymmetric);
150:   PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric);

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

170:     cols[0] = (i - 1 + n) % n;
171:     cols[1] = i % n;
172:     cols[2] = (i + 1) % n;
173:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
174:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
175:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
176:     ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
177:     MatSetValuesLocal(B, 1, &i, 3, cols, v, ADD_VALUES);
178:   }
179:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
180:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

182:   /* test MatView */
183:   PetscPrintf(PETSC_COMM_WORLD, "Test MatView\n");
184:   MatView(A, NULL);
185:   MatView(B, NULL);

187:   /* test CheckMat */
188:   PetscPrintf(PETSC_COMM_WORLD, "Test CheckMat\n");
189:   CheckMat(A, B, PETSC_FALSE, "CheckMat");

191:   /* test MatDuplicate and MatAXPY */
192:   PetscPrintf(PETSC_COMM_WORLD, "Test MatDuplicate and MatAXPY\n");
193:   MatDuplicate(A, MAT_COPY_VALUES, &A2);
194:   CheckMat(A, A2, PETSC_FALSE, "MatDuplicate and MatAXPY");

196:   /* test MatConvert */
197:   PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_IS_XAIJ\n");
198:   MatConvert(A2, MATAIJ, MAT_INITIAL_MATRIX, &B2);
199:   CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INITIAL_MATRIX");
200:   MatConvert(A2, MATAIJ, MAT_REUSE_MATRIX, &B2);
201:   CheckMat(B, B2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_REUSE_MATRIX");
202:   MatConvert(A2, MATAIJ, MAT_INPLACE_MATRIX, &A2);
203:   CheckMat(B, A2, PETSC_TRUE, "MatConvert_IS_XAIJ MAT_INPLACE_MATRIX");
204:   MatDestroy(&A2);
205:   MatDestroy(&B2);
206:   PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_XAIJ_IS\n");
207:   MatDuplicate(B, MAT_COPY_VALUES, &B2);
208:   MatConvert(B2, MATIS, MAT_INITIAL_MATRIX, &A2);
209:   CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INITIAL_MATRIX");
210:   MatConvert(B2, MATIS, MAT_REUSE_MATRIX, &A2);
211:   CheckMat(A, A2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_REUSE_MATRIX");
212:   MatConvert(B2, MATIS, MAT_INPLACE_MATRIX, &B2);
213:   CheckMat(A, B2, PETSC_TRUE, "MatConvert_XAIJ_IS MAT_INPLACE_MATRIX");
214:   MatDestroy(&A2);
215:   MatDestroy(&B2);
216:   PetscStrcmp(lmtype, MATSEQAIJ, &isaij);
217:   if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
218:     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};
219:     ISLocalToGlobalMapping tcmap, trmap;

221:     for (ri = 0; ri < 2; ri++) {
222:       PetscInt *r;

224:       r = (PetscInt *)(ri == 0 ? rr : rk);
225:       for (ci = 0; ci < 2; ci++) {
226:         PetscInt *c, rb, cb;

228:         c = (PetscInt *)(ci == 0 ? cr : ck);
229:         for (rb = 1; rb < 4; rb++) {
230:           ISCreateBlock(PETSC_COMM_SELF, rb, 3, r, PETSC_COPY_VALUES, &is);
231:           ISLocalToGlobalMappingCreateIS(is, &trmap);
232:           ISDestroy(&is);
233:           for (cb = 1; cb < 4; cb++) {
234:             Mat  T, lT, T2;
235:             char testname[256];

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

240:             ISCreateBlock(PETSC_COMM_SELF, cb, 4, c, PETSC_COPY_VALUES, &is);
241:             ISLocalToGlobalMappingCreateIS(is, &tcmap);
242:             ISDestroy(&is);

244:             MatCreate(PETSC_COMM_SELF, &T);
245:             MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, rb * 3, cb * 4);
246:             MatSetType(T, MATIS);
247:             MatSetLocalToGlobalMapping(T, trmap, tcmap);
248:             ISLocalToGlobalMappingDestroy(&tcmap);
249:             MatISGetLocalMat(T, &lT);
250:             MatSetType(lT, MATSEQAIJ);
251:             MatSeqAIJSetPreallocation(lT, cb * 4, NULL);
252:             MatSetRandom(lT, NULL);
253:             MatConvert(lT, lmtype, MAT_INPLACE_MATRIX, &lT);
254:             MatISRestoreLocalMat(T, &lT);
255:             MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
256:             MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);

258:             MatConvert(T, MATAIJ, MAT_INITIAL_MATRIX, &T2);
259:             CheckMat(T, T2, PETSC_TRUE, "MAT_INITIAL_MATRIX");
260:             MatConvert(T, MATAIJ, MAT_REUSE_MATRIX, &T2);
261:             CheckMat(T, T2, PETSC_TRUE, "MAT_REUSE_MATRIX");
262:             MatDestroy(&T2);
263:             MatDuplicate(T, MAT_COPY_VALUES, &T2);
264:             MatConvert(T2, MATAIJ, MAT_INPLACE_MATRIX, &T2);
265:             CheckMat(T, T2, PETSC_TRUE, "MAT_INPLACE_MATRIX");
266:             MatDestroy(&T);
267:             MatDestroy(&T2);
268:           }
269:           ISLocalToGlobalMappingDestroy(&trmap);
270:         }
271:       }
272:     }
273:   }

275:   /* test MatDiagonalScale */
276:   PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalScale\n");
277:   MatDuplicate(A, MAT_COPY_VALUES, &A2);
278:   MatDuplicate(B, MAT_COPY_VALUES, &B2);
279:   MatCreateVecs(A, &x, &y);
280:   VecSetRandom(x, NULL);
281:   if (issymmetric) {
282:     VecCopy(x, y);
283:   } else {
284:     VecSetRandom(y, NULL);
285:     VecScale(y, 8.);
286:   }
287:   MatDiagonalScale(A2, y, x);
288:   MatDiagonalScale(B2, y, x);
289:   CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalScale");
290:   MatDestroy(&A2);
291:   MatDestroy(&B2);
292:   VecDestroy(&x);
293:   VecDestroy(&y);

295:   /* test MatPtAP (A IS and B AIJ) */
296:   if (isaij && m == n) {
297:     PetscPrintf(PETSC_COMM_WORLD, "Test MatPtAP\n");
298:     MatISStoreL2L(A, PETSC_TRUE);
299:     MatPtAP(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &A2);
300:     MatPtAP(B, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B2);
301:     CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_INITIAL_MATRIX");
302:     MatPtAP(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &A2);
303:     CheckMat(A2, B2, PETSC_FALSE, "MatPtAP MAT_REUSE_MATRIX");
304:     MatDestroy(&A2);
305:     MatDestroy(&B2);
306:   }

308:   /* test MatGetLocalSubMatrix */
309:   PetscPrintf(PETSC_COMM_WORLD, "Test MatGetLocalSubMatrix\n");
310:   MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A2);
311:   ISCreateStride(PETSC_COMM_SELF, lm / 2 + lm % 2, 0, 2, &reven);
312:   ISComplement(reven, 0, lm, &rodd);
313:   ISCreateStride(PETSC_COMM_SELF, ln / 2 + ln % 2, 0, 2, &ceven);
314:   ISComplement(ceven, 0, ln, &codd);
315:   MatGetLocalSubMatrix(A2, reven, ceven, &Aee);
316:   MatGetLocalSubMatrix(A2, reven, codd, &Aeo);
317:   MatGetLocalSubMatrix(A2, rodd, ceven, &Aoe);
318:   MatGetLocalSubMatrix(A2, rodd, codd, &Aoo);
319:   for (i = 0; i < lm; i++) {
320:     PetscInt    j, je, jo, colse[3], colso[3];
321:     PetscScalar ve[3], vo[3];
322:     PetscScalar v[3];
323:     PetscInt    cols[3];
324:     PetscInt    row = i / 2;

326:     cols[0] = (i - 1 + n) % n;
327:     cols[1] = i % n;
328:     cols[2] = (i + 1) % n;
329:     v[0]    = -1. * (symmetric ? PetscMin(i + 1, cols[0] + 1) : i + 1);
330:     v[1]    = 2. * (symmetric ? PetscMin(i + 1, cols[1] + 1) : i + 1);
331:     v[2]    = -1. * (symmetric ? PetscMin(i + 1, cols[2] + 1) : i + 1);
332:     ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols);
333:     for (j = 0, je = 0, jo = 0; j < 3; j++) {
334:       if (cols[j] % 2) {
335:         vo[jo]      = v[j];
336:         colso[jo++] = cols[j] / 2;
337:       } else {
338:         ve[je]      = v[j];
339:         colse[je++] = cols[j] / 2;
340:       }
341:     }
342:     if (i % 2) {
343:       MatSetValuesLocal(Aoe, 1, &row, je, colse, ve, ADD_VALUES);
344:       MatSetValuesLocal(Aoo, 1, &row, jo, colso, vo, ADD_VALUES);
345:     } else {
346:       MatSetValuesLocal(Aee, 1, &row, je, colse, ve, ADD_VALUES);
347:       MatSetValuesLocal(Aeo, 1, &row, jo, colso, vo, ADD_VALUES);
348:     }
349:   }
350:   MatRestoreLocalSubMatrix(A2, reven, ceven, &Aee);
351:   MatRestoreLocalSubMatrix(A2, reven, codd, &Aeo);
352:   MatRestoreLocalSubMatrix(A2, rodd, ceven, &Aoe);
353:   MatRestoreLocalSubMatrix(A2, rodd, codd, &Aoo);
354:   ISDestroy(&reven);
355:   ISDestroy(&ceven);
356:   ISDestroy(&rodd);
357:   ISDestroy(&codd);
358:   MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
359:   MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
360:   MatAXPY(A2, -1., A, SAME_NONZERO_PATTERN);
361:   CheckMat(A2, NULL, PETSC_FALSE, "MatGetLocalSubMatrix");
362:   MatDestroy(&A2);

364:   /* test MatConvert_Nest_IS */
365:   testT = PETSC_FALSE;
366:   PetscOptionsGetBool(NULL, NULL, "-test_trans", &testT, NULL);

368:   PetscPrintf(PETSC_COMM_WORLD, "Test MatConvert_Nest_IS\n");
369:   nr = 2;
370:   nc = 2;
371:   PetscOptionsGetInt(NULL, NULL, "-nr", &nr, NULL);
372:   PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL);
373:   if (testT) {
374:     MatGetOwnershipRange(A, &cst, &cen);
375:     MatGetOwnershipRangeColumn(A, &rst, &ren);
376:   } else {
377:     MatGetOwnershipRange(A, &rst, &ren);
378:     MatGetOwnershipRangeColumn(A, &cst, &cen);
379:   }
380:   PetscMalloc3(nr, &rows, nc, &cols, 2 * nr * nc, &mats);
381:   for (i = 0; i < nr * nc; i++) {
382:     if (testT) {
383:       MatCreateTranspose(A, &mats[i]);
384:       MatTranspose(B, MAT_INITIAL_MATRIX, &mats[i + nr * nc]);
385:     } else {
386:       MatDuplicate(A, MAT_COPY_VALUES, &mats[i]);
387:       MatDuplicate(B, MAT_COPY_VALUES, &mats[i + nr * nc]);
388:     }
389:   }
390:   for (i = 0; i < nr; i++) ISCreateStride(PETSC_COMM_WORLD, ren - rst, i + rst, nr, &rows[i]);
391:   for (i = 0; i < nc; i++) ISCreateStride(PETSC_COMM_WORLD, cen - cst, i + cst, nc, &cols[i]);
392:   MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats, &A2);
393:   MatCreateNest(PETSC_COMM_WORLD, nr, rows, nc, cols, mats + nr * nc, &B2);
394:   for (i = 0; i < nr; i++) ISDestroy(&rows[i]);
395:   for (i = 0; i < nc; i++) ISDestroy(&cols[i]);
396:   for (i = 0; i < 2 * nr * nc; i++) MatDestroy(&mats[i]);
397:   PetscFree3(rows, cols, mats);
398:   MatConvert(B2, MATAIJ, MAT_INITIAL_MATRIX, &T);
399:   MatDestroy(&B2);
400:   MatConvert(A2, MATIS, MAT_INITIAL_MATRIX, &B2);
401:   CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INITIAL_MATRIX");
402:   MatConvert(A2, MATIS, MAT_REUSE_MATRIX, &B2);
403:   CheckMat(B2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_REUSE_MATRIX");
404:   MatDestroy(&B2);
405:   MatConvert(A2, MATIS, MAT_INPLACE_MATRIX, &A2);
406:   CheckMat(A2, T, PETSC_TRUE, "MatConvert_Nest_IS MAT_INPLACE_MATRIX");
407:   MatDestroy(&T);
408:   MatDestroy(&A2);

410:   /* test MatCreateSubMatrix */
411:   PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrix\n");
412:   if (rank == 0) {
413:     ISCreateStride(PETSC_COMM_WORLD, 1, 1, 1, &is);
414:     ISCreateStride(PETSC_COMM_WORLD, 2, 0, 1, &is2);
415:   } else if (rank == 1) {
416:     ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is);
417:     if (n > 3) {
418:       ISCreateStride(PETSC_COMM_WORLD, 1, 3, 1, &is2);
419:     } else {
420:       ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2);
421:     }
422:   } else if (rank == 2 && n > 4) {
423:     ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
424:     ISCreateStride(PETSC_COMM_WORLD, n - 4, 4, 1, &is2);
425:   } else {
426:     ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
427:     ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is2);
428:   }
429:   MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &A2);
430:   MatCreateSubMatrix(B, is, is, MAT_INITIAL_MATRIX, &B2);
431:   CheckMat(A2, B2, PETSC_TRUE, "first MatCreateSubMatrix");

433:   MatCreateSubMatrix(A, is, is, MAT_REUSE_MATRIX, &A2);
434:   MatCreateSubMatrix(B, is, is, MAT_REUSE_MATRIX, &B2);
435:   CheckMat(A2, B2, PETSC_FALSE, "reuse MatCreateSubMatrix");
436:   MatDestroy(&A2);
437:   MatDestroy(&B2);

439:   if (!issymmetric) {
440:     MatCreateSubMatrix(A, is, is2, MAT_INITIAL_MATRIX, &A2);
441:     MatCreateSubMatrix(B, is, is2, MAT_INITIAL_MATRIX, &B2);
442:     MatCreateSubMatrix(A, is, is2, MAT_REUSE_MATRIX, &A2);
443:     MatCreateSubMatrix(B, is, is2, MAT_REUSE_MATRIX, &B2);
444:     CheckMat(A2, B2, PETSC_FALSE, "second MatCreateSubMatrix");
445:   }

447:   MatDestroy(&A2);
448:   MatDestroy(&B2);
449:   ISDestroy(&is);
450:   ISDestroy(&is2);

452:   /* test MatCreateSubMatrices */
453:   PetscPrintf(PETSC_COMM_WORLD, "Test MatCreateSubMatrices\n");
454:   MatGetLayouts(A, &rlayout, &clayout);
455:   PetscLayoutGetRanges(rlayout, &rrange);
456:   PetscLayoutGetRanges(clayout, &crange);
457:   lrank = (size + rank - 1) % size;
458:   rrank = (rank + 1) % size;
459:   ISCreateStride(PETSC_COMM_SELF, (rrange[lrank + 1] - rrange[lrank]), rrange[lrank], 1, &irow[0]);
460:   ISCreateStride(PETSC_COMM_SELF, (crange[rrank + 1] - crange[rrank]), crange[rrank], 1, &icol[0]);
461:   ISCreateStride(PETSC_COMM_SELF, (rrange[rrank + 1] - rrange[rrank]), rrange[rrank], 1, &irow[1]);
462:   ISCreateStride(PETSC_COMM_SELF, (crange[lrank + 1] - crange[lrank]), crange[lrank], 1, &icol[1]);
463:   MatCreateSubMatrices(A, 2, irow, icol, MAT_INITIAL_MATRIX, &Asub);
464:   MatCreateSubMatrices(B, 2, irow, icol, MAT_INITIAL_MATRIX, &Bsub);
465:   CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]");
466:   CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]");
467:   MatCreateSubMatrices(A, 2, irow, icol, MAT_REUSE_MATRIX, &Asub);
468:   MatCreateSubMatrices(B, 2, irow, icol, MAT_REUSE_MATRIX, &Bsub);
469:   CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatCreateSubMatrices[0]");
470:   CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatCreateSubMatrices[1]");
471:   MatDestroySubMatrices(2, &Asub);
472:   MatDestroySubMatrices(2, &Bsub);
473:   ISDestroy(&irow[0]);
474:   ISDestroy(&irow[1]);
475:   ISDestroy(&icol[0]);
476:   ISDestroy(&icol[1]);

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

483:       st  = (m + 1) / 2;
484:       len = PetscMin(m / 2, PetscMax(m - (m + 1) / 2 - 1, 0));
485:       ISCreateStride(PETSC_COMM_WORLD, len, st, 1, &is);
486:     } else {
487:       ISCreateStride(PETSC_COMM_WORLD, 0, 0, 1, &is);
488:     }
489:   } else {
490:     ISCreateStride(PETSC_COMM_WORLD, 1, 0, 1, &is);
491:   }

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

497:     /* test MatDiagonalSet */
498:     PetscPrintf(PETSC_COMM_WORLD, "Test MatDiagonalSet\n");
499:     MatDuplicate(A, MAT_COPY_VALUES, &A2);
500:     MatDuplicate(B, MAT_COPY_VALUES, &B2);
501:     MatCreateVecs(A, NULL, &x);
502:     VecSetRandom(x, NULL);
503:     MatDiagonalSet(A2, x, INSERT_VALUES);
504:     MatDiagonalSet(B2, x, INSERT_VALUES);
505:     CheckMat(A2, B2, PETSC_FALSE, "MatDiagonalSet");
506:     VecDestroy(&x);
507:     MatDestroy(&A2);
508:     MatDestroy(&B2);

510:     /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
511:     PetscPrintf(PETSC_COMM_WORLD, "Test MatShift\n");
512:     MatDuplicate(A, MAT_COPY_VALUES, &A2);
513:     MatDuplicate(B, MAT_COPY_VALUES, &B2);
514:     MatShift(A2, 2.0);
515:     MatShift(B2, 2.0);
516:     CheckMat(A2, B2, PETSC_FALSE, "MatShift");
517:     MatDestroy(&A2);
518:     MatDestroy(&B2);

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

523:     /* test MatIncreaseOverlap */
524:     PetscPrintf(PETSC_COMM_WORLD, "Test MatIncreaseOverlap\n");
525:     MatGetOwnershipRange(A, &rst, &ren);
526:     n0 = (ren - rst) / 2;
527:     n1 = (ren - rst) / 3;
528:     PetscMalloc1(n0, &idx0);
529:     PetscMalloc1(n1, &idx1);
530:     for (i = 0; i < n0; i++) idx0[i] = ren - i - 1;
531:     for (i = 0; i < n1; i++) idx1[i] = rst + i;
532:     ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_OWN_POINTER, &Ais[0]);
533:     ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_OWN_POINTER, &Ais[1]);
534:     ISCreateGeneral(PETSC_COMM_WORLD, n0, idx0, PETSC_COPY_VALUES, &Bis[0]);
535:     ISCreateGeneral(PETSC_COMM_WORLD, n1, idx1, PETSC_COPY_VALUES, &Bis[1]);
536:     MatIncreaseOverlap(A, 2, Ais, 3);
537:     MatIncreaseOverlap(B, 2, Bis, 3);
538:     /* Non deterministic output! */
539:     ISSort(Ais[0]);
540:     ISSort(Ais[1]);
541:     ISSort(Bis[0]);
542:     ISSort(Bis[1]);
543:     ISView(Ais[0], NULL);
544:     ISView(Bis[0], NULL);
545:     ISView(Ais[1], NULL);
546:     ISView(Bis[1], NULL);
547:     MatCreateSubMatrices(A, 2, Ais, Ais, MAT_INITIAL_MATRIX, &Asub);
548:     MatCreateSubMatrices(B, 2, Bis, Ais, MAT_INITIAL_MATRIX, &Bsub);
549:     CheckMat(Asub[0], Bsub[0], PETSC_FALSE, "MatIncreaseOverlap[0]");
550:     CheckMat(Asub[1], Bsub[1], PETSC_FALSE, "MatIncreaseOverlap[1]");
551:     MatDestroySubMatrices(2, &Asub);
552:     MatDestroySubMatrices(2, &Bsub);
553:     ISDestroy(&Ais[0]);
554:     ISDestroy(&Ais[1]);
555:     ISDestroy(&Bis[0]);
556:     ISDestroy(&Bis[1]);
557:   }
558:   TestMatZeroRows(A, B, squaretest, is, 0.0);
559:   ISDestroy(&is);

561:   /* test MatTranspose */
562:   PetscPrintf(PETSC_COMM_WORLD, "Test MatTranspose\n");
563:   MatTranspose(A, MAT_INITIAL_MATRIX, &A2);
564:   MatTranspose(B, MAT_INITIAL_MATRIX, &B2);
565:   CheckMat(A2, B2, PETSC_FALSE, "initial matrix MatTranspose");

567:   MatTranspose(A, MAT_REUSE_MATRIX, &A2);
568:   CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (not in place) MatTranspose");
569:   MatDestroy(&A2);

571:   MatDuplicate(A, MAT_COPY_VALUES, &A2);
572:   MatTranspose(A2, MAT_INPLACE_MATRIX, &A2);
573:   CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (in place) MatTranspose");
574:   MatDestroy(&A2);

576:   MatTranspose(A, MAT_INITIAL_MATRIX, &A2);
577:   CheckMat(A2, B2, PETSC_FALSE, "reuse matrix (different type) MatTranspose");
578:   MatDestroy(&A2);
579:   MatDestroy(&B2);

581:   /* test MatISFixLocalEmpty */
582:   if (isaij) {
583:     PetscInt r[2];

585:     r[0] = 0;
586:     r[1] = PetscMin(m, n) - 1;
587:     PetscPrintf(PETSC_COMM_WORLD, "Test MatISFixLocalEmpty\n");
588:     MatDuplicate(A, MAT_COPY_VALUES, &A2);

590:     MatISFixLocalEmpty(A2, PETSC_TRUE);
591:     MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
592:     MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
593:     CheckMat(A2, B, PETSC_FALSE, "MatISFixLocalEmpty (null)");

595:     MatZeroRows(A2, 2, r, 0.0, NULL, NULL);
596:     MatViewFromOptions(A2, NULL, "-fixempty_view");
597:     MatDuplicate(B, MAT_COPY_VALUES, &B2);
598:     MatZeroRows(B2, 2, r, 0.0, NULL, NULL);
599:     MatISFixLocalEmpty(A2, PETSC_TRUE);
600:     MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
601:     MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
602:     MatViewFromOptions(A2, NULL, "-fixempty_view");
603:     CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows)");
604:     MatDestroy(&A2);

606:     MatDuplicate(A, MAT_COPY_VALUES, &A2);
607:     MatZeroRows(A2, 2, r, 0.0, NULL, NULL);
608:     MatTranspose(A2, MAT_INPLACE_MATRIX, &A2);
609:     MatTranspose(B2, MAT_INPLACE_MATRIX, &B2);
610:     MatViewFromOptions(A2, NULL, "-fixempty_view");
611:     MatISFixLocalEmpty(A2, PETSC_TRUE);
612:     MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
613:     MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
614:     MatViewFromOptions(A2, NULL, "-fixempty_view");
615:     CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (cols)");

617:     MatDestroy(&A2);
618:     MatDestroy(&B2);

620:     if (squaretest) {
621:       MatDuplicate(A, MAT_COPY_VALUES, &A2);
622:       MatDuplicate(B, MAT_COPY_VALUES, &B2);
623:       MatZeroRowsColumns(A2, 2, r, 0.0, NULL, NULL);
624:       MatZeroRowsColumns(B2, 2, r, 0.0, NULL, NULL);
625:       MatViewFromOptions(A2, NULL, "-fixempty_view");
626:       MatISFixLocalEmpty(A2, PETSC_TRUE);
627:       MatAssemblyBegin(A2, MAT_FINAL_ASSEMBLY);
628:       MatAssemblyEnd(A2, MAT_FINAL_ASSEMBLY);
629:       MatViewFromOptions(A2, NULL, "-fixempty_view");
630:       CheckMat(A2, B2, PETSC_FALSE, "MatISFixLocalEmpty (rows+cols)");
631:       MatDestroy(&A2);
632:       MatDestroy(&B2);
633:     }
634:   }

636:   /* test MatInvertBlockDiagonal
637:        special cases for block-diagonal matrices */
638:   if (m == n) {
639:     ISLocalToGlobalMapping map;
640:     Mat                    Abd, Bbd;
641:     IS                     is, bis;
642:     const PetscScalar     *isbd, *aijbd;
643:     PetscScalar           *vals;
644:     const PetscInt        *sts, *idxs;
645:     PetscInt              *idxs2, diff, perm, nl, bs, st, en, in;
646:     PetscBool              ok;

648:     for (diff = 0; diff < 3; diff++) {
649:       for (perm = 0; perm < 3; perm++) {
650:         for (bs = 1; bs < 4; bs++) {
651:           PetscPrintf(PETSC_COMM_WORLD, "Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", n, diff, perm, bs);
652:           PetscMalloc1(bs * bs, &vals);
653:           MatGetOwnershipRanges(A, &sts);
654:           switch (diff) {
655:           case 1: /* inverted layout by processes */
656:             in = 1;
657:             st = sts[size - rank - 1];
658:             en = sts[size - rank];
659:             nl = en - st;
660:             break;
661:           case 2: /* round-robin layout */
662:             in = size;
663:             st = rank;
664:             nl = n / size;
665:             if (rank < n % size) nl++;
666:             break;
667:           default: /* same layout */
668:             in = 1;
669:             st = sts[rank];
670:             en = sts[rank + 1];
671:             nl = en - st;
672:             break;
673:           }
674:           ISCreateStride(PETSC_COMM_WORLD, nl, st, in, &is);
675:           ISGetLocalSize(is, &nl);
676:           ISGetIndices(is, &idxs);
677:           PetscMalloc1(nl, &idxs2);
678:           for (i = 0; i < nl; i++) {
679:             switch (perm) { /* invert some of the indices */
680:             case 2:
681:               idxs2[i] = rank % 2 ? idxs[i] : idxs[nl - i - 1];
682:               break;
683:             case 1:
684:               idxs2[i] = rank % 2 ? idxs[nl - i - 1] : idxs[i];
685:               break;
686:             default:
687:               idxs2[i] = idxs[i];
688:               break;
689:             }
690:           }
691:           ISRestoreIndices(is, &idxs);
692:           ISCreateBlock(PETSC_COMM_WORLD, bs, nl, idxs2, PETSC_OWN_POINTER, &bis);
693:           ISLocalToGlobalMappingCreateIS(bis, &map);
694:           MatCreateIS(PETSC_COMM_WORLD, bs, PETSC_DECIDE, PETSC_DECIDE, bs * n, bs * n, map, map, &Abd);
695:           ISLocalToGlobalMappingDestroy(&map);
696:           MatISSetPreallocation(Abd, bs, NULL, 0, NULL);
697:           for (i = 0; i < nl; i++) {
698:             PetscInt b1, b2;

700:             for (b1 = 0; b1 < bs; b1++) {
701:               for (b2 = 0; b2 < bs; b2++) vals[b1 * bs + b2] = i * bs * bs + b1 * bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
702:             }
703:             MatSetValuesBlockedLocal(Abd, 1, &i, 1, &i, vals, INSERT_VALUES);
704:           }
705:           MatAssemblyBegin(Abd, MAT_FINAL_ASSEMBLY);
706:           MatAssemblyEnd(Abd, MAT_FINAL_ASSEMBLY);
707:           MatConvert(Abd, MATAIJ, MAT_INITIAL_MATRIX, &Bbd);
708:           MatInvertBlockDiagonal(Abd, &isbd);
709:           MatInvertBlockDiagonal(Bbd, &aijbd);
710:           MatGetLocalSize(Bbd, &nl, NULL);
711:           ok = PETSC_TRUE;
712:           for (i = 0; i < nl / bs; i++) {
713:             PetscInt b1, b2;

715:             for (b1 = 0; b1 < bs; b1++) {
716:               for (b2 = 0; b2 < bs; b2++) {
717:                 if (PetscAbsScalar(isbd[i * bs * bs + b1 * bs + b2] - aijbd[i * bs * bs + b1 * bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
718:                 if (!ok) {
719:                   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]));
720:                   break;
721:                 }
722:               }
723:               if (!ok) break;
724:             }
725:             if (!ok) break;
726:           }
727:           MatDestroy(&Abd);
728:           MatDestroy(&Bbd);
729:           PetscFree(vals);
730:           ISDestroy(&is);
731:           ISDestroy(&bis);
732:         }
733:       }
734:     }
735:   }

737:   /* test MatGetDiagonalBlock */
738:   PetscPrintf(PETSC_COMM_WORLD, "Test MatGetDiagonalBlock\n");
739:   MatGetDiagonalBlock(A, &A2);
740:   MatGetDiagonalBlock(B, &B2);
741:   CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock");
742:   MatScale(A, 2.0);
743:   MatScale(B, 2.0);
744:   MatGetDiagonalBlock(A, &A2);
745:   MatGetDiagonalBlock(B, &B2);
746:   CheckMat(A2, B2, PETSC_FALSE, "MatGetDiagonalBlock");

748:   /* free testing matrices */
749:   ISLocalToGlobalMappingDestroy(&cmap);
750:   ISLocalToGlobalMappingDestroy(&rmap);
751:   MatDestroy(&A);
752:   MatDestroy(&B);
753:   PetscFinalize();
754:   return 0;
755: }

757: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char *func)
758: {
759:   Mat       Bcheck;
760:   PetscReal error;

763:   if (!usemult && B) {
764:     PetscBool hasnorm;

766:     MatHasOperation(B, MATOP_NORM, &hasnorm);
767:     if (!hasnorm) usemult = PETSC_TRUE;
768:   }
769:   if (!usemult) {
770:     if (B) {
771:       MatType Btype;

773:       MatGetType(B, &Btype);
774:       MatConvert(A, Btype, MAT_INITIAL_MATRIX, &Bcheck);
775:     } else {
776:       MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck);
777:     }
778:     if (B) { /* if B is present, subtract it */
779:       MatAXPY(Bcheck, -1., B, DIFFERENT_NONZERO_PATTERN);
780:     }
781:     MatNorm(Bcheck, NORM_INFINITY, &error);
782:     if (error > PETSC_SQRT_MACHINE_EPSILON) {
783:       ISLocalToGlobalMapping rl2g, cl2g;

785:       PetscObjectSetName((PetscObject)Bcheck, "Bcheck");
786:       MatView(Bcheck, NULL);
787:       if (B) {
788:         PetscObjectSetName((PetscObject)B, "B");
789:         MatView(B, NULL);
790:         MatDestroy(&Bcheck);
791:         MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &Bcheck);
792:         PetscObjectSetName((PetscObject)Bcheck, "Assembled A");
793:         MatView(Bcheck, NULL);
794:       }
795:       MatDestroy(&Bcheck);
796:       PetscObjectSetName((PetscObject)A, "A");
797:       MatView(A, NULL);
798:       MatGetLocalToGlobalMapping(A, &rl2g, &cl2g);
799:       if (rl2g) ISLocalToGlobalMappingView(rl2g, NULL);
800:       if (cl2g) ISLocalToGlobalMappingView(cl2g, NULL);
801:       SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "ERROR ON %s: %g", func, (double)error);
802:     }
803:     MatDestroy(&Bcheck);
804:   } else {
805:     PetscBool ok, okt;

807:     MatMultEqual(A, B, 3, &ok);
808:     MatMultTransposeEqual(A, B, 3, &okt);
810:   }
811:   return 0;
812: }

814: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
815: {
816:   Mat                    B, Bcheck, B2 = NULL, lB;
817:   Vec                    x = NULL, b = NULL, b2 = NULL;
818:   ISLocalToGlobalMapping l2gr, l2gc;
819:   PetscReal              error;
820:   char                   diagstr[16];
821:   const PetscInt        *idxs;
822:   PetscInt               rst, ren, i, n, N, d;
823:   PetscMPIInt            rank;
824:   PetscBool              miss, haszerorows;

827:   if (diag == 0.) {
828:     PetscStrcpy(diagstr, "zero");
829:   } else {
830:     PetscStrcpy(diagstr, "nonzero");
831:   }
832:   ISView(is, NULL);
833:   MatGetLocalToGlobalMapping(A, &l2gr, &l2gc);
834:   /* tests MatDuplicate and MatCopy */
835:   if (diag == 0.) {
836:     MatDuplicate(A, MAT_COPY_VALUES, &B);
837:   } else {
838:     MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B);
839:     MatCopy(A, B, SAME_NONZERO_PATTERN);
840:   }
841:   MatISGetLocalMat(B, &lB);
842:   MatHasOperation(lB, MATOP_ZERO_ROWS, &haszerorows);
843:   if (squaretest && haszerorows) {
844:     MatCreateVecs(B, &x, &b);
845:     MatDuplicate(B, MAT_COPY_VALUES, &B2);
846:     VecSetLocalToGlobalMapping(b, l2gr);
847:     VecSetLocalToGlobalMapping(x, l2gc);
848:     VecSetRandom(x, NULL);
849:     VecSetRandom(b, NULL);
850:     /* mimic b[is] = x[is] */
851:     VecDuplicate(b, &b2);
852:     VecSetLocalToGlobalMapping(b2, l2gr);
853:     VecCopy(b, b2);
854:     ISGetLocalSize(is, &n);
855:     ISGetIndices(is, &idxs);
856:     VecGetSize(x, &N);
857:     for (i = 0; i < n; i++) {
858:       if (0 <= idxs[i] && idxs[i] < N) {
859:         VecSetValue(b2, idxs[i], diag, INSERT_VALUES);
860:         VecSetValue(x, idxs[i], 1., INSERT_VALUES);
861:       }
862:     }
863:     VecAssemblyBegin(b2);
864:     VecAssemblyEnd(b2);
865:     VecAssemblyBegin(x);
866:     VecAssemblyEnd(x);
867:     ISRestoreIndices(is, &idxs);
868:     /*  test ZeroRows on MATIS */
869:     PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr);
870:     MatZeroRowsIS(B, is, diag, x, b);
871:     PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRowsColumns (diag %s)\n", diagstr);
872:     MatZeroRowsColumnsIS(B2, is, diag, NULL, NULL);
873:   } else if (haszerorows) {
874:     /*  test ZeroRows on MATIS */
875:     PetscPrintf(PETSC_COMM_WORLD, "Test MatZeroRows (diag %s)\n", diagstr);
876:     MatZeroRowsIS(B, is, diag, NULL, NULL);
877:     b = b2 = x = NULL;
878:   } else {
879:     PetscPrintf(PETSC_COMM_WORLD, "Skipping MatZeroRows (diag %s)\n", diagstr);
880:     b = b2 = x = NULL;
881:   }

883:   if (squaretest && haszerorows) {
884:     VecAXPY(b2, -1., b);
885:     VecNorm(b2, NORM_INFINITY, &error);
887:   }

889:   /* test MatMissingDiagonal */
890:   PetscPrintf(PETSC_COMM_WORLD, "Test MatMissingDiagonal\n");
891:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
892:   MatMissingDiagonal(B, &miss, &d);
893:   MatGetOwnershipRange(B, &rst, &ren);
894:   PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
895:   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);
896:   PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
897:   PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);

899:   VecDestroy(&x);
900:   VecDestroy(&b);
901:   VecDestroy(&b2);

903:   /* check the result of ZeroRows with that from MPIAIJ routines
904:      assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
905:   if (haszerorows) {
906:     MatDuplicate(Afull, MAT_COPY_VALUES, &Bcheck);
907:     MatSetOption(Bcheck, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
908:     MatZeroRowsIS(Bcheck, is, diag, NULL, NULL);
909:     CheckMat(B, Bcheck, PETSC_FALSE, "Zerorows");
910:     MatDestroy(&Bcheck);
911:   }
912:   MatDestroy(&B);

914:   if (B2) { /* test MatZeroRowsColumns */
915:     MatDuplicate(Afull, MAT_COPY_VALUES, &B);
916:     MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
917:     MatZeroRowsColumnsIS(B, is, diag, NULL, NULL);
918:     CheckMat(B2, B, PETSC_FALSE, "MatZeroRowsColumns");
919:     MatDestroy(&B);
920:     MatDestroy(&B2);
921:   }
922:   return 0;
923: }

925: /*TEST

927:    test:
928:       args: -test_trans

930:    test:
931:       suffix: 2
932:       nsize: 4
933:       args: -matis_convert_local_nest -nr 3 -nc 4

935:    test:
936:       suffix: 3
937:       nsize: 5
938:       args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1

940:    test:
941:       suffix: 4
942:       nsize: 6
943:       args: -m 9 -n 12 -test_trans -nr 2 -nc 7

945:    test:
946:       suffix: 5
947:       nsize: 6
948:       args: -m 12 -n 12 -test_trans -nr 3 -nc 1

950:    test:
951:       suffix: 6
952:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap

954:    test:
955:       suffix: 7
956:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

958:    test:
959:       suffix: 8
960:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap

962:    test:
963:       suffix: 9
964:       nsize: 5
965:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap

967:    test:
968:       suffix: 10
969:       nsize: 5
970:       args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap

972:    test:
973:       suffix: vscat_default
974:       nsize: 5
975:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
976:       output_file: output/ex23_11.out

978:    test:
979:       suffix: 12
980:       nsize: 3
981:       args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3

983:    testset:
984:       output_file: output/ex23_13.out
985:       nsize: 3
986:       args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
987:       filter: grep -v "type:"
988:       test:
989:         suffix: baij
990:         args: -matis_localmat_type baij
991:       test:
992:         requires: viennacl
993:         suffix: viennacl
994:         args: -matis_localmat_type aijviennacl
995:       test:
996:         requires: cuda
997:         suffix: cusparse
998:         args: -matis_localmat_type aijcusparse

1000:    test:
1001:       suffix: negrep
1002:       nsize: {{1 3}separate output}
1003:       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}

1005: TEST*/