Actual source code: ex30.c
petsc-3.4.5 2014-06-29
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(NULL,"-m",&m,NULL);
34: PetscOptionsGetInt(NULL,"-n",&n,NULL);
35: PetscOptionsGetInt(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: 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(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,"-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,"-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: if (MATDSPL) {
110: printf("factored matrix:\n");
111: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF,PETSC_VIEWER_ASCII_INFO);
112: MatView(A,PETSC_VIEWER_STDOUT_SELF);
113: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
114: MatView(A,PETSC_VIEWER_STDOUT_SELF);
115: MatView(A,viewer2);
116: }
118: /* Solve A*y = b, then check the error */
119: MatSolve(A,b,y);
120: VecAXPY(y,-1.0,x);
121: VecNorm(y,NORM_2,&norm2);
122: MatDestroy(&A);
124: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
125: if (!LU && lf==0) {
126: MatDuplicate(C,MAT_COPY_VALUES,&A);
127: MatILUFactor(A,row,col,&info);
128: /*
129: printf("In-place factored matrix:\n");
130: MatView(C,PETSC_VIEWER_STDOUT_SELF);
131: */
132: MatSolve(A,b,y);
133: VecAXPY(y,-1.0,x);
134: VecNorm(y,NORM_2,&norm2_inplace);
135: 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);
136: MatDestroy(&A);
137: }
139: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
140: CHOLESKY = LU;
141: if (CHOLESKY) {
142: printf("Test Cholesky...\n");
143: lf = -1;
144: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&A);
145: MatCholeskyFactorSymbolic(A,C,row,&info);
146: } else {
147: printf("Test ICC...\n");
148: info.levels = lf;
149: info.fill = 1.0;
150: info.diagonal_fill = 0;
151: info.zeropivot = 0.0;
153: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_ICC,&A);
154: MatICCFactorSymbolic(A,C,row,&info);
155: }
156: MatCholeskyFactorNumeric(A,C,&info);
158: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
159: if (lf == -1) {
160: PetscOptionsHasName(NULL,"-triangular_solve",&TRIANGULAR);
161: if (TRIANGULAR) {
162: printf("Test MatForwardSolve...\n");
163: MatForwardSolve(A,b,ytmp);
164: printf("Test MatBackwardSolve...\n");
165: MatBackwardSolve(A,ytmp,y);
166: VecAXPY(y,-1.0,x);
167: VecNorm(y,NORM_2,&norm2);
168: if (norm2 > 1.e-14) {
169: PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%G\n",norm2);
170: }
171: }
172: }
174: MatSolve(A,b,y);
175: MatDestroy(&A);
176: VecAXPY(y,-1.0,x);
177: VecNorm(y,NORM_2,&norm2);
178: if (lf == -1 && norm2 > 1.e-14) {
179: PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %d, residual %g\n",lf,norm2);
180: }
182: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
183: if (!CHOLESKY && lf==0 && !matordering) {
184: MatConvert(C,MATSBAIJ,MAT_INITIAL_MATRIX,&A);
185: MatICCFactor(A,row,&info);
186: /*
187: printf("In-place factored matrix:\n");
188: MatView(A,PETSC_VIEWER_STDOUT_SELF);
189: */
190: MatSolve(A,b,y);
191: VecAXPY(y,-1.0,x);
192: VecNorm(y,NORM_2,&norm2_inplace);
193: 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);
194: MatDestroy(&A);
195: }
197: /* Free data structures */
198: ISDestroy(&row);
199: ISDestroy(&col);
200: MatDestroy(&C);
201: PetscViewerDestroy(&viewer1);
202: PetscViewerDestroy(&viewer2);
203: PetscRandomDestroy(&rdm);
204: VecDestroy(&x);
205: VecDestroy(&y);
206: VecDestroy(&ytmp);
207: VecDestroy(&b);
208: PetscFinalize();
209: return 0;
210: }