1 /*
2 Provides an interface to the KLUv1.2 sparse solver
3
4 When build with PETSC_USE_64BIT_INDICES this will use SuiteSparse_long as the
5 integer type in KLU, otherwise it will use int. This means
6 all integers in this file are simply declared as PetscInt. Also it means
7 that KLU SuiteSparse_long version MUST be built with 64-bit integers when used.
8
9 */
10 #include <../src/mat/impls/aij/seq/aij.h>
11
12 #if defined(PETSC_USE_64BIT_INDICES)
13 #define klu_K_defaults klu_l_defaults
14 #define klu_K_analyze(a, b, c, d) klu_l_analyze((SuiteSparse_long)a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, d)
15 #define klu_K_analyze_given(a, b, c, d, e, f) klu_l_analyze_given((SuiteSparse_long)a, (SuiteSparse_long *)b, (SuiteSparse_long *)c, (SuiteSparse_long *)d, (SuiteSparse_long *)e, f)
16 #define klu_K_free_symbolic klu_l_free_symbolic
17 #define klu_K_free_numeric klu_l_free_numeric
18 #define klu_K_common klu_l_common
19 #define klu_K_symbolic klu_l_symbolic
20 #define klu_K_numeric klu_l_numeric
21 #if defined(PETSC_USE_COMPLEX)
22 #define klu_K_factor(a, b, c, d, e) klu_zl_factor((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e);
23 #define klu_K_solve klu_zl_solve
24 #define klu_K_tsolve(a, b, c, d, e, f, g) klu_zl_tsolve((a), (b), (c), (d), (PetscReal *)(e), (f), (g))
25 #define klu_K_refactor klu_zl_refactor
26 #define klu_K_sort klu_zl_sort
27 #define klu_K_flops klu_zl_flops
28 #define klu_K_rgrowth klu_zl_rgrowth
29 #define klu_K_condest klu_zl_condest
30 #define klu_K_rcond klu_zl_rcond
31 #define klu_K_scale klu_zl_scale
32 #else
33 #define klu_K_factor(a, b, c, d, e) klu_l_factor((SuiteSparse_long *)a, (SuiteSparse_long *)b, c, d, e);
34 #define klu_K_solve klu_l_solve
35 #define klu_K_tsolve(a, b, c, d, e, f, g) klu_l_tsolve((a), (b), (c), (d), (e), (g))
36 #define klu_K_refactor klu_l_refactor
37 #define klu_K_sort klu_l_sort
38 #define klu_K_flops klu_l_flops
39 #define klu_K_rgrowth klu_l_rgrowth
40 #define klu_K_condest klu_l_condest
41 #define klu_K_rcond klu_l_rcond
42 #define klu_K_scale klu_l_scale
43 #endif
44 #else
45 #define klu_K_defaults klu_defaults
46 #define klu_K_analyze klu_analyze
47 #define klu_K_analyze_given klu_analyze_given
48 #define klu_K_free_symbolic klu_free_symbolic
49 #define klu_K_free_numeric klu_free_numeric
50 #define klu_K_common klu_common
51 #define klu_K_symbolic klu_symbolic
52 #define klu_K_numeric klu_numeric
53 #if defined(PETSC_USE_COMPLEX)
54 #define klu_K_factor klu_z_factor
55 #define klu_K_solve klu_z_solve
56 #define klu_K_tsolve(a, b, c, d, e, f, g) klu_z_tsolve((a), (b), (c), (d), (PetscReal *)(e), (f), (g))
57 #define klu_K_refactor klu_z_refactor
58 #define klu_K_sort klu_z_sort
59 #define klu_K_flops klu_z_flops
60 #define klu_K_rgrowth klu_z_rgrowth
61 #define klu_K_condest klu_z_condest
62 #define klu_K_rcond klu_z_rcond
63 #define klu_K_scale klu_z_scale
64 #else
65 #define klu_K_factor klu_factor
66 #define klu_K_solve klu_solve
67 #define klu_K_tsolve(a, b, c, d, e, f, g) klu_tsolve((a), (b), (c), (d), (e), (g))
68 #define klu_K_refactor klu_refactor
69 #define klu_K_sort klu_sort
70 #define klu_K_flops klu_flops
71 #define klu_K_rgrowth klu_rgrowth
72 #define klu_K_condest klu_condest
73 #define klu_K_rcond klu_rcond
74 #define klu_K_scale klu_scale
75 #endif
76 #endif
77
78 EXTERN_C_BEGIN
79 #include <klu.h>
80 EXTERN_C_END
81
82 static const char *KluOrderingTypes[] = {"AMD", "COLAMD"};
83 static const char *scale[] = {"NONE", "SUM", "MAX"};
84
85 typedef struct {
86 klu_K_common Common;
87 klu_K_symbolic *Symbolic;
88 klu_K_numeric *Numeric;
89 PetscInt *perm_c, *perm_r;
90 MatStructure flg;
91 PetscBool PetscMatOrdering;
92 PetscBool CleanUpKLU;
93 } Mat_KLU;
94
MatDestroy_KLU(Mat A)95 static PetscErrorCode MatDestroy_KLU(Mat A)
96 {
97 Mat_KLU *lu = (Mat_KLU *)A->data;
98
99 PetscFunctionBegin;
100 if (lu->CleanUpKLU) {
101 klu_K_free_symbolic(&lu->Symbolic, &lu->Common);
102 klu_K_free_numeric(&lu->Numeric, &lu->Common);
103 PetscCall(PetscFree2(lu->perm_r, lu->perm_c));
104 }
105 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
106 PetscCall(PetscFree(A->data));
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
MatSolveTranspose_KLU(Mat A,Vec b,Vec x)110 static PetscErrorCode MatSolveTranspose_KLU(Mat A, Vec b, Vec x)
111 {
112 Mat_KLU *lu = (Mat_KLU *)A->data;
113 PetscScalar *xa;
114 PetscInt status;
115
116 PetscFunctionBegin;
117 /* KLU uses a column major format, solve Ax = b by klu_*_solve */
118 PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */
119 PetscCall(VecGetArray(x, &xa));
120 status = klu_K_solve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, (PetscReal *)xa, &lu->Common);
121 PetscCheck(status == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed");
122 PetscCall(VecRestoreArray(x, &xa));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
MatSolve_KLU(Mat A,Vec b,Vec x)126 static PetscErrorCode MatSolve_KLU(Mat A, Vec b, Vec x)
127 {
128 Mat_KLU *lu = (Mat_KLU *)A->data;
129 PetscScalar *xa;
130
131 PetscFunctionBegin;
132 /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */
133 PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */
134 PetscCall(VecGetArray(x, &xa));
135 PetscCheck(klu_K_tsolve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, xa, 1, &lu->Common), PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed"); /* conjugate solve */
136 PetscCall(VecRestoreArray(x, &xa));
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
MatLUFactorNumeric_KLU(Mat F,Mat A,const MatFactorInfo * info)140 static PetscErrorCode MatLUFactorNumeric_KLU(Mat F, Mat A, const MatFactorInfo *info)
141 {
142 Mat_KLU *lu = (Mat_KLU *)F->data;
143 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
144 PetscInt *ai = a->i, *aj = a->j;
145 PetscScalar *av = a->a;
146
147 PetscFunctionBegin;
148 /* numeric factorization of A' */
149
150 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) klu_K_free_numeric(&lu->Numeric, &lu->Common);
151 lu->Numeric = klu_K_factor(ai, aj, (PetscReal *)av, lu->Symbolic, &lu->Common);
152 PetscCheck(lu->Numeric, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Numeric factorization failed");
153
154 lu->flg = SAME_NONZERO_PATTERN;
155 lu->CleanUpKLU = PETSC_TRUE;
156 F->ops->solve = MatSolve_KLU;
157 F->ops->solvetranspose = MatSolveTranspose_KLU;
158 PetscFunctionReturn(PETSC_SUCCESS);
159 }
160
MatLUFactorSymbolic_KLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)161 static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
162 {
163 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
164 Mat_KLU *lu = (Mat_KLU *)F->data;
165 PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n;
166 const PetscInt *ra, *ca;
167
168 PetscFunctionBegin;
169 if (lu->PetscMatOrdering) {
170 PetscCall(ISGetIndices(r, &ra));
171 PetscCall(ISGetIndices(c, &ca));
172 PetscCall(PetscMalloc2(m, &lu->perm_r, n, &lu->perm_c));
173 /* we cannot simply memcpy on 64-bit archs */
174 for (i = 0; i < m; i++) lu->perm_r[i] = ra[i];
175 for (i = 0; i < n; i++) lu->perm_c[i] = ca[i];
176 PetscCall(ISRestoreIndices(r, &ra));
177 PetscCall(ISRestoreIndices(c, &ca));
178 }
179
180 /* symbolic factorization of A' */
181 if (r) {
182 lu->PetscMatOrdering = PETSC_TRUE;
183 lu->Symbolic = klu_K_analyze_given(n, ai, aj, lu->perm_c, lu->perm_r, &lu->Common);
184 } else { /* use klu internal ordering */
185 lu->Symbolic = klu_K_analyze(n, ai, aj, &lu->Common);
186 }
187 PetscCheck(lu->Symbolic, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Symbolic Factorization failed");
188
189 lu->flg = DIFFERENT_NONZERO_PATTERN;
190 lu->CleanUpKLU = PETSC_TRUE;
191 F->ops->lufactornumeric = MatLUFactorNumeric_KLU;
192 PetscFunctionReturn(PETSC_SUCCESS);
193 }
194
MatView_Info_KLU(Mat A,PetscViewer viewer)195 static PetscErrorCode MatView_Info_KLU(Mat A, PetscViewer viewer)
196 {
197 Mat_KLU *lu = (Mat_KLU *)A->data;
198 klu_K_numeric *Numeric = (klu_K_numeric *)lu->Numeric;
199
200 PetscFunctionBegin;
201 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU stats:\n"));
202 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of diagonal blocks: %" PetscInt_FMT "\n", (PetscInt)Numeric->nblocks));
203 PetscCall(PetscViewerASCIIPrintf(viewer, " Total nonzeros=%" PetscInt_FMT "\n", (PetscInt)(Numeric->lnz + Numeric->unz)));
204 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU runtime parameters:\n"));
205 /* Control parameters used by numeric factorization */
206 PetscCall(PetscViewerASCIIPrintf(viewer, " Partial pivoting tolerance: %g\n", lu->Common.tol));
207 /* BTF preordering */
208 PetscCall(PetscViewerASCIIPrintf(viewer, " BTF preordering enabled: %" PetscInt_FMT "\n", (PetscInt)lu->Common.btf));
209 /* mat ordering */
210 if (!lu->PetscMatOrdering) PetscCall(PetscViewerASCIIPrintf(viewer, " Ordering: %s (not using the PETSc ordering)\n", KluOrderingTypes[(int)lu->Common.ordering]));
211 /* matrix row scaling */
212 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n", scale[(int)lu->Common.scale]));
213 PetscFunctionReturn(PETSC_SUCCESS);
214 }
215
MatView_KLU(Mat A,PetscViewer viewer)216 static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer)
217 {
218 PetscBool isascii;
219 PetscViewerFormat format;
220
221 PetscFunctionBegin;
222 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
223 if (isascii) {
224 PetscCall(PetscViewerGetFormat(viewer, &format));
225 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_KLU(A, viewer));
226 }
227 PetscFunctionReturn(PETSC_SUCCESS);
228 }
229
MatFactorGetSolverType_seqaij_klu(Mat A,MatSolverType * type)230 static PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type)
231 {
232 PetscFunctionBegin;
233 *type = MATSOLVERKLU;
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
237 /*MC
238 MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices
239 via the external package KLU.
240
241 `./configure --download-suitesparse` to install PETSc to use KLU
242
243 Use `-pc_type lu` `-pc_factor_mat_solver_type klu` to use this direct solver
244
245 Consult KLU documentation for more information on the options database keys below.
246
247 Options Database Keys:
248 + -mat_klu_pivot_tol <0.001> - Partial pivoting tolerance
249 . -mat_klu_use_btf <1> - Use BTF preordering
250 . -mat_klu_ordering <AMD> - KLU reordering scheme to reduce fill-in (choose one of) `AMD`, `COLAMD`, `PETSC`
251 - -mat_klu_row_scale <NONE> - Matrix row scaling (choose one of) `NONE`, `SUM`, `MAX`
252
253 Level: beginner
254
255 Note:
256 KLU is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html>
257
258 .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType`
259 M*/
260
MatGetFactor_seqaij_klu(Mat A,MatFactorType ftype,Mat * F)261 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F)
262 {
263 Mat B;
264 Mat_KLU *lu;
265 PetscInt m = A->rmap->n, n = A->cmap->n, idx = 0, status;
266 PetscBool flg;
267
268 PetscFunctionBegin;
269 /* Create the factorization matrix F */
270 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
271 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
272 PetscCall(PetscStrallocpy("klu", &((PetscObject)B)->type_name));
273 PetscCall(MatSetUp(B));
274
275 PetscCall(PetscNew(&lu));
276
277 B->data = lu;
278 B->ops->getinfo = MatGetInfo_External;
279 B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
280 B->ops->destroy = MatDestroy_KLU;
281 B->ops->view = MatView_KLU;
282
283 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_klu));
284
285 B->factortype = MAT_FACTOR_LU;
286 B->assembled = PETSC_TRUE; /* required by -ksp_view */
287 B->preallocated = PETSC_TRUE;
288
289 PetscCall(PetscFree(B->solvertype));
290 PetscCall(PetscStrallocpy(MATSOLVERKLU, &B->solvertype));
291 B->canuseordering = PETSC_TRUE;
292 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));
293
294 /* initializations */
295 /* get the default control parameters */
296 status = klu_K_defaults(&lu->Common);
297 PetscCheck(status > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Initialization failed");
298
299 lu->Common.scale = 0; /* No row scaling */
300
301 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "KLU Options", "Mat");
302 /* Partial pivoting tolerance */
303 PetscCall(PetscOptionsReal("-mat_klu_pivot_tol", "Partial pivoting tolerance", "None", lu->Common.tol, &lu->Common.tol, NULL));
304 /* BTF pre-ordering */
305 PetscCall(PetscOptionsInt("-mat_klu_use_btf", "Enable BTF preordering", "None", (PetscInt)lu->Common.btf, (PetscInt *)&lu->Common.btf, NULL));
306 /* Matrix reordering */
307 PetscCall(PetscOptionsEList("-mat_klu_ordering", "Internal ordering method", "None", KluOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(KluOrderingTypes), KluOrderingTypes[0], &idx, &flg));
308 lu->Common.ordering = (int)idx;
309 /* Matrix row scaling */
310 PetscCall(PetscOptionsEList("-mat_klu_row_scale", "Matrix row scaling", "None", scale, 3, scale[0], &idx, &flg));
311 PetscOptionsEnd();
312 *F = B;
313 PetscFunctionReturn(PETSC_SUCCESS);
314 }
315