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 PetscBool PetscMatOrdering; 88 89 /* Flag to clean up UMFPACK objects during Destroy */ 90 PetscBool CleanUpUMFPACK; 91 } Mat_UMFPACK; 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatDestroy_UMFPACK" 95 static PetscErrorCode MatDestroy_UMFPACK(Mat A) 96 { 97 PetscErrorCode ierr; 98 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 99 100 PetscFunctionBegin; 101 if (lu && lu->CleanUpUMFPACK) { 102 umfpack_UMF_free_symbolic(&lu->Symbolic); 103 umfpack_UMF_free_numeric(&lu->Numeric); 104 ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 105 ierr = PetscFree(lu->W);CHKERRQ(ierr); 106 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 107 } 108 ierr = MatDestroy(&lu->A);CHKERRQ(ierr); 109 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 110 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "MatSolve_UMFPACK_Private" 116 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag) 117 { 118 Mat_UMFPACK *lu = (Mat_UMFPACK*)A->spptr; 119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)lu->A->data; 120 PetscScalar *av = a->a,*xa; 121 const PetscScalar *ba; 122 PetscErrorCode ierr; 123 PetscInt *ai = a->i,*aj = a->j,status; 124 125 PetscFunctionBegin; 126 /* solve Ax = b by umfpack_*_wsolve */ 127 /* ----------------------------------*/ 128 129 if (!lu->Wi) { /* first time, allocate working space for wsolve */ 130 ierr = PetscMalloc1(A->rmap->n,&lu->Wi);CHKERRQ(ierr); 131 ierr = PetscMalloc1(5*A->rmap->n,&lu->W);CHKERRQ(ierr); 132 } 133 134 ierr = VecGetArrayRead(b,&ba); 135 ierr = VecGetArray(x,&xa); 136 #if defined(PETSC_USE_COMPLEX) 137 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); 138 #else 139 status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 140 #endif 141 umfpack_UMF_report_info(lu->Control, lu->Info); 142 if (status < 0) { 143 umfpack_UMF_report_status(lu->Control, status); 144 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed"); 145 } 146 147 ierr = VecRestoreArrayRead(b,&ba);CHKERRQ(ierr); 148 ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "MatSolve_UMFPACK" 154 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 155 { 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 160 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);CHKERRQ(ierr); 161 PetscFunctionReturn(0); 162 } 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "MatSolveTranspose_UMFPACK" 166 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x) 167 { 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */ 172 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 #undef __FUNCT__ 177 #define __FUNCT__ "MatLUFactorNumeric_UMFPACK" 178 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info) 179 { 180 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F)->spptr; 181 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 182 PetscInt *ai = a->i,*aj=a->j,status; 183 PetscScalar *av = a->a; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 /* numeric factorization of A' */ 188 /* ----------------------------*/ 189 190 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) { 191 umfpack_UMF_free_numeric(&lu->Numeric); 192 } 193 #if defined(PETSC_USE_COMPLEX) 194 status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 195 #else 196 status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 197 #endif 198 if (status < 0) { 199 umfpack_UMF_report_status(lu->Control, status); 200 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_numeric failed"); 201 } 202 /* report numeric factorization of A' when Control[PRL] > 3 */ 203 (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control); 204 205 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 206 ierr = MatDestroy(&lu->A);CHKERRQ(ierr); 207 208 lu->A = A; 209 lu->flg = SAME_NONZERO_PATTERN; 210 lu->CleanUpUMFPACK = PETSC_TRUE; 211 F->ops->solve = MatSolve_UMFPACK; 212 F->ops->solvetranspose = MatSolveTranspose_UMFPACK; 213 PetscFunctionReturn(0); 214 } 215 216 /* 217 Note the r permutation is ignored 218 */ 219 #undef __FUNCT__ 220 #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK" 221 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 222 { 223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 224 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F->spptr); 225 PetscErrorCode ierr; 226 PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n; 227 #if !defined(PETSC_USE_COMPLEX) 228 PetscScalar *av = a->a; 229 #endif 230 const PetscInt *ra; 231 PetscInt status; 232 233 PetscFunctionBegin; 234 if (lu->PetscMatOrdering) { 235 ierr = ISGetIndices(r,&ra);CHKERRQ(ierr); 236 ierr = PetscMalloc1(m,&lu->perm_c);CHKERRQ(ierr); 237 /* we cannot simply memcpy on 64 bit archs */ 238 for (i = 0; i < m; i++) lu->perm_c[i] = ra[i]; 239 ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr); 240 } 241 242 /* print the control parameters */ 243 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control); 244 245 /* symbolic factorization of A' */ 246 /* ---------------------------------------------------------------------- */ 247 if (lu->PetscMatOrdering) { /* use Petsc row ordering */ 248 #if !defined(PETSC_USE_COMPLEX) 249 status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 250 #else 251 status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 252 #endif 253 } else { /* use Umfpack col ordering */ 254 #if !defined(PETSC_USE_COMPLEX) 255 status = umfpack_UMF_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info); 256 #else 257 status = umfpack_UMF_symbolic(n,m,ai,aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info); 258 #endif 259 } 260 if (status < 0) { 261 umfpack_UMF_report_info(lu->Control, lu->Info); 262 umfpack_UMF_report_status(lu->Control, status); 263 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_symbolic failed"); 264 } 265 /* report sumbolic factorization of A' when Control[PRL] > 3 */ 266 (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control); 267 268 lu->flg = DIFFERENT_NONZERO_PATTERN; 269 lu->CleanUpUMFPACK = PETSC_TRUE; 270 (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNCT__ 275 #define __FUNCT__ "MatFactorInfo_UMFPACK" 276 static PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 277 { 278 Mat_UMFPACK *lu= (Mat_UMFPACK*)A->spptr; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 /* check if matrix is UMFPACK type */ 283 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 284 285 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 286 /* Control parameters used by reporting routiones */ 287 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 288 289 /* Control parameters used by symbolic factorization */ 290 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr); 291 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 292 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 293 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr); 294 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 295 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr); 296 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr); 297 298 /* Control parameters used by numeric factorization */ 299 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 300 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr); 301 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr); 302 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 303 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr); 304 305 /* Control parameters used by solve */ 306 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 307 308 /* mat ordering */ 309 if (!lu->PetscMatOrdering) { 310 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ORDERING]: %s (not using the PETSc ordering)\n",UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]);CHKERRQ(ierr); 311 } 312 PetscFunctionReturn(0); 313 } 314 315 #undef __FUNCT__ 316 #define __FUNCT__ "MatView_UMFPACK" 317 static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer) 318 { 319 PetscErrorCode ierr; 320 PetscBool iascii; 321 PetscViewerFormat format; 322 323 PetscFunctionBegin; 324 ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 325 326 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 327 if (iascii) { 328 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 329 if (format == PETSC_VIEWER_ASCII_INFO) { 330 ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 331 } 332 } 333 PetscFunctionReturn(0); 334 } 335 336 #undef __FUNCT__ 337 #define __FUNCT__ "MatFactorGetSolverPackage_seqaij_umfpack" 338 PetscErrorCode MatFactorGetSolverPackage_seqaij_umfpack(Mat A,const MatSolverPackage *type) 339 { 340 PetscFunctionBegin; 341 *type = MATSOLVERUMFPACK; 342 PetscFunctionReturn(0); 343 } 344 345 346 /*MC 347 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices 348 via the external package UMFPACK. 349 350 ./configure --download-suitesparse to install PETSc to use UMFPACK 351 352 Consult UMFPACK documentation for more information about the Control parameters 353 which correspond to the options database keys below. 354 355 Options Database Keys: 356 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE 357 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 358 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2 359 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 360 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW] 361 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE] 362 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 363 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE] 364 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ] 365 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE] 366 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 367 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE] 368 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX 369 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 370 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL] 371 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 372 373 Level: beginner 374 375 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 376 377 .seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage 378 M*/ 379 380 #undef __FUNCT__ 381 #define __FUNCT__ "MatGetFactor_seqaij_umfpack" 382 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F) 383 { 384 Mat B; 385 Mat_UMFPACK *lu; 386 PetscErrorCode ierr; 387 PetscInt m=A->rmap->n,n=A->cmap->n,idx; 388 389 const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC"}; 390 const char *scale[] ={"NONE","SUM","MAX"}; 391 PetscBool flg; 392 393 PetscFunctionBegin; 394 /* Create the factorization matrix F */ 395 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 396 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 397 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 398 ierr = MatSeqAIJSetPreallocation(B,0,NULL);CHKERRQ(ierr); 399 ierr = PetscNewLog(B,&lu);CHKERRQ(ierr); 400 401 B->spptr = lu; 402 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 403 B->ops->destroy = MatDestroy_UMFPACK; 404 B->ops->view = MatView_UMFPACK; 405 406 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_umfpack);CHKERRQ(ierr); 407 408 B->factortype = MAT_FACTOR_LU; 409 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 410 B->preallocated = PETSC_TRUE; 411 412 /* initializations */ 413 /* ------------------------------------------------*/ 414 /* get the default control parameters */ 415 umfpack_UMF_defaults(lu->Control); 416 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */ 417 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 418 419 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 420 /* Control parameters used by reporting routiones */ 421 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],NULL);CHKERRQ(ierr); 422 423 /* Control parameters for symbolic factorization */ 424 ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg);CHKERRQ(ierr); 425 if (flg) { 426 switch (idx) { 427 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 428 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 429 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 430 } 431 } 432 ierr = PetscOptionsEList("-mat_umfpack_ordering","Internal ordering method","None",UmfpackOrderingTypes,sizeof(UmfpackOrderingTypes)/sizeof(UmfpackOrderingTypes[0]),UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]],&idx,&flg);CHKERRQ(ierr); 433 if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx; 434 ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],NULL);CHKERRQ(ierr); 435 ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],NULL);CHKERRQ(ierr); 436 ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],NULL);CHKERRQ(ierr); 437 ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],NULL);CHKERRQ(ierr); 438 ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],NULL);CHKERRQ(ierr); 439 ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],NULL);CHKERRQ(ierr); 440 441 /* Control parameters used by numeric factorization */ 442 ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],NULL);CHKERRQ(ierr); 443 ierr = 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);CHKERRQ(ierr); 444 ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr); 445 if (flg) { 446 switch (idx) { 447 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 448 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 449 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 450 } 451 } 452 ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);CHKERRQ(ierr); 453 ierr = PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);CHKERRQ(ierr); 454 ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],NULL);CHKERRQ(ierr); 455 456 /* Control parameters used by solve */ 457 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],NULL);CHKERRQ(ierr); 458 459 /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */ 460 ierr = PetscOptionsHasName(NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOrdering);CHKERRQ(ierr); 461 PetscOptionsEnd(); 462 *F = B; 463 PetscFunctionReturn(0); 464 } 465 466 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 467 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 468 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*); 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "MatSolverPackageRegister_SuiteSparse" 472 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_SuiteSparse(void) 473 { 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 ierr = MatSolverPackageRegister(MATSOLVERUMFPACK,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 478 ierr = MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 479 ierr = MatSolverPackageRegister(MATSOLVERCHOLMOD,MATSEQSBAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 480 ierr = MatSolverPackageRegister(MATSOLVERKLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_klu);CHKERRQ(ierr); 481 PetscFunctionReturn(0); 482 } 483