Actual source code: ex23.c
petsc-3.8.4 2018-03-24
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] = { -1.*(i+1),2.*(i+1),-1.*(i+1)};
57: PetscInt cols[3] = {(i-1+n)%n,i%n,(i+1)%n};
58:
59: MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES);
60: }
61: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
64: /* test MatGetInfo */
65: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n");
66: MatISGetLocalMat(A,&B);
67: if (!PetscGlobalRank) {
68: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
70: }
71: MatGetInfo(A,MAT_LOCAL,&info);
72: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
73: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %D %D %D %D %D\n",PetscGlobalRank,(PetscInt)info.nz_used,
74: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
75: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
76: MatGetInfo(A,MAT_GLOBAL_MAX,&info);
77: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %D %D %D %D %D\n",(PetscInt)info.nz_used,
78: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
79: MatGetInfo(A,MAT_GLOBAL_SUM,&info);
80: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %D %D %D %D %D\n",(PetscInt)info.nz_used,
81: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
83: /* test MatView */
84: PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n");
85: MatView(A,NULL);
87: /* Create a MPIAIJ matrix, same as A */
88: MatCreate(PETSC_COMM_WORLD,&B);
89: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
90: MatSetType(B,MATAIJ);
91: MatSetUp(B);
92: MatSetLocalToGlobalMapping(B,rmap,cmap);
93: MatMPIAIJSetPreallocation(B,3,NULL,3,NULL);
94: for (i=0; i<m; i++) {
95: PetscScalar v[3] = { -1.*(i+1),2.*(i+1),-1.*(i+1)};
96: PetscInt cols[3] = {(i-1+n)%n,i%n,(i+1)%n};
97:
98: MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES);
99: }
100: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
102: ISLocalToGlobalMappingDestroy(&cmap);
103: ISLocalToGlobalMappingDestroy(&rmap);
105: /* test MatISGetMPIXAIJ */
106: PetscPrintf(PETSC_COMM_WORLD,"Test MatISGetMPIXAIJ\n");
107: CheckMat(A,B,PETSC_FALSE,"MatISGetMPIXAIJ");
109: /* test MatDiagonalScale */
110: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n");
111: MatDuplicate(A,MAT_COPY_VALUES,&A2);
112: MatDuplicate(B,MAT_COPY_VALUES,&B2);
113: MatCreateVecs(A,&x,&y);
114: VecSetRandom(x,NULL);
115: VecSetRandom(y,NULL);
116: VecScale(y,8.);
117: MatDiagonalScale(A2,y,x);
118: MatDiagonalScale(B2,y,x);
119: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale");
120: MatDestroy(&A2);
121: MatDestroy(&B2);
122: VecDestroy(&x);
123: VecDestroy(&y);
125: /* test MatGetLocalSubMatrix */
126: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n");
127: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2);
128: ISCreateStride(PETSC_COMM_WORLD,m/2+m%2,0,2,&reven);
129: ISCreateStride(PETSC_COMM_WORLD,m/2,1,2,&rodd);
130: ISCreateStride(PETSC_COMM_WORLD,n/2+n%2,0,2,&ceven);
131: ISCreateStride(PETSC_COMM_WORLD,n/2,1,2,&codd);
132: MatGetLocalSubMatrix(A2,reven,ceven,&Aee);
133: MatGetLocalSubMatrix(A2,reven,codd,&Aeo);
134: MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe);
135: MatGetLocalSubMatrix(A2,rodd,codd,&Aoo);
136: for (i=0; i<m; i++) {
137: PetscScalar v[3] = { -1.*(i+1),2.*(i+1),-1.*(i+1)};
138: PetscInt cols[3] = {(i-1+n)%n,i%n,(i+1)%n};
139: PetscInt j,je,jo,colse[3], colso[3];
140: PetscScalar ve[3], vo[3];
142: for (j=0,je=0,jo=0;j<3;j++) {
143: if (cols[j]%2) {
144: vo[jo] = v[j];
145: colso[jo++] = cols[j]/2;
146: } else {
147: ve[je] = v[j];
148: colse[je++] = cols[j]/2;
149: }
150: }
151: if (i%2) {
152: PetscInt row = i/2;
153: MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES);
154: MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES);
155: } else {
156: PetscInt row = i/2;
157: MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES);
158: MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES);
159: }
160: }
161: MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee);
162: MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo);
163: MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe);
164: MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo);
165: ISDestroy(&reven);
166: ISDestroy(&ceven);
167: ISDestroy(&rodd);
168: ISDestroy(&codd);
169: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
171: MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN);
172: CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix");
173: MatDestroy(&A2);
175: /* test MatConvert_Nest_IS */
176: testT = PETSC_FALSE;
177: PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL);
179: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n");
180: nr = 2;
181: nc = 2;
182: PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL);
183: PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
184: if (testT) {
185: MatGetOwnershipRange(A,&cst,&cen);
186: MatGetOwnershipRangeColumn(A,&rst,&ren);
187: } else {
188: MatGetOwnershipRange(A,&rst,&ren);
189: MatGetOwnershipRangeColumn(A,&cst,&cen);
190: }
191: PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats);
192: for (i=0;i<nr*nc;i++) {
193: if (testT) {
194: MatCreateTranspose(A,&mats[i]);
195: MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc]);
196: } else {
197: MatDuplicate(A,MAT_COPY_VALUES,&mats[i]);
198: MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc]);
199: }
200: }
201: for (i=0;i<nr;i++) {
202: ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i]);
203: }
204: for (i=0;i<nc;i++) {
205: ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i]);
206: }
207: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2);
208: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2);
209: for (i=0;i<nr;i++) {
210: ISDestroy(&rows[i]);
211: }
212: for (i=0;i<nc;i++) {
213: ISDestroy(&cols[i]);
214: }
215: for (i=0;i<2*nr*nc;i++) {
216: MatDestroy(&mats[i]);
217: }
218: PetscFree3(rows,cols,mats);
219: MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T);
220: MatDestroy(&B2);
221: MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2);
222: CheckMat(B2,T,testT,"MatConvert_Nest_IS MAT_INITIAL_MATRIX");
223: MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2);
224: CheckMat(B2,T,testT,"MatConvert_Nest_IS MAT_REUSE_MATRIX");
225: MatDestroy(&B2);
226: MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2);
227: CheckMat(A2,T,testT,"MatConvert_Nest_IS MAT_INPLACE_MATRIX");
228: MatDestroy(&T);
229: MatDestroy(&A2);
231: /* test MatCreateSubMatrix */
232: PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n");
233: if (!rank) {
234: ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is);
235: ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2);
236: } else if (rank == 1) {
237: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
238: ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2);
239: } else if (rank == 2 && n > 4) {
240: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
241: ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2);
242: } else {
243: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
244: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
245: }
246: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2);
247: MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2);
248: CheckMat(A2,B2,PETSC_FALSE,"first MatCreateSubMatrix");
250: MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2);
251: MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2);
252: CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix");
253: MatDestroy(&A2);
254: MatDestroy(&B2);
256: MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2);
257: MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2);
258: MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2);
259: MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2);
260: CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix");
262: MatDestroy(&A2);
263: MatDestroy(&B2);
264: ISDestroy(&is);
265: ISDestroy(&is2);
267: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
268: if (size > 1) {
269: if (!rank) {
270: PetscInt st = (m+1)/2;
271: PetscInt len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0));
272: ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is);
273: } else {
274: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
275: }
276: } else {
277: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
278: }
280: if (m == n) { /* tests for square matrices only */
281: /* test MatDiagonalSet */
282: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n");
283: MatDuplicate(A,MAT_COPY_VALUES,&A2);
284: MatDuplicate(B,MAT_COPY_VALUES,&B2);
285: MatCreateVecs(A,NULL,&x);
286: VecSetRandom(x,NULL);
287: MatDiagonalSet(A2,x,INSERT_VALUES);
288: MatDiagonalSet(B2,x,INSERT_VALUES);
289: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet");
290: VecDestroy(&x);
291: MatDestroy(&A2);
292: MatDestroy(&B2);
294: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
295: PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n");
296: MatDuplicate(A,MAT_COPY_VALUES,&A2);
297: MatDuplicate(B,MAT_COPY_VALUES,&B2);
298: MatShift(A2,2.0);
299: MatShift(B2,2.0);
300: CheckMat(A2,B2,PETSC_FALSE,"MatShift");
301: MatDestroy(&A2);
302: MatDestroy(&B2);
304: /* nonzero diag value is supported for square matrices only */
305: TestMatZeroRows(A,B,is,diag);
306: }
307: TestMatZeroRows(A,B,is,0.0);
308: ISDestroy(&is);
310: /* test MatTranspose */
311: PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n");
312: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
313: MatTranspose(B,MAT_INITIAL_MATRIX,&B2);
314: CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose");
316: MatTranspose(A,MAT_REUSE_MATRIX,&A2);
317: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose");
318: MatDestroy(&A2);
320: MatDuplicate(A,MAT_COPY_VALUES,&A2);
321: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
322: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose");
323: MatDestroy(&A2);
325: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
326: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose");
327: MatDestroy(&A2);
328: MatDestroy(&B2);
330: /* free testing matrices */
331: MatDestroy(&A);
332: MatDestroy(&B);
333: PetscFinalize();
334: return ierr;
335: }
337: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func)
338: {
339: Mat Bcheck;
340: PetscReal error;
344: if (!usemult) {
345: MatISGetMPIXAIJ(A,MAT_INITIAL_MATRIX,&Bcheck);
346: if (B) { /* if B is present, subtract it */
347: MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);
348: }
349: MatNorm(Bcheck,NORM_INFINITY,&error);
350: MatDestroy(&Bcheck);
351: if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,error);
352: } else {
353: Vec rv,rv2,lv,lv2;
354: PetscReal error1,error2;
356: MatCreateVecs(A,&rv,&lv);
357: MatCreateVecs(B,&rv2,&lv2);
359: VecSetRandom(rv,NULL);
360: MatMult(A,rv,lv);
361: MatMult(B,rv,lv2);
362: VecAXPY(lv,-1.,lv2);
363: VecNorm(lv,NORM_INFINITY,&error1);
365: VecSetRandom(lv,NULL);
366: MatMultTranspose(A,lv,rv);
367: MatMultTranspose(B,lv,rv2);
368: VecAXPY(rv,-1.,rv2);
369: VecNorm(rv,NORM_INFINITY,&error2);
371: VecDestroy(&lv);
372: VecDestroy(&lv2);
373: VecDestroy(&rv);
374: VecDestroy(&rv2);
376: 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);
377: }
378: return(0);
379: }
381: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, IS is, PetscScalar diag)
382: {
383: Mat B,Bcheck,B2 = NULL;
384: Vec x = NULL, b = NULL, b2 = NULL;
385: ISLocalToGlobalMapping l2gr,l2gc;
386: PetscReal error;
387: char diagstr[16];
388: const PetscInt *idxs;
389: PetscInt rst,ren,i,n,M,N,d;
390: PetscMPIInt rank;
391: PetscBool miss,square;
392: PetscErrorCode ierr;
395: if (diag == 0.) {
396: PetscStrcpy(diagstr,"zero");
397: } else {
398: PetscStrcpy(diagstr,"nonzero");
399: }
400: ISView(is,NULL);
401: MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);
402: /* tests MatDuplicate and MatCopy */
403: if (diag == 0.) {
404: MatDuplicate(A,MAT_COPY_VALUES,&B);
405: } else {
406: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
407: MatCopy(A,B,SAME_NONZERO_PATTERN);
408: }
409: MatGetSize(A,&M,&N);
410: square = M == N ? PETSC_TRUE : PETSC_FALSE;
411: if (square) {
412: MatCreateVecs(B,&x,&b);
413: MatDuplicate(B,MAT_COPY_VALUES,&B2);
414: VecSetLocalToGlobalMapping(b,l2gr);
415: VecSetLocalToGlobalMapping(x,l2gc);
416: VecSetRandom(x,NULL);
417: VecSetRandom(b,NULL);
418: /* mimic b[is] = x[is] */
419: VecDuplicate(b,&b2);
420: VecSetLocalToGlobalMapping(b2,l2gr);
421: VecCopy(b,b2);
422: ISGetLocalSize(is,&n);
423: ISGetIndices(is,&idxs);
424: VecGetSize(x,&N);
425: for (i=0;i<n;i++) {
426: if (0 <= idxs[i] && idxs[i] < N) {
427: VecSetValue(b2,idxs[i],diag,INSERT_VALUES);
428: VecSetValue(x,idxs[i],1.,INSERT_VALUES);
429: }
430: }
431: VecAssemblyBegin(b2);
432: VecAssemblyEnd(b2);
433: VecAssemblyBegin(x);
434: VecAssemblyEnd(x);
435: ISRestoreIndices(is,&idxs);
436: /* test ZeroRows on MATIS */
437: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
438: MatZeroRowsIS(B,is,diag,x,b);
439: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);
440: MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);
441: } else {
442: /* test ZeroRows on MATIS */
443: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
444: MatZeroRowsIS(B,is,diag,NULL,NULL);
445: b = b2 = x = NULL;
446: }
447: if (square) {
448: VecAXPY(b2,-1.,b);
449: VecNorm(b2,NORM_INFINITY,&error);
450: if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",error,diagstr);
451: }
452: /* test MatMissingDiagonal */
453: PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");
454: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
455: MatMissingDiagonal(B,&miss,&d);
456: MatGetOwnershipRange(B,&rst,&ren);
457: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
458: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%D,%D) Missing %d, row %D (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);
459: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
460: PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
462: VecDestroy(&x);
463: VecDestroy(&b);
464: VecDestroy(&b2);
466: /* check the result of ZeroRows with that from MPIAIJ routines
467: assuming that MatISGetMPIXAIJ and MatZeroRows_MPIAIJ work fine */
468: MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);
469: MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
470: MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);
471: CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");
472: MatDestroy(&Bcheck);
473: MatDestroy(&B);
475: if (B2) { /* test MatZeroRowsColumns */
476: MatDuplicate(Afull,MAT_COPY_VALUES,&B);
477: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
478: MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);
479: CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");
480: MatDestroy(&B);
481: MatDestroy(&B2);
482: }
483: return(0);
484: }