Actual source code: ex76.c


  2: static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Vec            x,y,b;
  9:   Mat            A;           /* linear system matrix */
 10:   Mat            sA,sC;       /* symmetric part of the matrices */
 11:   PetscInt       n,mbs=16,bs=1,nz=3,prob=1,i,j,col[3],block, row,Ii,J,n1,lvl;
 12:   PetscMPIInt    size;
 13:   PetscReal      norm2;
 14:   PetscScalar    neg_one = -1.0,four=4.0,value[3];
 15:   IS             perm,cperm;
 16:   PetscRandom    rdm;
 17:   PetscBool      reorder = PETSC_FALSE,displ = PETSC_FALSE;
 18:   MatFactorInfo  factinfo;
 19:   PetscBool      equal;
 20:   PetscBool      TestAIJ = PETSC_FALSE,TestBAIJ = PETSC_TRUE;
 21:   PetscInt       TestShift=0;

 23:   PetscInitialize(&argc,&args,(char*)0,help);
 24:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 26:   PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
 27:   PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);
 28:   PetscOptionsGetBool(NULL,NULL,"-reorder",&reorder,NULL);
 29:   PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL);
 30:   PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL);
 31:   PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);

 33:   n = mbs*bs;
 34:   if (TestAIJ) { /* A is in aij format */
 35:     MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A);
 36:     TestBAIJ = PETSC_FALSE;
 37:   } else { /* A is in baij format */
 38:     MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A);
 39:     TestAIJ = PETSC_FALSE;
 40:   }

 42:   /* Assemble matrix */
 43:   if (bs == 1) {
 44:     PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);
 45:     if (prob == 1) { /* tridiagonal matrix */
 46:       value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
 47:       for (i=1; i<n-1; i++) {
 48:         col[0] = i-1; col[1] = i; col[2] = i+1;
 49:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 50:       }
 51:       i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;

 53:       value[0]= 0.1; value[1]=-1; value[2]=2;
 54:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);

 56:       i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;

 58:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
 59:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 60:     } else if (prob ==2) { /* matrix for the five point stencil */
 61:       n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
 63:       for (i=0; i<n1; i++) {
 64:         for (j=0; j<n1; j++) {
 65:           Ii = j + n1*i;
 66:           if (i>0) {
 67:             J    = Ii - n1;
 68:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 69:           }
 70:           if (i<n1-1) {
 71:             J    = Ii + n1;
 72:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 73:           }
 74:           if (j>0) {
 75:             J    = Ii - 1;
 76:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 77:           }
 78:           if (j<n1-1) {
 79:             J    = Ii + 1;
 80:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 81:           }
 82:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 83:         }
 84:       }
 85:     }
 86:   } else { /* bs > 1 */
 87:     for (block=0; block<n/bs; block++) {
 88:       /* diagonal blocks */
 89:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 90:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 91:         col[0] = i-1; col[1] = i; col[2] = i+1;
 92:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 93:       }
 94:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;

 96:       value[0]=-1.0; value[1]=4.0;
 97:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);

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

101:       value[0]=4.0; value[1] = -1.0;
102:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
103:     }
104:     /* off-diagonal blocks */
105:     value[0]=-1.0;
106:     for (i=0; i<(n/bs-1)*bs; i++) {
107:       col[0]=i+bs;
108:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
109:       col[0]=i; row=i+bs;
110:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
111:     }
112:   }

114:   if (TestShift) {
115:     /* set diagonals in the 0-th block as 0 for testing shift numerical factor */
116:     for (i=0; i<bs; i++) {
117:       row  = i; col[0] = i; value[0] = 0.0;
118:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
119:     }
120:   }

122:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
123:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

125:   /* Test MatConvert */
126:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
127:   MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
128:   MatMultEqual(A,sA,20,&equal);

131:   /* Test MatGetOwnershipRange() */
132:   MatGetOwnershipRange(A,&Ii,&J);
133:   MatGetOwnershipRange(sA,&i,&j);

136:   /* Vectors */
137:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
138:   PetscRandomSetFromOptions(rdm);
139:   VecCreateSeq(PETSC_COMM_SELF,n,&x);
140:   VecDuplicate(x,&b);
141:   VecDuplicate(x,&y);
142:   VecSetRandom(x,rdm);

144:   /* Test MatReordering() - not work on sbaij matrix */
145:   if (reorder) {
146:     MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);
147:   } else {
148:     MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);
149:   }
150:   ISDestroy(&cperm);

152:   /* initialize factinfo */
153:   MatFactorInfoInitialize(&factinfo);
154:   if (TestShift == 1) {
155:     factinfo.shifttype   = (PetscReal)MAT_SHIFT_NONZERO;
156:     factinfo.shiftamount = 0.1;
157:   } else if (TestShift == 2) {
158:     factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
159:   }

161:   /* Test MatCholeskyFactor(), MatICCFactor() */
162:   /*------------------------------------------*/
163:   /* Test aij matrix A */
164:   if (TestAIJ) {
165:     if (displ) {
166:       PetscPrintf(PETSC_COMM_WORLD,"AIJ: \n");
167:     }
168:     i = 0;
169:     for (lvl=-1; lvl<10; lvl++) {
170:       if (lvl==-1) {  /* Cholesky factor */
171:         factinfo.fill = 5.0;

173:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
174:         MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
175:       } else {       /* incomplete Cholesky factor */
176:         factinfo.fill   = 5.0;
177:         factinfo.levels = lvl;

179:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
180:         MatICCFactorSymbolic(sC,A,perm,&factinfo);
181:       }
182:       MatCholeskyFactorNumeric(sC,A,&factinfo);

184:       MatMult(A,x,b);
185:       MatSolve(sC,b,y);
186:       MatDestroy(&sC);

188:       /* Check the residual */
189:       VecAXPY(y,neg_one,x);
190:       VecNorm(y,NORM_2,&norm2);

192:       if (displ) {
193:         PetscPrintf(PETSC_COMM_WORLD,"  lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
194:       }
195:     }
196:   }

198:   /* Test baij matrix A */
199:   if (TestBAIJ) {
200:     if (displ) {
201:       PetscPrintf(PETSC_COMM_WORLD,"BAIJ: \n");
202:     }
203:     i = 0;
204:     for (lvl=-1; lvl<10; lvl++) {
205:       if (lvl==-1) {  /* Cholesky factor */
206:         factinfo.fill = 5.0;

208:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
209:         MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
210:       } else {       /* incomplete Cholesky factor */
211:         factinfo.fill   = 5.0;
212:         factinfo.levels = lvl;

214:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
215:         MatICCFactorSymbolic(sC,A,perm,&factinfo);
216:       }
217:       MatCholeskyFactorNumeric(sC,A,&factinfo);

219:       MatMult(A,x,b);
220:       MatSolve(sC,b,y);
221:       MatDestroy(&sC);

223:       /* Check the residual */
224:       VecAXPY(y,neg_one,x);
225:       VecNorm(y,NORM_2,&norm2);
226:       if (displ) {
227:         PetscPrintf(PETSC_COMM_WORLD,"  lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
228:       }
229:     }
230:   }

232:   /* Test sbaij matrix sA */
233:   if (displ) {
234:     PetscPrintf(PETSC_COMM_WORLD,"SBAIJ: \n");
235:   }
236:   i = 0;
237:   for (lvl=-1; lvl<10; lvl++) {
238:     if (lvl==-1) {  /* Cholesky factor */
239:       factinfo.fill = 5.0;

241:       MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
242:       MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
243:     } else {       /* incomplete Cholesky factor */
244:       factinfo.fill   = 5.0;
245:       factinfo.levels = lvl;

247:       MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
248:       MatICCFactorSymbolic(sC,sA,perm,&factinfo);
249:     }
250:     MatCholeskyFactorNumeric(sC,sA,&factinfo);

252:     if (lvl==0 && bs==1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
253:       /*
254:         Mat B;
255:         MatDuplicate(sA,MAT_COPY_VALUES,&B);
256:         MatICCFactor(B,perm,&factinfo);
257:         MatEqual(sC,B,&equal);
258:         if (!equal) {
259:           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
260:         }
261:         MatDestroy(&B);
262:       */
263:     }

265:     MatMult(sA,x,b);
266:     MatSolve(sC,b,y);

268:     /* Test MatSolves() */
269:     if (bs == 1) {
270:       Vecs xx,bb;
271:       VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
272:       VecsDuplicate(xx,&bb);
273:       MatSolves(sC,bb,xx);
274:       VecsDestroy(xx);
275:       VecsDestroy(bb);
276:     }
277:     MatDestroy(&sC);

279:     /* Check the residual */
280:     VecAXPY(y,neg_one,x);
281:     VecNorm(y,NORM_2,&norm2);
282:     if (displ) {
283:       PetscPrintf(PETSC_COMM_WORLD,"  lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2);
284:     }
285:   }

287:   ISDestroy(&perm);
288:   MatDestroy(&A);
289:   MatDestroy(&sA);
290:   VecDestroy(&x);
291:   VecDestroy(&y);
292:   VecDestroy(&b);
293:   PetscRandomDestroy(&rdm);

295:   PetscFinalize();
296:   return 0;
297: }

299: /*TEST

301:    test:
302:       args: -bs {{1 2 3 4 5 6 7 8}}

304:    test:
305:       suffix: 3
306:       args: -testaij
307:       output_file: output/ex76_1.out

309: TEST*/