1 static char help[] = "Tests ILU and ICC factorization with and without matrix ordering on seqsbaij format. Modified from ex30.c\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, sC, sA;
15 PetscInt i, j, m = 5, n = 5, Ii, J, lf = 0;
16 PetscBool CHOLESKY = PETSC_FALSE, TRIANGULAR = PETSC_FALSE, flg;
17 PetscScalar v;
18 IS row, col;
19 MatFactorInfo info;
20 Vec x, y, b, ytmp;
21 PetscReal norm2, tol = 100 * PETSC_MACHINE_EPSILON;
22 PetscRandom rdm;
23 PetscMPIInt size;
24
25 PetscFunctionBeginUser;
26 PetscCall(PetscInitialize(&argc, &args, NULL, help));
27 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
28 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
29 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
30 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
31 PetscCall(PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL));
32
33 PetscCall(MatCreate(PETSC_COMM_SELF, &C));
34 PetscCall(MatSetSizes(C, m * n, m * n, m * n, m * n));
35 PetscCall(MatSetFromOptions(C));
36 PetscCall(MatSetUp(C));
37
38 /* Create matrix C in seqaij format and sC in seqsbaij. (This is five-point stencil with some extra elements) */
39 for (i = 0; i < m; i++) {
40 for (j = 0; j < n; j++) {
41 v = -1.0;
42 Ii = j + n * i;
43 J = Ii - n;
44 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
45 J = Ii + n;
46 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
47 J = Ii - 1;
48 if (J >= 0) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
49 J = Ii + 1;
50 if (J < m * n) PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
51 v = 4.0;
52 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
53 }
54 }
55 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
56 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
57
58 PetscCall(MatIsSymmetric(C, 0.0, &flg));
59 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "C is non-symmetric");
60 PetscCall(MatConvert(C, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sC));
61
62 /* Create vectors for error checking */
63 PetscCall(MatCreateVecs(C, &x, &b));
64 PetscCall(VecDuplicate(x, &y));
65 PetscCall(VecDuplicate(x, &ytmp));
66 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
67 PetscCall(PetscRandomSetFromOptions(rdm));
68 PetscCall(VecSetRandom(x, rdm));
69 PetscCall(MatMult(C, x, b));
70
71 PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col));
72
73 /* Compute CHOLESKY or ICC factor sA */
74 PetscCall(MatFactorInfoInitialize(&info));
75
76 info.fill = 1.0;
77 info.diagonal_fill = 0;
78 info.zeropivot = 0.0;
79
80 PetscCall(PetscOptionsHasName(NULL, NULL, "-cholesky", &CHOLESKY));
81 if (CHOLESKY) {
82 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test CHOLESKY...\n"));
83 PetscCall(MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sA));
84 PetscCall(MatCholeskyFactorSymbolic(sA, sC, row, &info));
85 } else {
86 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test ICC...\n"));
87 info.levels = lf;
88
89 PetscCall(MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_ICC, &sA));
90 PetscCall(MatICCFactorSymbolic(sA, sC, row, &info));
91 }
92 PetscCall(MatCholeskyFactorNumeric(sA, sC, &info));
93
94 /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
95 if (CHOLESKY) {
96 PetscCall(PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR));
97 if (TRIANGULAR) {
98 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test MatForwardSolve...\n"));
99 PetscCall(MatForwardSolve(sA, b, ytmp));
100 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Test MatBackwardSolve...\n"));
101 PetscCall(MatBackwardSolve(sA, ytmp, y));
102 PetscCall(VecAXPY(y, -1.0, x));
103 PetscCall(VecNorm(y, NORM_2, &norm2));
104 if (norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2));
105 }
106 }
107
108 PetscCall(MatSolve(sA, b, y));
109 PetscCall(MatDestroy(&sC));
110 PetscCall(MatDestroy(&sA));
111 PetscCall(VecAXPY(y, -1.0, x));
112 PetscCall(VecNorm(y, NORM_2, &norm2));
113 if (lf == -1 && norm2 > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2));
114
115 /* Free data structures */
116 PetscCall(MatDestroy(&C));
117 PetscCall(ISDestroy(&row));
118 PetscCall(ISDestroy(&col));
119 PetscCall(PetscRandomDestroy(&rdm));
120 PetscCall(VecDestroy(&x));
121 PetscCall(VecDestroy(&y));
122 PetscCall(VecDestroy(&ytmp));
123 PetscCall(VecDestroy(&b));
124 PetscCall(PetscFinalize());
125 return 0;
126 }
127
128 /*TEST
129
130 test:
131 output_file: output/ex128.out
132
133 test:
134 suffix: 2
135 args: -cholesky -triangular_solve
136
137 TEST*/
138