1 2 /* 3 Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1 4 5 When build with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the 6 integer type in UMFPACK, otherwise it will use int. This means 7 all integers in this file as simply declared as PetscInt. Also it means 8 that one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only] 9 10 */ 11 #include <../src/mat/impls/aij/seq/aij.h> 12 13 #if defined(PETSC_USE_64BIT_INDICES) 14 #if defined(PETSC_USE_COMPLEX) 15 #define umfpack_UMF_free_symbolic umfpack_zl_free_symbolic 16 #define umfpack_UMF_free_numeric umfpack_zl_free_numeric 17 /* the type casts are needed because PetscInt is long long while SuiteSparse_long is long and compilers warn even when they are identical */ 18 #define umfpack_UMF_wsolve(a,b,c,d,e,f,g,h,i,j,k,l,m,n) umfpack_zl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k,l,(SuiteSparse_long*)m,n) 19 #define umfpack_UMF_numeric(a,b,c,d,e,f,g,h) umfpack_zl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g,h) 20 #define umfpack_UMF_report_numeric umfpack_zl_report_numeric 21 #define umfpack_UMF_report_control umfpack_zl_report_control 22 #define umfpack_UMF_report_status umfpack_zl_report_status 23 #define umfpack_UMF_report_info umfpack_zl_report_info 24 #define umfpack_UMF_report_symbolic umfpack_zl_report_symbolic 25 #define umfpack_UMF_qsymbolic(a,b,c,d,e,f,g,h,i,j) umfpack_zl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,(SuiteSparse_long*)g,h,i,j) 26 #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h,i) umfpack_zl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h,i) 27 #define umfpack_UMF_defaults umfpack_zl_defaults 28 29 #else 30 #define umfpack_UMF_free_symbolic umfpack_dl_free_symbolic 31 #define umfpack_UMF_free_numeric umfpack_dl_free_numeric 32 #define umfpack_UMF_wsolve(a,b,c,d,e,f,g,h,i,j,k) umfpack_dl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k) 33 #define umfpack_UMF_numeric(a,b,c,d,e,f,g) umfpack_dl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g) 34 #define umfpack_UMF_report_numeric umfpack_dl_report_numeric 35 #define umfpack_UMF_report_control umfpack_dl_report_control 36 #define umfpack_UMF_report_status umfpack_dl_report_status 37 #define umfpack_UMF_report_info umfpack_dl_report_info 38 #define umfpack_UMF_report_symbolic umfpack_dl_report_symbolic 39 #define umfpack_UMF_qsymbolic(a,b,c,d,e,f,g,h,i) umfpack_dl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,(SuiteSparse_long*)f,g,h,i) 40 #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h) umfpack_dl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h) 41 #define umfpack_UMF_defaults umfpack_dl_defaults 42 #endif 43 44 #else 45 #if defined(PETSC_USE_COMPLEX) 46 #define umfpack_UMF_free_symbolic umfpack_zi_free_symbolic 47 #define umfpack_UMF_free_numeric umfpack_zi_free_numeric 48 #define umfpack_UMF_wsolve umfpack_zi_wsolve 49 #define umfpack_UMF_numeric umfpack_zi_numeric 50 #define umfpack_UMF_report_numeric umfpack_zi_report_numeric 51 #define umfpack_UMF_report_control umfpack_zi_report_control 52 #define umfpack_UMF_report_status umfpack_zi_report_status 53 #define umfpack_UMF_report_info umfpack_zi_report_info 54 #define umfpack_UMF_report_symbolic umfpack_zi_report_symbolic 55 #define umfpack_UMF_qsymbolic umfpack_zi_qsymbolic 56 #define umfpack_UMF_symbolic umfpack_zi_symbolic 57 #define umfpack_UMF_defaults umfpack_zi_defaults 58 59 #else 60 #define umfpack_UMF_free_symbolic umfpack_di_free_symbolic 61 #define umfpack_UMF_free_numeric umfpack_di_free_numeric 62 #define umfpack_UMF_wsolve umfpack_di_wsolve 63 #define umfpack_UMF_numeric umfpack_di_numeric 64 #define umfpack_UMF_report_numeric umfpack_di_report_numeric 65 #define umfpack_UMF_report_control umfpack_di_report_control 66 #define umfpack_UMF_report_status umfpack_di_report_status 67 #define umfpack_UMF_report_info umfpack_di_report_info 68 #define umfpack_UMF_report_symbolic umfpack_di_report_symbolic 69 #define umfpack_UMF_qsymbolic umfpack_di_qsymbolic 70 #define umfpack_UMF_symbolic umfpack_di_symbolic 71 #define umfpack_UMF_defaults umfpack_di_defaults 72 #endif 73 #endif 74 75 EXTERN_C_BEGIN 76 #include <umfpack.h> 77 EXTERN_C_END 78 79 static const char *const UmfpackOrderingTypes[] = {"CHOLMOD","AMD","GIVEN","METIS","BEST","NONE","USER","UmfpackOrderingTypes","UMFPACK_ORDERING_",0}; 80 81 typedef struct { 82 void *Symbolic, *Numeric; 83 double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W; 84 PetscInt *Wi,*perm_c; 85 Mat A; /* Matrix used for factorization */ 86 MatStructure flg; 87 88 /* Flag to clean up UMFPACK objects during Destroy */ 89 PetscBool CleanUpUMFPACK; 90 } Mat_UMFPACK; 91 92 static PetscErrorCode MatDestroy_UMFPACK(Mat A) 93 { 94 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->data; 95 96 PetscFunctionBegin; 97 if (lu->CleanUpUMFPACK) { 98 umfpack_UMF_free_symbolic(&lu->Symbolic); 99 umfpack_UMF_free_numeric(&lu->Numeric); 100 PetscCall(PetscFree(lu->Wi)); 101 PetscCall(PetscFree(lu->W)); 102 PetscCall(PetscFree(lu->perm_c)); 103 } 104 PetscCall(MatDestroy(&lu->A)); 105 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL)); 106 PetscCall(PetscFree(A->data)); 107 PetscFunctionReturn(0); 108 } 109 110 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag) 111 { 112 Mat_UMFPACK *lu = (Mat_UMFPACK*)A->data; 113 Mat_SeqAIJ *a = (Mat_SeqAIJ*)lu->A->data; 114 PetscScalar *av = a->a,*xa; 115 const PetscScalar *ba; 116 PetscInt *ai = a->i,*aj = a->j,status; 117 static PetscBool cite = PETSC_FALSE; 118 119 PetscFunctionBegin; 120 if (!A->rmap->n) PetscFunctionReturn(0); 121 PetscCall(PetscCitationsRegister("@article{davis2004algorithm,\n title={Algorithm 832: {UMFPACK} V4.3---An Unsymmetric-Pattern Multifrontal Method},\n author={Davis, Timothy A},\n journal={ACM Transactions on Mathematical Software (TOMS)},\n volume={30},\n number={2},\n pages={196--199},\n year={2004},\n publisher={ACM}\n}\n",&cite)); 122 /* solve Ax = b by umfpack_*_wsolve */ 123 /* ----------------------------------*/ 124 125 if (!lu->Wi) { /* first time, allocate working space for wsolve */ 126 PetscCall(PetscMalloc1(A->rmap->n,&lu->Wi)); 127 PetscCall(PetscMalloc1(5*A->rmap->n,&lu->W)); 128 } 129 130 PetscCall(VecGetArrayRead(b,&ba)); 131 PetscCall(VecGetArray(x,&xa)); 132 #if defined(PETSC_USE_COMPLEX) 133 status = umfpack_UMF_wsolve(uflag,ai,aj,(PetscReal*)av,NULL,(PetscReal*)xa,NULL,(PetscReal*)ba,NULL,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 134 #else 135 status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 136 #endif 137 umfpack_UMF_report_info(lu->Control, lu->Info); 138 if (status < 0) { 139 umfpack_UMF_report_status(lu->Control, status); 140 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed"); 141 } 142 143 PetscCall(VecRestoreArrayRead(b,&ba)); 144 PetscCall(VecRestoreArray(x,&xa)); 145 PetscFunctionReturn(0); 146 } 147 148 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 149 { 150 PetscFunctionBegin; 151 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 152 PetscCall(MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat)); 153 PetscFunctionReturn(0); 154 } 155 156 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x) 157 { 158 PetscFunctionBegin; 159 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 160 PetscCall(MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A)); 161 PetscFunctionReturn(0); 162 } 163 164 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info) 165 { 166 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F)->data; 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscInt *ai = a->i,*aj=a->j,status; 169 PetscScalar *av = a->a; 170 171 PetscFunctionBegin; 172 if (!A->rmap->n) PetscFunctionReturn(0); 173 /* numeric factorization of A' */ 174 /* ----------------------------*/ 175 176 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) { 177 umfpack_UMF_free_numeric(&lu->Numeric); 178 } 179 #if defined(PETSC_USE_COMPLEX) 180 status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 181 #else 182 status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 183 #endif 184 if (status < 0) { 185 umfpack_UMF_report_status(lu->Control, status); 186 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_numeric failed"); 187 } 188 /* report numeric factorization of A' when Control[PRL] > 3 */ 189 (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control); 190 191 PetscCall(PetscObjectReference((PetscObject)A)); 192 PetscCall(MatDestroy(&lu->A)); 193 194 lu->A = A; 195 lu->flg = SAME_NONZERO_PATTERN; 196 lu->CleanUpUMFPACK = PETSC_TRUE; 197 F->ops->solve = MatSolve_UMFPACK; 198 F->ops->solvetranspose = MatSolveTranspose_UMFPACK; 199 PetscFunctionReturn(0); 200 } 201 202 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 203 { 204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 205 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F->data); 206 PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n,status,idx; 207 #if !defined(PETSC_USE_COMPLEX) 208 PetscScalar *av = a->a; 209 #endif 210 const PetscInt *ra; 211 const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC"}; 212 const char *scale[] ={"NONE","SUM","MAX"}; 213 PetscBool flg; 214 215 PetscFunctionBegin; 216 (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 217 if (!n) PetscFunctionReturn(0); 218 219 /* Set options to F */ 220 PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"UMFPACK Options","Mat"); 221 /* Control parameters used by reporting routiones */ 222 PetscCall(PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],NULL)); 223 224 /* Control parameters for symbolic factorization */ 225 PetscCall(PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg)); 226 if (flg) { 227 switch (idx) { 228 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 229 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 230 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 231 } 232 } 233 PetscCall(PetscOptionsEList("-mat_umfpack_ordering","Internal ordering method","None",UmfpackOrderingTypes,PETSC_STATIC_ARRAY_LENGTH(UmfpackOrderingTypes),UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]],&idx,&flg)); 234 if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx; 235 PetscCall(PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],NULL)); 236 PetscCall(PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],NULL)); 237 PetscCall(PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],NULL)); 238 PetscCall(PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],NULL)); 239 PetscCall(PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],NULL)); 240 PetscCall(PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],NULL)); 241 242 /* Control parameters used by numeric factorization */ 243 PetscCall(PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],NULL)); 244 PetscCall(PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],NULL)); 245 PetscCall(PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg)); 246 if (flg) { 247 switch (idx) { 248 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 249 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 250 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 251 } 252 } 253 PetscCall(PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL)); 254 PetscCall(PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL)); 255 PetscCall(PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],NULL)); 256 257 /* Control parameters used by solve */ 258 PetscCall(PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],NULL)); 259 PetscOptionsEnd(); 260 261 if (r) { 262 PetscCall(ISGetIndices(r,&ra)); 263 PetscCall(PetscMalloc1(m,&lu->perm_c)); 264 /* we cannot simply memcpy on 64 bit archs */ 265 for (i = 0; i < m; i++) lu->perm_c[i] = ra[i]; 266 PetscCall(ISRestoreIndices(r,&ra)); 267 } 268 269 /* print the control parameters */ 270 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control); 271 272 /* symbolic factorization of A' */ 273 /* ---------------------------------------------------------------------- */ 274 if (r) { /* use Petsc row ordering */ 275 #if !defined(PETSC_USE_COMPLEX) 276 status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 277 #else 278 status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 279 #endif 280 } else { /* use Umfpack col ordering */ 281 #if !defined(PETSC_USE_COMPLEX) 282 status = umfpack_UMF_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info); 283 #else 284 status = umfpack_UMF_symbolic(n,m,ai,aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info); 285 #endif 286 } 287 if (status < 0) { 288 umfpack_UMF_report_info(lu->Control, lu->Info); 289 umfpack_UMF_report_status(lu->Control, status); 290 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_symbolic failed"); 291 } 292 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 293 (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control); 294 295 lu->flg = DIFFERENT_NONZERO_PATTERN; 296 lu->CleanUpUMFPACK = PETSC_TRUE; 297 PetscFunctionReturn(0); 298 } 299 300 static PetscErrorCode MatView_Info_UMFPACK(Mat A,PetscViewer viewer) 301 { 302 Mat_UMFPACK *lu= (Mat_UMFPACK*)A->data; 303 304 PetscFunctionBegin; 305 /* check if matrix is UMFPACK type */ 306 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 307 308 PetscCall(PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n")); 309 /* Control parameters used by reporting routiones */ 310 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL])); 311 312 /* Control parameters used by symbolic factorization */ 313 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY])); 314 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL])); 315 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW])); 316 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE])); 317 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE])); 318 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ])); 319 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE])); 320 321 /* Control parameters used by numeric factorization */ 322 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE])); 323 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE])); 324 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE])); 325 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT])); 326 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL])); 327 328 /* Control parameters used by solve */ 329 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP])); 330 331 /* mat ordering */ 332 if (!lu->perm_c) { 333 PetscCall(PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n",UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]])); 334 } 335 PetscFunctionReturn(0); 336 } 337 338 static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer) 339 { 340 PetscBool iascii; 341 PetscViewerFormat format; 342 343 PetscFunctionBegin; 344 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 345 if (iascii) { 346 PetscCall(PetscViewerGetFormat(viewer,&format)); 347 if (format == PETSC_VIEWER_ASCII_INFO) { 348 PetscCall(MatView_Info_UMFPACK(A,viewer)); 349 } 350 } 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A,MatSolverType *type) 355 { 356 PetscFunctionBegin; 357 *type = MATSOLVERUMFPACK; 358 PetscFunctionReturn(0); 359 } 360 361 /*MC 362 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices 363 via the external package UMFPACK. 364 365 Use ./configure --download-suitesparse to install PETSc to use UMFPACK 366 367 Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver 368 369 Consult UMFPACK documentation for more information about the Control parameters 370 which correspond to the options database keys below. 371 372 Options Database Keys: 373 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE 374 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 375 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2 376 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 377 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW] 378 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE] 379 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 380 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE] 381 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ] 382 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE] 383 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 384 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE] 385 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX 386 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 387 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL] 388 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 389 390 Level: beginner 391 392 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 393 394 .seealso: `PCLU`, `MATSOLVERSUPERLU`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType` 395 M*/ 396 397 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F) 398 { 399 Mat B; 400 Mat_UMFPACK *lu; 401 PetscInt m=A->rmap->n,n=A->cmap->n; 402 403 PetscFunctionBegin; 404 /* Create the factorization matrix F */ 405 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 406 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n)); 407 PetscCall(PetscStrallocpy("umfpack",&((PetscObject)B)->type_name)); 408 PetscCall(MatSetUp(B)); 409 410 PetscCall(PetscNewLog(B,&lu)); 411 412 B->data = lu; 413 B->ops->getinfo = MatGetInfo_External; 414 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 415 B->ops->destroy = MatDestroy_UMFPACK; 416 B->ops->view = MatView_UMFPACK; 417 B->ops->matsolve = NULL; 418 419 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_umfpack)); 420 421 B->factortype = MAT_FACTOR_LU; 422 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 423 B->preallocated = PETSC_TRUE; 424 425 PetscCall(PetscFree(B->solvertype)); 426 PetscCall(PetscStrallocpy(MATSOLVERUMFPACK,&B->solvertype)); 427 B->canuseordering = PETSC_TRUE; 428 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU])); 429 430 /* initializations */ 431 /* ------------------------------------------------*/ 432 /* get the default control parameters */ 433 umfpack_UMF_defaults(lu->Control); 434 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */ 435 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 436 437 *F = B; 438 PetscFunctionReturn(0); 439 } 440 441 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 442 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 443 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*); 444 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_spqr(Mat,MatFactorType,Mat*); 445 446 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void) 447 { 448 PetscFunctionBegin; 449 PetscCall(MatSolverTypeRegister(MATSOLVERUMFPACK,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_umfpack)); 450 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_cholmod)); 451 PetscCall(MatSolverTypeRegister(MATSOLVERCHOLMOD,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_cholmod)); 452 PetscCall(MatSolverTypeRegister(MATSOLVERKLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_klu)); 453 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR,MATSEQAIJ, MAT_FACTOR_QR,MatGetFactor_seqaij_spqr)); 454 if (!PetscDefined(USE_COMPLEX)) { 455 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR,MATNORMAL, MAT_FACTOR_QR,MatGetFactor_seqaij_spqr)); 456 } 457 PetscCall(MatSolverTypeRegister(MATSOLVERSPQR,MATNORMALHERMITIAN, MAT_FACTOR_QR,MatGetFactor_seqaij_spqr)); 458 PetscFunctionReturn(0); 459 } 460