Actual source code: ex75.c

  2: /* Program usage:  mpiexec -n <procs> ex75 [-help] [all PETSc options] */

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

 6:  #include petscmat.h

 10: int main(int argc,char **args)
 11: {
 12:   Vec               x,y,u,s1,s2;
 13:   Mat               A,sA,sB;
 14:   PetscRandom       rctx;
 15:   PetscReal         r1,r2,rnorm,tol=1.e-10;
 16:   PetscScalar       one=1.0, neg_one=-1.0, value[3], four=4.0,alpha=0.1;
 17:   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;
 18:   PetscErrorCode    ierr;
 19:   PetscMPIInt       size,rank;
 20:   PetscTruth        flg;
 21:   const MatType     type;

 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:   MatCreate(PETSC_COMM_WORLD,&sA);
 34:   MatSetSizes(sA,PETSC_DECIDE,PETSC_DECIDE,n,n);
 35:   MatSetType(sA,MATSBAIJ);
 36:   MatSetFromOptions(sA);
 37:   MatGetType(sA,&type);
 38:   /* printf(" mattype: %s\n",type); */
 39:   MatMPISBAIJSetPreallocation(sA,bs,d_nz,PETSC_NULL,o_nz,PETSC_NULL);

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

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

102:   /* Test MatView() */
103:   /*
104:   MatView(sA, PETSC_VIEWER_STDOUT_WORLD); 
105:   MatView(sA, PETSC_VIEWER_DRAW_WORLD);
106:   */
107:   /* Assemble MPIBAIJ matrix A */
108:   MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,n,n,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&A);

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

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

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

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

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

193:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
194:   PetscRandomSetFromOptions(rctx);
195:   VecSetRandom(x,rctx);
196:   VecSet(u,one);

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

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

227:   /* Test MatDiagonalScale() */
228:   MatDiagonalScale(A,x,x);
229:   MatDiagonalScale(sA,x,x);
230:   MatMultEqual(A,sA,10,&flg);
231:   if (!flg) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");
232: 
233:   /* Test MatGetDiagonal(), MatScale() */
234:   MatGetDiagonal(A,s1);
235:   MatGetDiagonal(sA,s2);
236:   VecNorm(s1,NORM_1,&r1);
237:   VecNorm(s2,NORM_1,&r2);
238:   r1 -= r2;
239:   if (r1<-tol || r1>tol) {
240:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDiagonalScale() or MatGetDiagonal(), r1=%G \n",rank,r1);
241:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
242:   }
243: 
244:   MatScale(A,alpha);
245:   MatScale(sA,alpha);

247:   /* Test MatGetRowMaxAbs() */
248:   MatGetRowMaxAbs(A,s1,PETSC_NULL);
249:   MatGetRowMaxAbs(sA,s2,PETSC_NULL);

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

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

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

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

300:   /* Test MatDuplicate() */
301:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
302:   MatEqual(sA,sB,&flg);
303:   if (!flg){
304:     PetscPrintf(PETSC_COMM_WORLD," Error in MatDuplicate(), sA != sB \n");
305:   }
306:   MatMultEqual(sA,sB,5,&flg);
307:   if (!flg){
308:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMult()\n",rank);
309:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
310:   }
311:   MatMultAddEqual(sA,sB,5,&flg);
312:   if (!flg){
313:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d], Error: MatDuplicate() or MatMultAdd(()\n",rank);
314:     PetscSynchronizedFlush(PETSC_COMM_WORLD);
315:   }
316:   MatDestroy(sB);
317: 
318:   VecDestroy(u);
319:   VecDestroy(x);
320:   VecDestroy(y);
321:   VecDestroy(s1);
322:   VecDestroy(s2);
323:   MatDestroy(sA);
324:   MatDestroy(A);
325:   PetscRandomDestroy(rctx);
326: 
327:   PetscFinalize();
328:   return 0;
329: }