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>

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

 26:   PetscInitialize(&argc,&args,(char *)0,help);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 28:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
 29:   PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 30:   PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
 31:   PetscOptionsGetInt(PETSC_NULL,"-reorder",&reorder,PETSC_NULL);
 32:   PetscOptionsGetBool(PETSC_NULL,"-testaij",&TestAIJ,PETSC_NULL);
 33:   PetscOptionsGetInt(PETSC_NULL,"-testShift",&TestShift,PETSC_NULL);
 34:   PetscOptionsGetInt(PETSC_NULL,"-displ",&displ,PETSC_NULL);

 36:   n = mbs*bs;
 37:   if (TestAIJ){ /* A is in aij format */
 38:     MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,PETSC_NULL,&A);
 39:     TestBAIJ = PETSC_FALSE;
 40:   } else { /* A is in baij format */
 41:     ierr=MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,PETSC_NULL,&A);
 42:     TestAIJ = PETSC_FALSE;
 43:   }

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

 58:       i = 0; col[0] = 0; col[1] = 1; col[2]=n-1;
 59:       value[0] = 2.0; value[1] = -1.0; value[2]=0.1;
 60:       MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 61:     } else if (prob ==2){ /* matrix for the five point stencil */
 62:       n1 = (int) (PetscSqrtReal((PetscReal)n) + 0.001);
 63:       if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
 64:       for (i=0; i<n1; i++) {
 65:         for (j=0; j<n1; j++) {
 66:           Ii = j + n1*i;
 67:           if (i>0)   {
 68:             J = Ii - n1;
 69:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 70:           }
 71:           if (i<n1-1) {
 72:             J = Ii + n1;
 73:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 74:           }
 75:           if (j>0)   {
 76:             J = Ii - 1;
 77:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 78:           }
 79:           if (j<n1-1) {
 80:             J = Ii + 1;
 81:             MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);
 82:           }
 83:           MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);
 84:         }
 85:       }
 86:     }
 87:   } else { /* bs > 1 */
 88:     for (block=0; block<n/bs; block++){
 89:       /* diagonal blocks */
 90:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 91:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 92:         col[0] = i-1; col[1] = i; col[2] = i+1;
 93:         MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 94:       }
 95:       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;
100:       value[0]=4.0; value[1] = -1.0;
101:       MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
102:     }
103:     /* off-diagonal blocks */
104:     value[0]=-1.0;
105:     for (i=0; i<(n/bs-1)*bs; i++){
106:       col[0]=i+bs;
107:       MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);
108:       col[0]=i; row=i+bs;
109:       MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);
110:     }
111:   }

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

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

124:   /* Test MatConvert */
125:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
126:   MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
127:   MatMultEqual(A,sA,20,&equal);
128:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA");

130:   /* Test MatGetOwnershipRange() */
131:   MatGetOwnershipRange(A,&Ii,&J);
132:   MatGetOwnershipRange(sA,&i,&j);
133:   if (i-Ii || j-J){
134:     PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");
135:   }

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

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

153:   /* initialize factinfo */
154:   MatFactorInfoInitialize(&factinfo);
155:   if (TestShift == 1){
156:     factinfo.shifttype   = (PetscReal)MAT_SHIFT_NONZERO;
157:     factinfo.shiftamount = 0.1;
158:   } else if (TestShift == 2){
159:     factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE;
160:   }
161: 
162:   /* Test MatCholeskyFactor(), MatICCFactor() */
163:   /*------------------------------------------*/
164:   /* Test aij matrix A */
165:   if (TestAIJ){
166:     if (displ>0){PetscPrintf(PETSC_COMM_SELF,"AIJ: \n");}
167:     i = 0;
168:     for (lvl=-1; lvl<10; lvl++){
169:       if (lvl==-1) {  /* Cholesky factor */
170:         factinfo.fill = 5.0;
171:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
172:         MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
173:       } else {       /* incomplete Cholesky factor */
174:         factinfo.fill   = 5.0;
175:         factinfo.levels = lvl;
176:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
177:         MatICCFactorSymbolic(sC,A,perm,&factinfo);
178:       }
179:       MatCholeskyFactorNumeric(sC,A,&factinfo);

181:       MatMult(A,x,b);
182:       MatSolve(sC,b,y);
183:       MatDestroy(&sC);

185:       /* Check the error */
186:       VecAXPY(y,neg_one,x);
187:       VecNorm(y,NORM_2,&norm2);
188:       if (displ>0){PetscPrintf(PETSC_COMM_SELF,"  lvl: %d, error: %G\n", lvl,norm2);}
189:       err[i++] = norm2;
190:     }
191:   }
192: 
193:   /* Test baij matrix A */
194:   if (TestBAIJ){
195:     if (displ>0){PetscPrintf(PETSC_COMM_SELF,"BAIJ: \n");}
196:     i = 0;
197:     for (lvl=-1; lvl<10; lvl++){
198:       if (lvl==-1) {  /* Cholesky factor */
199:         factinfo.fill = 5.0;
200:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
201:         MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
202:       } else {       /* incomplete Cholesky factor */
203:         factinfo.fill   = 5.0;
204:         factinfo.levels = lvl;
205:         MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
206:         MatICCFactorSymbolic(sC,A,perm,&factinfo);
207:       }
208:       MatCholeskyFactorNumeric(sC,A,&factinfo);

210:       MatMult(A,x,b);
211:       MatSolve(sC,b,y);
212:       MatDestroy(&sC);

214:       /* Check the error */
215:       VecAXPY(y,neg_one,x);
216:       VecNorm(y,NORM_2,&norm2);
217:       if (displ>0){PetscPrintf(PETSC_COMM_SELF,"  lvl: %d, error: %G\n", lvl,norm2);}
218:       err[i++] = norm2;
219:     }
220:   }

222:   /* Test sbaij matrix sA */
223:   if (displ>0){PetscPrintf(PETSC_COMM_SELF,"SBAIJ: \n");}
224:   i = 0;
225:   for (lvl=-1; lvl<10; lvl++){
226:     if (lvl==-1) {  /* Cholesky factor */
227:       factinfo.fill = 5.0;
228:       MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);
229:       MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);
230:     } else {       /* incomplete Cholesky factor */
231:       factinfo.fill   = 5.0;
232:       factinfo.levels = lvl;
233:       MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);
234:       MatICCFactorSymbolic(sC,sA,perm,&factinfo);
235:     }
236:     MatCholeskyFactorNumeric(sC,sA,&factinfo);

238:     if (lvl==0 && bs==1){ /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */
239:       /*
240:         Mat B;
241:         MatDuplicate(sA,MAT_COPY_VALUES,&B);
242:         MatICCFactor(B,perm,&factinfo);
243:         MatEqual(sC,B,&equal);
244:         if (!equal){
245:           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor");
246:         }
247:         MatDestroy(&B);
248:       */
249:     }


252:     MatMult(sA,x,b);
253:     MatSolve(sC,b,y);

255:     /* Test MatSolves() */
256:     if (bs == 1) {
257:       Vecs xx,bb;
258:       VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);
259:       VecsDuplicate(xx,&bb);
260:       MatSolves(sC,bb,xx);
261:       VecsDestroy(xx);
262:       VecsDestroy(bb);
263:     }
264:     MatDestroy(&sC);

266:     /* Check the error */
267:     VecAXPY(y,neg_one,x);
268:     VecNorm(y,NORM_2,&norm2);
269:     if (displ>0){PetscPrintf(PETSC_COMM_SELF,"  lvl: %d, error: %G\n", lvl,norm2); }
270:     err[i] -= norm2;
271:     if (err[i] > tol) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER," level: %d, err: %G\n", lvl,err[i]);
272:   }

274:   ISDestroy(&perm);
275:   MatDestroy(&A);
276:   MatDestroy(&sA);
277:   VecDestroy(&x);
278:   VecDestroy(&y);
279:   VecDestroy(&b);
280:   PetscRandomDestroy(&rdm);

282:   PetscFinalize();
283:   return 0;
284: }