Actual source code: ex30.c
petsc-3.3-p6 2013-02-11
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);
43: MatSetUp(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(PETSC_COMM_SELF,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: PetscOptionsHasName(PETSC_NULL,"-mat_ordering",&matordering);
72: if (matordering){
73: MatGetOrdering(C,MATORDERINGRCM,&row,&col);
74: } else {
75: MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);
76: }
78: PetscOptionsHasName(PETSC_NULL,"-display_matrices",&MATDSPL);
79: if (MATDSPL){
80: printf("original matrix:\n");
81: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
82: MatView(C,PETSC_VIEWER_STDOUT_SELF);
83: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
84: MatView(C,PETSC_VIEWER_STDOUT_SELF);
85: MatView(C,viewer1);
86: }
88: /* Compute LU or ILU factor A */
89: MatFactorInfoInitialize(&info);
90: info.fill = 1.0;
91: info.diagonal_fill = 0;
92: info.zeropivot = 0.0;
93: PetscOptionsHasName(PETSC_NULL,"-lu",&LU);
94: if (LU){
95: printf("Test LU...\n");
96: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A);
97: MatLUFactorSymbolic(A,C,row,col,&info);
98: } else {
99: printf("Test ILU...\n");
100: 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: if (MATDSPL){
107: printf("factored matrix:\n");
108: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
109: MatView(A,PETSC_VIEWER_STDOUT_SELF);
110: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
111: MatView(A,PETSC_VIEWER_STDOUT_SELF);
112: MatView(A,viewer2);
113: }
115: /* Solve A*y = b, then check the error */
116: MatSolve(A,b,y);
117: VecAXPY(y,-1.0,x);
118: VecNorm(y,NORM_2,&norm2);
119: MatDestroy(&A);
121: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
122: if (!LU && lf==0){
123: MatDuplicate(C,MAT_COPY_VALUES,&A);
124: MatILUFactor(A,row,col,&info);
125: /*
126: printf("In-place factored matrix:\n");
127: MatView(C,PETSC_VIEWER_STDOUT_SELF);
128: */
129: MatSolve(A,b,y);
130: VecAXPY(y,-1.0,x);
131: VecNorm(y,NORM_2,&norm2_inplace);
132: 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);
133: MatDestroy(&A);
134: }
136: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
137: CHOLESKY = LU;
138: if (CHOLESKY){
139: printf("Test Cholesky...\n");
140: lf = -1;
141: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
142: MatCholeskyFactorSymbolic(A,C,row,&info);
143: } else {
144: printf("Test ICC...\n");
145: info.levels = lf;
146: info.fill = 1.0;
147: info.diagonal_fill = 0;
148: info.zeropivot = 0.0;
149: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
150: MatICCFactorSymbolic(A,C,row,&info);
151: }
152: MatCholeskyFactorNumeric(A,C,&info);
154: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
155: if (lf == -1){
156: PetscOptionsHasName(PETSC_NULL,"-triangular_solve",&TRIANGULAR);
157: if (TRIANGULAR){
158: printf("Test MatForwardSolve...\n");
159: MatForwardSolve(A,b,ytmp);
160: printf("Test MatBackwardSolve...\n");
161: MatBackwardSolve(A,ytmp,y);
162: VecAXPY(y,-1.0,x);
163: VecNorm(y,NORM_2,&norm2);
164: if (norm2 > 1.e-14){
165: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
166: }
167: }
168: }
170: MatSolve(A,b,y);
171: MatDestroy(&A);
172: VecAXPY(y,-1.0,x);
173: VecNorm(y,NORM_2,&norm2);
174: if (lf == -1 && norm2 > 1.e-14){
175: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
176: }
178: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
179: if (!CHOLESKY && lf==0 && !matordering){
180: MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
181: MatICCFactor(A,row,&info);
182: /*
183: printf("In-place factored matrix:\n");
184: MatView(A,PETSC_VIEWER_STDOUT_SELF);
185: */
186: MatSolve(A,b,y);
187: VecAXPY(y,-1.0,x);
188: VecNorm(y,NORM_2,&norm2_inplace);
189: 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);
190: MatDestroy(&A);
191: }
193: /* Free data structures */
194: ISDestroy(&row);
195: ISDestroy(&col);
196: MatDestroy(&C);
197: PetscViewerDestroy(&viewer1);
198: PetscViewerDestroy(&viewer2);
199: PetscRandomDestroy(&rdm);
200: VecDestroy(&x);
201: VecDestroy(&y);
202: VecDestroy(&ytmp);
203: VecDestroy(&b);
204: PetscFinalize();
205: return 0;
206: }