Actual source code: ex23.c
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,lm,ln;
26: PetscInt rst,ren,cst,cen,nr,nc;
27: PetscMPIInt rank,size;
28: PetscBool testT,squaretest,isaij;
29: PetscBool permute = PETSC_FALSE, negmap = PETSC_FALSE, repmap = PETSC_FALSE;
30: PetscBool diffmap = PETSC_TRUE, symmetric = PETSC_FALSE, issymmetric;
31: PetscErrorCode ierr;
33: PetscInitialize(&argc,&args,(char*)0,help);
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: PetscOptionsGetBool(NULL,NULL,"-negmap",&negmap,NULL);
41: PetscOptionsGetBool(NULL,NULL,"-repmap",&repmap,NULL);
42: PetscOptionsGetBool(NULL,NULL,"-permmap",&permute,NULL);
43: PetscOptionsGetBool(NULL,NULL,"-diffmap",&diffmap,NULL);
48: /* create a MATIS matrix */
49: MatCreate(PETSC_COMM_WORLD,&A);
50: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
51: MatSetType(A,MATIS);
52: MatSetFromOptions(A);
53: if (!negmap && !repmap) {
54: /* This is not the proper setting for MATIS for finite elements, it is just used to test the routines
55: Here we use a one-to-one correspondence between local row/column spaces and global row/column spaces
56: Equivalent to passing NULL for the mapping */
57: ISCreateStride(PETSC_COMM_WORLD,n,0,1,&is);
58: } else if (negmap && !repmap) { /* non repeated but with negative indices */
59: ISCreateStride(PETSC_COMM_WORLD,n+2,-2,1,&is);
60: } else if (!negmap && repmap) { /* non negative but repeated indices */
61: IS isl[2];
63: ISCreateStride(PETSC_COMM_WORLD,n,0,1,&isl[0]);
64: ISCreateStride(PETSC_COMM_WORLD,n,n-1,-1,&isl[1]);
65: ISConcatenate(PETSC_COMM_WORLD,2,isl,&is);
66: ISDestroy(&isl[0]);
67: ISDestroy(&isl[1]);
68: } else { /* negative and repeated indices */
69: IS isl[2];
71: ISCreateStride(PETSC_COMM_WORLD,n+1,-1,1,&isl[0]);
72: ISCreateStride(PETSC_COMM_WORLD,n+1,n-1,-1,&isl[1]);
73: ISConcatenate(PETSC_COMM_WORLD,2,isl,&is);
74: ISDestroy(&isl[0]);
75: ISDestroy(&isl[1]);
76: }
77: ISLocalToGlobalMappingCreateIS(is,&cmap);
78: ISDestroy(&is);
80: if (m != n || diffmap) {
81: ISCreateStride(PETSC_COMM_WORLD,m,permute ? m-1 : 0,permute ? -1 : 1,&is);
82: ISLocalToGlobalMappingCreateIS(is,&rmap);
83: ISDestroy(&is);
84: } else {
85: PetscObjectReference((PetscObject)cmap);
86: rmap = cmap;
87: }
89: MatSetLocalToGlobalMapping(A,rmap,cmap);
90: MatISStoreL2L(A,PETSC_FALSE);
91: MatISSetPreallocation(A,3,NULL,3,NULL);
92: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap)); /* I do not want to precompute the pattern */
93: ISLocalToGlobalMappingGetSize(rmap,&lm);
94: ISLocalToGlobalMappingGetSize(cmap,&ln);
95: for (i=0; i<lm; i++) {
96: PetscScalar v[3];
97: PetscInt cols[3];
99: cols[0] = (i-1+n)%n;
100: cols[1] = i%n;
101: cols[2] = (i+1)%n;
102: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
103: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
104: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
105: ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols);
106: MatSetValuesLocal(A,1,&i,3,cols,v,ADD_VALUES);
107: }
108: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
111: /* activate tests for square matrices with same maps only */
112: MatHasCongruentLayouts(A,&squaretest);
113: if (squaretest && rmap != cmap) {
114: PetscInt nr, nc;
116: ISLocalToGlobalMappingGetSize(rmap,&nr);
117: ISLocalToGlobalMappingGetSize(cmap,&nc);
118: if (nr != nc) squaretest = PETSC_FALSE;
119: else {
120: const PetscInt *idxs1,*idxs2;
122: ISLocalToGlobalMappingGetIndices(rmap,&idxs1);
123: ISLocalToGlobalMappingGetIndices(cmap,&idxs2);
124: PetscArraycmp(idxs1,idxs2,nr,&squaretest);
125: ISLocalToGlobalMappingRestoreIndices(rmap,&idxs1);
126: ISLocalToGlobalMappingRestoreIndices(cmap,&idxs2);
127: }
128: MPIU_Allreduce(MPI_IN_PLACE,&squaretest,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
129: }
131: /* test MatISGetLocalMat */
132: MatISGetLocalMat(A,&B);
133: MatGetType(B,&lmtype);
135: /* test MatGetInfo */
136: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetInfo\n");
137: MatGetInfo(A,MAT_LOCAL,&info);
138: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
139: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD,"Process %2d: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",PetscGlobalRank,(PetscInt)info.nz_used,
140: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
141: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
142: MatGetInfo(A,MAT_GLOBAL_MAX,&info);
143: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalMax : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used,
144: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
145: MatGetInfo(A,MAT_GLOBAL_SUM,&info);
146: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"GlobalSum : %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",(PetscInt)info.nz_used,
147: (PetscInt)info.nz_allocated,(PetscInt)info.nz_unneeded,(PetscInt)info.assemblies,(PetscInt)info.mallocs);
149: /* test MatIsSymmetric */
150: MatIsSymmetric(A,0.0,&issymmetric);
151: PetscPrintf(PETSC_COMM_WORLD,"Test MatIsSymmetric: %d\n",issymmetric);
153: /* Create a MPIAIJ matrix, same as A */
154: MatCreate(PETSC_COMM_WORLD,&B);
155: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
156: MatSetType(B,MATAIJ);
157: MatSetFromOptions(B);
158: MatSetUp(B);
159: MatSetLocalToGlobalMapping(B,rmap,cmap);
160: MatMPIAIJSetPreallocation(B,3,NULL,3,NULL);
161: MatMPIBAIJSetPreallocation(B,1,3,NULL,3,NULL);
162: #if defined(PETSC_HAVE_HYPRE)
163: MatHYPRESetPreallocation(B,3,NULL,3,NULL);
164: #endif
165: MatISSetPreallocation(B,3,NULL,3,NULL);
166: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,(PetscBool)!(repmap || negmap)); /* I do not want to precompute the pattern */
167: for (i=0; i<lm; i++) {
168: PetscScalar v[3];
169: PetscInt cols[3];
171: cols[0] = (i-1+n)%n;
172: cols[1] = i%n;
173: cols[2] = (i+1)%n;
174: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
175: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
176: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
177: ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols);
178: MatSetValuesLocal(B,1,&i,3,cols,v,ADD_VALUES);
179: }
180: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
181: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
183: /* test MatView */
184: PetscPrintf(PETSC_COMM_WORLD,"Test MatView\n");
185: MatView(A,NULL);
186: MatView(B,NULL);
188: /* test CheckMat */
189: PetscPrintf(PETSC_COMM_WORLD,"Test CheckMat\n");
190: CheckMat(A,B,PETSC_FALSE,"CheckMat");
192: /* test MatDuplicate and MatAXPY */
193: PetscPrintf(PETSC_COMM_WORLD,"Test MatDuplicate and MatAXPY\n");
194: MatDuplicate(A,MAT_COPY_VALUES,&A2);
195: CheckMat(A,A2,PETSC_FALSE,"MatDuplicate and MatAXPY");
197: /* test MatConvert */
198: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_IS_XAIJ\n");
199: MatConvert(A2,MATAIJ,MAT_INITIAL_MATRIX,&B2);
200: CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INITIAL_MATRIX");
201: MatConvert(A2,MATAIJ,MAT_REUSE_MATRIX,&B2);
202: CheckMat(B,B2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_REUSE_MATRIX");
203: MatConvert(A2,MATAIJ,MAT_INPLACE_MATRIX,&A2);
204: CheckMat(B,A2,PETSC_TRUE,"MatConvert_IS_XAIJ MAT_INPLACE_MATRIX");
205: MatDestroy(&A2);
206: MatDestroy(&B2);
207: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_XAIJ_IS\n");
208: MatDuplicate(B,MAT_COPY_VALUES,&B2);
209: MatConvert(B2,MATIS,MAT_INITIAL_MATRIX,&A2);
210: CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INITIAL_MATRIX");
211: MatConvert(B2,MATIS,MAT_REUSE_MATRIX,&A2);
212: CheckMat(A,A2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_REUSE_MATRIX");
213: MatConvert(B2,MATIS,MAT_INPLACE_MATRIX,&B2);
214: CheckMat(A,B2,PETSC_TRUE,"MatConvert_XAIJ_IS MAT_INPLACE_MATRIX");
215: MatDestroy(&A2);
216: MatDestroy(&B2);
217: PetscStrcmp(lmtype,MATSEQAIJ,&isaij);
218: if (size == 1 && isaij) { /* tests special code paths in MatConvert_IS_XAIJ */
219: 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};
220: ISLocalToGlobalMapping tcmap,trmap;
222: for (ri = 0; ri < 2; ri++) {
223: PetscInt *r;
225: r = (PetscInt*)(ri == 0 ? rr : rk);
226: for (ci = 0; ci < 2; ci++) {
227: PetscInt *c,rb,cb;
229: c = (PetscInt*)(ci == 0 ? cr : ck);
230: for (rb = 1; rb < 4; rb++) {
231: ISCreateBlock(PETSC_COMM_SELF,rb,3,r,PETSC_COPY_VALUES,&is);
232: ISLocalToGlobalMappingCreateIS(is,&trmap);
233: ISDestroy(&is);
234: for (cb = 1; cb < 4; cb++) {
235: Mat T,lT,T2;
236: char testname[256];
238: PetscSNPrintf(testname,sizeof(testname),"MatConvert_IS_XAIJ special case (%" PetscInt_FMT " %" PetscInt_FMT ", bs %" PetscInt_FMT " %" PetscInt_FMT ")",ri,ci,rb,cb);
239: PetscPrintf(PETSC_COMM_WORLD,"Test %s\n",testname);
241: ISCreateBlock(PETSC_COMM_SELF,cb,4,c,PETSC_COPY_VALUES,&is);
242: ISLocalToGlobalMappingCreateIS(is,&tcmap);
243: ISDestroy(&is);
245: MatCreate(PETSC_COMM_SELF,&T);
246: MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,rb*3,cb*4);
247: MatSetType(T,MATIS);
248: MatSetLocalToGlobalMapping(T,trmap,tcmap);
249: ISLocalToGlobalMappingDestroy(&tcmap);
250: MatISGetLocalMat(T,&lT);
251: MatSetType(lT,MATSEQAIJ);
252: MatSeqAIJSetPreallocation(lT,cb*4,NULL);
253: MatSetRandom(lT,NULL);
254: MatConvert(lT,lmtype,MAT_INPLACE_MATRIX,&lT);
255: MatISRestoreLocalMat(T,&lT);
256: MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
257: MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);
259: MatConvert(T,MATAIJ,MAT_INITIAL_MATRIX,&T2);
260: CheckMat(T,T2,PETSC_TRUE,"MAT_INITIAL_MATRIX");
261: MatConvert(T,MATAIJ,MAT_REUSE_MATRIX,&T2);
262: CheckMat(T,T2,PETSC_TRUE,"MAT_REUSE_MATRIX");
263: MatDestroy(&T2);
264: MatDuplicate(T,MAT_COPY_VALUES,&T2);
265: MatConvert(T2,MATAIJ,MAT_INPLACE_MATRIX,&T2);
266: CheckMat(T,T2,PETSC_TRUE,"MAT_INPLACE_MATRIX");
267: MatDestroy(&T);
268: MatDestroy(&T2);
269: }
270: ISLocalToGlobalMappingDestroy(&trmap);
271: }
272: }
273: }
274: }
276: /* test MatDiagonalScale */
277: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalScale\n");
278: MatDuplicate(A,MAT_COPY_VALUES,&A2);
279: MatDuplicate(B,MAT_COPY_VALUES,&B2);
280: MatCreateVecs(A,&x,&y);
281: VecSetRandom(x,NULL);
282: if (issymmetric) {
283: VecCopy(x,y);
284: } else {
285: VecSetRandom(y,NULL);
286: VecScale(y,8.);
287: }
288: MatDiagonalScale(A2,y,x);
289: MatDiagonalScale(B2,y,x);
290: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalScale");
291: MatDestroy(&A2);
292: MatDestroy(&B2);
293: VecDestroy(&x);
294: VecDestroy(&y);
296: /* test MatPtAP (A IS and B AIJ) */
297: if (isaij && m == n) {
298: PetscPrintf(PETSC_COMM_WORLD,"Test MatPtAP\n");
299: MatISStoreL2L(A,PETSC_TRUE);
300: MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A2);
301: MatPtAP(B,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B2);
302: CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_INITIAL_MATRIX");
303: MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&A2);
304: CheckMat(A2,B2,PETSC_FALSE,"MatPtAP MAT_REUSE_MATRIX");
305: MatDestroy(&A2);
306: MatDestroy(&B2);
307: }
309: /* test MatGetLocalSubMatrix */
310: PetscPrintf(PETSC_COMM_WORLD,"Test MatGetLocalSubMatrix\n");
311: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&A2);
312: ISCreateStride(PETSC_COMM_SELF,lm/2+lm%2,0,2,&reven);
313: ISComplement(reven,0,lm,&rodd);
314: ISCreateStride(PETSC_COMM_SELF,ln/2+ln%2,0,2,&ceven);
315: ISComplement(ceven,0,ln,&codd);
316: MatGetLocalSubMatrix(A2,reven,ceven,&Aee);
317: MatGetLocalSubMatrix(A2,reven,codd,&Aeo);
318: MatGetLocalSubMatrix(A2,rodd,ceven,&Aoe);
319: MatGetLocalSubMatrix(A2,rodd,codd,&Aoo);
320: for (i=0; i<lm; i++) {
321: PetscInt j,je,jo,colse[3], colso[3];
322: PetscScalar ve[3], vo[3];
323: PetscScalar v[3];
324: PetscInt cols[3];
325: PetscInt row = i/2;
327: cols[0] = (i-1+n)%n;
328: cols[1] = i%n;
329: cols[2] = (i+1)%n;
330: v[0] = -1.*(symmetric ? PetscMin(i+1,cols[0]+1) : i+1);
331: v[1] = 2.*(symmetric ? PetscMin(i+1,cols[1]+1) : i+1);
332: v[2] = -1.*(symmetric ? PetscMin(i+1,cols[2]+1) : i+1);
333: ISGlobalToLocalMappingApply(cmap,IS_GTOLM_MASK,3,cols,NULL,cols);
334: for (j=0,je=0,jo=0;j<3;j++) {
335: if (cols[j]%2) {
336: vo[jo] = v[j];
337: colso[jo++] = cols[j]/2;
338: } else {
339: ve[je] = v[j];
340: colse[je++] = cols[j]/2;
341: }
342: }
343: if (i%2) {
344: MatSetValuesLocal(Aoe,1,&row,je,colse,ve,ADD_VALUES);
345: MatSetValuesLocal(Aoo,1,&row,jo,colso,vo,ADD_VALUES);
346: } else {
347: MatSetValuesLocal(Aee,1,&row,je,colse,ve,ADD_VALUES);
348: MatSetValuesLocal(Aeo,1,&row,jo,colso,vo,ADD_VALUES);
349: }
350: }
351: MatRestoreLocalSubMatrix(A2,reven,ceven,&Aee);
352: MatRestoreLocalSubMatrix(A2,reven,codd,&Aeo);
353: MatRestoreLocalSubMatrix(A2,rodd,ceven,&Aoe);
354: MatRestoreLocalSubMatrix(A2,rodd,codd,&Aoo);
355: ISDestroy(&reven);
356: ISDestroy(&ceven);
357: ISDestroy(&rodd);
358: ISDestroy(&codd);
359: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
360: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
361: MatAXPY(A2,-1.,A,SAME_NONZERO_PATTERN);
362: CheckMat(A2,NULL,PETSC_FALSE,"MatGetLocalSubMatrix");
363: MatDestroy(&A2);
365: /* test MatConvert_Nest_IS */
366: testT = PETSC_FALSE;
367: PetscOptionsGetBool(NULL,NULL,"-test_trans",&testT,NULL);
369: PetscPrintf(PETSC_COMM_WORLD,"Test MatConvert_Nest_IS\n");
370: nr = 2;
371: nc = 2;
372: PetscOptionsGetInt(NULL,NULL,"-nr",&nr,NULL);
373: PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL);
374: if (testT) {
375: MatGetOwnershipRange(A,&cst,&cen);
376: MatGetOwnershipRangeColumn(A,&rst,&ren);
377: } else {
378: MatGetOwnershipRange(A,&rst,&ren);
379: MatGetOwnershipRangeColumn(A,&cst,&cen);
380: }
381: PetscMalloc3(nr,&rows,nc,&cols,2*nr*nc,&mats);
382: for (i=0;i<nr*nc;i++) {
383: if (testT) {
384: MatCreateTranspose(A,&mats[i]);
385: MatTranspose(B,MAT_INITIAL_MATRIX,&mats[i+nr*nc]);
386: } else {
387: MatDuplicate(A,MAT_COPY_VALUES,&mats[i]);
388: MatDuplicate(B,MAT_COPY_VALUES,&mats[i+nr*nc]);
389: }
390: }
391: for (i=0;i<nr;i++) {
392: ISCreateStride(PETSC_COMM_WORLD,ren-rst,i+rst,nr,&rows[i]);
393: }
394: for (i=0;i<nc;i++) {
395: ISCreateStride(PETSC_COMM_WORLD,cen-cst,i+cst,nc,&cols[i]);
396: }
397: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats,&A2);
398: MatCreateNest(PETSC_COMM_WORLD,nr,rows,nc,cols,mats+nr*nc,&B2);
399: for (i=0;i<nr;i++) {
400: ISDestroy(&rows[i]);
401: }
402: for (i=0;i<nc;i++) {
403: ISDestroy(&cols[i]);
404: }
405: for (i=0;i<2*nr*nc;i++) {
406: MatDestroy(&mats[i]);
407: }
408: PetscFree3(rows,cols,mats);
409: MatConvert(B2,MATAIJ,MAT_INITIAL_MATRIX,&T);
410: MatDestroy(&B2);
411: MatConvert(A2,MATIS,MAT_INITIAL_MATRIX,&B2);
412: CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INITIAL_MATRIX");
413: MatConvert(A2,MATIS,MAT_REUSE_MATRIX,&B2);
414: CheckMat(B2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_REUSE_MATRIX");
415: MatDestroy(&B2);
416: MatConvert(A2,MATIS,MAT_INPLACE_MATRIX,&A2);
417: CheckMat(A2,T,PETSC_TRUE,"MatConvert_Nest_IS MAT_INPLACE_MATRIX");
418: MatDestroy(&T);
419: MatDestroy(&A2);
421: /* test MatCreateSubMatrix */
422: PetscPrintf(PETSC_COMM_WORLD,"Test MatCreateSubMatrix\n");
423: if (rank == 0) {
424: ISCreateStride(PETSC_COMM_WORLD,1,1,1,&is);
425: ISCreateStride(PETSC_COMM_WORLD,2,0,1,&is2);
426: } else if (rank == 1) {
427: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
428: if (n > 3) {
429: ISCreateStride(PETSC_COMM_WORLD,1,3,1,&is2);
430: } else {
431: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
432: }
433: } else if (rank == 2 && n > 4) {
434: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
435: ISCreateStride(PETSC_COMM_WORLD,n-4,4,1,&is2);
436: } else {
437: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
438: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is2);
439: }
440: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&A2);
441: MatCreateSubMatrix(B,is,is,MAT_INITIAL_MATRIX,&B2);
442: CheckMat(A2,B2,PETSC_TRUE,"first MatCreateSubMatrix");
444: MatCreateSubMatrix(A,is,is,MAT_REUSE_MATRIX,&A2);
445: MatCreateSubMatrix(B,is,is,MAT_REUSE_MATRIX,&B2);
446: CheckMat(A2,B2,PETSC_FALSE,"reuse MatCreateSubMatrix");
447: MatDestroy(&A2);
448: MatDestroy(&B2);
450: if (!issymmetric) {
451: MatCreateSubMatrix(A,is,is2,MAT_INITIAL_MATRIX,&A2);
452: MatCreateSubMatrix(B,is,is2,MAT_INITIAL_MATRIX,&B2);
453: MatCreateSubMatrix(A,is,is2,MAT_REUSE_MATRIX,&A2);
454: MatCreateSubMatrix(B,is,is2,MAT_REUSE_MATRIX,&B2);
455: CheckMat(A2,B2,PETSC_FALSE,"second MatCreateSubMatrix");
456: }
458: MatDestroy(&A2);
459: MatDestroy(&B2);
460: ISDestroy(&is);
461: ISDestroy(&is2);
463: /* Create an IS required by MatZeroRows(): just rank zero provides the rows to be eliminated */
464: if (size > 1) {
465: if (rank == 0) {
466: PetscInt st,len;
468: st = (m+1)/2;
469: len = PetscMin(m/2,PetscMax(m-(m+1)/2-1,0));
470: ISCreateStride(PETSC_COMM_WORLD,len,st,1,&is);
471: } else {
472: ISCreateStride(PETSC_COMM_WORLD,0,0,1,&is);
473: }
474: } else {
475: ISCreateStride(PETSC_COMM_WORLD,1,0,1,&is);
476: }
478: if (squaretest) { /* tests for square matrices only, with same maps for rows and columns */
479: /* test MatDiagonalSet */
480: PetscPrintf(PETSC_COMM_WORLD,"Test MatDiagonalSet\n");
481: MatDuplicate(A,MAT_COPY_VALUES,&A2);
482: MatDuplicate(B,MAT_COPY_VALUES,&B2);
483: MatCreateVecs(A,NULL,&x);
484: VecSetRandom(x,NULL);
485: MatDiagonalSet(A2,x,INSERT_VALUES);
486: MatDiagonalSet(B2,x,INSERT_VALUES);
487: CheckMat(A2,B2,PETSC_FALSE,"MatDiagonalSet");
488: VecDestroy(&x);
489: MatDestroy(&A2);
490: MatDestroy(&B2);
492: /* test MatShift (MatShift_IS internally uses MatDiagonalSet_IS with ADD_VALUES) */
493: PetscPrintf(PETSC_COMM_WORLD,"Test MatShift\n");
494: MatDuplicate(A,MAT_COPY_VALUES,&A2);
495: MatDuplicate(B,MAT_COPY_VALUES,&B2);
496: MatShift(A2,2.0);
497: MatShift(B2,2.0);
498: CheckMat(A2,B2,PETSC_FALSE,"MatShift");
499: MatDestroy(&A2);
500: MatDestroy(&B2);
502: /* nonzero diag value is supported for square matrices only */
503: TestMatZeroRows(A,B,PETSC_TRUE,is,diag);
504: }
505: TestMatZeroRows(A,B,squaretest,is,0.0);
506: ISDestroy(&is);
508: /* test MatTranspose */
509: PetscPrintf(PETSC_COMM_WORLD,"Test MatTranspose\n");
510: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
511: MatTranspose(B,MAT_INITIAL_MATRIX,&B2);
512: CheckMat(A2,B2,PETSC_FALSE,"initial matrix MatTranspose");
514: MatTranspose(A,MAT_REUSE_MATRIX,&A2);
515: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (not in place) MatTranspose");
516: MatDestroy(&A2);
518: MatDuplicate(A,MAT_COPY_VALUES,&A2);
519: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
520: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (in place) MatTranspose");
521: MatDestroy(&A2);
523: MatTranspose(A,MAT_INITIAL_MATRIX,&A2);
524: CheckMat(A2,B2,PETSC_FALSE,"reuse matrix (different type) MatTranspose");
525: MatDestroy(&A2);
526: MatDestroy(&B2);
528: /* test MatISFixLocalEmpty */
529: if (isaij) {
530: PetscInt r[2];
532: r[0] = 0;
533: r[1] = PetscMin(m,n)-1;
534: PetscPrintf(PETSC_COMM_WORLD,"Test MatISFixLocalEmpty\n");
535: MatDuplicate(A,MAT_COPY_VALUES,&A2);
537: MatISFixLocalEmpty(A2,PETSC_TRUE);
538: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
539: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
540: CheckMat(A2,B,PETSC_FALSE,"MatISFixLocalEmpty (null)");
542: MatZeroRows(A2,2,r,0.0,NULL,NULL);
543: MatViewFromOptions(A2,NULL,"-fixempty_view");
544: MatDuplicate(B,MAT_COPY_VALUES,&B2);
545: MatZeroRows(B2,2,r,0.0,NULL,NULL);
546: MatISFixLocalEmpty(A2,PETSC_TRUE);
547: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
548: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
549: MatViewFromOptions(A2,NULL,"-fixempty_view");
550: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows)");
551: MatDestroy(&A2);
553: MatDuplicate(A,MAT_COPY_VALUES,&A2);
554: MatZeroRows(A2,2,r,0.0,NULL,NULL);
555: MatTranspose(A2,MAT_INPLACE_MATRIX,&A2);
556: MatTranspose(B2,MAT_INPLACE_MATRIX,&B2);
557: MatViewFromOptions(A2,NULL,"-fixempty_view");
558: MatISFixLocalEmpty(A2,PETSC_TRUE);
559: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
560: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
561: MatViewFromOptions(A2,NULL,"-fixempty_view");
562: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (cols)");
564: MatDestroy(&A2);
565: MatDestroy(&B2);
567: if (squaretest) {
568: MatDuplicate(A,MAT_COPY_VALUES,&A2);
569: MatDuplicate(B,MAT_COPY_VALUES,&B2);
570: MatZeroRowsColumns(A2,2,r,0.0,NULL,NULL);
571: MatZeroRowsColumns(B2,2,r,0.0,NULL,NULL);
572: MatViewFromOptions(A2,NULL,"-fixempty_view");
573: MatISFixLocalEmpty(A2,PETSC_TRUE);
574: MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY);
575: MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY);
576: MatViewFromOptions(A2,NULL,"-fixempty_view");
577: CheckMat(A2,B2,PETSC_FALSE,"MatISFixLocalEmpty (rows+cols)");
578: MatDestroy(&A2);
579: MatDestroy(&B2);
580: }
582: }
584: /* test MatInvertBlockDiagonal
585: special cases for block-diagonal matrices */
586: if (m == n) {
587: ISLocalToGlobalMapping map;
588: Mat Abd,Bbd;
589: IS is,bis;
590: const PetscScalar *isbd,*aijbd;
591: PetscScalar *vals;
592: const PetscInt *sts,*idxs;
593: PetscInt *idxs2,diff,perm,nl,bs,st,en,in;
594: PetscBool ok;
596: for (diff = 0; diff < 3; diff++) {
597: for (perm = 0; perm < 3; perm++) {
598: for (bs = 1; bs < 4; bs++) {
599: PetscPrintf(PETSC_COMM_WORLD,"Test MatInvertBlockDiagonal blockdiag %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",n,diff,perm,bs);
600: PetscMalloc1(bs*bs,&vals);
601: MatGetOwnershipRanges(A,&sts);
602: switch (diff) {
603: case 1: /* inverted layout by processes */
604: in = 1;
605: st = sts[size - rank - 1];
606: en = sts[size - rank];
607: nl = en - st;
608: break;
609: case 2: /* round-robin layout */
610: in = size;
611: st = rank;
612: nl = n/size;
613: if (rank < n%size) nl++;
614: break;
615: default: /* same layout */
616: in = 1;
617: st = sts[rank];
618: en = sts[rank + 1];
619: nl = en - st;
620: break;
621: }
622: ISCreateStride(PETSC_COMM_WORLD,nl,st,in,&is);
623: ISGetLocalSize(is,&nl);
624: ISGetIndices(is,&idxs);
625: PetscMalloc1(nl,&idxs2);
626: for (i=0;i<nl;i++) {
627: switch (perm) { /* invert some of the indices */
628: case 2:
629: idxs2[i] = rank%2 ? idxs[i] : idxs[nl-i-1];
630: break;
631: case 1:
632: idxs2[i] = rank%2 ? idxs[nl-i-1] : idxs[i];
633: break;
634: default:
635: idxs2[i] = idxs[i];
636: break;
637: }
638: }
639: ISRestoreIndices(is,&idxs);
640: ISCreateBlock(PETSC_COMM_WORLD,bs,nl,idxs2,PETSC_OWN_POINTER,&bis);
641: ISLocalToGlobalMappingCreateIS(bis,&map);
642: MatCreateIS(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,bs*n,bs*n,map,map,&Abd);
643: ISLocalToGlobalMappingDestroy(&map);
644: MatISSetPreallocation(Abd,bs,NULL,0,NULL);
645: for (i=0;i<nl;i++) {
646: PetscInt b1,b2;
648: for (b1=0;b1<bs;b1++) {
649: for (b2=0;b2<bs;b2++) {
650: vals[b1*bs + b2] = i*bs*bs + b1*bs + b2 + 1 + (b1 == b2 ? 1.0 : 0);
651: }
652: }
653: MatSetValuesBlockedLocal(Abd,1,&i,1,&i,vals,INSERT_VALUES);
654: }
655: MatAssemblyBegin(Abd,MAT_FINAL_ASSEMBLY);
656: MatAssemblyEnd(Abd,MAT_FINAL_ASSEMBLY);
657: MatConvert(Abd,MATAIJ,MAT_INITIAL_MATRIX,&Bbd);
658: MatInvertBlockDiagonal(Abd,&isbd);
659: MatInvertBlockDiagonal(Bbd,&aijbd);
660: MatGetLocalSize(Bbd,&nl,NULL);
661: ok = PETSC_TRUE;
662: for (i=0;i<nl/bs;i++) {
663: PetscInt b1,b2;
665: for (b1=0;b1<bs;b1++) {
666: for (b2=0;b2<bs;b2++) {
667: if (PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]-aijbd[i*bs*bs+b1*bs + b2]) > PETSC_SMALL) ok = PETSC_FALSE;
668: if (!ok) {
669: PetscPrintf(PETSC_COMM_SELF,"[%d] ERROR block %" PetscInt_FMT ", entry %" PetscInt_FMT " %" PetscInt_FMT ": %g %g\n",rank,i,b1,b2,(double)PetscAbsScalar(isbd[i*bs*bs+b1*bs + b2]),(double)PetscAbsScalar(aijbd[i*bs*bs+b1*bs + b2]));
670: break;
671: }
672: }
673: if (!ok) break;
674: }
675: if (!ok) break;
676: }
677: MatDestroy(&Abd);
678: MatDestroy(&Bbd);
679: PetscFree(vals);
680: ISDestroy(&is);
681: ISDestroy(&bis);
682: }
683: }
684: }
685: }
686: /* free testing matrices */
687: ISLocalToGlobalMappingDestroy(&cmap);
688: ISLocalToGlobalMappingDestroy(&rmap);
689: MatDestroy(&A);
690: MatDestroy(&B);
691: PetscFinalize();
692: return 0;
693: }
695: PetscErrorCode CheckMat(Mat A, Mat B, PetscBool usemult, const char* func)
696: {
697: Mat Bcheck;
698: PetscReal error;
701: if (!usemult && B) {
702: PetscBool hasnorm;
704: MatHasOperation(B,MATOP_NORM,&hasnorm);
705: if (!hasnorm) usemult = PETSC_TRUE;
706: }
707: if (!usemult) {
708: if (B) {
709: MatType Btype;
711: MatGetType(B,&Btype);
712: MatConvert(A,Btype,MAT_INITIAL_MATRIX,&Bcheck);
713: } else {
714: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
715: }
716: if (B) { /* if B is present, subtract it */
717: MatAXPY(Bcheck,-1.,B,DIFFERENT_NONZERO_PATTERN);
718: }
719: MatNorm(Bcheck,NORM_INFINITY,&error);
720: if (error > PETSC_SQRT_MACHINE_EPSILON) {
721: ISLocalToGlobalMapping rl2g,cl2g;
723: PetscObjectSetName((PetscObject)Bcheck,"Assembled Bcheck");
724: MatView(Bcheck,NULL);
725: if (B) {
726: PetscObjectSetName((PetscObject)B,"Assembled AIJ");
727: MatView(B,NULL);
728: MatDestroy(&Bcheck);
729: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&Bcheck);
730: PetscObjectSetName((PetscObject)Bcheck,"Assembled IS");
731: MatView(Bcheck,NULL);
732: }
733: MatDestroy(&Bcheck);
734: PetscObjectSetName((PetscObject)A,"MatIS");
735: MatView(A,NULL);
736: MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);
737: ISLocalToGlobalMappingView(rl2g,NULL);
738: ISLocalToGlobalMappingView(cl2g,NULL);
739: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"ERROR ON %s: %g",func,(double)error);
740: }
741: MatDestroy(&Bcheck);
742: } else {
743: PetscBool ok,okt;
745: MatMultEqual(A,B,3,&ok);
746: MatMultTransposeEqual(A,B,3,&okt);
748: }
749: return 0;
750: }
752: PetscErrorCode TestMatZeroRows(Mat A, Mat Afull, PetscBool squaretest, IS is, PetscScalar diag)
753: {
754: Mat B,Bcheck,B2 = NULL,lB;
755: Vec x = NULL, b = NULL, b2 = NULL;
756: ISLocalToGlobalMapping l2gr,l2gc;
757: PetscReal error;
758: char diagstr[16];
759: const PetscInt *idxs;
760: PetscInt rst,ren,i,n,N,d;
761: PetscMPIInt rank;
762: PetscBool miss,haszerorows;
765: if (diag == 0.) {
766: PetscStrcpy(diagstr,"zero");
767: } else {
768: PetscStrcpy(diagstr,"nonzero");
769: }
770: ISView(is,NULL);
771: MatGetLocalToGlobalMapping(A,&l2gr,&l2gc);
772: /* tests MatDuplicate and MatCopy */
773: if (diag == 0.) {
774: MatDuplicate(A,MAT_COPY_VALUES,&B);
775: } else {
776: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&B);
777: MatCopy(A,B,SAME_NONZERO_PATTERN);
778: }
779: MatISGetLocalMat(B,&lB);
780: MatHasOperation(lB,MATOP_ZERO_ROWS,&haszerorows);
781: if (squaretest && haszerorows) {
783: MatCreateVecs(B,&x,&b);
784: MatDuplicate(B,MAT_COPY_VALUES,&B2);
785: VecSetLocalToGlobalMapping(b,l2gr);
786: VecSetLocalToGlobalMapping(x,l2gc);
787: VecSetRandom(x,NULL);
788: VecSetRandom(b,NULL);
789: /* mimic b[is] = x[is] */
790: VecDuplicate(b,&b2);
791: VecSetLocalToGlobalMapping(b2,l2gr);
792: VecCopy(b,b2);
793: ISGetLocalSize(is,&n);
794: ISGetIndices(is,&idxs);
795: VecGetSize(x,&N);
796: for (i=0;i<n;i++) {
797: if (0 <= idxs[i] && idxs[i] < N) {
798: VecSetValue(b2,idxs[i],diag,INSERT_VALUES);
799: VecSetValue(x,idxs[i],1.,INSERT_VALUES);
800: }
801: }
802: VecAssemblyBegin(b2);
803: VecAssemblyEnd(b2);
804: VecAssemblyBegin(x);
805: VecAssemblyEnd(x);
806: ISRestoreIndices(is,&idxs);
807: /* test ZeroRows on MATIS */
808: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
809: MatZeroRowsIS(B,is,diag,x,b);
810: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRowsColumns (diag %s)\n",diagstr);
811: MatZeroRowsColumnsIS(B2,is,diag,NULL,NULL);
812: } else if (haszerorows) {
813: /* test ZeroRows on MATIS */
814: PetscPrintf(PETSC_COMM_WORLD,"Test MatZeroRows (diag %s)\n",diagstr);
815: MatZeroRowsIS(B,is,diag,NULL,NULL);
816: b = b2 = x = NULL;
817: } else {
818: PetscPrintf(PETSC_COMM_WORLD,"Skipping MatZeroRows (diag %s)\n",diagstr);
819: b = b2 = x = NULL;
820: }
822: if (squaretest && haszerorows) {
823: VecAXPY(b2,-1.,b);
824: VecNorm(b2,NORM_INFINITY,&error);
826: }
828: /* test MatMissingDiagonal */
829: PetscPrintf(PETSC_COMM_WORLD,"Test MatMissingDiagonal\n");
830: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
831: MatMissingDiagonal(B,&miss,&d);
832: MatGetOwnershipRange(B,&rst,&ren);
833: PetscViewerASCIIPushSynchronized(PETSC_VIEWER_STDOUT_WORLD);
834: PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_WORLD, "[%d] [%" PetscInt_FMT ",%" PetscInt_FMT ") Missing %d, row %" PetscInt_FMT " (diag %s)\n",rank,rst,ren,(int)miss,d,diagstr);
835: PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD);
836: PetscViewerASCIIPopSynchronized(PETSC_VIEWER_STDOUT_WORLD);
838: VecDestroy(&x);
839: VecDestroy(&b);
840: VecDestroy(&b2);
842: /* check the result of ZeroRows with that from MPIAIJ routines
843: assuming that MatConvert_IS_XAIJ and MatZeroRows_MPIAIJ work fine */
844: if (haszerorows) {
845: MatDuplicate(Afull,MAT_COPY_VALUES,&Bcheck);
846: MatSetOption(Bcheck,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
847: MatZeroRowsIS(Bcheck,is,diag,NULL,NULL);
848: CheckMat(B,Bcheck,PETSC_FALSE,"Zerorows");
849: MatDestroy(&Bcheck);
850: }
851: MatDestroy(&B);
853: if (B2) { /* test MatZeroRowsColumns */
854: MatDuplicate(Afull,MAT_COPY_VALUES,&B);
855: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
856: MatZeroRowsColumnsIS(B,is,diag,NULL,NULL);
857: CheckMat(B2,B,PETSC_FALSE,"MatZeroRowsColumns");
858: MatDestroy(&B);
859: MatDestroy(&B2);
860: }
861: return 0;
862: }
864: /*TEST
866: test:
867: args: -test_trans
869: test:
870: suffix: 2
871: nsize: 4
872: args: -matis_convert_local_nest -nr 3 -nc 4
874: test:
875: suffix: 3
876: nsize: 5
877: args: -m 11 -n 10 -matis_convert_local_nest -nr 2 -nc 1
879: test:
880: suffix: 4
881: nsize: 6
882: args: -m 9 -n 12 -test_trans -nr 2 -nc 7
884: test:
885: suffix: 5
886: nsize: 6
887: args: -m 12 -n 12 -test_trans -nr 3 -nc 1
889: test:
890: suffix: 6
891: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
893: test:
894: suffix: 7
895: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
897: test:
898: suffix: 8
899: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
901: test:
902: suffix: 9
903: nsize: 5
904: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap
906: test:
907: suffix: 10
908: nsize: 5
909: args: -m 12 -n 12 -test_trans -nr 2 -nc 3 -diffmap -permmap
911: test:
912: suffix: vscat_default
913: nsize: 5
914: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -permmap
915: output_file: output/ex23_11.out
917: test:
918: suffix: 12
919: nsize: 3
920: args: -m 12 -n 12 -symmetric -matis_localmat_type sbaij -test_trans -nr 2 -nc 3
922: testset:
923: output_file: output/ex23_13.out
924: nsize: 3
925: args: -m 12 -n 17 -test_trans -nr 2 -nc 3 -diffmap -permmap
926: filter: grep -v "type:"
927: test:
928: suffix: baij
929: args: -matis_localmat_type baij
930: test:
931: requires: viennacl
932: suffix: viennacl
933: args: -matis_localmat_type aijviennacl
934: test:
935: requires: cuda
936: suffix: cusparse
937: args: -matis_localmat_type aijcusparse
939: test:
940: suffix: negrep
941: nsize: {{1 3}separate output}
942: args: -m {{5 7}separate output} -n {{5 7}separate output} -test_trans -nr 2 -nc 3 -negmap {{0 1}separate output} -repmap {{0 1}separate output} -permmap -diffmap {{0 1}separate output}
944: TEST*/