Actual source code: ex30.c
petsc-3.7.3 2016-08-01
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, tol = 100.*PETSC_MACHINE_EPSILON;
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(NULL,NULL,"-m",&m,NULL);
34: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
35: PetscOptionsGetInt(NULL,NULL,"-lf",&lf,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: MatCreateVecs(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(NULL,NULL,"-mat_ordering",&matordering);
72: if (matordering) {
73: MatGetOrdering(C,MATORDERINGRCM,&row,&col);
74: } else {
75: MatGetOrdering(C,MATORDERINGNATURAL,&row,&col);
76: }
78: PetscOptionsHasName(NULL,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);
91: info.fill = 1.0;
92: info.diagonal_fill = 0;
93: info.zeropivot = 0.0;
95: PetscOptionsHasName(NULL,NULL,"-lu",&LU);
96: if (LU) {
97: printf("Test LU...\n");
98: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A);
99: MatLUFactorSymbolic(A,C,row,col,&info);
100: } else {
101: printf("Test ILU...\n");
102: info.levels = lf;
104: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ILU,&A);
105: MatILUFactorSymbolic(A,C,row,col,&info);
106: }
107: MatLUFactorNumeric(A,C,&info);
109: /* Solve A*y = b, then check the error */
110: MatSolve(A,b,y);
111: VecAXPY(y,-1.0,x);
112: VecNorm(y,NORM_2,&norm2);
113: MatDestroy(&A);
115: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
116: if (!LU && lf==0) {
117: MatDuplicate(C,MAT_COPY_VALUES,&A);
118: MatILUFactor(A,row,col,&info);
119: /*
120: printf("In-place factored matrix:\n");
121: MatView(C,PETSC_VIEWER_STDOUT_SELF);
122: */
123: MatSolve(A,b,y);
124: VecAXPY(y,-1.0,x);
125: VecNorm(y,NORM_2,&norm2_inplace);
126: 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);
127: MatDestroy(&A);
128: }
130: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
131: CHOLESKY = LU;
132: if (CHOLESKY) {
133: printf("Test Cholesky...\n");
134: lf = -1;
135: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
136: MatCholeskyFactorSymbolic(A,C,row,&info);
137: } else {
138: printf("Test ICC...\n");
139: info.levels = lf;
140: info.fill = 1.0;
141: info.diagonal_fill = 0;
142: info.zeropivot = 0.0;
144: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
145: MatICCFactorSymbolic(A,C,row,&info);
146: }
147: MatCholeskyFactorNumeric(A,C,&info);
149: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
150: if (lf == -1) {
151: PetscOptionsHasName(NULL,NULL,"-triangular_solve",&TRIANGULAR);
152: if (TRIANGULAR) {
153: printf("Test MatForwardSolve...\n");
154: MatForwardSolve(A,b,ytmp);
155: printf("Test MatBackwardSolve...\n");
156: MatBackwardSolve(A,ytmp,y);
157: VecAXPY(y,-1.0,x);
158: VecNorm(y,NORM_2,&norm2);
159: if (norm2 > tol) {
160: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g\n",(double)norm2);
161: }
162: }
163: }
165: MatSolve(A,b,y);
166: MatDestroy(&A);
167: VecAXPY(y,-1.0,x);
168: VecNorm(y,NORM_2,&norm2);
169: if (lf == -1 && norm2 > tol) {
170: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
171: }
173: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
174: if (!CHOLESKY && lf==0 && !matordering) {
175: MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
176: MatICCFactor(A,row,&info);
177: /*
178: printf("In-place factored matrix:\n");
179: MatView(A,PETSC_VIEWER_STDOUT_SELF);
180: */
181: MatSolve(A,b,y);
182: VecAXPY(y,-1.0,x);
183: VecNorm(y,NORM_2,&norm2_inplace);
184: 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);
185: MatDestroy(&A);
186: }
188: /* Free data structures */
189: ISDestroy(&row);
190: ISDestroy(&col);
191: MatDestroy(&C);
192: PetscViewerDestroy(&viewer1);
193: PetscViewerDestroy(&viewer2);
194: PetscRandomDestroy(&rdm);
195: VecDestroy(&x);
196: VecDestroy(&y);
197: VecDestroy(&ytmp);
198: VecDestroy(&b);
199: PetscFinalize();
200: return 0;
201: }