Actual source code: ex74.c

  2: static char help[] = "Tests the various sequential routines in MatSBAIJ format.\n";

  4: #include <petscmat.h>

  8: int main(int argc,char **args)
  9: {
 10:   PetscMPIInt        size;
 11:   PetscErrorCode     ierr;
 12:   Vec                x,y,b,s1,s2;
 13:   Mat                A;                /* linear system matrix */
 14:   Mat                sA,sB,sC;         /* symmetric part of the matrices */
 15:   PetscInt           n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc;
 16:   PetscReal          norm1,norm2,rnorm,tol=1.e-10;
 17:   PetscScalar        neg_one = -1.0,four=4.0,value[3];
 18:   IS                 perm, iscol;
 19:   PetscRandom        rdm;
 20:   PetscBool          doIcc=PETSC_TRUE,equal;
 21:   MatInfo            minfo1,minfo2;
 22:   MatFactorInfo      factinfo;
 23:   const MatType      type;

 25:   PetscInitialize(&argc,&args,(char *)0,help);
 26:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 27:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 28:   PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 29:   PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);

 31:   n = mbs*bs;
 32:   MatCreate(PETSC_COMM_SELF,&A);
 33:   MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 34:   MatSetType(A,MATSEQBAIJ);
 35:   MatSetFromOptions(A);
 36:   MatSeqBAIJSetPreallocation(A,bs,nz,PETSC_NULL);

 38:   MatCreate(PETSC_COMM_SELF,&sA);
 39:   MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 40:   MatSetType(sA,MATSEQSBAIJ);
 41:   MatSetFromOptions(sA);
 42:   MatGetType(sA,&type);
 43:   PetscTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
 44:   MatSeqSBAIJSetPreallocation(sA,bs,nz,PETSC_NULL);
 45:   MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);

 47:   /* Test MatGetOwnershipRange() */
 48:   MatGetOwnershipRange(A,&Ii,&J);
 49:   MatGetOwnershipRange(sA,&i,&j);
 50:   if (i-Ii || j-J){
 51:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
 52:   }

 54:   /* Assemble matrix */
 55:   if (bs == 1){
 56:     PetscOptionsGetInt(PETSC_NULL,"-test_problem",&prob,PETSC_NULL);
 57:     if (prob == 1){ /* tridiagonal matrix */
 58:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 59:       for (i=1; i<n-1; i++) {
 60:         col[0] = i-1; col[1] = i; col[2] = i+1;
 61:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 62:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 63:       }
 64:       i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
 65:       value[0]= 0.1; value[1]=-1; value[2]=2;
 66:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 67:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

 69:       i = 0;
 70:       col[0] = n-1;  col[1] = 1; col[2]=0;
 71:       value[0] = 0.1; value[1] = -1.0; value[2]=2;
 72:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 73:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

 75:     } else if (prob ==2){ /* matrix for the five point stencil */
 76:       n1 = (int) (PetscSqrtReal((PetscReal)n) + 0.001);
 77:       if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
 78:       for (i=0; i<n1; i++) {
 79:         for (j=0; j<n1; j++) {
 80:           Ii = j + n1*i;
 81:           if (i>0)   {
 82:             J = Ii - n1;
 83:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 84:             MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 85:           }
 86:           if (i<n1-1) {
 87:             J = Ii + n1;
 88:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 89:             MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 90:           }
 91:           if (j>0)   {
 92:             J = Ii - 1;
 93:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 94:             MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 95:           }
 96:           if (j<n1-1) {
 97:             J = Ii + 1;
 98:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 99:             MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
100:           }
101:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
102:           MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
103:         }
104:       }
105:     }

107:   } else { /* bs > 1 */
108:     for (block=0; block<n/bs; block++){
109:       /* diagonal blocks */
110:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
111:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
112:         col[0] = i-1; col[1] = i; col[2] = i+1;
113:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
114:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
115:       }
116:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
117:       value[0]=-1.0; value[1]=4.0;
118:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
119:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

121:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
122:       value[0]=4.0; value[1] = -1.0;
123:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
124:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
125:     }
126:     /* off-diagonal blocks */
127:     value[0]=-1.0;
128:     for (i=0; i<(n/bs-1)*bs; i++){
129:       col[0]=i+bs;
130:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
131:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
132:       col[0]=i; row=i+bs;
133:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
134:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
135:     }
136:   }
137:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
138:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
139: 
140:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
141:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);

143:   /* Test MatGetInfo() of A and sA */
144:   MatGetInfo(A,MAT_LOCAL,&minfo1);
145:   MatGetInfo(sA,MAT_LOCAL,&minfo2);
146:   /*
147:   printf("A matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated); 
148:   printf("sA matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated); 
149:   */
150:   i = (int) (minfo1.nz_used - minfo2.nz_used);
151:   j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
152:   k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
153:   k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
154:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
155:     PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n");
156:   }

158:   /* Test MatDuplicate() */
159:   MatNorm(A,NORM_FROBENIUS,&norm1);
160:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
161:   MatEqual(sA,sB,&equal);
162:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");

164:   /* Test MatNorm() */
165:   MatNorm(A,NORM_FROBENIUS,&norm1);
166:   MatNorm(sB,NORM_FROBENIUS,&norm2);
167:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
168:   if (rnorm > tol){
169:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
170:   }
171:   MatNorm(A,NORM_INFINITY,&norm1);
172:   MatNorm(sB,NORM_INFINITY,&norm2);
173:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
174:   if (rnorm > tol){
175:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
176:   }
177:   MatNorm(A,NORM_1,&norm1);
178:   MatNorm(sB,NORM_1,&norm2);
179:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
180:   if (rnorm > tol){
181:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
182:   }

184:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
185:   MatGetInfo(A,MAT_LOCAL,&minfo1);
186:   MatGetInfo(sB,MAT_LOCAL,&minfo2);
187:   /*
188:   printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated); 
189:   printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated); 
190:   */
191:   i = (int) (minfo1.nz_used - minfo2.nz_used);
192:   j = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
193:   k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
194:   k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
195:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
196:     PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n");
197:   }

199:   MatGetSize(A,&Ii,&J);
200:   MatGetSize(sB,&i,&j);
201:   if (i-Ii || j-J) {
202:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
203:   }
204: 
205:   MatGetBlockSize(A, &Ii);
206:   MatGetBlockSize(sB, &i);
207:   if (i-Ii){
208:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
209:   }

211:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
212:   PetscRandomSetFromOptions(rdm);
213:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
214:   VecDuplicate(x,&s1);
215:   VecDuplicate(x,&s2);
216:   VecDuplicate(x,&y);
217:   VecDuplicate(x,&b);
218:   VecSetRandom(x,rdm);

220:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
221: #if !defined(PETSC_USE_COMPLEX)
222:   /* Scaling matrix with complex numbers results non-spd matrix, 
223:      causing crash of MatForwardSolve() and MatBackwardSolve() */
224:   MatDiagonalScale(A,x,x);
225:   MatDiagonalScale(sB,x,x);
226:   MatMultEqual(A,sB,10,&equal);
227:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");

229:   MatGetDiagonal(A,s1);
230:   MatGetDiagonal(sB,s2);
231:   VecAXPY(s2,neg_one,s1);
232:   VecNorm(s2,NORM_1,&norm1);
233:   if ( norm1>tol) {
234:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%G\n",norm1);
235:   }

237:   {
238:     PetscScalar alpha=0.1;
239:     MatScale(A,alpha);
240:     MatScale(sB,alpha);
241:   }
242: #endif

244:   /* Test MatGetRowMaxAbs() */
245:   MatGetRowMaxAbs(A,s1,PETSC_NULL);
246:   MatGetRowMaxAbs(sB,s2,PETSC_NULL);
247:   VecNorm(s1,NORM_1,&norm1);
248:   VecNorm(s2,NORM_1,&norm2);
249:   norm1 -= norm2;
250:   if (norm1<-tol || norm1>tol) {
251:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n");
252:   }

254:   /* Test MatMult() */
255:   for (i=0; i<40; i++) {
256:     VecSetRandom(x,rdm);
257:     MatMult(A,x,s1);
258:     MatMult(sB,x,s2);
259:     VecNorm(s1,NORM_1,&norm1);
260:     VecNorm(s2,NORM_1,&norm2);
261:     norm1 -= norm2;
262:     if (norm1<-tol || norm1>tol) {
263:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %G\n",norm1);
264:     }
265:   }

267:   /* MatMultAdd() */
268:   for (i=0; i<40; i++) {
269:     VecSetRandom(x,rdm);
270:     VecSetRandom(y,rdm);
271:     MatMultAdd(A,x,y,s1);
272:     MatMultAdd(sB,x,y,s2);
273:     VecNorm(s1,NORM_1,&norm1);
274:     VecNorm(s2,NORM_1,&norm2);
275:     norm1 -= norm2;
276:     if (norm1<-tol || norm1>tol) {
277:       PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(),  norm1-norm2: %G\n",norm1);
278:     }
279:   }

281:   /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
282:   MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol);
283:   ISDestroy(&iscol);
284:   norm1 = tol;
285:   inc   = bs;

287:   /* initialize factinfo */
288:   PetscMemzero(&factinfo,sizeof(MatFactorInfo));

290:   for (lf=-1; lf<10; lf += inc){
291:     if (lf==-1) {  /* Cholesky factor of sB (duplicate sA) */
292:       factinfo.fill = 5.0;
293:       MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
294:       MatCholeskyFactorSymbolic(sC,sB,perm,&factinfo);
295:     } else if (!doIcc){
296:       break;
297:     } else {       /* incomplete Cholesky factor */
298:       factinfo.fill   = 5.0;
299:       factinfo.levels = lf;
300:       MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
301:       MatICCFactorSymbolic(sC,sB,perm,&factinfo);
302:     }
303:     MatCholeskyFactorNumeric(sC,sB,&factinfo);
304:     /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */

306:     /* test MatGetDiagonal on numeric factor */
307:     /*
308:     if (lf == -1) {
309:       MatGetDiagonal(sC,s1);  
310:       printf(" in ex74.c, diag: \n");
311:       VecView(s1,PETSC_VIEWER_STDOUT_SELF);
312:     }
313:     */

315:     MatMult(sB,x,b);

317:     /* test MatForwardSolve() and MatBackwardSolve() */
318:     if (lf == -1){
319:       MatForwardSolve(sC,b,s1);
320:       MatBackwardSolve(sC,s1,s2);
321:       VecAXPY(s2,neg_one,x);
322:       VecNorm(s2,NORM_2,&norm2);
323:       if (10*norm1 < norm2){
324:         PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G, bs=%d\n",norm2,bs);
325:       }
326:     }

328:     /* test MatSolve() */
329:     MatSolve(sC,b,y);
330:     MatDestroy(&sC);
331:     /* Check the error */
332:     VecAXPY(y,neg_one,x);
333:     VecNorm(y,NORM_2,&norm2);
334:     /* printf("lf: %d, error: %G\n", lf,norm2); */
335:     if (10*norm1 < norm2 && lf-inc != -1){
336:       PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%G, %G\n",lf-inc,lf,norm1,norm2);
337:     }
338:     norm1 = norm2;
339:     if (norm2 < tol && lf != -1) break;
340:   }

342:   ISDestroy(&perm);

344:   MatDestroy(&A);
345:   MatDestroy(&sB);
346:   MatDestroy(&sA);
347:   VecDestroy(&x);
348:   VecDestroy(&y);
349:   VecDestroy(&s1);
350:   VecDestroy(&s2);
351:   VecDestroy(&b);
352:   PetscRandomDestroy(&rdm);

354:   PetscFinalize();
355:   return 0;
356: }