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