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