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: PetscBool LU=PETSC_FALSE,CHOLESKY,TRIANGULAR=PETSC_FALSE,MATDSPL=PETSC_FALSE,flg,matordering;
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: PetscMPIInt size;
30: PetscInitialize(&argc,&args,(char *)0,help);
31: MPI_Comm_size(PETSC_COMM_WORLD,&size);
32: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
33: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
34: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
35: PetscOptionsGetInt(PETSC_NULL,"-lf",&lf,PETSC_NULL);
37: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&viewer1);
38: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,400,0,400,400,&viewer2);
40: MatCreate(PETSC_COMM_SELF,&C);
41: MatSetSizes(C,m*n,m*n,m*n,m*n);
42: MatSetFromOptions(C);
44: /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
45: for (i=0; i<m; i++) {
46: for (j=0; j<n; j++) {
47: v = -1.0; Ii = j + n*i;
48: J = Ii - n; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
49: J = Ii + n; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
50: J = Ii - 1; if (J>=0) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
51: J = Ii + 1; if (J<m*n) {MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
52: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
53: }
54: }
55: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
58: MatIsSymmetric(C,0.0,&flg);
59: if (!flg) SETERRQ(PETSC_COMM_SELF,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: PetscOptionsHasName(PETSC_NULL,"-mat_ordering",&matordering);
71: if (matordering){
72: MatGetOrdering(C,MATORDERINGRCM,&row,&col);
73: } else {
74: MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);
75: }
77: PetscOptionsHasName(PETSC_NULL,"-display_matrices",&MATDSPL);
78: if (MATDSPL){
79: printf("original matrix:\n");
80: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
81: MatView(C,PETSC_VIEWER_STDOUT_SELF);
82: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
83: MatView(C,PETSC_VIEWER_STDOUT_SELF);
84: MatView(C,viewer1);
85: }
87: /* Compute LU or ILU factor A */
88: MatFactorInfoInitialize(&info);
89: info.fill = 1.0;
90: info.diagonal_fill = 0;
91: info.zeropivot = 0.0;
92: PetscOptionsHasName(PETSC_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;
100: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ILU,&A);
101: MatILUFactorSymbolic(A,C,row,col,&info);
102: }
103: MatLUFactorNumeric(A,C,&info);
105: if (MATDSPL){
106: printf("factored matrix:\n");
107: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
108: MatView(A,PETSC_VIEWER_STDOUT_SELF);
109: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
110: MatView(A,PETSC_VIEWER_STDOUT_SELF);
111: MatView(A,viewer2);
112: }
114: /* Solve A*y = b, then check the error */
115: MatSolve(A,b,y);
116: VecAXPY(y,-1.0,x);
117: VecNorm(y,NORM_2,&norm2);
118: MatDestroy(&A);
120: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
121: if (!LU && lf==0){
122: MatDuplicate(C,MAT_COPY_VALUES,&A);
123: MatILUFactor(A,row,col,&info);
124: /*
125: printf("In-place factored matrix:\n");
126: MatView(C,PETSC_VIEWER_STDOUT_SELF);
127: */
128: MatSolve(A,b,y);
129: VecAXPY(y,-1.0,x);
130: VecNorm(y,NORM_2,&norm2_inplace);
131: if (PetscAbs(norm2 - norm2_inplace) > 1.e-14) SETERRQ2(PETSC_COMM_SELF,1,"ILU(0) %G and in-place ILU(0) %G give different residuals",norm2,norm2_inplace);
132: MatDestroy(&A);
133: }
135: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
136: CHOLESKY = LU;
137: if (CHOLESKY){
138: printf("Test Cholesky...\n");
139: lf = -1;
140: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
141: MatCholeskyFactorSymbolic(A,C,row,&info);
142: } else {
143: printf("Test ICC...\n");
144: info.levels = lf;
145: info.fill = 1.0;
146: info.diagonal_fill = 0;
147: info.zeropivot = 0.0;
148: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
149: MatICCFactorSymbolic(A,C,row,&info);
150: }
151: MatCholeskyFactorNumeric(A,C,&info);
153: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
154: if (lf == -1){
155: PetscOptionsHasName(PETSC_NULL,"-triangular_solve",&TRIANGULAR);
156: if (TRIANGULAR){
157: printf("Test MatForwardSolve...\n");
158: MatForwardSolve(A,b,ytmp);
159: printf("Test MatBackwardSolve...\n");
160: MatBackwardSolve(A,ytmp,y);
161: VecAXPY(y,-1.0,x);
162: VecNorm(y,NORM_2,&norm2);
163: if (norm2 > 1.e-14){
164: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
165: }
166: }
167: }
169: MatSolve(A,b,y);
170: MatDestroy(&A);
171: VecAXPY(y,-1.0,x);
172: VecNorm(y,NORM_2,&norm2);
173: if (lf == -1 && norm2 > 1.e-14){
174: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
175: }
177: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
178: if (!CHOLESKY && lf==0 && !matordering){
179: MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
180: MatICCFactor(A,row,&info);
181: /*
182: printf("In-place factored matrix:\n");
183: MatView(A,PETSC_VIEWER_STDOUT_SELF);
184: */
185: MatSolve(A,b,y);
186: VecAXPY(y,-1.0,x);
187: VecNorm(y,NORM_2,&norm2_inplace);
188: if (PetscAbs(norm2 - norm2_inplace) > 1.e-14) SETERRQ2(PETSC_COMM_SELF,1,"ICC(0) %G and in-place ICC(0) %G give different residuals",norm2,norm2_inplace);
189: MatDestroy(&A);
190: }
192: /* Free data structures */
193: ISDestroy(&row);
194: ISDestroy(&col);
195: MatDestroy(&C);
196: PetscViewerDestroy(&viewer1);
197: PetscViewerDestroy(&viewer2);
198: PetscRandomDestroy(&rdm);
199: VecDestroy(&x);
200: VecDestroy(&y);
201: VecDestroy(&ytmp);
202: VecDestroy(&b);
203: PetscFinalize();
204: return 0;
205: }