1 2 /* 3 Provides an interface to the KLUv1.2 sparse solver 4 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. 9 10 */ 11 #include <../src/mat/impls/aij/seq/aij.h> 12 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 78 79 EXTERN_C_BEGIN 80 #include <klu.h> 81 EXTERN_C_END 82 83 static const char *KluOrderingTypes[] = {"AMD", "COLAMD"}; 84 static const char *scale[] = {"NONE", "SUM", "MAX"}; 85 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; 95 96 static PetscErrorCode MatDestroy_KLU(Mat A) 97 { 98 Mat_KLU *lu = (Mat_KLU *)A->data; 99 100 PetscFunctionBegin; 101 if (lu->CleanUpKLU) { 102 klu_K_free_symbolic(&lu->Symbolic, &lu->Common); 103 klu_K_free_numeric(&lu->Numeric, &lu->Common); 104 PetscCall(PetscFree2(lu->perm_r, lu->perm_c)); 105 } 106 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL)); 107 PetscCall(PetscFree(A->data)); 108 PetscFunctionReturn(PETSC_SUCCESS); 109 } 110 111 static PetscErrorCode MatSolveTranspose_KLU(Mat A, Vec b, Vec x) 112 { 113 Mat_KLU *lu = (Mat_KLU *)A->data; 114 PetscScalar *xa; 115 PetscInt status; 116 117 PetscFunctionBegin; 118 /* KLU uses a column major format, solve Ax = b by klu_*_solve */ 119 PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */ 120 PetscCall(VecGetArray(x, &xa)); 121 status = klu_K_solve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, (PetscReal *)xa, &lu->Common); 122 PetscCheck(status == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed"); 123 PetscCall(VecRestoreArray(x, &xa)); 124 PetscFunctionReturn(PETSC_SUCCESS); 125 } 126 127 static PetscErrorCode MatSolve_KLU(Mat A, Vec b, Vec x) 128 { 129 Mat_KLU *lu = (Mat_KLU *)A->data; 130 PetscScalar *xa; 131 PetscInt status; 132 133 PetscFunctionBegin; 134 /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */ 135 PetscCall(VecCopy(b, x)); /* klu_solve stores the solution in rhs */ 136 PetscCall(VecGetArray(x, &xa)); 137 #if defined(PETSC_USE_COMPLEX) 138 PetscInt conj_solve = 1; 139 status = klu_K_tsolve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, (PetscReal *)xa, conj_solve, &lu->Common); /* conjugate solve */ 140 #else 141 status = klu_K_tsolve(lu->Symbolic, lu->Numeric, A->rmap->n, 1, xa, &lu->Common); 142 #endif 143 PetscCheck(status == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Solve failed"); 144 PetscCall(VecRestoreArray(x, &xa)); 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 static PetscErrorCode MatLUFactorNumeric_KLU(Mat F, Mat A, const MatFactorInfo *info) 149 { 150 Mat_KLU *lu = (Mat_KLU *)(F)->data; 151 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 152 PetscInt *ai = a->i, *aj = a->j; 153 PetscScalar *av = a->a; 154 155 PetscFunctionBegin; 156 /* numeric factorization of A' */ 157 158 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) klu_K_free_numeric(&lu->Numeric, &lu->Common); 159 lu->Numeric = klu_K_factor(ai, aj, (PetscReal *)av, lu->Symbolic, &lu->Common); 160 PetscCheck(lu->Numeric, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Numeric factorization failed"); 161 162 lu->flg = SAME_NONZERO_PATTERN; 163 lu->CleanUpKLU = PETSC_TRUE; 164 F->ops->solve = MatSolve_KLU; 165 F->ops->solvetranspose = MatSolveTranspose_KLU; 166 PetscFunctionReturn(PETSC_SUCCESS); 167 } 168 169 static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) 170 { 171 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 172 Mat_KLU *lu = (Mat_KLU *)(F->data); 173 PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n; 174 const PetscInt *ra, *ca; 175 176 PetscFunctionBegin; 177 if (lu->PetscMatOrdering) { 178 PetscCall(ISGetIndices(r, &ra)); 179 PetscCall(ISGetIndices(c, &ca)); 180 PetscCall(PetscMalloc2(m, &lu->perm_r, n, &lu->perm_c)); 181 /* we cannot simply memcpy on 64-bit archs */ 182 for (i = 0; i < m; i++) lu->perm_r[i] = ra[i]; 183 for (i = 0; i < n; i++) lu->perm_c[i] = ca[i]; 184 PetscCall(ISRestoreIndices(r, &ra)); 185 PetscCall(ISRestoreIndices(c, &ca)); 186 } 187 188 /* symbolic factorization of A' */ 189 if (r) { 190 lu->PetscMatOrdering = PETSC_TRUE; 191 lu->Symbolic = klu_K_analyze_given(n, ai, aj, lu->perm_c, lu->perm_r, &lu->Common); 192 } else { /* use klu internal ordering */ 193 lu->Symbolic = klu_K_analyze(n, ai, aj, &lu->Common); 194 } 195 PetscCheck(lu->Symbolic, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Symbolic Factorization failed"); 196 197 lu->flg = DIFFERENT_NONZERO_PATTERN; 198 lu->CleanUpKLU = PETSC_TRUE; 199 (F)->ops->lufactornumeric = MatLUFactorNumeric_KLU; 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 203 static PetscErrorCode MatView_Info_KLU(Mat A, PetscViewer viewer) 204 { 205 Mat_KLU *lu = (Mat_KLU *)A->data; 206 klu_K_numeric *Numeric = (klu_K_numeric *)lu->Numeric; 207 208 PetscFunctionBegin; 209 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU stats:\n")); 210 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of diagonal blocks: %" PetscInt_FMT "\n", (PetscInt)(Numeric->nblocks))); 211 PetscCall(PetscViewerASCIIPrintf(viewer, " Total nonzeros=%" PetscInt_FMT "\n", (PetscInt)(Numeric->lnz + Numeric->unz))); 212 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU runtime parameters:\n")); 213 /* Control parameters used by numeric factorization */ 214 PetscCall(PetscViewerASCIIPrintf(viewer, " Partial pivoting tolerance: %g\n", lu->Common.tol)); 215 /* BTF preordering */ 216 PetscCall(PetscViewerASCIIPrintf(viewer, " BTF preordering enabled: %" PetscInt_FMT "\n", (PetscInt)(lu->Common.btf))); 217 /* mat ordering */ 218 if (!lu->PetscMatOrdering) PetscCall(PetscViewerASCIIPrintf(viewer, " Ordering: %s (not using the PETSc ordering)\n", KluOrderingTypes[(int)lu->Common.ordering])); 219 /* matrix row scaling */ 220 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n", scale[(int)lu->Common.scale])); 221 PetscFunctionReturn(PETSC_SUCCESS); 222 } 223 224 static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer) 225 { 226 PetscBool iascii; 227 PetscViewerFormat format; 228 229 PetscFunctionBegin; 230 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 231 if (iascii) { 232 PetscCall(PetscViewerGetFormat(viewer, &format)); 233 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_KLU(A, viewer)); 234 } 235 PetscFunctionReturn(PETSC_SUCCESS); 236 } 237 238 PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type) 239 { 240 PetscFunctionBegin; 241 *type = MATSOLVERKLU; 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 /*MC 246 MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices 247 via the external package KLU. 248 249 `./configure --download-suitesparse` to install PETSc to use KLU 250 251 Use `-pc_type lu` `-pc_factor_mat_solver_type klu` to use this direct solver 252 253 Consult KLU documentation for more information on the options database keys below. 254 255 Options Database Keys: 256 + -mat_klu_pivot_tol <0.001> - Partial pivoting tolerance 257 . -mat_klu_use_btf <1> - Use BTF preordering 258 . -mat_klu_ordering <AMD> - KLU reordering scheme to reduce fill-in (choose one of) `AMD`, `COLAMD`, `PETSC` 259 - -mat_klu_row_scale <NONE> - Matrix row scaling (choose one of) `NONE`, `SUM`, `MAX` 260 261 Level: beginner 262 263 Note: 264 KLU is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 265 266 .seealso: [](ch_matrices), `Mat`, `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType` 267 M*/ 268 269 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F) 270 { 271 Mat B; 272 Mat_KLU *lu; 273 PetscInt m = A->rmap->n, n = A->cmap->n, idx = 0, status; 274 PetscBool flg; 275 276 PetscFunctionBegin; 277 /* Create the factorization matrix F */ 278 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 279 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 280 PetscCall(PetscStrallocpy("klu", &((PetscObject)B)->type_name)); 281 PetscCall(MatSetUp(B)); 282 283 PetscCall(PetscNew(&lu)); 284 285 B->data = lu; 286 B->ops->getinfo = MatGetInfo_External; 287 B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU; 288 B->ops->destroy = MatDestroy_KLU; 289 B->ops->view = MatView_KLU; 290 291 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_klu)); 292 293 B->factortype = MAT_FACTOR_LU; 294 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 295 B->preallocated = PETSC_TRUE; 296 297 PetscCall(PetscFree(B->solvertype)); 298 PetscCall(PetscStrallocpy(MATSOLVERKLU, &B->solvertype)); 299 B->canuseordering = PETSC_TRUE; 300 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU])); 301 302 /* initializations */ 303 /* get the default control parameters */ 304 status = klu_K_defaults(&lu->Common); 305 PetscCheck(status > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Initialization failed"); 306 307 lu->Common.scale = 0; /* No row scaling */ 308 309 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "KLU Options", "Mat"); 310 /* Partial pivoting tolerance */ 311 PetscCall(PetscOptionsReal("-mat_klu_pivot_tol", "Partial pivoting tolerance", "None", lu->Common.tol, &lu->Common.tol, NULL)); 312 /* BTF pre-ordering */ 313 PetscCall(PetscOptionsInt("-mat_klu_use_btf", "Enable BTF preordering", "None", (PetscInt)lu->Common.btf, (PetscInt *)&lu->Common.btf, NULL)); 314 /* Matrix reordering */ 315 PetscCall(PetscOptionsEList("-mat_klu_ordering", "Internal ordering method", "None", KluOrderingTypes, PETSC_STATIC_ARRAY_LENGTH(KluOrderingTypes), KluOrderingTypes[0], &idx, &flg)); 316 lu->Common.ordering = (int)idx; 317 /* Matrix row scaling */ 318 PetscCall(PetscOptionsEList("-mat_klu_row_scale", "Matrix row scaling", "None", scale, 3, scale[0], &idx, &flg)); 319 PetscOptionsEnd(); 320 *F = B; 321 PetscFunctionReturn(PETSC_SUCCESS); 322 } 323