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:   PetscTruth     equal;
 23:   PetscTruth     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(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:   PetscOptionsGetTruth(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) (sqrt((PetscReal)n) + 0.001);
 63:       if (n1*n1 - n) SETERRQ(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:   MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
126:   MatMultEqual(A,sA,20,&equal);
127:   if (!equal) SETERRQ(PETSC_ERR_USER,"A != sA");

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

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,MATORDERING_RCM,&perm,&cperm);
147:   } else {
148:     MatGetOrdering(A,MATORDERING_NATURAL,&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:   }
160: 
161:   /* Test MatCholeskyFactor(), MatICCFactor() */
162:   /*------------------------------------------*/
163:   /* Test aij matrix A */
164:   if (TestAIJ){
165:     if (displ>0){PetscPrintf(PETSC_COMM_SELF,"AIJ: \n");}
166:     i = 0;
167:     for (lvl=-1; lvl<10; lvl++){
168:       if (lvl==-1) {  /* Cholesky factor */
169:         factinfo.fill = 5.0;
170:         MatGetFactor(A,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&sC);
171:         MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);
172:       } else {       /* incomplete Cholesky factor */
173:         factinfo.fill   = 5.0;
174:         factinfo.levels = lvl;
175:         MatGetFactor(A,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&sC);
176:         MatICCFactorSymbolic(sC,A,perm,&factinfo);
177:       }
178:       MatCholeskyFactorNumeric(sC,A,&factinfo);

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

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

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

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

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

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


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

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

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

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

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