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 klu_l_analyze 16 #define klu_K_analyze_given klu_l_analyze_given 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 klu_zl_factor 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 klu_l_factor 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 80 #define SuiteSparse_long long long 81 #define SuiteSparse_long_max LONG_LONG_MAX 82 #define SuiteSparse_long_id "%lld" 83 84 EXTERN_C_BEGIN 85 #include <klu.h> 86 EXTERN_C_END 87 88 static const char *KluOrderingTypes[] = {"AMD","COLAMD","PETSC"}; 89 static const char *scale[] ={"NONE","SUM","MAX"}; 90 91 typedef struct { 92 klu_K_common Common; 93 klu_K_symbolic *Symbolic; 94 klu_K_numeric *Numeric; 95 PetscInt *perm_c,*perm_r; 96 MatStructure flg; 97 PetscBool PetscMatOrdering; 98 99 /* Flag to clean up KLU objects during Destroy */ 100 PetscBool CleanUpKLU; 101 } Mat_KLU; 102 103 #undef __FUNCT__ 104 #define __FUNCT__ "MatDestroy_KLU" 105 static PetscErrorCode MatDestroy_KLU(Mat A) 106 { 107 PetscErrorCode ierr; 108 Mat_KLU *lu=(Mat_KLU*)A->spptr; 109 110 PetscFunctionBegin; 111 if (lu && lu->CleanUpKLU) { 112 klu_K_free_symbolic(&lu->Symbolic,&lu->Common); 113 klu_K_free_numeric(&lu->Numeric,&lu->Common); 114 ierr = PetscFree2(lu->perm_r,lu->perm_c);CHKERRQ(ierr); 115 ierr = PetscFree(lu->perm_r);CHKERRQ(ierr); 116 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 117 } 118 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 119 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 #undef __FUNCT__ 124 #define __FUNCT__ "MatSolveTranspose_KLU" 125 static PetscErrorCode MatSolveTranspose_KLU(Mat A,Vec b,Vec x) 126 { 127 Mat_KLU *lu = (Mat_KLU*)A->spptr; 128 PetscScalar *xa; 129 PetscErrorCode ierr; 130 PetscInt status; 131 132 PetscFunctionBegin; 133 /* KLU uses a column major format, solve Ax = b by klu_*_solve */ 134 /* ----------------------------------*/ 135 ierr = VecCopy(b,x); /* klu_solve stores the solution in rhs */ 136 ierr = VecGetArray(x,&xa); 137 status = klu_K_solve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,&lu->Common); 138 if (status != 1) { 139 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed"); 140 } 141 ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatSolve_KLU" 147 static PetscErrorCode MatSolve_KLU(Mat A,Vec b,Vec x) 148 { 149 Mat_KLU *lu = (Mat_KLU*)A->spptr; 150 PetscScalar *xa; 151 PetscErrorCode ierr; 152 PetscInt status; 153 154 PetscFunctionBegin; 155 /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */ 156 /* ----------------------------------*/ 157 ierr = VecCopy(b,x); /* klu_solve stores the solution in rhs */ 158 ierr = VecGetArray(x,&xa); 159 #if defined(PETSC_USE_COMPLEX) 160 PetscInt conj_solve=1; 161 status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,conj_solve,&lu->Common); /* conjugate solve */ 162 #else 163 status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,xa,&lu->Common); 164 #endif 165 if (status != 1) { 166 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed"); 167 } 168 ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatLUFactorNumeric_KLU" 174 static PetscErrorCode MatLUFactorNumeric_KLU(Mat F,Mat A,const MatFactorInfo *info) 175 { 176 Mat_KLU *lu = (Mat_KLU*)(F)->spptr; 177 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 178 PetscInt *ai = a->i,*aj=a->j; 179 PetscScalar *av = a->a; 180 181 PetscFunctionBegin; 182 /* numeric factorization of A' */ 183 /* ----------------------------*/ 184 185 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) { 186 klu_K_free_numeric(&lu->Numeric,&lu->Common); 187 } 188 lu->Numeric = klu_K_factor(ai,aj,(PetscReal*)av,lu->Symbolic,&lu->Common); 189 if(!lu->Numeric) { 190 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Numeric factorization failed"); 191 } 192 193 lu->flg = SAME_NONZERO_PATTERN; 194 lu->CleanUpKLU = PETSC_TRUE; 195 F->ops->solve = MatSolve_KLU; 196 F->ops->solvetranspose = MatSolveTranspose_KLU; 197 PetscFunctionReturn(0); 198 } 199 200 #undef __FUNCT__ 201 #define __FUNCT__ "MatLUFactorSymbolic_KLU" 202 static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 203 { 204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 205 Mat_KLU *lu = (Mat_KLU*)(F->spptr); 206 PetscErrorCode ierr; 207 PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n; 208 const PetscInt *ra,*ca; 209 210 PetscFunctionBegin; 211 if (lu->PetscMatOrdering) { 212 ierr = ISGetIndices(r,&ra);CHKERRQ(ierr); 213 ierr = ISGetIndices(c,&ca);CHKERRQ(ierr); 214 ierr = PetscMalloc2(m,&lu->perm_c,n,&lu->perm_r);CHKERRQ(ierr); 215 /* we cannot simply memcpy on 64 bit archs */ 216 for (i = 0; i < m; i++) lu->perm_r[i] = ra[i]; 217 for (i = 0; i < n; i++) lu->perm_c[i] = ca[i]; 218 ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr); 219 ierr = ISRestoreIndices(c,&ca);CHKERRQ(ierr); 220 } 221 222 /* symbolic factorization of A' */ 223 /* ---------------------------------------------------------------------- */ 224 if (lu->PetscMatOrdering) { /* use Petsc ordering */ 225 lu->Symbolic = klu_K_analyze_given(n,ai,aj,lu->perm_c,lu->perm_r,&lu->Common); 226 } else { /* use klu internal ordering */ 227 lu->Symbolic = klu_K_analyze(n,ai,aj,&lu->Common); 228 } 229 if (!lu->Symbolic) { 230 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Symbolic Factorization failed"); 231 } 232 233 lu->flg = DIFFERENT_NONZERO_PATTERN; 234 lu->CleanUpKLU = PETSC_TRUE; 235 (F)->ops->lufactornumeric = MatLUFactorNumeric_KLU; 236 PetscFunctionReturn(0); 237 } 238 239 #undef __FUNCT__ 240 #define __FUNCT__ "MatFactorInfo_KLU" 241 static PetscErrorCode MatFactorInfo_KLU(Mat A,PetscViewer viewer) 242 { 243 Mat_KLU *lu= (Mat_KLU*)A->spptr; 244 klu_K_numeric *Numeric=(klu_K_numeric*)lu->Numeric; 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 /* check if matrix is KLU type */ 249 if (A->ops->solve != MatSolve_KLU) PetscFunctionReturn(0); 250 251 ierr = PetscViewerASCIIPrintf(viewer,"KLU stats:\n");CHKERRQ(ierr); 252 ierr = PetscViewerASCIIPrintf(viewer," Number of diagonal blocks: %d\n",Numeric->nblocks); 253 ierr = PetscViewerASCIIPrintf(viewer," Total nonzeros=%d\n",Numeric->lnz+Numeric->unz);CHKERRQ(ierr); 254 255 ierr = PetscViewerASCIIPrintf(viewer,"KLU runtime parameters:\n");CHKERRQ(ierr); 256 257 /* Control parameters used by numeric factorization */ 258 ierr = PetscViewerASCIIPrintf(viewer," Partial pivoting tolerance: %g\n",lu->Common.tol);CHKERRQ(ierr); 259 /* BTF preordering */ 260 ierr = PetscViewerASCIIPrintf(viewer," BTF preordering enabled: %d\n",lu->Common.btf);CHKERRQ(ierr); 261 /* mat ordering */ 262 if (!lu->PetscMatOrdering) { 263 ierr = PetscViewerASCIIPrintf(viewer," Ordering: %s (not using the PETSc ordering)\n",KluOrderingTypes[(int)lu->Common.ordering]);CHKERRQ(ierr); 264 } else { 265 ierr = PetscViewerASCIIPrintf(viewer," Using PETSc ordering\n");CHKERRQ(ierr); 266 } 267 /* matrix row scaling */ 268 ierr = PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n",scale[(int)lu->Common.scale]);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "MatView_KLU" 274 static PetscErrorCode MatView_KLU(Mat A,PetscViewer viewer) 275 { 276 PetscErrorCode ierr; 277 PetscBool iascii; 278 PetscViewerFormat format; 279 280 PetscFunctionBegin; 281 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 282 283 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 284 if (iascii) { 285 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 286 if (format == PETSC_VIEWER_ASCII_INFO) { 287 ierr = MatFactorInfo_KLU(A,viewer);CHKERRQ(ierr); 288 } 289 } 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_klu" 295 PetscErrorCode MatFactorGetSolverPackage_seqaij_klu(Mat A,const MatSolverPackage *type) 296 { 297 PetscFunctionBegin; 298 *type = MATSOLVERKLU; 299 PetscFunctionReturn(0); 300 } 301 302 303 /*MC 304 MATSOLVERKLU = "klu" - A matrix type providing direct solvers (LU) for sequential matrices 305 via the external package KLU. 306 307 ./configure --download-SuiteSparse --SuiteSparse-enable-klu to install PETSc to use KLU 308 309 Consult KLU documentation for more information on the options database keys below. 310 311 Options Database Keys: 312 + -mat_klu_pivot_tol <0.001> - Partial pivoting tolerance 313 . -mat_klu_use_btf <1> - Use BTF preordering 314 . -mat_klu_ordering <AMD> - KLU reordering scheme to reduce fill-in (choose one of) AMD COLAMD PETSC 315 - -mat_klu_row_scale <NONE> - Matrix row scaling (choose one of) NONE SUM MAX 316 317 Level: beginner 318 319 .seealso: PCLU, MATSOLVERUMFPACK, MATSOLVERCHOLMOD, PCFactorSetMatSolverPackage(), MatSolverPackage 320 M*/ 321 322 #undef __FUNCT__ 323 #define __FUNCT__ "MatGetFactor_seqaij_klu" 324 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A,MatFactorType ftype,Mat *F) 325 { 326 Mat B; 327 Mat_KLU *lu; 328 PetscErrorCode ierr; 329 PetscInt m=A->rmap->n,n=A->cmap->n,idx,status; 330 PetscBool flg; 331 332 PetscFunctionBegin; 333 /* Create the factorization matrix F */ 334 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 335 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 336 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 337 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 338 ierr = PetscNewLog(B,&lu);CHKERRQ(ierr); 339 340 B->spptr = lu; 341 B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU; 342 B->ops->destroy = MatDestroy_KLU; 343 B->ops->view = MatView_KLU; 344 345 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_klu);CHKERRQ(ierr); 346 347 B->factortype = MAT_FACTOR_LU; 348 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 349 B->preallocated = PETSC_TRUE; 350 351 /* initializations */ 352 /* ------------------------------------------------*/ 353 /* get the default control parameters */ 354 status = klu_K_defaults(&lu->Common); 355 if(status <= 0) { 356 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Initialization failed"); 357 } 358 lu->Common.scale = 0; /* No row scaling */ 359 360 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"KLU Options","Mat");CHKERRQ(ierr); 361 /* Partial pivoting tolerance */ 362 ierr = PetscOptionsReal("-mat_klu_pivot_tol","Partial pivoting tolerance","None",lu->Common.tol,&lu->Common.tol,NULL);CHKERRQ(ierr); 363 /* BTF pre-ordering */ 364 ierr = PetscOptionsInt("-mat_klu_use_btf","Enable BTF preordering","None",lu->Common.btf,&lu->Common.btf,NULL);CHKERRQ(ierr); 365 /* Matrix reordering */ 366 ierr = PetscOptionsEList("-mat_klu_ordering","Internal ordering method","None",KluOrderingTypes,sizeof(KluOrderingTypes)/sizeof(KluOrderingTypes[0]),KluOrderingTypes[0],&idx,&flg);CHKERRQ(ierr); 367 if (flg) { 368 if ((int)idx == 2) lu->PetscMatOrdering = PETSC_TRUE; /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Klu perm_c) */ 369 else lu->Common.ordering = (int)idx; 370 } 371 /* Matrix row scaling */ 372 ierr = PetscOptionsEList("-mat_klu_row_scale","Matrix row scaling","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr); 373 PetscOptionsEnd(); 374 *F = B; 375 PetscFunctionReturn(0); 376 } 377