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";
9
10 #include <petscmat.h>
11
main(int argc,char ** args)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];
26
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));
34
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));
37
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));
42
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));
62
63 PetscCall(MatIsSymmetric(C, 0.0, &flg));
64 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric");
65
66 PetscCall(MatSetOption(C, MAT_SPD, PETSC_TRUE));
67
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));
76
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 }
83
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 }
93
94 /* Compute LU or ILU factor A */
95 PetscCall(MatFactorInfoInitialize(&info));
96
97 info.fill = 1.0;
98 info.diagonal_fill = 0;
99 info.zeropivot = 0.0;
100
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
105
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;
115
116 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ILU, &A));
117 PetscCall(MatILUFactorSymbolic(A, C, row, col, &info));
118 }
119 PetscCall(MatLUFactorNumeric(A, C, &info));
120
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 }
134
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));
140
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 }
155
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;
170
171 PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_ICC, &A));
172 PetscCall(MatICCFactorSymbolic(A, C, row, &info));
173 }
174 PetscCall(MatCholeskyFactorNumeric(A, C, &info));
175
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 }
189
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));
195
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 }
210
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 }
225
226 /*TEST
227
228 test:
229 args: -mat_ordering -display_matrices -nox
230 filter: grep -v " MPI process"
231
232 test:
233 suffix: 2
234 args: -mat_ordering -display_matrices -nox -lu -chol
235
236 test:
237 suffix: 3
238 args: -mat_ordering -lu -chol -triangular_solve
239
240 test:
241 suffix: 4
242
243 test:
244 suffix: 5
245 args: -lu -chol
246
247 test:
248 suffix: 6
249 args: -lu -chol -triangular_solve
250 output_file: output/ex30_3.out
251
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
257
258 test:
259 suffix: 8
260 requires: mkl_pardiso
261 args: -lu -mat_solver_type mkl_pardiso -triangular_solve
262
263 TEST*/
264