Actual source code: ex30.c

  2: static char help[] = "Tests ILU and ICC factorization with matrix ordering, and illustrates drawing of matrix sparsity structure with MatView().\n\
  3:   Input parameters are:\n\
  4:   -lf <level> : level of fill for ILU (default is 0)\n\
  5:   -lu : use full LU or Cholesky factorization\n\
  6:   -m <value>,-n <value> : grid dimensions\n\
  7: Note that most users should employ the KSP interface to the\n\
  8: linear solvers instead of using the factorization routines\n\
  9: directly.\n\n";

 11:  #include petscmat.h

 15: int main(int argc,char **args)
 16: {
 17:   Mat            C,A,sC,sA;
 18:   PetscInt       i,j,m = 5,n = 5,Ii,J,lf = 0;
 20:   PetscTruth     LU=PETSC_FALSE,flg;
 21:   PetscScalar    v;
 22:   IS             row,col;
 23:   PetscViewer    viewer1,viewer2;
 24:   MatFactorInfo  info;
 25:   Vec            x,y,b,ytmp;
 26:   PetscReal      norm2,norm2_inplace;
 27:   PetscRandom    rdm;
 28:   PetscInt       *ii;

 30:   PetscInitialize(&argc,&args,(char *)0,help);
 31:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
 32:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 33:   PetscOptionsGetInt(PETSC_NULL,"-lf",&lf,PETSC_NULL);

 35:   PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&viewer1);
 36:   PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,400,0,400,400,&viewer2);

 38:   MatCreate(PETSC_COMM_SELF,&C);
 39:   MatSetSizes(C,m*n,m*n,m*n,m*n);
 40:   MatSetFromOptions(C);

 42:   /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
 43:   for (i=0; i<m; i++) {
 44:     for (j=0; j<n; j++) {
 45:       v = -1.0;  Ii = j + n*i;
 46:       J = Ii - n; if (J>=0)  {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 47:       J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 48:       J = Ii - 1; if (J>=0)  {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 49:       J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 50:       v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 51:     }
 52:   }
 53:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 54:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 56:   MatConvert(C,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sC);

 58:   MatIsSymmetric(C,0.0,&flg);
 59:   if (!flg) SETERRQ(1,"C is non-symmetric");

 61:   /* Create vectors for error checking */
 62:   MatGetVecs(C,&x,&b);
 63:   VecDuplicate(x,&y);
 64:   VecDuplicate(x,&ytmp);
 65:   PetscRandomCreate(PETSC_COMM_SELF,&rdm);
 66:   PetscRandomSetFromOptions(rdm);
 67:   VecSetRandom(x,rdm);
 68:   MatMult(C,x,b);

 70:   MatGetOrdering(C,MATORDERING_RCM,&row,&col);
 71:   /* replace row or col with natural ordering for testing */
 72:   PetscOptionsHasName(PETSC_NULL,"-no_rowperm",&flg);
 73:   if (flg){
 74:     ISDestroy(row);
 75:     PetscMalloc(m*n*sizeof(PetscInt),&ii);
 76:     for (i=0; i<m*n; i++) ii[i] = i;
 77:     ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&row);
 78:     PetscFree(ii);
 79:     ISSetIdentity(row);
 80:     ISSetPermutation(row);
 81:   }
 82:   PetscOptionsHasName(PETSC_NULL,"-no_colperm",&flg);
 83:   if (flg){
 84:     ISDestroy(col);
 85:     PetscMalloc(m*n*sizeof(PetscInt),&ii);
 86:     for (i=0; i<m*n; i++) ii[i] = i;
 87:     ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&col);
 88:     PetscFree(ii);
 89:     ISSetIdentity(col);
 90:     ISSetPermutation(col);
 91:   }

 93:   printf("original matrix:\n");
 94:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
 95:   MatView(C,PETSC_VIEWER_STDOUT_SELF);
 96:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
 97:   MatView(C,PETSC_VIEWER_STDOUT_SELF);
 98:   MatView(C,viewer1);

100:   /* Compute LU or ILU factor A */
101:   MatFactorInfoInitialize(&info);
102:   info.fill          = 1.0;
103:   info.diagonal_fill = 0;
104:   info.shiftnz       = 0;
105:   info.zeropivot     = 0.0;
106:   PetscOptionsHasName(PETSC_NULL,"-lu",&LU);
107:   if (LU){
108:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&A);
109:     MatLUFactorSymbolic(A,C,row,col,&info);
110:   } else {
111:     info.levels = lf;
112:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_ILU,&A);
113:     MatILUFactorSymbolic(A,C,row,col,&info);
114:   }
115:   MatLUFactorNumeric(A,C,&info);

117:   printf("factored matrix:\n");
118:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
119:   MatView(A,PETSC_VIEWER_STDOUT_SELF);
120:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
121:   MatView(A,PETSC_VIEWER_STDOUT_SELF);
122:   MatView(A,viewer2);

124:   /* Solve A*y = b, then check the error */
125:   MatSolve(A,b,y);
126:   VecAXPY(y,-1.0,x);
127:   VecNorm(y,NORM_2,&norm2);
128:   MatDestroy(A);

130:   /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
131:   if (!LU && lf==0){
132:     MatDuplicate(C,MAT_COPY_VALUES,&A);
133:     MatILUFactor(A,row,col,&info);
134:     /*
135:     printf("In-place factored matrix:\n");
136:     MatView(C,PETSC_VIEWER_STDOUT_SELF);
137:     */
138:     MatSolve(A,b,y);
139:     VecAXPY(y,-1.0,x);
140:     VecNorm(y,NORM_2,&norm2_inplace);
141:     if (PetscAbs(norm2 - norm2_inplace) > 1.e-16) SETERRQ2(1,"ILU(0) %G and in-place ILU(0) %G give different residuals",norm2,norm2_inplace);
142:     MatDestroy(A);
143:   }

145:   /* Test Cholesky and ICC on seqaij matrix with matrix reordering */
146:   if (LU){
147:     lf = -1;
148:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&A);
149:     MatCholeskyFactorSymbolic(A,C,row,&info);
150:   } else {
151:     info.levels        = lf;
152:     info.fill          = 1.0;
153:     info.diagonal_fill = 0;
154:     info.shiftnz       = 0;
155:     info.zeropivot     = 0.0;
156:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&A);
157:     MatICCFactorSymbolic(A,C,row,&info);
158:   }
159:   MatCholeskyFactorNumeric(A,C,&info);

161:   /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering */
162:   if (lf == -1){
163:     MatForwardSolve(A,b,ytmp);
164:     MatBackwardSolve(A,ytmp,y);
165:     VecAXPY(y,-1.0,x);
166:     VecNorm(y,NORM_2,&norm2);
167:     if (norm2 > 1.e-14){
168:       PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
169:     }
170:   }

172:   MatSolve(A,b,y);
173:   MatDestroy(A);
174:   VecAXPY(y,-1.0,x);
175:   VecNorm(y,NORM_2,&norm2);
176:   if (lf == -1 && norm2 > 1.e-14){
177:     PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ:   Cholesky/ICC levels %d, residual %g\n",lf,norm2);
178:   }
179: 
180:   /* Test Cholesky and ICC on seqaij matrix without matrix reordering */
181:   ISDestroy(row);
182:   ISDestroy(col);
183:   MatGetOrdering(C,MATORDERING_NATURAL,&row,&col);
184:   if (LU){
185:     lf = -1;
186:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&A);
187:     MatCholeskyFactorSymbolic(A,C,row,&info);
188:   } else {
189:     info.levels        = lf;
190:     info.fill          = 1.0;
191:     info.diagonal_fill = 0;
192:     info.shiftnz       = 0;
193:     info.zeropivot     = 0.0;
194:     MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&A);
195:     MatICCFactorSymbolic(A,C,row,&info);
196:   }
197:   MatCholeskyFactorNumeric(A,C,&info);

199:   /* test MatForwardSolve() and MatBackwardSolve() */
200:   if (lf == -1){
201:     MatForwardSolve(A,b,ytmp);
202:     MatBackwardSolve(A,ytmp,y);
203:     VecAXPY(y,-1.0,x);
204:     VecNorm(y,NORM_2,&norm2);
205:     if (norm2 > 1.e-14){
206:       PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
207:     }
208:   }

210:   /* Test MatSolve() */
211:   MatSolve(A,b,y);
212:   VecAXPY(y,-1.0,x);
213:   VecNorm(y,NORM_2,&norm2);
214:   if (lf == -1 && norm2 > 1.e-14){
215:     printf(" SEQAIJ:   Cholesky/ICC levels %d, residual %g\n",lf,norm2);
216:   }

218:   /* Test Cholesky and ICC on seqsbaij matrix without matrix reordering */
219:   if (LU){
220:     MatGetFactor(sC,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&sA);
221:     MatCholeskyFactorSymbolic(sA,sC,row,&info);
222:   } else {
223:     MatGetFactor(sC,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&sA);
224:     MatICCFactorSymbolic(sA,sC,row,&info);
225:   }
226:   MatCholeskyFactorNumeric(sA,sC,&info);
227:   MatEqual(A,sA,&flg);
228:   if (!flg) SETERRQ(1,"CholeskyFactors for aij and sbaij matrices are different");
229:   MatDestroy(sC);
230:   MatDestroy(sA);
231:   MatDestroy(A);

233:   /* Free data structures */
234:   MatDestroy(C);
235:   ISDestroy(row);
236:   ISDestroy(col);
237:   PetscViewerDestroy(viewer1);
238:   PetscViewerDestroy(viewer2);
239:   PetscRandomDestroy(rdm);
240:   VecDestroy(x);
241:   VecDestroy(y);
242:   VecDestroy(ytmp);
243:   VecDestroy(b);
244:   PetscFinalize();
245:   return 0;
246: }