Actual source code: ex128.c
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";
10: #include <petscmat.h>
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;
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));
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));
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));
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));
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));
71: PetscCall(MatGetOrdering(C, MATORDERINGNATURAL, &row, &col));
73: /* Compute CHOLESKY or ICC factor sA */
74: PetscCall(MatFactorInfoInitialize(&info));
76: info.fill = 1.0;
77: info.diagonal_fill = 0;
78: info.zeropivot = 0.0;
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;
89: PetscCall(MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_ICC, &sA));
90: PetscCall(MatICCFactorSymbolic(sA, sC, row, &info));
91: }
92: PetscCall(MatCholeskyFactorNumeric(sA, sC, &info));
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: }
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));
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: }
128: /*TEST
130: test:
131: output_file: output/ex128.out
133: test:
134: suffix: 2
135: args: -cholesky -triangular_solve
137: TEST*/