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