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,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:   PetscTruth         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_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:   ierr=MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,n,n,nz,PETSC_NULL, &A);
 33:   MatCreate(PETSC_COMM_SELF,&sA);
 34:   MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 35:   MatSetType(sA,MATSEQSBAIJ);
 36:   MatSetFromOptions(sA);
 37:   MatGetType(sA,&type);
 38:   PetscTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
 39:   MatSeqSBAIJSetPreallocation(sA,bs,nz,PETSC_NULL);
 40:   MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);

 42:   /* Test MatGetOwnershipRange() */
 43:   MatGetOwnershipRange(A,&Ii,&J);
 44:   MatGetOwnershipRange(sA,&i,&j);
 45:   if (i-Ii || j-J){
 46:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
 47:   }

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

 64:       i = 0;
 65:       col[0] = n-1;  col[1] = 1; col[2]=0;
 66:       value[0] = 0.1; value[1] = -1.0; value[2]=2;
 67:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 68:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
 69:     }
 70:     else if (prob ==2){ /* matrix for the five point stencil */
 71:       n1 = (int) (sqrt((PetscReal)n) + 0.001);
 72:       if (n1*n1 - n) SETERRQ(PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
 73:       for (i=0; i<n1; i++) {
 74:         for (j=0; j<n1; j++) {
 75:           Ii = j + n1*i;
 76:           if (i>0)   {
 77:             J = Ii - n1;
 78:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 79:             MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 80:           }
 81:           if (i<n1-1) {
 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 (j>0)   {
 87:             J = Ii - 1;
 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<n1-1) {
 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:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 97:           MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 98:         }
 99:       }
100:     }
101:   }
102:   else { /* bs > 1 */
103:     for (block=0; block<n/bs; block++){
104:       /* diagonal blocks */
105:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
106:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
107:         col[0] = i-1; col[1] = i; col[2] = i+1;
108:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
109:         MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);
110:       }
111:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
112:       value[0]=-1.0; value[1]=4.0;
113:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
114:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

116:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
117:       value[0]=4.0; value[1] = -1.0;
118:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
119:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
120:     }
121:     /* off-diagonal blocks */
122:     value[0]=-1.0;
123:     for (i=0; i<(n/bs-1)*bs; i++){
124:       col[0]=i+bs;
125:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
126:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);
127:       col[0]=i; row=i+bs;
128:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
129:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
130:     }
131:   }
132:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
133:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
134:   /* PetscPrintf(PETSC_COMM_SELF,"\n The Matrix: \n");
135:   MatView(A, PETSC_VIEWER_DRAW_WORLD);
136:   MatView(A, PETSC_VIEWER_STDOUT_WORLD); */

138:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
139:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);
140:   /* PetscPrintf(PETSC_COMM_SELF,"\n Symmetric Part of Matrix: \n");
141:   MatView(sA, PETSC_VIEWER_DRAW_WORLD); 
142:   MatView(sA, PETSC_VIEWER_STDOUT_WORLD); 
143:   */

145:   /* Test MatDuplicate() */
146:   MatNorm(A,NORM_FROBENIUS,&norm1);
147:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
148:   MatEqual(sA,sB,&equal);
149:   if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()");

151:   /* Test MatNorm() */
152:   MatNorm(A,NORM_FROBENIUS,&norm1);
153:   MatNorm(sB,NORM_FROBENIUS,&norm2);
154:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
155:   if (rnorm > tol){
156:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
157:   }
158:   MatNorm(A,NORM_INFINITY,&norm1);
159:   MatNorm(sB,NORM_INFINITY,&norm2);
160:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
161:   if (rnorm > tol){
162:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
163:   }
164:   MatNorm(A,NORM_1,&norm1);
165:   MatNorm(sB,NORM_1,&norm2);
166:   rnorm = PetscAbsScalar(norm1-norm2)/norm2;
167:   if (rnorm > tol){
168:     PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",norm1,norm2);
169:   }

171:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
172:   MatGetInfo(A,MAT_LOCAL,&minfo1);
173:   MatGetInfo(sB,MAT_LOCAL,&minfo2);
174:   /*
175:   printf("matrix nonzeros (BAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo1.nz_used,(int)minfo1.nz_allocated); 
176:   printf("matrix nonzeros(SBAIJ format) = %d, allocated nonzeros= %d\n", (int)minfo2.nz_used,(int)minfo2.nz_allocated); 
177:   */
178:   i = (int) (minfo1.nz_used - minfo2.nz_used);
179:   j = (int) (minfo2.nz_allocated - minfo2.nz_used);
180:   if (i<0 || j<0) {
181:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetInfo()\n");
182:   }

184:   MatGetSize(A,&Ii,&J);
185:   MatGetSize(sB,&i,&j);
186:   if (i-Ii || j-J) {
187:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
188:   }
189: 
190:   MatGetBlockSize(A, &Ii);
191:   MatGetBlockSize(sB, &i);
192:   if (i-Ii){
193:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
194:   }

196:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
197:   PetscRandomSetFromOptions(rdm);
198:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
199:   VecDuplicate(x,&s1);
200:   VecDuplicate(x,&s2);
201:   VecDuplicate(x,&y);
202:   VecDuplicate(x,&b);
203:   VecSetRandom(x,rdm);

205:   /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
206: #if !defined(PETSC_USE_COMPLEX)
207:   /* Scaling matrix with complex numbers results non-spd matrix, 
208:      causing crash of MatForwardSolve() and MatBackwardSolve() */
209:   MatDiagonalScale(A,x,x);
210:   MatDiagonalScale(sB,x,x);
211:   MatMultEqual(A,sB,10,&equal);
212:   if (!equal) SETERRQ(PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale");

214:   MatGetDiagonal(A,s1);
215:   MatGetDiagonal(sB,s2);
216:   VecAXPY(s2,neg_one,s1);
217:   VecNorm(s2,NORM_1,&norm1);
218:   if ( norm1>tol) {
219:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%G\n",norm1);
220:   }

222:   {
223:     PetscScalar alpha=0.1;
224:     MatScale(A,alpha);
225:     MatScale(sB,alpha);
226:   }
227: #endif

229:   /* Test MatGetRowMaxAbs() */
230:   MatGetRowMaxAbs(A,s1,PETSC_NULL);
231:   MatGetRowMaxAbs(sB,s2,PETSC_NULL);
232:   VecNorm(s1,NORM_1,&norm1);
233:   VecNorm(s2,NORM_1,&norm2);
234:   norm1 -= norm2;
235:   if (norm1<-tol || norm1>tol) {
236:     PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n");
237:   }

239:   /* Test MatMult() */
240:   for (i=0; i<40; i++) {
241:     VecSetRandom(x,rdm);
242:     MatMult(A,x,s1);
243:     MatMult(sB,x,s2);
244:     VecNorm(s1,NORM_1,&norm1);
245:     VecNorm(s2,NORM_1,&norm2);
246:     norm1 -= norm2;
247:     if (norm1<-tol || norm1>tol) {
248:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %G\n",norm1);
249:     }
250:   }

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

266:   /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
267:   MatGetOrdering(A,MATORDERING_NATURAL,&perm,&iscol);
268:   ISDestroy(iscol);
269:   norm1 = tol;
270:   inc   = bs;

272:   /* initialize factinfo */
273:   PetscMemzero(&factinfo,sizeof(MatFactorInfo));

275:   for (lf=-1; lf<10; lf += inc){
276:     if (lf==-1) {  /* Cholesky factor of sB (duplicate sA) */
277:       factinfo.fill = 5.0;
278:       MatGetFactor(sB,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&sC);
279:       MatCholeskyFactorSymbolic(sC,sB,perm,&factinfo);
280:     } else if (!doIcc){
281:       break;
282:     } else {       /* incomplete Cholesky factor */
283:       factinfo.fill   = 5.0;
284:       factinfo.levels = lf;
285:       MatGetFactor(sB,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&sC);
286:       MatICCFactorSymbolic(sC,sB,perm,&factinfo);
287:     }
288:     MatCholeskyFactorNumeric(sC,sB,&factinfo);
289:     /* MatView(sC, PETSC_VIEWER_DRAW_WORLD); */

291:     /* test MatGetDiagonal on numeric factor */
292:     /*
293:     if (lf == -1) {
294:       MatGetDiagonal(sC,s1);  
295:       printf(" in ex74.c, diag: \n");
296:       VecView(s1,PETSC_VIEWER_STDOUT_SELF);
297:     }
298:     */

300:     MatMult(sB,x,b);

302:     /* test MatForwardSolve() and MatBackwardSolve() */
303:     if (lf == -1){
304:       MatForwardSolve(sC,b,s1);
305:       MatBackwardSolve(sC,s1,s2);
306:       VecAXPY(s2,neg_one,x);
307:       VecNorm(s2,NORM_2,&norm2);
308:       if (10*norm1 < norm2){
309:         PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G, bs=%d\n",norm2,bs);
310:       }
311:     }

313:     /* test MatSolve() */
314:     MatSolve(sC,b,y);
315:     MatDestroy(sC);
316:     /* Check the error */
317:     VecAXPY(y,neg_one,x);
318:     VecNorm(y,NORM_2,&norm2);
319:     /* printf("lf: %d, error: %G\n", lf,norm2); */
320:     if (10*norm1 < norm2 && lf-inc != -1){
321:       PetscPrintf(PETSC_COMM_SELF,"lf=%D, %D, Norm of error=%G, %G\n",lf-inc,lf,norm1,norm2);
322:     }
323:     norm1 = norm2;
324:     if (norm2 < tol && lf != -1) break;
325:   }

327:   ISDestroy(perm);

329:   MatDestroy(A);
330:   MatDestroy(sB);
331:   MatDestroy(sA);
332:   VecDestroy(x);
333:   VecDestroy(y);
334:   VecDestroy(s1);
335:   VecDestroy(s2);
336:   VecDestroy(b);
337:   PetscRandomDestroy(rdm);

339:   PetscFinalize();
340:   return 0;
341: }