Actual source code: ex75.c

  1: /*$Id: ex75.c,v 1.28 2001/09/11 16:32:50 bsmith Exp $*/

  3: /* Program usage:  mpirun -np <procs> ex75 [-help] [all PETSc options] */

  5: static char help[] = "Tests the vatious routines in MatMPISBAIJ format.\n";

 7:  #include petscmat.h

 11: int main(int argc,char **args)
 12: {
 13:   int         bs=1, mbs=16, d_nz=3, o_nz=3, prob=2;
 14:   Vec         x,y,u,s1,s2;
 15:   Mat         A,sA,sB;
 16:   PetscRandom rctx;
 17:   double      r1,r2,tol=1.e-10;
 18:   int         i,j,i2,j2,I,J,ierr;
 19:   PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1,*vr;
 20:   int         n,rank,size,col[3],n1,block,row;
 21:   int         ncols,*cols,rstart,rend;
 22:   PetscTruth  flg;

 24:   PetscInitialize(&argc,&args,(char *)0,help);
 25:   PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
 26:   PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 27: 
 28:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 29:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 30: 
 31:   n = mbs*bs;
 32: 
 33:   /* Assemble MPISBAIJ matrix sA */
 34:   MatCreateMPISBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&sA);

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

 48:       i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
 49:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
 50:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 51:     }
 52:     else if (prob ==2){ /* matrix for the five point stencil */
 53:       n1 =  (int) sqrt((double)n);
 54:       if (n1*n1 != n){
 55:         SETERRQ(PETSC_ERR_ARG_SIZ,"n must be a perfect square of n1");
 56:       }
 57: 
 58:       for (i=0; i<n1; i++) {
 59:         for (j=0; j<n1; j++) {
 60:           I = j + n1*i;
 61:           if (i>0)   {J = I - n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
 62:           if (i<n1-1) {J = I + n1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
 63:           if (j>0)   {J = I - 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
 64:           if (j<n1-1) {J = I + 1; MatSetValues(sA,1,&I,1,&J,&neg_one,INSERT_VALUES);}
 65:           MatSetValues(sA,1,&I,1,&I,&four,INSERT_VALUES);
 66:         }
 67:       }
 68:     }
 69:   } /* end of if (bs == 1) */
 70:   else {  /* bs > 1 */
 71:   for (block=0; block<n/bs; block++){
 72:     /* diagonal blocks */
 73:     value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 74:     for (i=1+block*bs; i<bs-1+block*bs; i++) {
 75:       col[0] = i-1; col[1] = i; col[2] = i+1;
 76:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 77:     }
 78:     i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
 79:     value[0]=-1.0; value[1]=4.0;
 80:     MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

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

 98:   /* Test MatView() */
 99:   /*
100:   MatView(sA, PETSC_VIEWER_STDOUT_WORLD); 
101:   MatView(sA, PETSC_VIEWER_DRAW_WORLD);
102:   */
103:   /* Assemble MPIBAIJ matrix A */
104:   MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);

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

118:       i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
119:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
120:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
121:     }
122:     else if (prob ==2){ /* matrix for the five point stencil */
123:       n1 = (int) sqrt((double)n);
124:       for (i=0; i<n1; i++) {
125:         for (j=0; j<n1; j++) {
126:           I = j + n1*i;
127:           if (i>0)   {J = I - n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
128:           if (i<n1-1) {J = I + n1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
129:           if (j>0)   {J = I - 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
130:           if (j<n1-1) {J = I + 1; MatSetValues(A,1,&I,1,&J,&neg_one,INSERT_VALUES);}
131:           MatSetValues(A,1,&I,1,&I,&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);
170:   }
171: 
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);
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,RANDOM_DEFAULT,&rctx);
191:   VecSetRandom(rctx,x);
192:   VecSet(&one,u);

194:   /* Test MatNorm() */
195:   MatNorm(A,NORM_FROBENIUS,&r1);
196:   MatNorm(sA,NORM_FROBENIUS,&r2);
197:   r1 -= r2;
198:   if (r1<-tol || r1>tol){
199:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatNorm(), A_fnorm - sA_fnorm = %16.14e\n",rank,r1);
200:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
201:   }

203:   /* Test MatGetOwnershipRange() */
204:   MatGetOwnershipRange(sA,&rstart,&rend);
205:   MatGetOwnershipRange(A,&i2,&j2);
206:   i2 -= rstart; j2 -= rend;
207:   if (i2 || j2) {
208:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MaGetOwnershipRange()\n",rank);
209:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
210:   }

212:   /* Test MatGetRow(): can only obtain rows associated with the given processor */
213:   for (i=rstart; i<rstart+1; i++) {
214:     MatGetRow(sA,i,&ncols,&cols,&vr);
215:     MatRestoreRow(sA,i,&ncols,&cols,&vr);
216:   }

218:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
219:   MatDiagonalScale(A,x,x);
220:   MatDiagonalScale(sA,x,x);
221:   MatGetDiagonal(A,s1);
222:   MatGetDiagonal(sA,s2);
223:   VecNorm(s1,NORM_1,&r1);
224:   VecNorm(s2,NORM_1,&r2);
225:   r1 -= r2;
226:   if (r1<-tol || r1>tol) {
227:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%g \n",rank,r1);
228:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
229:   }
230: 
231:   MatScale(&alpha,A);
232:   MatScale(&alpha,sA);

234:   /* Test MatGetRowMax() */
235:   MatGetRowMax(A,s1);
236:   MatGetRowMax(sA,s2);

238:   VecNorm(s1,NORM_1,&r1);
239:   VecNorm(s2,NORM_1,&r2);
240:   r1 -= r2;
241:   if (r1<-tol || r1>tol) {
242:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMax() \n");
243:   }

245:   /* Test MatMult(), MatMultAdd() */
246:   for (i=0; i<10; i++) {
247:     VecSetRandom(rctx,x);
248:     MatMult(A,x,s1);
249:     MatMult(sA,x,s2);
250:     VecNorm(s1,NORM_1,&r1);
251:     VecNorm(s2,NORM_1,&r2);
252:     r1 -= r2;
253:     if (r1<-tol || r1>tol) {
254:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMult() or MatScale(), err=%g\n",rank,r1);
255:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
256:     }
257:   }

259:   for (i=0; i<10; i++) {
260:     VecSetRandom(rctx,x);
261:     VecSetRandom(rctx,y);
262:     MatMultAdd(A,x,y,s1);
263:     MatMultAdd(sA,x,y,s2);
264:     VecNorm(s1,NORM_1,&r1);
265:     VecNorm(s2,NORM_1,&r2);
266:     r1 -= r2;
267:     if (r1<-tol || r1>tol) {
268:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,r1);
269:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
270:     }
271:   }

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

302:   /* Test MatDuplicate() */
303:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
304:   MatEqual(sA,sB,&flg);
305:   if (!flg){
306:     PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
307:   }
308:   for (i=0; i<5; i++) {
309:     VecSetRandom(rctx,x);
310:     MatMult(A,x,s1);
311:     MatMult(sB,x,s2);
312:     VecNorm(s1,NORM_1,&r1);
313:     VecNorm(s2,NORM_1,&r2);
314:     r1 -= r2;
315:     if (r1<-tol || r1>tol) {
316:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult(), err=%g\n",rank,r1);
317:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
318:     }
319:   }

321:   for (i=0; i<5; i++) {
322:     VecSetRandom(rctx,x);
323:     VecSetRandom(rctx,y);
324:     MatMultAdd(A,x,y,s1);
325:     MatMultAdd(sB,x,y,s2);
326:     VecNorm(s1,NORM_1,&r1);
327:     VecNorm(s2,NORM_1,&r2);
328:     r1 -= r2;
329:     if (r1<-tol || r1>tol) {
330:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(), err=%g \n",rank,r1);
331:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
332:     }
333:   }
334:   MatDestroy(sB);
335: 
336:   VecDestroy(u);
337:   VecDestroy(x);
338:   VecDestroy(y);
339:   VecDestroy(s1);
340:   VecDestroy(s2);
341:   MatDestroy(sA);
342:   MatDestroy(A);
343:   PetscRandomDestroy(rctx);
344: 
345:   PetscFinalize();
346:   return 0;
347: }