Actual source code: ex30.c


  2: static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqaij format, 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>

 13: int main(int argc,char **args)
 14: {
 15:   Mat            C,A;
 16:   PetscInt       i,j,m = 5,n = 5,Ii,J,lf = 0;
 17:   PetscBool      LU=PETSC_FALSE,CHOLESKY,TRIANGULAR=PETSC_FALSE,MATDSPL=PETSC_FALSE,flg,matordering;
 18:   PetscScalar    v;
 19:   IS             row,col;
 20:   PetscViewer    viewer1,viewer2;
 21:   MatFactorInfo  info;
 22:   Vec            x,y,b,ytmp;
 23:   PetscReal      norm2,norm2_inplace, tol = 100.*PETSC_MACHINE_EPSILON;
 24:   PetscRandom    rdm;
 25:   PetscMPIInt    size;

 27:   PetscInitialize(&argc,&args,(char*)0,help);
 28:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 30:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 31:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 32:   PetscOptionsGetInt(NULL,NULL,"-lf",&lf,NULL);

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

 37:   MatCreate(PETSC_COMM_SELF,&C);
 38:   MatSetSizes(C,m*n,m*n,m*n,m*n);
 39:   MatSetFromOptions(C);
 40:   MatSetUp(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:   MatIsSymmetric(C,0.0,&flg);

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

 68:   PetscOptionsHasName(NULL,NULL,"-mat_ordering",&matordering);
 69:   if (matordering) {
 70:     MatGetOrdering(C,MATORDERINGRCM,&row,&col);
 71:   } else {
 72:     MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);
 73:   }

 75:   PetscOptionsHasName(NULL,NULL,"-display_matrices",&MATDSPL);
 76:   if (MATDSPL) {
 77:     printf("original matrix:\n");
 78:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
 79:     MatView(C,PETSC_VIEWER_STDOUT_SELF);
 80:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
 81:     MatView(C,PETSC_VIEWER_STDOUT_SELF);
 82:     MatView(C,viewer1);
 83:   }

 85:   /* Compute LU or ILU factor A */
 86:   MatFactorInfoInitialize(&info);

 88:   info.fill          = 1.0;
 89:   info.diagonal_fill = 0;
 90:   info.zeropivot     = 0.0;

 92:   PetscOptionsHasName(NULL,NULL,"-lu",&LU);
 93:   if (LU) {
 94:     printf("Test LU...\n");
 95:     MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A);
 96:     MatLUFactorSymbolic(A,C,row,col,&info);
 97:   } else {
 98:     printf("Test ILU...\n");
 99:     info.levels = lf;

101:     MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ILU,&A);
102:     MatILUFactorSymbolic(A,C,row,col,&info);
103:   }
104:   MatLUFactorNumeric(A,C,&info);

106:   /* Solve A*y = b, then check the error */
107:   MatSolve(A,b,y);
108:   VecAXPY(y,-1.0,x);
109:   VecNorm(y,NORM_2,&norm2);
110:   MatDestroy(&A);

112:   /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
113:   if (!LU && lf==0) {
114:     MatDuplicate(C,MAT_COPY_VALUES,&A);
115:     MatILUFactor(A,row,col,&info);
116:     /*
117:     printf("In-place factored matrix:\n");
118:     MatView(C,PETSC_VIEWER_STDOUT_SELF);
119:     */
120:     MatSolve(A,b,y);
121:     VecAXPY(y,-1.0,x);
122:     VecNorm(y,NORM_2,&norm2_inplace);
124:     MatDestroy(&A);
125:   }

127:   /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
128:   CHOLESKY = LU;
129:   if (CHOLESKY) {
130:     printf("Test Cholesky...\n");
131:     lf   = -1;
132:     MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
133:     MatCholeskyFactorSymbolic(A,C,row,&info);
134:   } else {
135:     printf("Test ICC...\n");
136:     info.levels        = lf;
137:     info.fill          = 1.0;
138:     info.diagonal_fill = 0;
139:     info.zeropivot     = 0.0;

141:     MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
142:     MatICCFactorSymbolic(A,C,row,&info);
143:   }
144:   MatCholeskyFactorNumeric(A,C,&info);

146:   /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
147:   if (lf == -1) {
148:     PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR);
149:     if (TRIANGULAR) {
150:       printf("Test MatForwardSolve...\n");
151:       MatForwardSolve(A,b,ytmp);
152:       printf("Test MatBackwardSolve...\n");
153:       MatBackwardSolve(A,ytmp,y);
154:       VecAXPY(y,-1.0,x);
155:       VecNorm(y,NORM_2,&norm2);
156:       if (norm2 > tol) {
157:         PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2);
158:       }
159:     }
160:   }

162:   MatSolve(A,b,y);
163:   MatDestroy(&A);
164:   VecAXPY(y,-1.0,x);
165:   VecNorm(y,NORM_2,&norm2);
166:   if (lf == -1 && norm2 > tol) {
167:     PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ:   Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n",lf,(double)norm2);
168:   }

170:   /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
171:   if (!CHOLESKY && lf==0 && !matordering) {
172:     MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
173:     MatICCFactor(A,row,&info);
174:     /*
175:     printf("In-place factored matrix:\n");
176:     MatView(A,PETSC_VIEWER_STDOUT_SELF);
177:     */
178:     MatSolve(A,b,y);
179:     VecAXPY(y,-1.0,x);
180:     VecNorm(y,NORM_2,&norm2_inplace);
182:     MatDestroy(&A);
183:   }

185:   /* Free data structures */
186:   ISDestroy(&row);
187:   ISDestroy(&col);
188:   MatDestroy(&C);
189:   PetscViewerDestroy(&viewer1);
190:   PetscViewerDestroy(&viewer2);
191:   PetscRandomDestroy(&rdm);
192:   VecDestroy(&x);
193:   VecDestroy(&y);
194:   VecDestroy(&ytmp);
195:   VecDestroy(&b);
196:   PetscFinalize();
197:   return 0;
198: }

200: /*TEST

202:    test:
203:       args: -mat_ordering -display_matrices -nox
204:       filter: grep -v "MPI processes"

206:    test:
207:       suffix: 2
208:       args: -mat_ordering -display_matrices -nox -lu

210:    test:
211:       suffix: 3
212:       args: -mat_ordering -lu -triangular_solve

214:    test:
215:       suffix: 4

217:    test:
218:       suffix: 5
219:       args: -lu

221:    test:
222:       suffix: 6
223:       args: -lu -triangular_solve
224:       output_file: output/ex30_3.out

226: TEST*/