Actual source code: ex74.c


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

  4: #include <petscmat.h>

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

 22:   PetscInitialize(&argc,&args,(char*)0,help);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 25:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 26:   PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);

 28:   n    = mbs*bs;
 29:   MatCreate(PETSC_COMM_SELF,&A);
 30:   MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 31:   MatSetType(A,MATSEQBAIJ);
 32:   MatSetFromOptions(A);
 33:   MatSeqBAIJSetPreallocation(A,bs,nz,NULL);

 35:   MatCreate(PETSC_COMM_SELF,&sA);
 36:   MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
 37:   MatSetType(sA,MATSEQSBAIJ);
 38:   MatSetFromOptions(sA);
 39:   MatGetType(sA,&type);
 40:   PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);
 41:   MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL);
 42:   MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);

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

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

 63:       value[0]= 0.1; value[1]=-1; value[2]=2;

 65:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 66:       MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);

 68:       i        = 0;
 69:       col[0]   = n-1;   col[1] = 1;      col[2] = 0;
 70:       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 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
 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;

118:       value[0]=-1.0; value[1]=4.0;

120:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
121:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);

123:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;

125:       value[0]=4.0; value[1] = -1.0;

127:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
128:       MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES);
129:     }
130:     /* off-diagonal blocks */
131:     value[0]=-1.0;
132:     for (i=0; i<(n/bs-1)*bs; i++) {
133:       col[0]=i+bs;

135:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
136:       MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES);

138:       col[0]=i; row=i+bs;

140:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
141:       MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES);
142:     }
143:   }
144:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
145:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

147:   MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY);
148:   MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY);

150:   /* Test MatGetInfo() of A and sA */
151:   MatGetInfo(A,MAT_LOCAL,&minfo1);
152:   MatGetInfo(sA,MAT_LOCAL,&minfo2);
153:   i  = (int) (minfo1.nz_used - minfo2.nz_used);
154:   j  = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
155:   k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
156:   k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
157:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
158:     PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n");
159:   }

161:   /* Test MatDuplicate() */
162:   MatNorm(A,NORM_FROBENIUS,&norm1);
163:   MatDuplicate(sA,MAT_COPY_VALUES,&sB);
164:   MatEqual(sA,sB,&equal);

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

187:   /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
188:   MatGetInfo(A,MAT_LOCAL,&minfo1);
189:   MatGetInfo(sB,MAT_LOCAL,&minfo2);
190:   i  = (int) (minfo1.nz_used - minfo2.nz_used);
191:   j  = (int) (minfo1.nz_allocated - minfo2.nz_allocated);
192:   k1 = (int) (minfo1.nz_allocated - minfo1.nz_used);
193:   k2 = (int) (minfo2.nz_allocated - minfo2.nz_used);
194:   if (i < 0 || j < 0 || k1 < 0 || k2 < 0) {
195:     PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n");
196:   }

198:   MatGetSize(A,&Ii,&J);
199:   MatGetSize(sB,&i,&j);
200:   if (i-Ii || j-J) {
201:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n");
202:   }

204:   MatGetBlockSize(A, &Ii);
205:   MatGetBlockSize(sB, &i);
206:   if (i-Ii) {
207:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n");
208:   }

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

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

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

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

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

253:   /* Test MatMult() */
254:   for (i=0; i<40; i++) {
255:     VecSetRandom(x,rdm);
256:     MatMult(A,x,s1);
257:     MatMult(sB,x,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: MatMult(), norm1-norm2: %g\n",(double)norm1);
263:     }
264:   }

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

280:   /* Test MatMatMult() for sbaij and dense matrices */
281:   MatCreateSeqDense(PETSC_COMM_SELF,n,5*n,NULL,&B);
282:   MatSetRandom(B,rdm);
283:   MatMatMult(sA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
284:   MatMatMultEqual(sA,B,C,5*n,&equal);
286:   MatDestroy(&C);
287:   MatDestroy(&B);

289:   /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */
290:   MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol);
291:   ISDestroy(&iscol);
292:   norm1 = tol;
293:   inc   = bs;

295:   /* initialize factinfo */
296:   PetscMemzero(&factinfo,sizeof(MatFactorInfo));

298:   for (lf=-1; lf<10; lf += inc) {
299:     if (lf==-1) {  /* Cholesky factor of sB (duplicate sA) */
300:       factinfo.fill = 5.0;

302:       MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sFactor);
303:       MatCholeskyFactorSymbolic(sFactor,sB,perm,&factinfo);
304:     } else if (!doIcc) break;
305:     else {       /* incomplete Cholesky factor */
306:       factinfo.fill   = 5.0;
307:       factinfo.levels = lf;

309:       MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sFactor);
310:       MatICCFactorSymbolic(sFactor,sB,perm,&factinfo);
311:     }
312:     MatCholeskyFactorNumeric(sFactor,sB,&factinfo);
313:     /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */

315:     /* test MatGetDiagonal on numeric factor */
316:     /*
317:     if (lf == -1) {
318:       MatGetDiagonal(sFactor,s1);
319:       printf(" in ex74.c, diag: \n");
320:       VecView(s1,PETSC_VIEWER_STDOUT_SELF);
321:     }
322:     */

324:     MatMult(sB,x,b);

326:     /* test MatForwardSolve() and MatBackwardSolve() */
327:     if (lf == -1) {
328:       MatForwardSolve(sFactor,b,s1);
329:       MatBackwardSolve(sFactor,s1,s2);
330:       VecAXPY(s2,neg_one,x);
331:       VecNorm(s2,NORM_2,&norm2);
332:       if (10*norm1 < norm2) {
333:         PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n",(double)norm2,bs);
334:       }
335:     }

337:     /* test MatSolve() */
338:     MatSolve(sFactor,b,y);
339:     MatDestroy(&sFactor);
340:     /* Check the error */
341:     VecAXPY(y,neg_one,x);
342:     VecNorm(y,NORM_2,&norm2);
343:     if (10*norm1 < norm2 && lf-inc != -1) {
344:       PetscPrintf(PETSC_COMM_SELF,"lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n",lf-inc,lf,(double)norm1,(double)norm2);
345:     }
346:     norm1 = norm2;
347:     if (norm2 < tol && lf != -1) break;
348:   }

350: #if defined(PETSC_HAVE_MUMPS)
351:   MatGetFactor(sA,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&sFactor);
352:   MatCholeskyFactorSymbolic(sFactor,sA,NULL,NULL);
353:   MatCholeskyFactorNumeric(sFactor,sA,NULL);
354:   for (i=0; i<10; i++) {
355:     VecSetRandom(b,rdm);
356:     MatSolve(sFactor,b,y);
357:     /* Check the error */
358:     MatMult(sA,y,x);
359:     VecAXPY(x,neg_one,b);
360:     VecNorm(x,NORM_2,&norm2);
361:     if (norm2>tol) {
362:       PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve(), norm2: %g\n",(double)norm2);
363:     }
364:   }
365:   MatDestroy(&sFactor);
366: #endif

368:   ISDestroy(&perm);

370:   MatDestroy(&A);
371:   MatDestroy(&sB);
372:   MatDestroy(&sA);
373:   VecDestroy(&x);
374:   VecDestroy(&y);
375:   VecDestroy(&s1);
376:   VecDestroy(&s2);
377:   VecDestroy(&b);
378:   PetscRandomDestroy(&rdm);

380:   PetscFinalize();
381:   return 0;
382: }

384: /*TEST

386:    test:
387:       args: -bs {{1 2 3 4 5 6 7 8}}

389: TEST*/