Actual source code: ex28.c
1: static char help[] = "Illustrate how to do one symbolic factorization and multiple numeric factorizations using same matrix 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];
20: PetscInitialize(&argc, &args, (char *)0, help);
21: 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: MatCreate(PETSC_COMM_WORLD, &A[k]);
26: MatSetSizes(A[k], PETSC_DECIDE, PETSC_DECIDE, N, N);
27: MatSetFromOptions(A[k]);
28: MatSetUp(A[k]);
29: 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: MatSetValues(A[k], 1, &i, 2, col + 1, value + 1, INSERT_VALUES);
40: } else if (i == N - 1) {
41: MatSetValues(A[k], 1, &i, 2, col, value, INSERT_VALUES);
42: } else {
43: MatSetValues(A[k], 1, &i, 3, col, value, INSERT_VALUES);
44: }
45: }
46: MatAssemblyBegin(A[k], MAT_FINAL_ASSEMBLY);
47: MatAssemblyEnd(A[k], MAT_FINAL_ASSEMBLY);
48: MatSetOption(A[k], MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
49: }
51: /* Create vectors */
52: MatCreateVecs(A[0], &x, &b);
53: VecDuplicate(x, &u);
55: /* Set rhs vector b */
56: VecSet(b, 1.0);
58: /* Get a symbolic factor F from A[0] */
59: PetscStrncpy(solvertype, "petsc", sizeof(solvertype));
60: PetscOptionsGetString(NULL, NULL, "-mat_solver_type", solvertype, sizeof(solvertype), NULL);
61: PetscOptionsGetEnum(NULL, NULL, "-mat_factor_type", MatFactorTypes, (PetscEnum *)&facttype, NULL);
63: MatGetFactor(A[0], solvertype, facttype, &F);
64: /* test mumps options */
65: #if defined(PETSC_HAVE_MUMPS)
66: MatMumpsSetIcntl(F, 7, 5);
67: #endif
68: PetscStrncpy(factortype, MatFactorTypes[facttype], sizeof(factortype));
69: PetscStrtoupper(solvertype);
70: PetscStrtoupper(factortype);
71: PetscPrintf(PETSC_COMM_WORLD, " %s %s:\n", solvertype, factortype);
73: MatFactorInfoInitialize(&info);
74: info.fill = 5.0;
75: MatGetOrdering(A[0], MATORDERINGNATURAL, &perm, &iperm);
76: switch (facttype) {
77: case MAT_FACTOR_LU:
78: MatLUFactorSymbolic(F, A[0], perm, iperm, &info);
79: break;
80: case MAT_FACTOR_ILU:
81: MatILUFactorSymbolic(F, A[0], perm, iperm, &info);
82: break;
83: case MAT_FACTOR_ICC:
84: MatICCFactorSymbolic(F, A[0], perm, &info);
85: break;
86: case MAT_FACTOR_CHOLESKY:
87: 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: MatLUFactorNumeric(F, A[k], &info);
99: break;
100: case MAT_FACTOR_ICC:
101: case MAT_FACTOR_CHOLESKY:
102: 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: MatSolve(F, b, x);
111: /* Check the residual */
112: MatMult(A[k], x, u);
113: VecAXPY(u, -1.0, b);
114: VecNorm(u, NORM_INFINITY, &norm);
115: if (norm > tol) 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++) MatDestroy(&A[k]);
120: MatDestroy(&F);
121: ISDestroy(&perm);
122: ISDestroy(&iperm);
123: VecDestroy(&x);
124: VecDestroy(&b);
125: VecDestroy(&u);
126: 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*/