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
15: int main(int argc,char **args)
16: {
17: Mat C,A;
18: PetscInt i,j,m = 5,n = 5,Ii,J,lf = 0;
20: PetscTruth LU=PETSC_FALSE,CHOLESKY,TRIANGULAR=PETSC_FALSE,MATDSPL=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;
29: PetscMPIInt size;
31: PetscInitialize(&argc,&args,(char *)0,help);
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");
34: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
35: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
36: PetscOptionsGetInt(PETSC_NULL,"-lf",&lf,PETSC_NULL);
38: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&viewer1);
39: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,400,0,400,400,&viewer2);
41: MatCreate(PETSC_COMM_SELF,&C);
42: MatSetSizes(C,m*n,m*n,m*n,m*n);
43: MatSetFromOptions(C);
45: /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
46: for (i=0; i<m; i++) {
47: for (j=0; j<n; j++) {
48: v = -1.0; Ii = j + n*i;
49: J = Ii - n; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
50: J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
51: J = Ii - 1; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
52: J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
53: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
54: }
55: }
56: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
59: MatIsSymmetric(C,0.0,&flg);
60: if (!flg) SETERRQ(1,"C is non-symmetric");
62: /* Create vectors for error checking */
63: MatGetVecs(C,&x,&b);
64: VecDuplicate(x,&y);
65: VecDuplicate(x,&ytmp);
66: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
67: PetscRandomSetFromOptions(rdm);
68: VecSetRandom(x,rdm);
69: MatMult(C,x,b);
71: MatGetOrdering(C,MATORDERING_RCM,&row,&col);
72: /* replace row or col with natural ordering for testing */
73: PetscOptionsHasName(PETSC_NULL,"-no_rowperm",&flg);
74: if (flg){
75: ISDestroy(row);
76: PetscMalloc(m*n*sizeof(PetscInt),&ii);
77: for (i=0; i<m*n; i++) ii[i] = i;
78: ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&row);
79: PetscFree(ii);
80: ISSetIdentity(row);
81: ISSetPermutation(row);
82: }
83: PetscOptionsHasName(PETSC_NULL,"-no_colperm",&flg);
84: if (flg){
85: ISDestroy(col);
86: PetscMalloc(m*n*sizeof(PetscInt),&ii);
87: for (i=0; i<m*n; i++) ii[i] = i;
88: ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&col);
89: PetscFree(ii);
90: ISSetIdentity(col);
91: ISSetPermutation(col);
92: }
94: PetscOptionsHasName(PETSC_NULL,"-display_matrices",&MATDSPL);
95: if (MATDSPL){
96: printf("original matrix:\n");
97: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
98: MatView(C,PETSC_VIEWER_STDOUT_SELF);
99: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
100: MatView(C,PETSC_VIEWER_STDOUT_SELF);
101: MatView(C,viewer1);
102: }
104: /* Compute LU or ILU factor A */
105: MatFactorInfoInitialize(&info);
106: info.fill = 1.0;
107: info.diagonal_fill = 0;
108: info.zeropivot = 0.0;
109: PetscOptionsHasName(PETSC_NULL,"-lu",&LU);
110: if (LU){
111: printf("Test LU...\n");
112: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&A);
113: MatLUFactorSymbolic(A,C,row,col,&info);
114: } else {
115: printf("Test ILU...\n");
116: info.levels = lf;
117: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_ILU,&A);
118: MatILUFactorSymbolic(A,C,row,col,&info);
119: }
120: MatLUFactorNumeric(A,C,&info);
122: if (MATDSPL){
123: printf("factored matrix:\n");
124: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
125: MatView(A,PETSC_VIEWER_STDOUT_SELF);
126: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
127: MatView(A,PETSC_VIEWER_STDOUT_SELF);
128: MatView(A,viewer2);
129: }
131: /* Solve A*y = b, then check the error */
132: MatSolve(A,b,y);
133: VecAXPY(y,-1.0,x);
134: VecNorm(y,NORM_2,&norm2);
135: MatDestroy(A);
137: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
138: if (!LU && lf==0){
139: MatDuplicate(C,MAT_COPY_VALUES,&A);
140: MatILUFactor(A,row,col,&info);
141: /*
142: printf("In-place factored matrix:\n");
143: MatView(C,PETSC_VIEWER_STDOUT_SELF);
144: */
145: MatSolve(A,b,y);
146: VecAXPY(y,-1.0,x);
147: VecNorm(y,NORM_2,&norm2_inplace);
148: 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);
149: MatDestroy(A);
150: }
152: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
153: CHOLESKY = LU;
154: if (CHOLESKY){
155: printf("Test Cholesky...\n");
156: lf = -1;
157: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_CHOLESKY,&A);
158: MatCholeskyFactorSymbolic(A,C,row,&info);
159: } else {
160: printf("Test ICC...\n");
161: info.levels = lf;
162: info.fill = 1.0;
163: info.diagonal_fill = 0;
164: info.zeropivot = 0.0;
165: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_ICC,&A);
166: MatICCFactorSymbolic(A,C,row,&info);
167: }
168: MatCholeskyFactorNumeric(A,C,&info);
170: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
171: if (lf == -1){
172: PetscOptionsHasName(PETSC_NULL,"-triangular_solve",&TRIANGULAR);
173: if (TRIANGULAR){
174: printf("Test MatForwardSolve...\n");
175: MatForwardSolve(A,b,ytmp);
176: printf("Test MatBackwardSolve...\n");
177: MatBackwardSolve(A,ytmp,y);
178: VecAXPY(y,-1.0,x);
179: VecNorm(y,NORM_2,&norm2);
180: if (norm2 > 1.e-14){
181: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
182: }
183: }
184: }
186: MatSolve(A,b,y);
187: MatDestroy(A);
188: VecAXPY(y,-1.0,x);
189: VecNorm(y,NORM_2,&norm2);
190: if (lf == -1 && norm2 > 1.e-14){
191: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
192: }
193: ISDestroy(row);
194: ISDestroy(col);
196: /* Free data structures */
197: MatDestroy(C);
198: PetscViewerDestroy(viewer1);
199: PetscViewerDestroy(viewer2);
200: PetscRandomDestroy(rdm);
201: VecDestroy(x);
202: VecDestroy(y);
203: VecDestroy(ytmp);
204: VecDestroy(b);
205: PetscFinalize();
206: return 0;
207: }