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 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(0); 108 } 109 110 static PetscErrorCode MatSolveTranspose_KLU(Mat A, Vec b, Vec x) { 111 Mat_KLU *lu = (Mat_KLU *)A->data; 112 PetscScalar *xa; 113 PetscInt status; 114 115 PetscFunctionBegin; 116 /* KLU uses a column major format, solve Ax = b by klu_*_solve */ 117 /* ----------------------------------*/ 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(0); 124 } 125 126 static PetscErrorCode MatSolve_KLU(Mat A, Vec b, Vec x) { 127 Mat_KLU *lu = (Mat_KLU *)A->data; 128 PetscScalar *xa; 129 PetscInt status; 130 131 PetscFunctionBegin; 132 /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */ 133 /* ----------------------------------*/ 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(0); 145 } 146 147 static PetscErrorCode MatLUFactorNumeric_KLU(Mat F, Mat A, const MatFactorInfo *info) { 148 Mat_KLU *lu = (Mat_KLU *)(F)->data; 149 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 150 PetscInt *ai = a->i, *aj = a->j; 151 PetscScalar *av = a->a; 152 153 PetscFunctionBegin; 154 /* numeric factorization of A' */ 155 /* ----------------------------*/ 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(0); 166 } 167 168 static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info) { 169 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 170 Mat_KLU *lu = (Mat_KLU *)(F->data); 171 PetscInt i, *ai = a->i, *aj = a->j, m = A->rmap->n, n = A->cmap->n; 172 const PetscInt *ra, *ca; 173 174 PetscFunctionBegin; 175 if (lu->PetscMatOrdering) { 176 PetscCall(ISGetIndices(r, &ra)); 177 PetscCall(ISGetIndices(c, &ca)); 178 PetscCall(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 PetscCall(ISRestoreIndices(r, &ra)); 183 PetscCall(ISRestoreIndices(c, &ca)); 184 } 185 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 } 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(0); 200 } 201 202 static PetscErrorCode MatView_Info_KLU(Mat A, PetscViewer viewer) { 203 Mat_KLU *lu = (Mat_KLU *)A->data; 204 klu_K_numeric *Numeric = (klu_K_numeric *)lu->Numeric; 205 206 PetscFunctionBegin; 207 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU stats:\n")); 208 PetscCall(PetscViewerASCIIPrintf(viewer, " Number of diagonal blocks: %" PetscInt_FMT "\n", (PetscInt)(Numeric->nblocks))); 209 PetscCall(PetscViewerASCIIPrintf(viewer, " Total nonzeros=%" PetscInt_FMT "\n", (PetscInt)(Numeric->lnz + Numeric->unz))); 210 PetscCall(PetscViewerASCIIPrintf(viewer, "KLU runtime parameters:\n")); 211 /* Control parameters used by numeric factorization */ 212 PetscCall(PetscViewerASCIIPrintf(viewer, " Partial pivoting tolerance: %g\n", lu->Common.tol)); 213 /* BTF preordering */ 214 PetscCall(PetscViewerASCIIPrintf(viewer, " BTF preordering enabled: %" PetscInt_FMT "\n", (PetscInt)(lu->Common.btf))); 215 /* mat ordering */ 216 if (!lu->PetscMatOrdering) PetscCall(PetscViewerASCIIPrintf(viewer, " Ordering: %s (not using the PETSc ordering)\n", KluOrderingTypes[(int)lu->Common.ordering])); 217 /* matrix row scaling */ 218 PetscCall(PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n", scale[(int)lu->Common.scale])); 219 PetscFunctionReturn(0); 220 } 221 222 static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer) { 223 PetscBool iascii; 224 PetscViewerFormat format; 225 226 PetscFunctionBegin; 227 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 228 if (iascii) { 229 PetscCall(PetscViewerGetFormat(viewer, &format)); 230 if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_KLU(A, viewer)); 231 } 232 PetscFunctionReturn(0); 233 } 234 235 PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type) { 236 PetscFunctionBegin; 237 *type = MATSOLVERKLU; 238 PetscFunctionReturn(0); 239 } 240 241 /*MC 242 MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices 243 via the external package KLU. 244 245 ./configure --download-suitesparse to install PETSc to use KLU 246 247 Use -pc_type lu -pc_factor_mat_solver_type klu to use this direct solver 248 249 Consult KLU documentation for more information on the options database keys below. 250 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 256 257 Note: KLU is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 258 259 Level: beginner 260 261 .seealso: `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType` 262 M*/ 263 264 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F) { 265 Mat B; 266 Mat_KLU *lu; 267 PetscInt m = A->rmap->n, n = A->cmap->n, idx = 0, status; 268 PetscBool flg; 269 270 PetscFunctionBegin; 271 /* Create the factorization matrix F */ 272 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 273 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n)); 274 PetscCall(PetscStrallocpy("klu", &((PetscObject)B)->type_name)); 275 PetscCall(MatSetUp(B)); 276 277 PetscCall(PetscNewLog(B, &lu)); 278 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; 284 285 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_klu)); 286 287 B->factortype = MAT_FACTOR_LU; 288 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 289 B->preallocated = PETSC_TRUE; 290 291 PetscCall(PetscFree(B->solvertype)); 292 PetscCall(PetscStrallocpy(MATSOLVERKLU, &B->solvertype)); 293 B->canuseordering = PETSC_TRUE; 294 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU])); 295 296 /* initializations */ 297 /* ------------------------------------------------*/ 298 /* get the default control parameters */ 299 status = klu_K_defaults(&lu->Common); 300 PetscCheck(status > 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "KLU Initialization failed"); 301 302 lu->Common.scale = 0; /* No row scaling */ 303 304 PetscOptionsBegin(PetscObjectComm((PetscObject)B), ((PetscObject)B)->prefix, "KLU Options", "Mat"); 305 /* Partial pivoting tolerance */ 306 PetscCall(PetscOptionsReal("-mat_klu_pivot_tol", "Partial pivoting tolerance", "None", lu->Common.tol, &lu->Common.tol, NULL)); 307 /* BTF pre-ordering */ 308 PetscCall(PetscOptionsInt("-mat_klu_use_btf", "Enable BTF preordering", "None", (PetscInt)lu->Common.btf, (PetscInt *)&lu->Common.btf, NULL)); 309 /* Matrix reordering */ 310 PetscCall(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 PetscCall(PetscOptionsEList("-mat_klu_row_scale", "Matrix row scaling", "None", scale, 3, scale[0], &idx, &flg)); 314 PetscOptionsEnd(); 315 *F = B; 316 PetscFunctionReturn(0); 317 } 318