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 klu_zl_tsolve 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 klu_l_tsolve 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 klu_z_tsolve 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 klu_tsolve 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 PetscInt status; 131 132 PetscFunctionBegin; 133 /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */ 134 PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */ 135 PetscCall(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 142 PetscCheck(status == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed"); 143 PetscCall(VecRestoreArray(x, &xa)); 144 PetscFunctionReturn(PETSC_SUCCESS); 145 } 146 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; 153 154 PetscFunctionBegin; 155 /* numeric factorization of A' */ 156 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); 159 PetscCheck(lu->Numeric, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Numeric factorization failed"); 160 161 lu->flg = SAME_NONZERO_PATTERN; 162 lu->CleanUpKLU = PETSC_TRUE; 163 F->ops->solve = MatSolve_KLU; 164 F->ops->solvetranspose = MatSolveTranspose_KLU; 165 PetscFunctionReturn(PETSC_SUCCESS); 166 } 167 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; 174 175 PetscFunctionBegin; 176 if (lu->PetscMatOrdering) { 177 PetscCall(ISGetIndices(r, &ra)); 178 PetscCall(ISGetIndices(c, &ca)); 179 PetscCall(PetscMalloc2(m, &lu->perm_r, n, &lu->perm_c)); 180 /* we cannot simply memcpy on 64-bit archs */ 181 for (i = 0; i < m; i++) lu->perm_r[i] = ra[i]; 182 for (i = 0; i < n; i++) lu->perm_c[i] = ca[i]; 183 PetscCall(ISRestoreIndices(r, &ra)); 184 PetscCall(ISRestoreIndices(c, &ca)); 185 } 186 187 /* symbolic factorization of A' */ 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 } 194 PetscCheck(lu->Symbolic, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Symbolic Factorization failed"); 195 196 lu->flg = DIFFERENT_NONZERO_PATTERN; 197 lu->CleanUpKLU = PETSC_TRUE; 198 F->ops->lufactornumeric = MatLUFactorNumeric_KLU; 199 PetscFunctionReturn(PETSC_SUCCESS); 200 } 201 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; 206 207 PetscFunctionBegin; 208 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU stats:\n")); 209 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of diagonal blocks: %" PetscInt_FMT "\n", (PetscInt)Numeric->nblocks)); 210 PetscCall(PetscViewerASCIIPrintf(viewer, " Total nonzeros=%" PetscInt_FMT "\n", (PetscInt)(Numeric->lnz + Numeric->unz))); 211 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU runtime parameters:\n")); 212 /* Control parameters used by numeric factorization */ 213 PetscCall(PetscViewerASCIIPrintf(viewer, " Partial pivoting tolerance: %g\n", lu->Common.tol)); 214 /* BTF preordering */ 215 PetscCall(PetscViewerASCIIPrintf(viewer, " BTF preordering enabled: %" PetscInt_FMT "\n", (PetscInt)lu->Common.btf)); 216 /* mat ordering */ 217 if (!lu->PetscMatOrdering) PetscCall(PetscViewerASCIIPrintf(viewer, " Ordering: %s (not using the PETSc ordering)\n", KluOrderingTypes[(int)lu->Common.ordering])); 218 /* matrix row scaling */ 219 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n", scale[(int)lu->Common.scale])); 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer) 224 { 225 PetscBool isascii; 226 PetscViewerFormat format; 227 228 PetscFunctionBegin; 229 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 230 if (isascii) { 231 PetscCall(PetscViewerGetFormat(viewer, &format)); 232 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_KLU(A, viewer)); 233 } 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 static PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type) 238 { 239 PetscFunctionBegin; 240 *type = MATSOLVERKLU; 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /*MC 245 MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices 246 via the external package KLU. 247 248 `./configure --download-suitesparse` to install PETSc to use KLU 249 250 Use `-pc_type lu` `-pc_factor_mat_solver_type klu` to use this direct solver 251 252 Consult KLU documentation for more information on the options database keys below. 253 254 Options Database Keys: 255 + -mat_klu_pivot_tol <0.001> - Partial pivoting tolerance 256 . -mat_klu_use_btf <1> - Use BTF preordering 257 . -mat_klu_ordering <AMD> - KLU reordering scheme to reduce fill-in (choose one of) `AMD`, `COLAMD`, `PETSC` 258 - -mat_klu_row_scale <NONE> - Matrix row scaling (choose one of) `NONE`, `SUM`, `MAX` 259 260 Level: beginner 261 262 Note: 263 KLU is part of SuiteSparse <http://faculty.cse.tamu.edu/davis/suitesparse.html> 264 265 .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType` 266 M*/ 267 268 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F) 269 { 270 Mat B; 271 Mat_KLU *lu; 272 PetscInt m = A->rmap->n, n = A->cmap->n, idx = 0, status; 273 PetscBool flg; 274 275 PetscFunctionBegin; 276 /* Create the factorization matrix F */ 277 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 278 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 279 PetscCall(PetscStrallocpy("klu", &((PetscObject)B)->type_name)); 280 PetscCall(MatSetUp(B)); 281 282 PetscCall(PetscNew(&lu)); 283 284 B->data = lu; 285 B->ops->getinfo = MatGetInfo_External; 286 B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU; 287 B->ops->destroy = MatDestroy_KLU; 288 B->ops->view = MatView_KLU; 289 290 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_klu)); 291 292 B->factortype = MAT_FACTOR_LU; 293 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 294 B->preallocated = PETSC_TRUE; 295 296 PetscCall(PetscFree(B->solvertype)); 297 PetscCall(PetscStrallocpy(MATSOLVERKLU, &B->solvertype)); 298 B->canuseordering = PETSC_TRUE; 299 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU])); 300 301 /* initializations */ 302 /* get the default control parameters */ 303 status = klu_K_defaults(&lu->Common); 304 PetscCheck(status > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Initialization failed"); 305 306 lu->Common.scale = 0; /* No row scaling */ 307 308 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "KLU Options", "Mat"); 309 /* Partial pivoting tolerance */ 310 PetscCall(PetscOptionsReal("-mat_klu_pivot_tol", "Partial pivoting tolerance", "None", lu->Common.tol, &lu->Common.tol, NULL)); 311 /* BTF pre-ordering */ 312 PetscCall(PetscOptionsInt("-mat_klu_use_btf", "Enable BTF preordering", "None", (PetscInt)lu->Common.btf, (PetscInt *)&lu->Common.btf, NULL)); 313 /* Matrix reordering */ 314 PetscCall(PetscOptionsEList("-mat_klu_ordering", "Internal ordering method", "None", KluOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(KluOrderingTypes), KluOrderingTypes[0], &idx, &flg)); 315 lu->Common.ordering = (int)idx; 316 /* Matrix row scaling */ 317 PetscCall(PetscOptionsEList("-mat_klu_row_scale", "Matrix row scaling", "None", scale, 3, scale[0], &idx, &flg)); 318 PetscOptionsEnd(); 319 *F = B; 320 PetscFunctionReturn(PETSC_SUCCESS); 321 } 322