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*/