Actual source code: ex23.c
petsc-3.14.6 2021-03-30
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 permute = PETSC_FALSE;
30: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE;
31: PetscErrorCode ierr;
33: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
34: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
35: MPI_Comm_size(PETSC_COMM_WORLD,&size);
36: m = n = 2*size;
37: PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);
38: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
39: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
40: 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");
41: 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");
42: if (n < 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of cols should be larger or equal 2");
43: if (symmetric) m = n = PetscMax(m,n);
45: /* create a MATIS matrix */
46: MatCreate(PETSC_COMM_WORLD,&A);
47: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
48: MatSetType(A,MATIS);
49: MatSetFromOptions(A);
50: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines */
51: ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);
52: ISLocalToGlobalMappingCreateIS(is,&cmap);
53: ISDestroy(&is);
55: PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL);
56: PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL);
57: if (!symmetric && (diffmap || m != n)) {
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: }
539: /* test MatInvertBlockDiagonal
540: special cases for block-diagonal matrices */
541: if (m == n) {
542: ISLocalToGlobalMapping map;
543: Mat Abd,Bbd;
544: IS is,bis;
545: const PetscScalar *isbd,*aijbd;
546: PetscScalar *vals;
547: const PetscInt *sts,*idxs;
548: PetscInt *idxs2,diff,perm,nl,bs,st,en,in;
549: PetscBool ok;
551: for (diff = 0; diff < 3; diff++) {
552: for (perm = 0; perm < 3; perm++) {
553: for (bs = 1; bs < 4; bs++) {
554: PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %D %D %D %D\n",n,diff,perm,bs);
555: PetscMalloc1(bs*bs,&vals);
556: MatGetOwnershipRanges(A,&sts);
557: switch (diff) {
558: case 1: /* inverted layout by processes */
559: in = 1;
560: st = sts[size - rank - 1];
561: en = sts[size - rank];
562: nl = en - st;
563: break;
564: case 2: /* round-robin layout */
565: in = size;
566: st = rank;
567: nl = n/size;
568: if (rank < n%size) nl++;
569: break;
570: default: /* same layout */
571: in = 1;
572: st = sts[rank];
573: en = sts[rank + 1];
574: nl = en - st;
575: break;
576: }
577: ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is);
578: ISGetLocalSize(is,&nl);
579: ISGetIndices(is,&idxs);
580: PetscMalloc1(nl,&idxs2);
581: for (i=0;i<nl;i++) {
582: switch (perm) { /* invert some of the indices */
583: case 2:
584: idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1];
585: break;
586: case 1:
587: idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i];
588: break;
589: default:
590: idxs2[i] = idxs[i];
591: break;
592: }
593: }
594: ISRestoreIndices(is,&idxs);
595: ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis);
596: ISLocalToGlobalMappingCreateIS(bis,&map);
597: MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,NULL,&Abd);
598: ISLocalToGlobalMappingDestroy(&map);
599: MatISSetPreallocation(Abd,bs,NULL,0,NULL);
600: for (i=0;i<nl;i++) {
601: PetscInt b1,b2;
603: for (b1=0;b1<bs;b1++) {
604: for (b2=0;b2<bs;b2++) {
605: vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
606: }
607: }
608: MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES);
609: }
610: MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY);
611: MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY);
612: MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd);
613: MatInvertBlockDiagonal(Abd,&isbd);
614: MatInvertBlockDiagonal(Bbd,&aijbd);
615: MatGetLocalSize(Bbd,&nl,NULL);
616: ok = PETSC_TRUE;
617: for (i=0;i<nl/bs;i++) {
618: PetscInt b1,b2;
620: for (b1=0;b1<bs;b1++) {
621: for (b2=0;b2<bs;b2++) {
622: if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
623: if (!ok) {
624: PetscPrintf(PETSC_COMM_SELF,"[%d] ERROR block %d, entry %d %d: %g %g\n",rank,i,b1,b2,isbd[i*bs*bs+b1*bs + b2],aijbd[i*bs*bs+b1*bs + b2]);
625: break;
626: }
627: }
628: if (!ok) break;
629: }
630: if (!ok) break;
631: }
632: MatDestroy(&Abd);
633: MatDestroy(&Bbd);
634: PetscFree(vals);
635: ISDestroy(&is);
636: ISDestroy(&bis);
637: }
638: }
639: }
640: }
641: /* free testing matrices */
642: MatDestroy(&A);
643: MatDestroy(&B);
644: PetscFinalize();
645: return ierr;
646: }
648: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func)
649: {
650: Mat Bcheck;
651: PetscReal error;
655: if (!usemult && B) {
656: PetscBool hasnorm;
658: MatHasOperation(B,MATOP_NORM,&hasnorm);
659: if (!hasnorm) usemult = PETSC_TRUE;
660: }
661: if (!usemult) {
662: if (B) {
663: MatType Btype;
665: MatGetType(B,&Btype);
666: MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck);
667: } else {
668: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
669: }
670: if (B) { /* if B is present, subtract it */
671: MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);
672: }
673: MatNorm(Bcheck,NORM_INFINITY,&error);
674: if (error > PETSC_SQRT_MACHINE_EPSILON) {
675: ISLocalToGlobalMapping rl2g,cl2g;
677: PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck");
678: MatView(Bcheck,NULL);
679: if (B) {
680: PetscObjectSetName((PetscObject)B,"Assembled AIJ");
681: MatView(B,NULL);
682: MatDestroy(&Bcheck);
683: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
684: PetscObjectSetName((PetscObject)Bcheck,"Assembled IS");
685: MatView(Bcheck,NULL);
686: }
687: MatDestroy(&Bcheck);
688: PetscObjectSetName((PetscObject)A,"MatIS");
689: MatView(A,NULL);
690: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
691: ISLocalToGlobalMappingView(rl2g,NULL);
692: ISLocalToGlobalMappingView(cl2g,NULL);
693: SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,error);
694: }
695: MatDestroy(&Bcheck);
696: } else {
697: PetscBool ok,okt;
699: MatMultEqual(A,B,3,&ok);
700: MatMultTransposeEqual(A,B,3,&okt);
701: if (!ok || !okt) SETERRQ3(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: mult ok ? %d, multtranspose ok ? %d",func,ok,okt);
702: }
703: return(0);
704: }
706: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
707: {
708: Mat B,Bcheck,B2 = NULL,lB;
709: Vec x = NULL, b = NULL, b2 = NULL;
710: ISLocalToGlobalMapping l2gr,l2gc;
711: PetscReal error;
712: char diagstr[16];
713: const PetscInt *idxs;
714: PetscInt rst,ren,i,n,N,d;
715: PetscMPIInt rank;
716: PetscBool miss,haszerorows;
717: PetscErrorCode ierr;
720: if (diag == 0.) {
721: PetscStrcpy(diagstr,"zero");
722: } else {
723: PetscStrcpy(diagstr,"nonzero");
724: }
725: ISView(is,NULL);
726: MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);
727: /* tests MatDuplicate and MatCopy */
728: if (diag == 0.) {
729: MatDuplicate(A,MAT_COPY_VALUES,&B);
730: } else {
731: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
732: MatCopy(A,B,SAME_NONZERO_PATTERN);
733: }
734: MatISGetLocalMat(B,&lB);
735: MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows);
736: if (squaretest && haszerorows) {
738: MatCreateVecs(B,&x,&b);
739: MatDuplicate(B,MAT_COPY_VALUES,&B2);
740: VecSetLocalToGlobalMapping(b,l2gr);
741: VecSetLocalToGlobalMapping(x,l2gc);
742: VecSetRandom(x,NULL);
743: VecSetRandom(b,NULL);
744: /* mimic b[is] = x[is] */
745: VecDuplicate(b,&b2);
746: VecSetLocalToGlobalMapping(b2,l2gr);
747: VecCopy(b,b2);
748: ISGetLocalSize(is,&n);
749: ISGetIndices(is,&idxs);
750: VecGetSize(x,&N);
751: for (i=0;i<n;i++) {
752: if (0 <= idxs[i] && idxs[i] < N) {
753: VecSetValue(b2,idxs[i],diag,INSERT_VALUES);
754: VecSetValue(x,idxs[i],1.,INSERT_VALUES);
755: }
756: }
757: VecAssemblyBegin(b2);
758: VecAssemblyEnd(b2);
759: VecAssemblyBegin(x);
760: VecAssemblyEnd(x);
761: ISRestoreIndices(is,&idxs);
762: /* test ZeroRows on MATIS */
763: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
764: MatZeroRowsIS(B,is,diag,x,b);
765: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);
766: MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);
767: } else if (haszerorows) {
768: /* test ZeroRows on MATIS */
769: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
770: MatZeroRowsIS(B,is,diag,NULL,NULL);
771: b = b2 = x = NULL;
772: } else {
773: PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr);
774: b = b2 = x = NULL;
775: }
777: if (squaretest && haszerorows) {
778: VecAXPY(b2,-1.,b);
779: VecNorm(b2,NORM_INFINITY,&error);
780: if (error > PETSC_SQRT_MACHINE_EPSILON) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR IN ZEROROWS ON B %g (diag %s)",error,diagstr);
781: }
783: /* test MatMissingDiagonal */
784: PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");
785: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
786: MatMissingDiagonal(B,&miss,&d);
787: MatGetOwnershipRange(B,&rst,&ren);
788: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
789: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%D,%D) Missing %d, row %D (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);
790: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
791: PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
793: VecDestroy(&x);
794: VecDestroy(&b);
795: VecDestroy(&b2);
797: /* check the result of ZeroRows with that from MPIAIJ routines
798: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
799: if (haszerorows) {
800: MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);
801: MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
802: MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);
803: CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");
804: MatDestroy(&Bcheck);
805: }
806: MatDestroy(&B);
808: if (B2) { /* test MatZeroRowsColumns */
809: MatDuplicate(Afull,MAT_COPY_VALUES,&B);
810: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
811: MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);
812: CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");
813: MatDestroy(&B);
814: MatDestroy(&B2);
815: }
816: return(0);
817: }
820: /*TEST
822: test:
823: args: -test_trans
825: test:
826: suffix: 2
827: nsize: 4
828: args: -matis_convert_local_nest -nr 3 -nc 4
830: test:
831: suffix: 3
832: nsize: 5
833: args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1
835: test:
836: suffix: 4
837: nsize: 6
838: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
840: test:
841: suffix: 5
842: nsize: 6
843: args: -m 12 -n 12 -test_trans -nr 3 -nc 1
845: test:
846: suffix: 6
847: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
849: test:
850: suffix: 7
851: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
853: test:
854: suffix: 8
855: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
857: test:
858: suffix: 9
859: nsize: 5
860: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
862: test:
863: suffix: 10
864: nsize: 5
865: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
867: test:
868: suffix: vscat_default
869: nsize: 5
870: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
871: output_file: output/ex23_11.out
873: test:
874: suffix: 12
875: nsize: 3
876: args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3
878: testset:
879: output_file: output/ex23_13.out
880: nsize: 3
881: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
882: filter: grep -v "type:"
883: test:
884: suffix: baij
885: args: -matis_localmat_type baij
886: test:
887: requires: viennacl
888: suffix: viennacl
889: args: -matis_localmat_type aijviennacl
890: test:
891: requires: cuda
892: suffix: cusparse
893: args: -matis_localmat_type aijcusparse
895: TEST*/