Actual source code: ex28.c
1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix nonzero structure. \n\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: PetscInt i, rstart, rend, N = 10, num_numfac = 5, col[3], k;
8: Mat A[5], F;
9: Vec u, x, b;
10: PetscMPIInt rank;
11: PetscScalar value[3];
12: PetscReal norm, tol = 100 * PETSC_MACHINE_EPSILON;
13: IS perm, iperm;
14: MatFactorInfo info;
15: MatFactorType facttype = MAT_FACTOR_LU;
16: char solvertype[64];
17: char factortype[64];
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, NULL, help));
21: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
23: /* Create and assemble matrices, all have same data structure */
24: for (k = 0; k < num_numfac; k++) {
25: PetscCall(MatCreate(PETSC_COMM_WORLD, &A[k]));
26: PetscCall(MatSetSizes(A[k], PETSC_DECIDE, PETSC_DECIDE, N, N));
27: PetscCall(MatSetFromOptions(A[k]));
28: PetscCall(MatSetUp(A[k]));
29: PetscCall(MatGetOwnershipRange(A[k], &rstart, &rend));
31: value[0] = -1.0 * (k + 1);
32: value[1] = 2.0 * (k + 1);
33: value[2] = -1.0 * (k + 1);
34: for (i = rstart; i < rend; i++) {
35: col[0] = i - 1;
36: col[1] = i;
37: col[2] = i + 1;
38: if (i == 0) {
39: PetscCall(MatSetValues(A[k], 1, &i, 2, col + 1, value + 1, INSERT_VALUES));
40: } else if (i == N - 1) {
41: PetscCall(MatSetValues(A[k], 1, &i, 2, col, value, INSERT_VALUES));
42: } else {
43: PetscCall(MatSetValues(A[k], 1, &i, 3, col, value, INSERT_VALUES));
44: }
45: }
46: PetscCall(MatAssemblyBegin(A[k], MAT_FINAL_ASSEMBLY));
47: PetscCall(MatAssemblyEnd(A[k], MAT_FINAL_ASSEMBLY));
48: PetscCall(MatSetOption(A[k], MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
49: }
51: /* Create vectors */
52: PetscCall(MatCreateVecs(A[0], &x, &b));
53: PetscCall(VecDuplicate(x, &u));
55: /* Set rhs vector b */
56: PetscCall(VecSet(b, 1.0));
58: /* Get a symbolic factor F from A[0] */
59: PetscCall(PetscStrncpy(solvertype, "petsc", sizeof(solvertype)));
60: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solvertype, sizeof(solvertype), NULL));
61: PetscCall(PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&facttype, NULL));
63: PetscCall(MatGetFactor(A[0], solvertype, facttype, &F));
64: /* test mumps options */
65: #if defined(PETSC_HAVE_MUMPS)
66: PetscCall(MatMumpsSetIcntl(F, 7, 5));
67: #endif
68: PetscCall(PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype)));
69: PetscCall(PetscStrtoupper(solvertype));
70: PetscCall(PetscStrtoupper(factortype));
71: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype));
73: PetscCall(MatFactorInfoInitialize(&info));
74: info.fill = 5.0;
75: PetscCall(MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm));
76: switch (facttype) {
77: case MAT_FACTOR_LU:
78: PetscCall(MatLUFactorSymbolic(F, A[0], perm, iperm, &info));
79: break;
80: case MAT_FACTOR_ILU:
81: PetscCall(MatILUFactorSymbolic(F, A[0], perm, iperm, &info));
82: break;
83: case MAT_FACTOR_ICC:
84: PetscCall(MatICCFactorSymbolic(F, A[0], perm, &info));
85: break;
86: case MAT_FACTOR_CHOLESKY:
87: PetscCall(MatCholeskyFactorSymbolic(F, A[0], perm, &info));
88: break;
89: default:
90: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype);
91: }
93: /* Compute numeric factors using same F, then solve */
94: for (k = 0; k < num_numfac; k++) {
95: switch (facttype) {
96: case MAT_FACTOR_LU:
97: case MAT_FACTOR_ILU:
98: PetscCall(MatLUFactorNumeric(F, A[k], &info));
99: break;
100: case MAT_FACTOR_ICC:
101: case MAT_FACTOR_CHOLESKY:
102: PetscCall(MatCholeskyFactorNumeric(F, A[k], &info));
103: break;
104: default:
105: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not for factor type %s", factortype);
106: }
108: /* Solve A[k] * x = b */
109: PetscCall(MatSolve(F, b, x));
111: /* Check the residual */
112: PetscCall(MatMult(A[k], x, u));
113: PetscCall(VecAXPY(u, -1.0, b));
114: PetscCall(VecNorm(u, NORM_INFINITY, &norm));
115: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the %s numfact and solve: residual %g\n", k, factortype, (double)norm));
116: }
118: /* Free data structures */
119: for (k = 0; k < num_numfac; k++) PetscCall(MatDestroy(&A[k]));
120: PetscCall(MatDestroy(&F));
121: PetscCall(ISDestroy(&perm));
122: PetscCall(ISDestroy(&iperm));
123: PetscCall(VecDestroy(&x));
124: PetscCall(VecDestroy(&b));
125: PetscCall(VecDestroy(&u));
126: PetscCall(PetscFinalize());
127: return 0;
128: }
130: /*TEST
132: test:
134: test:
135: suffix: 2
136: args: -mat_solver_type superlu
137: requires: superlu
139: test:
140: suffix: 3
141: nsize: 2
142: requires: mumps
143: args: -mat_solver_type mumps
145: test:
146: suffix: 4
147: args: -mat_solver_type cusparse -mat_type aijcusparse -mat_factor_type {{lu cholesky ilu icc}separate output}
148: requires: cuda
150: TEST*/