Actual source code: ex30.c
1: 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\
2: Input parameters are:\n\
3: -lf <level> : level of fill for ILU (default is 0)\n\
4: -lu : use full LU or Cholesky factorization\n\
5: -m <value>,-n <value> : grid dimensions\n\
6: Note that most users should employ the KSP interface to the\n\
7: linear solvers instead of using the factorization routines\n\
8: directly.\n\n";
10: #include <petscmat.h>
12: int main(int argc, char **args)
13: {
14: Mat C, A;
15: PetscInt i, j, m = 5, n = 5, Ii, J, lf = 0;
16: PetscBool LU = PETSC_FALSE, CHOLESKY = PETSC_FALSE, TRIANGULAR = PETSC_FALSE, MATDSPL = PETSC_FALSE, flg, matordering, use_mkl_pardiso = PETSC_FALSE;
17: PetscScalar v;
18: IS row, col;
19: PetscViewer viewer1, viewer2;
20: MatFactorInfo info;
21: Vec x, y, b, ytmp;
22: PetscReal norm2, norm2_inplace, tol = 100. * PETSC_MACHINE_EPSILON;
23: PetscRandom rdm;
24: PetscMPIInt size;
25: char pack[PETSC_MAX_PATH_LEN];
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &args, NULL, help));
29: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
30: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
31: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
32: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
33: PetscCall(PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL));
35: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 0, 0, 400, 400, &viewer1));
36: PetscCall(PetscViewerDrawOpen(PETSC_COMM_SELF, 0, 0, 400, 0, 400, 400, &viewer2));
38: PetscCall(MatCreate(PETSC_COMM_SELF, &C));
39: PetscCall(MatSetSizes(C, m * n, m * n, m * n, m * n));
40: PetscCall(MatSetFromOptions(C));
41: PetscCall(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) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
50: J = Ii + n;
51: if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
52: J = Ii - 1;
53: if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
54: J = Ii + 1;
55: if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
56: v = 4.0;
57: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
58: }
59: }
60: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
61: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
63: PetscCall(MatIsSymmetric(C, 0.0, &flg));
64: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric");
66: PetscCall(MatSetOption(C, MAT_SPD, PETSC_TRUE));
68: /* Create vectors for error checking */
69: PetscCall(MatCreateVecs(C, &x, &b));
70: PetscCall(VecDuplicate(x, &y));
71: PetscCall(VecDuplicate(x, &ytmp));
72: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
73: PetscCall(PetscRandomSetFromOptions(rdm));
74: PetscCall(VecSetRandom(x, rdm));
75: PetscCall(MatMult(C, x, b));
77: PetscCall(PetscOptionsHasName(NULL, NULL, "-mat_ordering", &matordering));
78: if (matordering) {
79: PetscCall(MatGetOrdering(C, MATORDERINGRCM, &row, &col));
80: } else {
81: PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col));
82: }
84: PetscCall(PetscOptionsHasName(NULL, NULL, "-display_matrices", &MATDSPL));
85: if (MATDSPL) {
86: printf("original matrix:\n");
87: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO));
88: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF));
89: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF));
90: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_SELF));
91: PetscCall(MatView(C, viewer1));
92: }
94: /* Compute LU or ILU factor A */
95: PetscCall(MatFactorInfoInitialize(&info));
97: info.fill = 1.0;
98: info.diagonal_fill = 0;
99: info.zeropivot = 0.0;
101: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
102: #if defined(PETSC_HAVE_MKL_PARDISO)
103: PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &use_mkl_pardiso));
104: #endif
106: PetscCall(PetscOptionsHasName(NULL, NULL, "-lu", &LU));
107: if (LU) {
108: printf("Test LU...\n");
109: if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &A));
110: else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A));
111: PetscCall(MatLUFactorSymbolic(A, C, row, col, &info));
112: } else {
113: printf("Test ILU...\n");
114: info.levels = lf;
116: PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A));
117: PetscCall(MatILUFactorSymbolic(A, C, row, col, &info));
118: }
119: PetscCall(MatLUFactorNumeric(A, C, &info));
121: /* test MatForwardSolve() and MatBackwardSolve() with MKL Pardiso*/
122: if (LU && use_mkl_pardiso) {
123: PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR));
124: if (TRIANGULAR) {
125: printf("Test MatForwardSolve...\n");
126: PetscCall(MatForwardSolve(A, b, ytmp));
127: printf("Test MatBackwardSolve...\n");
128: PetscCall(MatBackwardSolve(A, ytmp, y));
129: PetscCall(VecAXPY(y, -1.0, x));
130: PetscCall(VecNorm(y, NORM_2, &norm2));
131: if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
132: }
133: }
135: /* Solve A*y = b, then check the error */
136: PetscCall(MatSolve(A, b, y));
137: PetscCall(VecAXPY(y, -1.0, x));
138: PetscCall(VecNorm(y, NORM_2, &norm2));
139: PetscCall(MatDestroy(&A));
141: /* Test in-place ILU(0) and compare it with the out-place ILU(0) */
142: if (!LU && lf == 0) {
143: PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &A));
144: PetscCall(MatILUFactor(A, row, col, &info));
145: /*
146: printf("In-place factored matrix:\n");
147: PetscCall(MatView(C,PETSC_VIEWER_STDOUT_SELF));
148: */
149: PetscCall(MatSolve(A, b, y));
150: PetscCall(VecAXPY(y, -1.0, x));
151: PetscCall(VecNorm(y, NORM_2, &norm2_inplace));
152: PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ILU(0) %g and in-place ILU(0) %g give different residuals", (double)norm2, (double)norm2_inplace);
153: PetscCall(MatDestroy(&A));
154: }
156: /* Test Cholesky and ICC on seqaij matrix with matrix reordering on aij matrix C */
157: PetscCall(PetscOptionsHasName(NULL, NULL, "-chol", &CHOLESKY));
158: if (CHOLESKY) {
159: printf("Test Cholesky...\n");
160: lf = -1;
161: if (use_mkl_pardiso) PetscCall(MatGetFactor(C, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &A));
162: else PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &A));
163: PetscCall(MatCholeskyFactorSymbolic(A, C, row, &info));
164: } else {
165: printf("Test ICC...\n");
166: info.levels = lf;
167: info.fill = 1.0;
168: info.diagonal_fill = 0;
169: info.zeropivot = 0.0;
171: PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A));
172: PetscCall(MatICCFactorSymbolic(A, C, row, &info));
173: }
174: PetscCall(MatCholeskyFactorNumeric(A, C, &info));
176: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
177: if (CHOLESKY) {
178: PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR));
179: if (TRIANGULAR) {
180: printf("Test MatForwardSolve...\n");
181: PetscCall(MatForwardSolve(A, b, ytmp));
182: printf("Test MatBackwardSolve...\n");
183: PetscCall(MatBackwardSolve(A, ytmp, y));
184: PetscCall(VecAXPY(y, -1.0, x));
185: PetscCall(VecNorm(y, NORM_2, &norm2));
186: if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
187: }
188: }
190: PetscCall(MatSolve(A, b, y));
191: PetscCall(MatDestroy(&A));
192: PetscCall(VecAXPY(y, -1.0, x));
193: PetscCall(VecNorm(y, NORM_2, &norm2));
194: if (lf == -1 && norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2));
196: /* Test in-place ICC(0) and compare it with the out-place ICC(0) */
197: if (!CHOLESKY && lf == 0 && !matordering) {
198: PetscCall(MatConvert(C, MATSBAIJ, MAT_INITIAL_MATRIX, &A));
199: PetscCall(MatICCFactor(A, row, &info));
200: /*
201: printf("In-place factored matrix:\n");
202: PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF));
203: */
204: PetscCall(MatSolve(A, b, y));
205: PetscCall(VecAXPY(y, -1.0, x));
206: PetscCall(VecNorm(y, NORM_2, &norm2_inplace));
207: PetscCheck(PetscAbs(norm2 - norm2_inplace) <= tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "ICC(0) %g and in-place ICC(0) %g give different residuals", (double)norm2, (double)norm2_inplace);
208: PetscCall(MatDestroy(&A));
209: }
211: /* Free data structures */
212: PetscCall(ISDestroy(&row));
213: PetscCall(ISDestroy(&col));
214: PetscCall(MatDestroy(&C));
215: PetscCall(PetscViewerDestroy(&viewer1));
216: PetscCall(PetscViewerDestroy(&viewer2));
217: PetscCall(PetscRandomDestroy(&rdm));
218: PetscCall(VecDestroy(&x));
219: PetscCall(VecDestroy(&y));
220: PetscCall(VecDestroy(&ytmp));
221: PetscCall(VecDestroy(&b));
222: PetscCall(PetscFinalize());
223: return 0;
224: }
226: /*TEST
228: test:
229: args: -mat_ordering -display_matrices -nox
230: filter: grep -v " MPI process"
232: test:
233: suffix: 2
234: args: -mat_ordering -display_matrices -nox -lu -chol
236: test:
237: suffix: 3
238: args: -mat_ordering -lu -chol -triangular_solve
240: test:
241: suffix: 4
243: test:
244: suffix: 5
245: args: -lu -chol
247: test:
248: suffix: 6
249: args: -lu -chol -triangular_solve
250: output_file: output/ex30_3.out
252: test:
253: suffix: 7
254: requires: mkl_pardiso
255: args: -lu -chol -mat_solver_type mkl_pardiso
256: output_file: output/ex30_5.out
258: test:
259: suffix: 8
260: requires: mkl_pardiso
261: args: -lu -mat_solver_type mkl_pardiso -triangular_solve
263: TEST*/