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