Actual source code: ex30.c
petsc-3.10.5 2019-03-28
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;
18: PetscBool LU=PETSC_FALSE,CHOLESKY,TRIANGULAR=PETSC_FALSE,MATDSPL=PETSC_FALSE,flg,matordering;
19: PetscScalar v;
20: IS row,col;
21: PetscViewer viewer1,viewer2;
22: MatFactorInfo info;
23: Vec x,y,b,ytmp;
24: PetscReal norm2,norm2_inplace, tol = 100.*PETSC_MACHINE_EPSILON;
25: PetscRandom rdm;
26: PetscMPIInt size;
28: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
29: MPI_Comm_size(PETSC_COMM_WORLD,&size);
30: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
31: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-lf",&lf,NULL);
35: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,0,0,400,400,&viewer1);
36: PetscViewerDrawOpen(PETSC_COMM_SELF,0,0,400,0,400,400,&viewer2);
38: MatCreate(PETSC_COMM_SELF,&C);
39: MatSetSizes(C,m*n,m*n,m*n,m*n);
40: MatSetFromOptions(C);
41: MatSetUp(C);
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: MatIsSymmetric(C,0.0,&flg);
58: if (!flg) SETERRQ(PETSC_COMM_SELF,1,"C is non-symmetric");
60: /* Create vectors for error checking */
61: MatCreateVecs(C,&x,&b);
62: VecDuplicate(x,&y);
63: VecDuplicate(x,&ytmp);
64: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
65: PetscRandomSetFromOptions(rdm);
66: VecSetRandom(x,rdm);
67: MatMult(C,x,b);
69: PetscOptionsHasName(NULL,NULL,"-mat_ordering",&matordering);
70: if (matordering) {
71: MatGetOrdering(C,MATORDERINGRCM,&row,&col);
72: } else {
73: MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);
74: }
76: PetscOptionsHasName(NULL,NULL,"-display_matrices",&MATDSPL);
77: if (MATDSPL) {
78: printf("original matrix:\n");
79: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
80: MatView(C,PETSC_VIEWER_STDOUT_SELF);
81: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
82: MatView(C,PETSC_VIEWER_STDOUT_SELF);
83: MatView(C,viewer1);
84: }
86: /* Compute LU or ILU factor A */
87: MatFactorInfoInitialize(&info);
89: info.fill = 1.0;
90: info.diagonal_fill = 0;
91: info.zeropivot = 0.0;
93: PetscOptionsHasName(NULL,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;
102: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ILU,&A);
103: MatILUFactorSymbolic(A,C,row,col,&info);
104: }
105: MatLUFactorNumeric(A,C,&info);
107: /* Solve A*y = b, then check the error */
108: MatSolve(A,b,y);
109: VecAXPY(y,-1.0,x);
110: VecNorm(y,NORM_2,&norm2);
111: MatDestroy(&A);
113: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
114: if (!LU && lf==0) {
115: MatDuplicate(C,MAT_COPY_VALUES,&A);
116: MatILUFactor(A,row,col,&info);
117: /*
118: printf("In-place factored matrix:\n");
119: MatView(C,PETSC_VIEWER_STDOUT_SELF);
120: */
121: MatSolve(A,b,y);
122: VecAXPY(y,-1.0,x);
123: VecNorm(y,NORM_2,&norm2_inplace);
124: if (PetscAbs(norm2 - norm2_inplace) > tol) SETERRQ2(PETSC_COMM_SELF,1,"ILU(0) %g and in-place ILU(0) %g give different residuals",(double)norm2,(double)norm2_inplace);
125: MatDestroy(&A);
126: }
128: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
129: CHOLESKY = LU;
130: if (CHOLESKY) {
131: printf("Test Cholesky...\n");
132: lf = -1;
133: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
134: MatCholeskyFactorSymbolic(A,C,row,&info);
135: } else {
136: printf("Test ICC...\n");
137: info.levels = lf;
138: info.fill = 1.0;
139: info.diagonal_fill = 0;
140: info.zeropivot = 0.0;
142: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
143: MatICCFactorSymbolic(A,C,row,&info);
144: }
145: MatCholeskyFactorNumeric(A,C,&info);
147: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
148: if (lf == -1) {
149: PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR);
150: if (TRIANGULAR) {
151: printf("Test MatForwardSolve...\n");
152: MatForwardSolve(A,b,ytmp);
153: printf("Test MatBackwardSolve...\n");
154: MatBackwardSolve(A,ytmp,y);
155: VecAXPY(y,-1.0,x);
156: VecNorm(y,NORM_2,&norm2);
157: if (norm2 > tol) {
158: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2);
159: }
160: }
161: }
163: MatSolve(A,b,y);
164: MatDestroy(&A);
165: VecAXPY(y,-1.0,x);
166: VecNorm(y,NORM_2,&norm2);
167: if (lf == -1 && norm2 > tol) {
168: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
169: }
171: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
172: if (!CHOLESKY && lf==0 && !matordering) {
173: MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
174: MatICCFactor(A,row,&info);
175: /*
176: printf("In-place factored matrix:\n");
177: MatView(A,PETSC_VIEWER_STDOUT_SELF);
178: */
179: MatSolve(A,b,y);
180: VecAXPY(y,-1.0,x);
181: VecNorm(y,NORM_2,&norm2_inplace);
182: if (PetscAbs(norm2 - norm2_inplace) > tol) SETERRQ2(PETSC_COMM_SELF,1,"ICC(0) %g and in-place ICC(0) %g give different residuals",(double)norm2,(double)norm2_inplace);
183: MatDestroy(&A);
184: }
186: /* Free data structures */
187: ISDestroy(&row);
188: ISDestroy(&col);
189: MatDestroy(&C);
190: PetscViewerDestroy(&viewer1);
191: PetscViewerDestroy(&viewer2);
192: PetscRandomDestroy(&rdm);
193: VecDestroy(&x);
194: VecDestroy(&y);
195: VecDestroy(&ytmp);
196: VecDestroy(&b);
197: PetscFinalize();
198: return ierr;
199: }
202: /*TEST
204: test:
205: args: -mat_ordering -display_matrices -nox
206: filter: grep -v "MPI processes"
208: test:
209: suffix: 2
210: args: -mat_ordering -display_matrices -nox -lu
212: test:
213: suffix: 3
214: args: -mat_ordering -lu -triangular_solve
216: test:
217: suffix: 4
219: test:
220: suffix: 5
221: args: -lu
223: test:
224: suffix: 6
225: args: -lu -triangular_solve
226: output_file: output/ex30_3.out
228: TEST*/