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;
29: PetscInitialize(&argc,&args,(char *)0,help);
30: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
31: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
32: PetscOptionsGetInt(PETSC_NULL,"-lf",&lf,PETSC_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: MatSeqBDiagSetPreallocation(C,0,1,PETSC_NULL,PETSC_NULL);
41: MatSeqAIJSetPreallocation(C,5,PETSC_NULL);
43: /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
44: for (i=0; i<m; i++) {
45: for (j=0; j<n; j++) {
46: v = -1.0; Ii = j + n*i;
47: J = Ii - n; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
48: J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
49: J = Ii - 1; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
50: J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
51: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
52: }
53: }
54: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
57: MatConvert(C,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sC);
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: PetscInt *ii;
77: PetscMalloc(m*n*sizeof(PetscInt),&ii);
78: for (i=0; i<m*n; i++) ii[i] = i;
79: ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&row);
80: PetscFree(ii);
81: ISSetIdentity(row);
82: ISSetPermutation(row);
83: }
84: PetscOptionsHasName(PETSC_NULL,"-no_colperm",&flg);
85: if (flg){
86: ISDestroy(col);
87: PetscInt *ii;
88: PetscMalloc(m*n*sizeof(PetscInt),&ii);
89: for (i=0; i<m*n; i++) ii[i] = i;
90: ISCreateGeneral(PETSC_COMM_SELF,m*n,ii,&col);
91: PetscFree(ii);
92: ISSetIdentity(col);
93: ISSetPermutation(col);
94: }
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);
103: /* Compute LU or ILU factor A */
104: MatFactorInfoInitialize(&info);
105: info.fill = 1.0;
106: info.diagonal_fill = 0;
107: info.shiftnz = 0;
108: info.zeropivot = 0.0;
109: PetscOptionsHasName(PETSC_NULL,"-lu",&LU);
110: if (LU){
111: MatLUFactorSymbolic(C,row,col,&info,&A);
112: } else {
113: info.levels = lf;
114: MatILUFactorSymbolic(C,row,col,&info,&A);
115: }
116: MatLUFactorNumeric(C,&info,&A);
118: printf("factored matrix:\n");
119: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
120: MatView(A,PETSC_VIEWER_STDOUT_SELF);
121: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
122: MatView(A,PETSC_VIEWER_STDOUT_SELF);
123: MatView(A,viewer2);
125: /* Solve A*y = b, then check the error */
126: MatSolve(A,b,y);
127: VecAXPY(y,-1.0,x);
128: VecNorm(y,NORM_2,&norm2);
129: MatDestroy(A);
131: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
132: if (!LU && lf==0){
133: MatDuplicate(C,MAT_COPY_VALUES,&A);
134: MatILUFactor(A,row,col,&info);
135: /*
136: printf("In-place factored matrix:\n");
137: MatView(C,PETSC_VIEWER_STDOUT_SELF);
138: */
139: MatSolve(A,b,y);
140: VecAXPY(y,-1.0,x);
141: VecNorm(y,NORM_2,&norm2_inplace);
142: 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);
143: MatDestroy(A);
144: }
146: /* Test Cholesky and ICC on seqaij matrix with matrix reordering */
147: if (LU){
148: lf = -1;
149: MatCholeskyFactorSymbolic(C,row,&info,&A);
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: MatICCFactorSymbolic(C,row,&info,&A);
157: }
158: MatCholeskyFactorNumeric(C,&info,&A);
160: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering */
161: if (lf == -1){
162: MatForwardSolve(A,b,ytmp);
163: MatBackwardSolve(A,ytmp,y);
164: VecAXPY(y,-1.0,x);
165: VecNorm(y,NORM_2,&norm2);
166: if (norm2 > 1.e-14){
167: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
168: }
169: }
171: MatSolve(A,b,y);
172: MatDestroy(A);
173: VecAXPY(y,-1.0,x);
174: VecNorm(y,NORM_2,&norm2);
175: if (lf == -1 && norm2 > 1.e-14){
176: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
177: }
178:
179: /* Test Cholesky and ICC on seqaij matrix without matrix reordering */
180: ISDestroy(row);
181: ISDestroy(col);
182: MatGetOrdering(C,MATORDERING_NATURAL,&row,&col);
183: if (LU){
184: lf = -1;
185: MatCholeskyFactorSymbolic(C,row,&info,&A);
186: } else {
187: info.levels = lf;
188: info.fill = 1.0;
189: info.diagonal_fill = 0;
190: info.shiftnz = 0;
191: info.zeropivot = 0.0;
192: MatICCFactorSymbolic(C,row,&info,&A);
193: }
194: MatCholeskyFactorNumeric(C,&info,&A);
196: /* test MatForwardSolve() and MatBackwardSolve() */
197: if (lf == -1){
198: MatForwardSolve(A,b,ytmp);
199: MatBackwardSolve(A,ytmp,y);
200: VecAXPY(y,-1.0,x);
201: VecNorm(y,NORM_2,&norm2);
202: if (norm2 > 1.e-14){
203: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
204: }
205: }
207: /* Test MatSolve() */
208: MatSolve(A,b,y);
209: VecAXPY(y,-1.0,x);
210: VecNorm(y,NORM_2,&norm2);
211: if (lf == -1 && norm2 > 1.e-14){
212: printf(" SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
213: }
215: /* Test Cholesky and ICC on seqsbaij matrix without matrix reordering */
216: if (LU){
217: MatCholeskyFactorSymbolic(sC,row,&info,&sA);
218: } else {
219: MatICCFactorSymbolic(sC,row,&info,&sA);
220: }
221: MatCholeskyFactorNumeric(sC,&info,&sA);
222: MatEqual(A,sA,&flg);
223: if (!flg) SETERRQ(1,"CholeskyFactors for aij and sbaij matrices are different");
224: MatDestroy(sC);
225: MatDestroy(sA);
226: MatDestroy(A);
228: /* Free data structures */
229: MatDestroy(C);
230: ISDestroy(row);
231: ISDestroy(col);
232: PetscViewerDestroy(viewer1);
233: PetscViewerDestroy(viewer2);
234: PetscRandomDestroy(rdm);
235: VecDestroy(x);
236: VecDestroy(y);
237: VecDestroy(ytmp);
238: VecDestroy(b);
239: PetscFinalize();
240: return 0;
241: }