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;
27: PetscInitialize(&argc, &args, (char *)0, help);
28: MPI_Comm_size(PETSC_COMM_WORLD, &size);
30: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
31: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
32: PetscOptionsGetInt(NULL, NULL, "-lf", &lf, NULL);
34: MatCreate(PETSC_COMM_SELF, &C);
35: MatSetSizes(C, m * n, m * n, m * n, m * n);
36: MatSetFromOptions(C);
37: 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) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
46: J = Ii + n;
47: if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
48: J = Ii - 1;
49: if (J >= 0) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
50: J = Ii + 1;
51: if (J < m * n) MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
52: v = 4.0;
53: MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
54: }
55: }
56: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
59: MatIsSymmetric(C, 0.0, &flg);
61: MatConvert(C, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sC);
63: /* Create vectors for error checking */
64: MatCreateVecs(C, &x, &b);
65: VecDuplicate(x, &y);
66: VecDuplicate(x, &ytmp);
67: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
68: PetscRandomSetFromOptions(rdm);
69: VecSetRandom(x, rdm);
70: MatMult(C, x, b);
72: MatGetOrdering(C, MATORDERINGNATURAL, &row, &col);
74: /* Compute CHOLESKY or ICC factor sA */
75: MatFactorInfoInitialize(&info);
77: info.fill = 1.0;
78: info.diagonal_fill = 0;
79: info.zeropivot = 0.0;
81: PetscOptionsHasName(NULL, NULL, "-cholesky", &CHOLESKY);
82: if (CHOLESKY) {
83: PetscPrintf(PETSC_COMM_SELF, "Test CHOLESKY...\n");
84: MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &sA);
85: MatCholeskyFactorSymbolic(sA, sC, row, &info);
86: } else {
87: PetscPrintf(PETSC_COMM_SELF, "Test ICC...\n");
88: info.levels = lf;
90: MatGetFactor(sC, MATSOLVERPETSC, MAT_FACTOR_ICC, &sA);
91: MatICCFactorSymbolic(sA, sC, row, &info);
92: }
93: MatCholeskyFactorNumeric(sA, sC, &info);
95: /* test MatForwardSolve() and MatBackwardSolve() with matrix reordering on aij matrix C */
96: if (CHOLESKY) {
97: PetscOptionsHasName(NULL, NULL, "-triangular_solve", &TRIANGULAR);
98: if (TRIANGULAR) {
99: PetscPrintf(PETSC_COMM_SELF, "Test MatForwardSolve...\n");
100: MatForwardSolve(sA, b, ytmp);
101: PetscPrintf(PETSC_COMM_SELF, "Test MatBackwardSolve...\n");
102: MatBackwardSolve(sA, ytmp, y);
103: VecAXPY(y, -1.0, x);
104: VecNorm(y, NORM_2, &norm2);
105: if (norm2 > tol) PetscPrintf(PETSC_COMM_SELF, "MatForwardSolve and BackwardSolve: Norm of error=%g\n", (double)norm2);
106: }
107: }
109: MatSolve(sA, b, y);
110: MatDestroy(&sC);
111: MatDestroy(&sA);
112: VecAXPY(y, -1.0, x);
113: VecNorm(y, NORM_2, &norm2);
114: if (lf == -1 && norm2 > tol) PetscPrintf(PETSC_COMM_SELF, " reordered SEQAIJ: Cholesky/ICC levels %" PetscInt_FMT ", residual %g\n", lf, (double)norm2);
116: /* Free data structures */
117: MatDestroy(&C);
118: ISDestroy(&row);
119: ISDestroy(&col);
120: PetscRandomDestroy(&rdm);
121: VecDestroy(&x);
122: VecDestroy(&y);
123: VecDestroy(&ytmp);
124: VecDestroy(&b);
125: 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*/