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;
 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;

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

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

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

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

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

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

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

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

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

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

189:   PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
190:   VecSetRandom(rctx,x);
191:   VecSet(&one,u);

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

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

211:   /* Test MatGetRow(): can only obtain rows associated with the given processor */
212:   for (i=rstart; i<rstart+1; i++) {
213:     MatGetRow(sA,i,&ncols,&cols,&vr);
214:     /*
215:     PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"[%d] get row %d: ",rank,i);
216:     for (j=0; j<ncols; j++) {
217:       PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"%d %g  ",cols[j],vr[j]);
218:     }
219:     PetscSynchronizedFPrintf(PETSC_COMM_WORLD,stdout,"\n");
220:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
221:     */
222:     MatRestoreRow(sA,i,&ncols,&cols,&vr);
223:   }

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

241:   /* Test MatGetRowMax() */
242:   MatGetRowMax(A,s1);
243:   MatGetRowMax(sA,s2);

245:   VecNorm(s1,NORM_1,&r1);
246:   VecNorm(s2,NORM_1,&r2);
247:   /* MatView(A,PETSC_VIEWER_STDOUT_WORLD);
248:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], r1: %g, r2: %g\n",rank,r1,r2);
249:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
250:       */
251:   r1 -= r2;
252:   if (r1<-tol || r1>tol) {
253:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetRowMax() \n");
254:   }

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

270:   for (i=0; i<10; i++) {
271:     VecSetRandom(rctx,x);
272:     VecSetRandom(rctx,y);
273:     MatMultAdd(A,x,y,s1);
274:     MatMultAdd(sA,x,y,s2);
275:     VecNorm(s1,NORM_1,&r1);
276:     VecNorm(s2,NORM_1,&r2);
277:     r1 -= r2;
278:     if (r1<-tol || r1>tol) {
279:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatMultAdd(), err=%g \n",rank,r1);
280:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
281:     }
282:   }

284:   /* MatView(sA, PETSC_VIEWER_STDOUT_WORLD);  */
285:   /* MatView(sA, PETSC_VIEWER_DRAW_WORLD);  */
286: 
287:   VecDestroy(u);
288:   VecDestroy(x);
289:   VecDestroy(y);
290:   VecDestroy(s1);
291:   VecDestroy(s2);
292:   MatDestroy(sA);
293:   MatDestroy(A);
294:   PetscRandomDestroy(rctx);
295: 
296:   PetscFinalize();
297:   return 0;
298: }