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 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 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 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 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 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 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 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 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 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