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;

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

 47:   /* create a MATIS matrix */
 48:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 49:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
 50:   PetscCall(MatSetType(A, MATIS));
 51:   PetscCall(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:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n, 0, 1, &is));
 57:   } else if (negmap && !repmap) { /* non repeated but with negative indices */
 58:     PetscCall(ISCreateStride(PETSC_COMM_WORLD, n + 2, -2, 1, &is));
 59:   } else if (!negmap && repmap) { /* non negative but repeated indices */
 60:     IS isl[2];

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

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

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

 88:   PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
 89:   PetscCall(MatISStoreL2L(A, PETSC_FALSE));
 90:   PetscCall(MatISSetPreallocation(A, 3, NULL, 3, NULL));
 91:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, (PetscBool) !(repmap || negmap))); /* I do not want to precompute the pattern */
 92:   PetscCall(ISLocalToGlobalMappingGetSize(rmap, &lm));
 93:   PetscCall(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:     PetscCall(ISGlobalToLocalMappingApply(cmap, IS_GTOLM_MASK, 3, cols, NULL, cols));
105:     PetscCall(MatSetValuesLocal(A, 1, &i, 3, cols, v, ADD_VALUES));
106:   }
107:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
108:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

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

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

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

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

134:   /* test MatGetInfo */
135:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatGetInfo\n"));
136:   PetscCall(MatGetInfo(A, MAT_LOCAL, &info));
137:   PetscCall(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:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
141:   PetscCall(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:   PetscCall(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:   PetscCall(MatIsSymmetric(A, 0.0, &issymmetric));
150:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Test MatIsSymmetric: %d\n", issymmetric));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

616:     PetscCall(MatDestroy(&A2));
617:     PetscCall(MatDestroy(&B2));

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

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

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

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

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

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

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

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

761:   PetscFunctionBeginUser;
762:   if (!usemult && B) {
763:     PetscBool hasnorm;

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

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

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

806:     PetscCall(MatMultEqual(A, B, 3, &ok));
807:     PetscCall(MatMultTransposeEqual(A, B, 3, &okt));
808:     PetscCheck(ok && okt, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR ON %s: mult ok ?  %d, multtranspose ok ? %d", func, ok, okt);
809:   }
810:   PetscFunctionReturn(PETSC_SUCCESS);
811: }

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

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

882:   if (squaretest && haszerorows) {
883:     PetscCall(VecAXPY(b2, -1., b));
884:     PetscCall(VecNorm(b2, NORM_INFINITY, &error));
885:     PetscCheck(error <= PETSC_SQRT_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "ERROR IN ZEROROWS ON B %g (diag %s)", (double)error, diagstr);
886:   }

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

898:   PetscCall(VecDestroy(&x));
899:   PetscCall(VecDestroy(&b));
900:   PetscCall(VecDestroy(&b2));

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

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

924: /*TEST

926:    test:
927:       args: -test_trans

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

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

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

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

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

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

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

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

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

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

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

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

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

1004: TEST*/