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;
28: PetscInitialize(&argc, &args, (char *)0, help);
29: MPI_Comm_size(PETSC_COMM_WORLD, &size);
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;
47: Ii = j + n * i;
48: J = Ii - n;
49: if (J >= 0) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
50: J = Ii + n;
51: if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
52: J = Ii - 1;
53: if (J >= 0) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
54: J = Ii + 1;
55: if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
56: v = 4.0;
57: MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
58: }
59: }
60: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
61: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
63: MatIsSymmetric(C, 0.0, &flg);
66: /* Create vectors for error checking */
67: MatCreateVecs(C, &x, &b);
68: VecDuplicate(x, &y);
69: VecDuplicate(x, &ytmp);
70: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
71: PetscRandomSetFromOptions(rdm);
72: VecSetRandom(x, rdm);
73: MatMult(C, x, b);
75: PetscOptionsHasName(NULL, NULL, "-mat_ordering", &matordering);
76: if (matordering) {
77: MatGetOrdering(C, MATORDERINGRCM, &row, &col);
78: } else {
79: MatGetOrdering(C, MATORDERINGNATURAL, &row, &col);
80: }
82: PetscOptionsHasName(NULL, NULL, "-display_matrices", &MATDSPL);
83: if (MATDSPL) {
84: printf("original matrix:\n");
85: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO);
86: MatView(C, PETSC_VIEWER_STDOUT_SELF);
87: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF);
88: MatView(C, PETSC_VIEWER_STDOUT_SELF);
89: MatView(C, viewer1);
90: }
92: /* Compute LU or ILU factor A */
93: MatFactorInfoInitialize(&info);
95: info.fill = 1.0;
96: info.diagonal_fill = 0;
97: info.zeropivot = 0.0;
99: PetscOptionsHasName(NULL, NULL, "-lu", &LU);
100: if (LU) {
101: printf("Test LU...\n");
102: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A);
103: MatLUFactorSymbolic(A, C, row, col, &info);
104: } else {
105: printf("Test ILU...\n");
106: info.levels = lf;
108: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A);
109: MatILUFactorSymbolic(A, C, row, col, &info);
110: }
111: MatLUFactorNumeric(A, C, &info);
113: /* Solve A*y = b, then check the error */
114: MatSolve(A, b, y);
115: VecAXPY(y, -1.0, x);
116: VecNorm(y, NORM_2, &norm2);
117: MatDestroy(&A);
119: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
120: if (!LU && lf == 0) {
121: MatDuplicate(C, MAT_COPY_VALUES, &A);
122: MatILUFactor(A, row, col, &info);
123: /*
124: printf("In-place factored matrix:\n");
125: MatView(C,PETSC_VIEWER_STDOUT_SELF);
126: */
127: MatSolve(A, b, y);
128: VecAXPY(y, -1.0, x);
129: VecNorm(y, NORM_2, &norm2_inplace);
131: MatDestroy(&A);
132: }
134: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
135: CHOLESKY = LU;
136: if (CHOLESKY) {
137: printf("Test Cholesky...\n");
138: lf = -1;
139: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &A);
140: MatCholeskyFactorSymbolic(A, C, row, &info);
141: } else {
142: printf("Test ICC...\n");
143: info.levels = lf;
144: info.fill = 1.0;
145: info.diagonal_fill = 0;
146: info.zeropivot = 0.0;
148: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A);
149: MatICCFactorSymbolic(A, C, row, &info);
150: }
151: MatCholeskyFactorNumeric(A, C, &info);
153: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
154: if (lf == -1) {
155: PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR);
156: if (TRIANGULAR) {
157: printf("Test MatForwardSolve...\n");
158: MatForwardSolve(A, b, ytmp);
159: printf("Test MatBackwardSolve...\n");
160: MatBackwardSolve(A, ytmp, y);
161: VecAXPY(y, -1.0, x);
162: VecNorm(y, NORM_2, &norm2);
163: if (norm2 > tol) PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2);
164: }
165: }
167: MatSolve(A, b, y);
168: MatDestroy(&A);
169: VecAXPY(y, -1.0, x);
170: VecNorm(y, NORM_2, &norm2);
171: if (lf == -1 && norm2 > tol) PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2);
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);
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: }
203: /*TEST
205: test:
206: args: -mat_ordering -display_matrices -nox
207: filter: grep -v " MPI process"
209: test:
210: suffix: 2
211: args: -mat_ordering -display_matrices -nox -lu
213: test:
214: suffix: 3
215: args: -mat_ordering -lu -triangular_solve
217: test:
218: suffix: 4
220: test:
221: suffix: 5
222: args: -lu
224: test:
225: suffix: 6
226: args: -lu -triangular_solve
227: output_file: output/ex30_3.out
229: TEST*/