Actual source code: ex23.c
petsc-3.10.5 2019-03-28
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 MatPtAP(). It also tests some\n\
6: conversion routines.\n";
8: #include <petscmat.h>
10: PetscErrorCode TestMatZeroRows(Mat,Mat,PetscBool,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: MatType lmtype;
24: PetscScalar diag = 2.;
25: PetscInt n,m,i;
26: PetscInt rst,ren,cst,cen,nr,nc;
27: PetscMPIInt rank,size;
28: PetscBool testT,squaretest,isaij;
29: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE;
30: PetscErrorCode ierr;
32: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
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: if (size > 1 && m < 4) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 4 for parallel runs");
40: if (size == 1 && m < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of rows should be larger or equal 2 for uniprocessor runs");
41: if (n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be larger or equal 2");
42: if (symmetric) m = n = PetscMax(m,n);
44: /* create a MATIS matrix */
45: MatCreate(PETSC_COMM_WORLD,&A);
46: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
47: MatSetType(A,MATIS);
48: MatSetFromOptions(A);
49: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines */
50: ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);
51: ISLocalToGlobalMappingCreateIS(is,&cmap);
52: ISDestroy(&is);
54: PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL);
55: if (!symmetric && (diffmap || m != n)) {
56: PetscBool permute = PETSC_FALSE;
58: PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL);
59: ISCreateStride(PETSC_COMM_WORLD,m,permute ? m -1 : 0,permute ? -1 : 1,&is);
60: ISLocalToGlobalMappingCreateIS(is,&rmap);
61: ISDestroy(&is);
62: if (m==n && !permute) squaretest = PETSC_TRUE;
63: else squaretest = PETSC_FALSE;
64: } else {
65: PetscObjectReference((PetscObject)cmap);
66: rmap = cmap;
67: squaretest = PETSC_TRUE;
68: }
69: MatSetLocalToGlobalMapping(A,rmap,cmap);
70: MatISStoreL2L(A,PETSC_FALSE);
71: MatISSetPreallocation(A,3,NULL,0,NULL);
72: for (i=0; i<m; i++) {
73: PetscScalar v[3];
74: PetscInt cols[3];
76: cols[0] = (i-1+n)%n;
77: cols[1] = i%n;
78: cols[2] = (i+1)%n;
79: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
80: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
81: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
82: MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES);
83: }
84: if (symmetric) {
85: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
86: }
87: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
90: MatISGetLocalMat(A,&B);
91: MatGetType(B,&lmtype);
93: /* test MatGetInfo */
94: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n");
95: if (!PetscGlobalRank) {
96: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
97: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
98: }
99: MatGetInfo(A,MAT_LOCAL,&info);
100: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
101: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %D %D %D %D %D\n",PetscGlobalRank,(PetscInt)info.nz_used,
102: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
103: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
104: MatGetInfo(A,MAT_GLOBAL_MAX,&info);
105: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %D %D %D %D %D\n",(PetscInt)info.nz_used,
106: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
107: MatGetInfo(A,MAT_GLOBAL_SUM,&info);
108: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %D %D %D %D %D\n",(PetscInt)info.nz_used,
109: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
111: /* test MatView */
112: PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n");
113: MatView(A,NULL);
115: /* Create a MPIAIJ matrix, same as A */
116: MatCreate(PETSC_COMM_WORLD,&B);
117: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
118: MatSetType(B,MATAIJ);
119: MatSetFromOptions(B);
120: MatSetUp(B);
121: MatSetLocalToGlobalMapping(B,rmap,cmap);
122: MatMPIAIJSetPreallocation(B,3,NULL,3,NULL);
123: MatMPIBAIJSetPreallocation(B,1,3,NULL,3,NULL);
124: #if defined(PETSC_HAVE_HYPRE)
125: MatHYPRESetPreallocation(B,3,NULL,3,NULL);
126: #endif
127: MatISSetPreallocation(B,3,NULL,3,NULL);
128: for (i=0; i<m; i++) {
129: PetscScalar v[3];
130: PetscInt cols[3];
132: cols[0] = (i-1+n)%n;
133: cols[1] = i%n;
134: cols[2] = (i+1)%n;
135: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
136: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
137: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
138: MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES);
139: }
140: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
142: ISLocalToGlobalMappingDestroy(&cmap);
143: ISLocalToGlobalMappingDestroy(&rmap);
145: /* test CheckMat */
146: PetscPrintf(PETSC_COMM_WORLD,"Test CheckMat\n");
147: CheckMat(A,B,PETSC_FALSE,"CheckMat");
149: /* test MatDuplicate and MatAXPY */
150: PetscPrintf(PETSC_COMM_WORLD,"Test MatDuplicate and MatAXPY\n");
151: MatDuplicate(A,MAT_COPY_VALUES,&A2);
152: CheckMat(A,A2,PETSC_FALSE,"MatDuplicate and MatAXPY");
154: /* test MatConvert */
155: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_IS_XAIJ\n");
156: MatConvert(A2,MATAIJ,MAT_INITIAL_MATRIX,&B2);
157: CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INITIAL_MATRIX");
158: MatConvert(A2,MATAIJ,MAT_REUSE_MATRIX,&B2);
159: CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_REUSE_MATRIX");
160: MatConvert(A2,MATAIJ,MAT_INPLACE_MATRIX,&A2);
161: CheckMat(B,A2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INPLACE_MATRIX");
162: MatDestroy(&A2);
163: MatDestroy(&B2);
164: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_XAIJ_IS\n");
165: MatDuplicate(B,MAT_COPY_VALUES,&B2);
166: MatConvert(B2,MATIS,MAT_INITIAL_MATRIX,&A2);
167: CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INITIAL_MATRIX");
168: MatConvert(B2,MATIS,MAT_REUSE_MATRIX,&A2);
169: CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_REUSE_MATRIX");
170: MatConvert(B2,MATIS,MAT_INPLACE_MATRIX,&B2);
171: CheckMat(A,B2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INPLACE_MATRIX");
172: MatDestroy(&A2);
173: MatDestroy(&B2);
174: PetscStrcmp(lmtype,MATSEQAIJ,&isaij);
175: if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
176: 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};
178: for (ri = 0; ri < 2; ri++) {
179: PetscInt *r;
181: r = (PetscInt*)(ri == 0 ? rr : rk);
182: for (ci = 0; ci < 2; ci++) {
183: PetscInt *c,rb,cb;
185: c = (PetscInt*)(ci == 0 ? cr : ck);
186: for (rb = 1; rb < 4; rb++) {
187: ISCreateBlock(PETSC_COMM_SELF,rb,3,r,PETSC_COPY_VALUES,&is);
188: ISLocalToGlobalMappingCreateIS(is,&rmap);
189: ISDestroy(&is);
190: for (cb = 1; cb < 4; cb++) {
191: Mat T,lT,T2;
192: char testname[256];
194: PetscSNPrintf(testname,sizeof(testname),"MatConvert_IS_XAIJ special case (%D %D, bs %D %D)",ri,ci,rb,cb);
195: PetscPrintf(PETSC_COMM_WORLD,"Test %s\n",testname);
197: ISCreateBlock(PETSC_COMM_SELF,cb,4,c,PETSC_COPY_VALUES,&is);
198: ISLocalToGlobalMappingCreateIS(is,&cmap);
199: ISDestroy(&is);
201: MatCreate(PETSC_COMM_SELF,&T);
202: MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,rb*3,cb*4);
203: MatSetType(T,MATIS);
204: MatSetLocalToGlobalMapping(T,rmap,cmap);
205: ISLocalToGlobalMappingDestroy(&cmap);
206: MatISGetLocalMat(T,&lT);
207: MatSetType(lT,MATSEQAIJ);
208: MatSeqAIJSetPreallocation(lT,cb*4,NULL);
209: MatSetRandom(lT,NULL);
210: MatConvert(lT,lmtype,MAT_INPLACE_MATRIX,&lT);
211: MatISRestoreLocalMat(T,&lT);
212: MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
213: MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);
215: MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&T2);
216: CheckMat(T,T2,PETSC_TRUE,"MAT_INITIAL_MATRIX");
217: MatConvert(T,MATAIJ,MAT_REUSE_MATRIX,&T2);
218: CheckMat(T,T2,PETSC_TRUE,"MAT_REUSE_MATRIX");
219: MatDestroy(&T2);
220: MatDuplicate(T,MAT_COPY_VALUES,&T2);
221: MatConvert(T2,MATAIJ,MAT_INPLACE_MATRIX,&T2);
222: CheckMat(T,T2,PETSC_TRUE,"MAT_INPLACE_MATRIX");
223: MatDestroy(&T);
224: MatDestroy(&T2);
225: }
226: ISLocalToGlobalMappingDestroy(&rmap);
227: }
228: }
229: }
230: }
232: /* test MatDiagonalScale */
233: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n");
234: MatDuplicate(A,MAT_COPY_VALUES,&A2);
235: MatDuplicate(B,MAT_COPY_VALUES,&B2);
236: MatCreateVecs(A,&x,&y);
237: VecSetRandom(x,NULL);
238: if (symmetric) {
239: VecCopy(x,y);
240: } else {
241: VecSetRandom(y,NULL);
242: VecScale(y,8.);
243: }
244: MatDiagonalScale(A2,y,x);
245: MatDiagonalScale(B2,y,x);
246: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale");
247: MatDestroy(&A2);
248: MatDestroy(&B2);
249: VecDestroy(&x);
250: VecDestroy(&y);
252: /* test MatPtAP (A IS and B AIJ) */
253: if (isaij && m == n) {
254: PetscPrintf(PETSC_COMM_WORLD,"Test MatPtAP\n");
255: MatISStoreL2L(A,PETSC_TRUE);
256: MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A2);
257: MatPtAP(B,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B2);
258: CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_INITIAL_MATRIX");
259: MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&A2);
260: CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_REUSE_MATRIX");
261: MatDestroy(&A2);
262: MatDestroy(&B2);
263: }
265: /* test MatGetLocalSubMatrix */
266: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n");
267: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2);
268: ISCreateStride(PETSC_COMM_WORLD,m/2+m%2,0,2,&reven);
269: ISCreateStride(PETSC_COMM_WORLD,m/2,1,2,&rodd);
270: ISCreateStride(PETSC_COMM_WORLD,n/2+n%2,0,2,&ceven);
271: ISCreateStride(PETSC_COMM_WORLD,n/2,1,2,&codd);
272: MatGetLocalSubMatrix(A2,reven,ceven,&Aee);
273: MatGetLocalSubMatrix(A2,reven,codd,&Aeo);
274: MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe);
275: MatGetLocalSubMatrix(A2,rodd,codd,&Aoo);
276: for (i=0; i<m; i++) {
277: PetscInt j,je,jo,colse[3], colso[3];
278: PetscScalar ve[3], vo[3];
279: PetscScalar v[3];
280: PetscInt cols[3];
282: cols[0] = (i-1+n)%n;
283: cols[1] = i%n;
284: cols[2] = (i+1)%n;
285: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
286: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
287: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
288: for (j=0,je=0,jo=0;j<3;j++) {
289: if (cols[j]%2) {
290: vo[jo] = v[j];
291: colso[jo++] = cols[j]/2;
292: } else {
293: ve[je] = v[j];
294: colse[je++] = cols[j]/2;
295: }
296: }
297: if (i%2) {
298: PetscInt row = i/2;
299: MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES);
300: MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES);
301: } else {
302: PetscInt row = i/2;
303: MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES);
304: MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES);
305: }
306: }
307: MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee);
308: MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo);
309: MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe);
310: MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo);
311: ISDestroy(&reven);
312: ISDestroy(&ceven);
313: ISDestroy(&rodd);
314: ISDestroy(&codd);
315: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
316: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
317: MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN);
318: CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix");
319: MatDestroy(&A2);
321: /* test MatConvert_Nest_IS */
322: testT = PETSC_FALSE;
323: PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL);
325: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n");
326: nr = 2;
327: nc = 2;
328: PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL);
329: PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
330: if (testT) {
331: MatGetOwnershipRange(A,&cst,&cen);
332: MatGetOwnershipRangeColumn(A,&rst,&ren);
333: } else {
334: MatGetOwnershipRange(A,&rst,&ren);
335: MatGetOwnershipRangeColumn(A,&cst,&cen);
336: }
337: PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats);
338: for (i=0;i<nr*nc;i++) {
339: if (testT) {
340: MatCreateTranspose(A,&mats[i]);
341: MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc]);
342: } else {
343: MatDuplicate(A,MAT_COPY_VALUES,&mats[i]);
344: MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc]);
345: }
346: }
347: for (i=0;i<nr;i++) {
348: ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i]);
349: }
350: for (i=0;i<nc;i++) {
351: ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i]);
352: }
353: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2);
354: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2);
355: for (i=0;i<nr;i++) {
356: ISDestroy(&rows[i]);
357: }
358: for (i=0;i<nc;i++) {
359: ISDestroy(&cols[i]);
360: }
361: for (i=0;i<2*nr*nc;i++) {
362: MatDestroy(&mats[i]);
363: }
364: PetscFree3(rows,cols,mats);
365: MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T);
366: MatDestroy(&B2);
367: MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2);
368: CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INITIAL_MATRIX");
369: MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2);
370: CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_REUSE_MATRIX");
371: MatDestroy(&B2);
372: MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2);
373: CheckMat(A2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INPLACE_MATRIX");
374: MatDestroy(&T);
375: MatDestroy(&A2);
377: /* test MatCreateSubMatrix */
378: PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n");
379: if (!rank) {
380: ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is);
381: ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2);
382: } else if (rank == 1) {
383: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
384: if (n > 3) {
385: ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2);
386: } else {
387: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
388: }
389: } else if (rank == 2 && n > 4) {
390: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
391: ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2);
392: } else {
393: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
394: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
395: }
396: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2);
397: MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2);
398: CheckMat(A2,B2,PETSC_TRUE,"first MatCreateSubMatrix");
400: MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2);
401: MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2);
402: CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix");
403: MatDestroy(&A2);
404: MatDestroy(&B2);
406: if (!symmetric) {
407: MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2);
408: MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2);
409: MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2);
410: MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2);
411: CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix");
412: }
414: MatDestroy(&A2);
415: MatDestroy(&B2);
416: ISDestroy(&is);
417: ISDestroy(&is2);
419: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
420: if (size > 1) {
421: if (!rank) {
422: PetscInt st,len;
424: st = (m+1)/2;
425: len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0));
426: ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is);
427: } else {
428: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
429: }
430: } else {
431: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
432: }
434: if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
435: /* test MatDiagonalSet */
436: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n");
437: MatDuplicate(A,MAT_COPY_VALUES,&A2);
438: MatDuplicate(B,MAT_COPY_VALUES,&B2);
439: MatCreateVecs(A,NULL,&x);
440: VecSetRandom(x,NULL);
441: MatDiagonalSet(A2,x,INSERT_VALUES);
442: MatDiagonalSet(B2,x,INSERT_VALUES);
443: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet");
444: VecDestroy(&x);
445: MatDestroy(&A2);
446: MatDestroy(&B2);
448: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
449: PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n");
450: MatDuplicate(A,MAT_COPY_VALUES,&A2);
451: MatDuplicate(B,MAT_COPY_VALUES,&B2);
452: MatShift(A2,2.0);
453: MatShift(B2,2.0);
454: CheckMat(A2,B2,PETSC_FALSE,"MatShift");
455: MatDestroy(&A2);
456: MatDestroy(&B2);
458: /* nonzero diag value is supported for square matrices only */
459: TestMatZeroRows(A,B,PETSC_TRUE,is,diag);
460: }
461: TestMatZeroRows(A,B,squaretest,is,0.0);
462: ISDestroy(&is);
464: /* test MatTranspose */
465: PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n");
466: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
467: MatTranspose(B,MAT_INITIAL_MATRIX,&B2);
468: CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose");
470: MatTranspose(A,MAT_REUSE_MATRIX,&A2);
471: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose");
472: MatDestroy(&A2);
474: MatDuplicate(A,MAT_COPY_VALUES,&A2);
475: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
476: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose");
477: MatDestroy(&A2);
479: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
480: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose");
481: MatDestroy(&A2);
482: MatDestroy(&B2);
484: /* test MatISFixLocalEmpty */
485: if (isaij) {
486: PetscInt r[2];
488: r[0] = 0;
489: r[1] = PetscMin(m,n)-1;
490: PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n");
491: MatDuplicate(A,MAT_COPY_VALUES,&A2);
492: MatISFixLocalEmpty(A2,PETSC_TRUE);
493: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
494: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
495: CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)");
497: MatZeroRows(A2,2,r,0.0,NULL,NULL);
498: MatViewFromOptions(A2,NULL,"-fixempty_view");
499: MatDuplicate(B,MAT_COPY_VALUES,&B2);
500: MatZeroRows(B2,2,r,0.0,NULL,NULL);
501: MatISFixLocalEmpty(A2,PETSC_TRUE);
502: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
503: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
504: MatViewFromOptions(A2,NULL,"-fixempty_view");
505: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)");
506: MatDestroy(&A2);
508: MatDuplicate(A,MAT_COPY_VALUES,&A2);
509: MatZeroRows(A2,2,r,0.0,NULL,NULL);
510: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
511: MatTranspose(B2,MAT_INPLACE_MATRIX,&B2);
512: MatViewFromOptions(A2,NULL,"-fixempty_view");
513: MatISFixLocalEmpty(A2,PETSC_TRUE);
514: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
515: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
516: MatViewFromOptions(A2,NULL,"-fixempty_view");
517: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)");
519: MatDestroy(&A2);
520: MatDestroy(&B2);
522: if (squaretest) {
523: MatDuplicate(A,MAT_COPY_VALUES,&A2);
524: MatDuplicate(B,MAT_COPY_VALUES,&B2);
525: MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL);
526: MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL);
527: MatViewFromOptions(A2,NULL,"-fixempty_view");
528: MatISFixLocalEmpty(A2,PETSC_TRUE);
529: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
530: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
531: MatViewFromOptions(A2,NULL,"-fixempty_view");
532: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)");
533: MatDestroy(&A2);
534: MatDestroy(&B2);
535: }
537: }
538: /* free testing matrices */
539: MatDestroy(&A);
540: MatDestroy(&B);
541: PetscFinalize();
542: return ierr;
543: }
545: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func)
546: {
547: Mat Bcheck;
548: PetscReal error;
552: if (!usemult) {
553: if (B) {
554: MatType Btype;
556: MatGetType(B,&Btype);
557: MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck);
558: } else {
559: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
560: }
561: if (B) { /* if B is present, subtract it */
562: MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);
563: }
564: MatNorm(Bcheck,NORM_INFINITY,&error);
565: if (error > PETSC_SQRT_MACHINE_EPSILON) {
566: ISLocalToGlobalMapping rl2g,cl2g;
568: PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck");
569: MatView(Bcheck,NULL);
570: if (B) {
571: PetscObjectSetName((PetscObject)B,"Assembled AIJ");
572: MatView(B,NULL);
573: MatDestroy(&Bcheck);
574: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
575: PetscObjectSetName((PetscObject)Bcheck,"Assembled IS");
576: MatView(Bcheck,NULL);
577: }
578: MatDestroy(&Bcheck);
579: PetscObjectSetName((PetscObject)A,"MatIS");
580: MatView(A,NULL);
581: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
582: ISLocalToGlobalMappingView(rl2g,NULL);
583: ISLocalToGlobalMappingView(cl2g,NULL);
584: SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,error);
585: }
586: MatDestroy(&Bcheck);
587: } else {
588: PetscBool ok,okt;
590: MatMultEqual(A,B,3,&ok);
591: MatMultTransposeEqual(A,B,3,&okt);
592: if (!ok || !okt) SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt);
593: }
594: return(0);
595: }
597: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
598: {
599: Mat B,Bcheck,B2 = NULL,lB;
600: Vec x = NULL, b = NULL, b2 = NULL;
601: ISLocalToGlobalMapping l2gr,l2gc;
602: PetscReal error;
603: char diagstr[16];
604: const PetscInt *idxs;
605: PetscInt rst,ren,i,n,N,d;
606: PetscMPIInt rank;
607: PetscBool miss,haszerorows;
608: PetscErrorCode ierr;
611: if (diag == 0.) {
612: PetscStrcpy(diagstr,"zero");
613: } else {
614: PetscStrcpy(diagstr,"nonzero");
615: }
616: ISView(is,NULL);
617: MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);
618: /* tests MatDuplicate and MatCopy */
619: if (diag == 0.) {
620: MatDuplicate(A,MAT_COPY_VALUES,&B);
621: } else {
622: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
623: MatCopy(A,B,SAME_NONZERO_PATTERN);
624: }
625: MatISGetLocalMat(B,&lB);
626: MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows);
627: if (squaretest && haszerorows) {
629: MatCreateVecs(B,&x,&b);
630: MatDuplicate(B,MAT_COPY_VALUES,&B2);
631: VecSetLocalToGlobalMapping(b,l2gr);
632: VecSetLocalToGlobalMapping(x,l2gc);
633: VecSetRandom(x,NULL);
634: VecSetRandom(b,NULL);
635: /* mimic b[is] = x[is] */
636: VecDuplicate(b,&b2);
637: VecSetLocalToGlobalMapping(b2,l2gr);
638: VecCopy(b,b2);
639: ISGetLocalSize(is,&n);
640: ISGetIndices(is,&idxs);
641: VecGetSize(x,&N);
642: for (i=0;i<n;i++) {
643: if (0 <= idxs[i] && idxs[i] < N) {
644: VecSetValue(b2,idxs[i],diag,INSERT_VALUES);
645: VecSetValue(x,idxs[i],1.,INSERT_VALUES);
646: }
647: }
648: VecAssemblyBegin(b2);
649: VecAssemblyEnd(b2);
650: VecAssemblyBegin(x);
651: VecAssemblyEnd(x);
652: ISRestoreIndices(is,&idxs);
653: /* test ZeroRows on MATIS */
654: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
655: MatZeroRowsIS(B,is,diag,x,b);
656: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);
657: MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);
658: } else if (haszerorows) {
659: /* test ZeroRows on MATIS */
660: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
661: MatZeroRowsIS(B,is,diag,NULL,NULL);
662: b = b2 = x = NULL;
663: } else {
664: PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr);
665: b = b2 = x = NULL;
666: }
668: if (squaretest && haszerorows) {
669: VecAXPY(b2,-1.,b);
670: VecNorm(b2,NORM_INFINITY,&error);
671: if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",error,diagstr);
672: }
674: /* test MatMissingDiagonal */
675: PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");
676: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
677: MatMissingDiagonal(B,&miss,&d);
678: MatGetOwnershipRange(B,&rst,&ren);
679: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
680: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%D,%D) Missing %d, row %D (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);
681: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
682: PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
684: VecDestroy(&x);
685: VecDestroy(&b);
686: VecDestroy(&b2);
688: /* check the result of ZeroRows with that from MPIAIJ routines
689: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
690: if (haszerorows) {
691: MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);
692: MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
693: MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);
694: CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");
695: MatDestroy(&Bcheck);
696: }
697: MatDestroy(&B);
699: if (B2) { /* test MatZeroRowsColumns */
700: MatDuplicate(Afull,MAT_COPY_VALUES,&B);
701: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
702: MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);
703: CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");
704: MatDestroy(&B);
705: MatDestroy(&B2);
706: }
707: return(0);
708: }
711: /*TEST
713: test:
714: args: -test_trans
716: test:
717: suffix: 2
718: nsize: 4
719: args: -matis_convert_local_nest -nr 3 -nc 4
721: test:
722: suffix: 3
723: nsize: 5
724: args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1
726: test:
727: suffix: 4
728: nsize: 6
729: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
731: test:
732: suffix: 5
733: nsize: 6
734: args: -m 12 -n 12 -test_trans -nr 3 -nc 1
736: test:
737: suffix: 6
738: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
740: test:
741: suffix: 7
742: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
744: test:
745: suffix: 8
746: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
748: test:
749: suffix: 9
750: nsize: 5
751: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
753: test:
754: suffix: 10
755: nsize: 5
756: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
758: test:
759: suffix: 11
760: nsize: 5
761: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
763: test:
764: suffix: 12
765: nsize: 3
766: args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3
768: testset:
769: output_file: output/ex23_13.out
770: nsize: 3
771: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
772: filter: grep -v "type:"
773: test:
774: suffix: baij
775: args: -matis_localmat_type baij
776: test:
777: requires: viennacl
778: suffix: viennacl
779: args: -matis_localmat_type aijviennacl
780: test:
781: requires: veccuda
782: suffix: cusparse
783: args: -matis_localmat_type aijcusparse
785: TEST*/