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