Actual source code: ex23.c

petsc-3.9.3 2018-07-02
Report Typos and Errors

  2: static char help[] = "Tests the use of interface functions for MATIS matrices.\n\
  3: This example tests: MatZeroRows(), MatZeroRowsLocal(), MatView(), MatDuplicate(),\n\
  4: MatCopy(), MatCreateSubMatrix(), MatGetLocalSubMatrix(), MatAXPY(), MatShift()\n\
  5: MatDiagonalSet(), MatTranspose() and MatISGetMPIXAIJ(). It also tests some\n\
  6: conversion routines.\n";

  8:  #include <petscmat.h>

 10: PetscErrorCode TestMatZeroRows(Mat,Mat,IS,PetscScalar);
 11: PetscErrorCode CheckMat(Mat,Mat,PetscBool,const char*);

 13: int main(int argc,char **args)
 14: {
 15:   Mat                    A,B,A2,B2,T;
 16:   Mat                    Aee,Aeo,Aoe,Aoo;
 17:   Mat                    *mats;
 18:   Vec                    x,y;
 19:   MatInfo                info;
 20:   ISLocalToGlobalMapping cmap,rmap;
 21:   IS                     is,is2,reven,rodd,ceven,codd;
 22:   IS                     *rows,*cols;
 23:   PetscScalar            diag = 2.;
 24:   PetscInt               n,m,i;
 25:   PetscInt               rst,ren,cst,cen,nr,nc;
 26:   PetscMPIInt            rank,size;
 27:   PetscBool              testT;
 28:   PetscErrorCode         ierr;

 30:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 31:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 32:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 33:   m = n = 2*size;
 34:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 35:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 36:   if (size > 1 && m < 4) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be more than 4 for parallel runs");
 37:   if (size == 1 && m < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be more than 2 for uniprocessor runs");
 38:   if (n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be more than 2");

 40:   /* create a MATIS matrix */
 41:   MatCreate(PETSC_COMM_WORLD,&A);
 42:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
 43:   MatSetType(A,MATIS);
 44:   /* Each process has the same l2gmap for rows and cols
 45:      This is not the proper setting for MATIS for finite elements, it is just used to test the routines */
 46:   ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);
 47:   ISLocalToGlobalMappingCreateIS(is,&cmap);
 48:   ISDestroy(&is);
 49:   ISCreateStride(PETSC_COMM_WORLD,m,0,1,&is);
 50:   ISLocalToGlobalMappingCreateIS(is,&rmap);
 51:   ISDestroy(&is);
 52:   MatSetLocalToGlobalMapping(A,rmap,cmap);
 53:   MatSetUp(A);
 54:   MatISSetPreallocation(A,3,NULL,0,NULL);
 55:   for (i=0; i<m; i++) {
 56:     PetscScalar v[3];
 57:     PetscInt    cols[3];

 59:     v[0]    = -1.*(i+1);
 60:     v[1]    = 2.*(i+1);
 61:     v[2]    = -1.*(i+1);
 62:     cols[0] = (i-1+n)%n;
 63:     cols[1] = i%n;
 64:     cols[2] = (i+1)%n;
 65:     MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES);
 66:   }
 67:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 68:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 70:   /* test MatGetInfo */
 71:   PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n");
 72:   MatISGetLocalMat(A,&B);
 73:   if (!PetscGlobalRank) {
 74:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 75:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 76:   }
 77:   MatGetInfo(A,MAT_LOCAL,&info);
 78:   PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
 79:   PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process  %2d: %D %D %D %D %D\n",PetscGlobalRank,(PetscInt)info.nz_used,
 80:                                             (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
 81:   PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
 82:   MatGetInfo(A,MAT_GLOBAL_MAX,&info);
 83:   PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax  : %D %D %D %D %D\n",(PetscInt)info.nz_used,
 84:                                 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
 85:   MatGetInfo(A,MAT_GLOBAL_SUM,&info);
 86:   PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum  : %D %D %D %D %D\n",(PetscInt)info.nz_used,
 87:                                 (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);

 89:   /* test MatView */
 90:   PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n");
 91:   MatView(A,NULL);

 93:   /* Create a MPIAIJ matrix, same as A */
 94:   MatCreate(PETSC_COMM_WORLD,&B);
 95:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
 96:   MatSetType(B,MATAIJ);
 97:   MatSetUp(B);
 98:   MatSetLocalToGlobalMapping(B,rmap,cmap);
 99:   MatMPIAIJSetPreallocation(B,3,NULL,3,NULL);
100:   for (i=0; i<m; i++) {
101:     PetscScalar v[3];
102:     PetscInt    cols[3];

104:     v[0]    = -1.*(i+1);
105:     v[1]    = 2.*(i+1);
106:     v[2]    = -1.*(i+1);
107:     cols[0] = (i-1+n)%n;
108:     cols[1] = i%n;
109:     cols[2] = (i+1)%n;
110:     MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES);
111:   }
112:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
113:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
114:   ISLocalToGlobalMappingDestroy(&cmap);
115:   ISLocalToGlobalMappingDestroy(&rmap);

117:   /* test MatISGetMPIXAIJ */
118:   PetscPrintf(PETSC_COMM_WORLD,"Test MatISGetMPIXAIJ\n");
119:   CheckMat(A,B,PETSC_FALSE,"MatISGetMPIXAIJ");

121:   /* test MatDiagonalScale */
122:   PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n");
123:   MatDuplicate(A,MAT_COPY_VALUES,&A2);
124:   MatDuplicate(B,MAT_COPY_VALUES,&B2);
125:   MatCreateVecs(A,&x,&y);
126:   VecSetRandom(x,NULL);
127:   VecSetRandom(y,NULL);
128:   VecScale(y,8.);
129:   MatDiagonalScale(A2,y,x);
130:   MatDiagonalScale(B2,y,x);
131:   CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale");
132:   MatDestroy(&A2);
133:   MatDestroy(&B2);
134:   VecDestroy(&x);
135:   VecDestroy(&y);

137:   /* test MatGetLocalSubMatrix */
138:   PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n");
139:   MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2);
140:   ISCreateStride(PETSC_COMM_WORLD,m/2+m%2,0,2,&reven);
141:   ISCreateStride(PETSC_COMM_WORLD,m/2,1,2,&rodd);
142:   ISCreateStride(PETSC_COMM_WORLD,n/2+n%2,0,2,&ceven);
143:   ISCreateStride(PETSC_COMM_WORLD,n/2,1,2,&codd);
144:   MatGetLocalSubMatrix(A2,reven,ceven,&Aee);
145:   MatGetLocalSubMatrix(A2,reven,codd,&Aeo);
146:   MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe);
147:   MatGetLocalSubMatrix(A2,rodd,codd,&Aoo);
148:   for (i=0; i<m; i++) {
149:     PetscInt    j,je,jo,colse[3], colso[3];
150:     PetscScalar ve[3], vo[3];
151:     PetscScalar v[3];
152:     PetscInt    cols[3];

154:     v[0]    = -1.*(i+1);
155:     v[1]    = 2.*(i+1);
156:     v[2]    = -1.*(i+1);
157:     cols[0] = (i-1+n)%n;
158:     cols[1] = i%n;
159:     cols[2] = (i+1)%n;
160:     for (j=0,je=0,jo=0;j<3;j++) {
161:       if (cols[j]%2) {
162:         vo[jo] = v[j];
163:         colso[jo++] = cols[j]/2;
164:       } else {
165:         ve[je] = v[j];
166:         colse[je++] = cols[j]/2;
167:       }
168:     }
169:     if (i%2) {
170:       PetscInt row = i/2;
171:       MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES);
172:       MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES);
173:     } else {
174:       PetscInt row = i/2;
175:       MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES);
176:       MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES);
177:     }
178:   }
179:   MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee);
180:   MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo);
181:   MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe);
182:   MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo);
183:   ISDestroy(&reven);
184:   ISDestroy(&ceven);
185:   ISDestroy(&rodd);
186:   ISDestroy(&codd);
187:   MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
188:   MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
189:   MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN);
190:   CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix");
191:   MatDestroy(&A2);

193:   /* test MatConvert_Nest_IS */
194:   testT = PETSC_FALSE;
195:   PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL);

197:   PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n");
198:   nr   = 2;
199:   nc   = 2;
200:   PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL);
201:   PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
202:   if (testT) {
203:     MatGetOwnershipRange(A,&cst,&cen);
204:     MatGetOwnershipRangeColumn(A,&rst,&ren);
205:   } else {
206:     MatGetOwnershipRange(A,&rst,&ren);
207:     MatGetOwnershipRangeColumn(A,&cst,&cen);
208:   }
209:   PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats);
210:   for (i=0;i<nr*nc;i++) {
211:     if (testT) {
212:       MatCreateTranspose(A,&mats[i]);
213:       MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc]);
214:     } else {
215:       MatDuplicate(A,MAT_COPY_VALUES,&mats[i]);
216:       MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc]);
217:     }
218:   }
219:   for (i=0;i<nr;i++) {
220:     ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i]);
221:   }
222:   for (i=0;i<nc;i++) {
223:     ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i]);
224:   }
225:   MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2);
226:   MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2);
227:   for (i=0;i<nr;i++) {
228:     ISDestroy(&rows[i]);
229:   }
230:   for (i=0;i<nc;i++) {
231:     ISDestroy(&cols[i]);
232:   }
233:   for (i=0;i<2*nr*nc;i++) {
234:     MatDestroy(&mats[i]);
235:   }
236:   PetscFree3(rows,cols,mats);
237:   MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T);
238:   MatDestroy(&B2);
239:   MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2);
240:   CheckMat(B2,T,testT,"MatConvert_Nest_IS MAT_INITIAL_MATRIX");
241:   MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2);
242:   CheckMat(B2,T,testT,"MatConvert_Nest_IS MAT_REUSE_MATRIX");
243:   MatDestroy(&B2);
244:   MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2);
245:   CheckMat(A2,T,testT,"MatConvert_Nest_IS MAT_INPLACE_MATRIX");
246:   MatDestroy(&T);
247:   MatDestroy(&A2);

249:   /* test MatCreateSubMatrix */
250:   PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n");
251:   if (!rank) {
252:     ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is);
253:     ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2);
254:   } else if (rank == 1) {
255:     ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
256:     ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2);
257:   } else if (rank == 2 && n > 4) {
258:     ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
259:     ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2);
260:   } else {
261:     ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
262:     ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
263:   }
264:   MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2);
265:   MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2);
266:   CheckMat(A2,B2,PETSC_FALSE,"first MatCreateSubMatrix");

268:   MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2);
269:   MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2);
270:   CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix");
271:   MatDestroy(&A2);
272:   MatDestroy(&B2);

274:   MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2);
275:   MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2);
276:   MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2);
277:   MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2);
278:   CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix");

280:   MatDestroy(&A2);
281:   MatDestroy(&B2);
282:   ISDestroy(&is);
283:   ISDestroy(&is2);

285:   /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
286:   if (size > 1) {
287:     if (!rank) {
288:       PetscInt st,len;

290:       st   = (m+1)/2;
291:       len  = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0));
292:       ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is);
293:     } else {
294:       ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
295:     }
296:   } else {
297:     ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
298:   }

300:   if (m == n) { /* tests for square matrices only */
301:     /* test MatDiagonalSet */
302:     PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n");
303:     MatDuplicate(A,MAT_COPY_VALUES,&A2);
304:     MatDuplicate(B,MAT_COPY_VALUES,&B2);
305:     MatCreateVecs(A,NULL,&x);
306:     VecSetRandom(x,NULL);
307:     MatDiagonalSet(A2,x,INSERT_VALUES);
308:     MatDiagonalSet(B2,x,INSERT_VALUES);
309:     CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet");
310:     VecDestroy(&x);
311:     MatDestroy(&A2);
312:     MatDestroy(&B2);

314:     /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
315:     PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n");
316:     MatDuplicate(A,MAT_COPY_VALUES,&A2);
317:     MatDuplicate(B,MAT_COPY_VALUES,&B2);
318:     MatShift(A2,2.0);
319:     MatShift(B2,2.0);
320:     CheckMat(A2,B2,PETSC_FALSE,"MatShift");
321:     MatDestroy(&A2);
322:     MatDestroy(&B2);

324:     /* nonzero diag value is supported for square matrices only */
325:     TestMatZeroRows(A,B,is,diag);
326:   }
327:   TestMatZeroRows(A,B,is,0.0);
328:   ISDestroy(&is);

330:   /* test MatTranspose */
331:   PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n");
332:   MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
333:   MatTranspose(B,MAT_INITIAL_MATRIX,&B2);
334:   CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose");

336:   MatTranspose(A,MAT_REUSE_MATRIX,&A2);
337:   CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose");
338:   MatDestroy(&A2);

340:   MatDuplicate(A,MAT_COPY_VALUES,&A2);
341:   MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
342:   CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose");
343:   MatDestroy(&A2);

345:   MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
346:   CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose");
347:   MatDestroy(&A2);
348:   MatDestroy(&B2);

350:   /* free testing matrices */
351:   MatDestroy(&A);
352:   MatDestroy(&B);
353:   PetscFinalize();
354:   return ierr;
355: }

357: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func)
358: {
359:   Mat            Bcheck;
360:   PetscReal      error;

364:   if (!usemult) {
365:     MatISGetMPIXAIJ(A,MAT_INITIAL_MATRIX,&Bcheck);
366:     if (B) { /* if B is present, subtract it */
367:       MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);
368:     }
369:     MatNorm(Bcheck,NORM_INFINITY,&error);
370:     MatDestroy(&Bcheck);
371:     if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,error);
372:   } else {
373:     Vec       rv,rv2,lv,lv2;
374:     PetscReal error1,error2;

376:     MatCreateVecs(A,&rv,&lv);
377:     MatCreateVecs(B,&rv2,&lv2);

379:     VecSetRandom(rv,NULL);
380:     MatMult(A,rv,lv);
381:     MatMult(B,rv,lv2);
382:     VecAXPY(lv,-1.,lv2);
383:     VecNorm(lv,NORM_INFINITY,&error1);

385:     VecSetRandom(lv,NULL);
386:     MatMultTranspose(A,lv,rv);
387:     MatMultTranspose(B,lv,rv2);
388:     VecAXPY(rv,-1.,rv2);
389:     VecNorm(rv,NORM_INFINITY,&error2);

391:     VecDestroy(&lv);
392:     VecDestroy(&lv2);
393:     VecDestroy(&rv);
394:     VecDestroy(&rv2);

396:     if (error1 > PETSC_SQRT_MACHINE_EPSILON || error2 > PETSC_SQRT_MACHINE_EPSILON) SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: (mult) %g, (multt) %g",func,error1,error2);
397:   }
398:   return(0);
399: }

401: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, IS is, PetscScalar diag)
402: {
403:   Mat                    B,Bcheck,B2 = NULL;
404:   Vec                    x = NULL, b = NULL, b2 = NULL;
405:   ISLocalToGlobalMapping l2gr,l2gc;
406:   PetscReal              error;
407:   char                   diagstr[16];
408:   const PetscInt         *idxs;
409:   PetscInt               rst,ren,i,n,M,N,d;
410:   PetscMPIInt            rank;
411:   PetscBool              miss,square;
412:   PetscErrorCode         ierr;

415:   if (diag == 0.) {
416:     PetscStrcpy(diagstr,"zero");
417:   } else {
418:     PetscStrcpy(diagstr,"nonzero");
419:   }
420:   ISView(is,NULL);
421:   MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);
422:   /* tests MatDuplicate and MatCopy */
423:   if (diag == 0.) {
424:     MatDuplicate(A,MAT_COPY_VALUES,&B);
425:   } else {
426:     MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
427:     MatCopy(A,B,SAME_NONZERO_PATTERN);
428:   }
429:   MatGetSize(A,&M,&N);
430:   square = M == N ? PETSC_TRUE : PETSC_FALSE;
431:   if (square) {
432:     MatCreateVecs(B,&x,&b);
433:     MatDuplicate(B,MAT_COPY_VALUES,&B2);
434:     VecSetLocalToGlobalMapping(b,l2gr);
435:     VecSetLocalToGlobalMapping(x,l2gc);
436:     VecSetRandom(x,NULL);
437:     VecSetRandom(b,NULL);
438:     /* mimic b[is] = x[is] */
439:     VecDuplicate(b,&b2);
440:     VecSetLocalToGlobalMapping(b2,l2gr);
441:     VecCopy(b,b2);
442:     ISGetLocalSize(is,&n);
443:     ISGetIndices(is,&idxs);
444:     VecGetSize(x,&N);
445:     for (i=0;i<n;i++) {
446:       if (0 <= idxs[i] && idxs[i] < N) {
447:         VecSetValue(b2,idxs[i],diag,INSERT_VALUES);
448:         VecSetValue(x,idxs[i],1.,INSERT_VALUES);
449:       }
450:     }
451:     VecAssemblyBegin(b2);
452:     VecAssemblyEnd(b2);
453:     VecAssemblyBegin(x);
454:     VecAssemblyEnd(x);
455:     ISRestoreIndices(is,&idxs);
456:     /*  test ZeroRows on MATIS */
457:     PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
458:     MatZeroRowsIS(B,is,diag,x,b);
459:     PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);
460:     MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);
461:   } else {
462:     /*  test ZeroRows on MATIS */
463:     PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
464:     MatZeroRowsIS(B,is,diag,NULL,NULL);
465:     b = b2 = x = NULL;
466:   }
467:   if (square) {
468:     VecAXPY(b2,-1.,b);
469:     VecNorm(b2,NORM_INFINITY,&error);
470:     if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",error,diagstr);
471:   }
472:   /* test MatMissingDiagonal */
473:   PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");
474:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
475:   MatMissingDiagonal(B,&miss,&d);
476:   MatGetOwnershipRange(B,&rst,&ren);
477:   PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
478:   PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%D,%D) Missing %d, row %D (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);
479:   PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
480:   PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);

482:   VecDestroy(&x);
483:   VecDestroy(&b);
484:   VecDestroy(&b2);

486:   /* check the result of ZeroRows with that from MPIAIJ routines
487:      assuming that MatISGetMPIXAIJ and MatZeroRows_MPIAIJ work fine */
488:   MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);
489:   MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
490:   MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);
491:   CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");
492:   MatDestroy(&Bcheck);
493:   MatDestroy(&B);

495:   if (B2) { /* test MatZeroRowsColumns */
496:     MatDuplicate(Afull,MAT_COPY_VALUES,&B);
497:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
498:     MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);
499:     CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");
500:     MatDestroy(&B);
501:     MatDestroy(&B2);
502:   }
503:   return(0);
504: }


507: /*TEST

509:    test:
510:       args: -test_trans

512:    test:
513:       suffix: 2
514:       nsize: 4
515:       args: -matis_convert_local_nest -nr 3 -nc 4

517:    test:
518:       suffix: 3
519:       nsize: 5
520:       args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1

522:    test:
523:       suffix: 4
524:       nsize: 6
525:       args: -m 9 -n 12 -test_trans -nr 2 -nc 7

527: TEST*/