Actual source code: ex75.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  2: static char help[] = "Tests the vatious routines in MatMPISBAIJ format.\n";

  4: #include <petscmat.h>

  8: int main(int argc,char **args)
  9: {
 10:   Vec            x,y,u,s1,s2;
 11:   Mat            A,sA,sB;
 12:   PetscRandom    rctx;
 13:   PetscReal      r1,r2,rnorm,tol=1.e-10;
 14:   PetscScalar    one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
 15:   PetscInt       n,col[3],n1,block,row,i,j,i2,j2,Ii,J,rstart,rend,bs=1,mbs=16,d_nz=3,o_nz=3,prob=2;
 17:   PetscMPIInt    size,rank;
 18:   PetscBool      flg;
 19:   MatType        type;

 21:   PetscInitialize(&argc,&args,(char*)0,help);
 22:   PetscOptionsGetInt(NULL,"-mbs",&mbs,NULL);
 23:   PetscOptionsGetInt(NULL,"-bs",&bs,NULL);

 25:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 26:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 28:   n = mbs*bs;

 30:   /* Assemble MPISBAIJ matrix sA */
 31:   MatCreate(PETSC_COMM_WORLD,&sA);
 32:   MatSetSizes(sA,PETSC_DECIDE,PETSC_DECIDE,n,n);
 33:   MatSetType(sA,MATSBAIJ);
 34:   MatSetFromOptions(sA);
 35:   MatGetType(sA,&type);
 36:   /* printf(" mattype: %s\n",type); */
 37:   MatMPISBAIJSetPreallocation(sA,bs,d_nz,NULL,o_nz,NULL);
 38:   MatSeqSBAIJSetPreallocation(sA,bs,d_nz,NULL);
 39:   MatSetOption(sA,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 41:   if (bs == 1) {
 42:     if (prob == 1) { /* tridiagonal matrix */
 43:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 44:       for (i=1; i<n-1; i++) {
 45:         col[0] = i-1; col[1] = i; col[2] = i+1;
 46:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 47:       }
 48:       i       = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
 49:       value[0]= 0.1; value[1]=-1; value[2]=2;
 50:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

 52:       i        = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
 53:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
 54:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 55:     } else if (prob ==2) { /* matrix for the five point stencil */
 56:       n1 =  (int) PetscSqrtReal((PetscReal)n);
 57:       if (n1*n1 != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");

 59:       for (i=0; i<n1; i++) {
 60:         for (j=0; j<n1; j++) {
 61:           Ii = j + n1*i;
 62:           if (i>0)    {J = Ii - n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 63:           if (i<n1-1) {J = Ii + n1; MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 64:           if (j>0)    {J = Ii - 1;  MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 65:           if (j<n1-1) {J = Ii + 1;  MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
 66:           MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 67:         }
 68:       }
 69:     }
 70:     /* end of if (bs == 1) */
 71:   } else {  /* bs > 1 */
 72:     for (block=0; block<n/bs; block++) {
 73:       /* diagonal blocks */
 74:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 75:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 76:         col[0] = i-1; col[1] = i; col[2] = i+1;
 77:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 78:       }
 79:       i       = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
 80:       value[0]=-1.0; value[1]=4.0;
 81:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

 83:       i       = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
 84:       value[0]=4.0; value[1] = -1.0;
 85:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
 86:     }
 87:     /* off-diagonal blocks */
 88:     value[0]=-1.0;
 89:     for (i=0; i<(n/bs-1)*bs; i++) {
 90:       col[0]=i+bs;
 91:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
 92:       col[0]=i; row=i+bs;
 93:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
 94:     }
 95:   }
 96:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
 97:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);

 99:   /* Test MatView() */
100:   /*
101:   MatView(sA, PETSC_VIEWER_STDOUT_WORLD);
102:   MatView(sA, PETSC_VIEWER_DRAW_WORLD);
103:   */
104:   /* Assemble MPIBAIJ matrix A */
105:   MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,NULL,o_nz,NULL,&A);
106:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

108:   if (bs == 1) {
109:     if (prob == 1) { /* tridiagonal matrix */
110:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
111:       for (i=1; i<n-1; i++) {
112:         col[0] = i-1; col[1] = i; col[2] = i+1;
113:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
114:       }
115:       i       = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
116:       value[0]= 0.1; value[1]=-1; value[2]=2;
117:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);

119:       i        = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
120:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
121:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
122:     } else if (prob ==2) { /* matrix for the five point stencil */
123:       n1 = (int) PetscSqrtReal((PetscReal)n);
124:       for (i=0; i<n1; i++) {
125:         for (j=0; j<n1; j++) {
126:           Ii = j + n1*i;
127:           if (i>0)    {J = Ii - n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
128:           if (i<n1-1) {J = Ii + n1; MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
129:           if (j>0)    {J = Ii - 1;  MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
130:           if (j<n1-1) {J = Ii + 1;  MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);}
131:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
132:         }
133:       }
134:     }
135:     /* end of if (bs == 1) */
136:   } else {  /* bs > 1 */
137:     for (block=0; block<n/bs; block++) {
138:       /* diagonal blocks */
139:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
140:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
141:         col[0] = i-1; col[1] = i; col[2] = i+1;
142:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
143:       }
144:       i       = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
145:       value[0]=-1.0; value[1]=4.0;
146:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);

148:       i       = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
149:       value[0]=4.0; value[1] = -1.0;
150:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
151:     }
152:     /* off-diagonal blocks */
153:     value[0]=-1.0;
154:     for (i=0; i<(n/bs-1)*bs; i++) {
155:       col[0]=i+bs;
156:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
157:       col[0]=i; row=i+bs;
158:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
159:     }
160:   }
161:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
162:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

164:   /* Test MatGetSize(), MatGetLocalSize() */
165:   MatGetSize(sA, &i,&j); MatGetSize(A, &i2,&j2);
166:   i   -= i2; j -= j2;
167:   if (i || j) {
168:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetSize()\n",rank);
169:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
170:   }

172:   MatGetLocalSize(sA, &i,&j); MatGetLocalSize(A, &i2,&j2);
173:   i2  -= i; j2 -= j;
174:   if (i2 || j2) {
175:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatGetLocalSize()\n",rank);
176:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
177:   }

179:   /* vectors */
180:   /*--------------------*/
181:   /* i is obtained from MatGetLocalSize() */
182:   VecCreate(PETSC_COMM_WORLD,&x);
183:   VecSetSizes(x,i,PETSC_DECIDE);
184:   VecSetFromOptions(x);
185:   VecDuplicate(x,&y);
186:   VecDuplicate(x,&u);
187:   VecDuplicate(x,&s1);
188:   VecDuplicate(x,&s2);

190:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
191:   PetscRandomSetFromOptions(rctx);
192:   VecSetRandom(x,rctx);
193:   VecSet(u,one);

195:   /* Test MatNorm() */
196:   MatNorm(A,NORM_FROBENIUS,&r1);
197:   MatNorm(sA,NORM_FROBENIUS,&r2);
198:   rnorm = PetscAbsScalar(r1-r2)/r2;
199:   if (rnorm > tol && !rank) {
200:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
201:   }
202:   MatNorm(A,NORM_INFINITY,&r1);
203:   MatNorm(sA,NORM_INFINITY,&r2);
204:   rnorm = PetscAbsScalar(r1-r2)/r2;
205:   if (rnorm > tol && !rank) {
206:     PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_INFINITY(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
207:   }
208:   MatNorm(A,NORM_1,&r1);
209:   MatNorm(sA,NORM_1,&r2);
210:   rnorm = PetscAbsScalar(r1-r2)/r2;
211:   if (rnorm > tol && !rank) {
212:     PetscPrintf(PETSC_COMM_WORLD,"Error: MatNorm_1(), Anorm=%16.14e, sAnorm=%16.14e bs=%D\n",r1,r2,bs);
213:   }

215:   /* Test MatGetOwnershipRange() */
216:   MatGetOwnershipRange(sA,&rstart,&rend);
217:   MatGetOwnershipRange(A,&i2,&j2);
218:   i2  -= rstart; j2 -= rend;
219:   if (i2 || j2) {
220:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
221:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
222:   }

224:   /* Test MatDiagonalScale() */
225:   MatDiagonalScale(A,x,x);
226:   MatDiagonalScale(sA,x,x);
227:   MatMultEqual(A,sA,10,&flg);
228:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");

230:   /* Test MatGetDiagonal(), MatScale() */
231:   MatGetDiagonal(A,s1);
232:   MatGetDiagonal(sA,s2);
233:   VecNorm(s1,NORM_1,&r1);
234:   VecNorm(s2,NORM_1,&r2);
235:   r1  -= r2;
236:   if (r1<-tol || r1>tol) {
237:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,(double)r1);
238:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
239:   }

241:   MatScale(A,alpha);
242:   MatScale(sA,alpha);

244:   /* Test MatGetRowMaxAbs() */
245:   MatGetRowMaxAbs(A,s1,NULL);
246:   MatGetRowMaxAbs(sA,s2,NULL);

248:   VecNorm(s1,NORM_1,&r1);
249:   VecNorm(s2,NORM_1,&r2);
250:   r1  -= r2;
251:   if (r1<-tol || r1>tol) {
252:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMaxAbs() \n");
253:   }

255:   /* Test MatMult(), MatMultAdd() */
256:   MatMultEqual(A,sA,10,&flg);
257:   if (!flg) {
258:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale()\n",rank);
259:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
260:   }

262:   MatMultAddEqual(A,sA,10,&flg);
263:   if (!flg) {
264:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd()\n",rank);
265:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
266:   }

268:   /* Test MatMultTranspose(), MatMultTransposeAdd() */
269:   for (i=0; i<10; i++) {
270:     VecSetRandom(x,rctx);
271:     MatMultTranspose(A,x,s1);
272:     MatMultTranspose(sA,x,s2);
273:     VecNorm(s1,NORM_1,&r1);
274:     VecNorm(s2,NORM_1,&r2);
275:     r1  -= r2;
276:     if (r1<-tol || r1>tol) {
277:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,(double)r1);
278:       PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
279:     }
280:   }
281:   for (i=0; i<10; i++) {
282:     VecSetRandom(x,rctx);
283:     VecSetRandom(y,rctx);
284:     MatMultTransposeAdd(A,x,y,s1);
285:     MatMultTransposeAdd(sA,x,y,s2);
286:     VecNorm(s1,NORM_1,&r1);
287:     VecNorm(s2,NORM_1,&r2);
288:     r1  -= r2;
289:     if (r1<-tol || r1>tol) {
290:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,(double)r1);
291:       PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
292:     }
293:   }
294:   /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD);  */
295:   /* MatView(sA, PETSC_VIEWER_DRAW_WORLD);  */

297:   /* Test MatDuplicate() */
298:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
299:   MatEqual(sA,sB,&flg);
300:   if (!flg) {
301:     PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
302:   }
303:   MatMultEqual(sA,sB,5,&flg);
304:   if (!flg) {
305:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
306:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
307:   }
308:   MatMultAddEqual(sA,sB,5,&flg);
309:   if (!flg) {
310:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
311:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
312:   }
313:   MatDestroy(&sB);

315:   VecDestroy(&u);
316:   VecDestroy(&x);
317:   VecDestroy(&y);
318:   VecDestroy(&s1);
319:   VecDestroy(&s2);
320:   MatDestroy(&sA);
321:   MatDestroy(&A);
322:   PetscRandomDestroy(&rctx);

324:   PetscFinalize();
325:   return 0;
326: }