1be1d678aSKris Buschelman #define PETSCMAT_DLL 21677d0b8SKris Buschelman 31677d0b8SKris Buschelman /* 42d4e2982SHong Zhang Provides an interface to the UMFPACKv5.1 sparse solver 52d4e2982SHong Zhang 62d4e2982SHong Zhang This interface uses the "UF_long version" of the UMFPACK API 72d4e2982SHong Zhang (*_dl_* and *_zl_* routines, instead of *_di_* and *_zi_* routines) 82d4e2982SHong Zhang so that UMFPACK can address more than 2Gb of memory on 64 bit 92d4e2982SHong Zhang machines. 102d4e2982SHong Zhang 112d4e2982SHong Zhang If sizeof(UF_long) == 32 the interface only allocates the memory 122d4e2982SHong Zhang necessary for UMFPACK's working arrays (W, Wi and possibly 132d4e2982SHong Zhang perm_c). If sizeof(UF_long) == 64, in addition to allocating the 142d4e2982SHong Zhang working arrays, the interface also has to re-allocate the matrix 152d4e2982SHong Zhang index arrays (ai and aj, which must be stored as UF_long). 162d4e2982SHong Zhang 172d4e2982SHong Zhang The interface is implemented for both real and complex 182d4e2982SHong Zhang arithmetic. Complex numbers are assumed to be in packed format, 192d4e2982SHong Zhang which requires UMFPACK >= 4.4. 202d4e2982SHong Zhang 212d4e2982SHong Zhang We thank Christophe Geuzaine <geuzaine@acm.caltech.edu> for upgrading this interface to the UMFPACKv5.1 221677d0b8SKris Buschelman */ 231677d0b8SKris Buschelman 241677d0b8SKris Buschelman #include "src/mat/impls/aij/seq/aij.h" 251677d0b8SKris Buschelman 261677d0b8SKris Buschelman EXTERN_C_BEGIN 271677d0b8SKris Buschelman #include "umfpack.h" 281677d0b8SKris Buschelman EXTERN_C_END 291677d0b8SKris Buschelman 301677d0b8SKris Buschelman typedef struct { 311677d0b8SKris Buschelman void *Symbolic, *Numeric; 321677d0b8SKris Buschelman double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W; 332d4e2982SHong Zhang UF_long *Wi,*ai,*aj,*perm_c; 341677d0b8SKris Buschelman PetscScalar *av; 351677d0b8SKris Buschelman MatStructure flg; 361677d0b8SKris Buschelman PetscTruth PetscMatOdering; 371677d0b8SKris Buschelman 381677d0b8SKris Buschelman /* Flag to clean up UMFPACK objects during Destroy */ 391677d0b8SKris Buschelman PetscTruth CleanUpUMFPACK; 40f0c56d0fSKris Buschelman } Mat_UMFPACK; 41f0c56d0fSKris Buschelman 421677d0b8SKris Buschelman #undef __FUNCT__ 43f0c56d0fSKris Buschelman #define __FUNCT__ "MatDestroy_UMFPACK" 44dfbe8321SBarry Smith PetscErrorCode MatDestroy_UMFPACK(Mat A) 45dfbe8321SBarry Smith { 46dfbe8321SBarry Smith PetscErrorCode ierr; 47f0c56d0fSKris Buschelman Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr; 481677d0b8SKris Buschelman 491677d0b8SKris Buschelman PetscFunctionBegin; 50fb731535SKris Buschelman if (lu->CleanUpUMFPACK) { 512d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 522d4e2982SHong Zhang umfpack_zl_free_symbolic(&lu->Symbolic); 532d4e2982SHong Zhang umfpack_zl_free_numeric(&lu->Numeric); 542d4e2982SHong Zhang #else 552d4e2982SHong Zhang umfpack_dl_free_symbolic(&lu->Symbolic); 562d4e2982SHong Zhang umfpack_dl_free_numeric(&lu->Numeric); 572d4e2982SHong Zhang #endif 581677d0b8SKris Buschelman ierr = PetscFree(lu->Wi);CHKERRQ(ierr); 591677d0b8SKris Buschelman ierr = PetscFree(lu->W);CHKERRQ(ierr); 602d4e2982SHong Zhang if(sizeof(UF_long) != sizeof(int)){ 612d4e2982SHong Zhang ierr = PetscFree(lu->ai);CHKERRQ(ierr); 622d4e2982SHong Zhang ierr = PetscFree(lu->aj);CHKERRQ(ierr); 632d4e2982SHong Zhang } 641677d0b8SKris Buschelman if (lu->PetscMatOdering) { 651677d0b8SKris Buschelman ierr = PetscFree(lu->perm_c);CHKERRQ(ierr); 661677d0b8SKris Buschelman } 671677d0b8SKris Buschelman } 68b24902e0SBarry Smith ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 691677d0b8SKris Buschelman PetscFunctionReturn(0); 701677d0b8SKris Buschelman } 711677d0b8SKris Buschelman 721677d0b8SKris Buschelman #undef __FUNCT__ 73f0c56d0fSKris Buschelman #define __FUNCT__ "MatSolve_UMFPACK" 746849ba73SBarry Smith PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x) 756849ba73SBarry Smith { 76f0c56d0fSKris Buschelman Mat_UMFPACK *lu = (Mat_UMFPACK*)A->spptr; 771677d0b8SKris Buschelman PetscScalar *av=lu->av,*ba,*xa; 78dfbe8321SBarry Smith PetscErrorCode ierr; 792d4e2982SHong Zhang UF_long *ai=lu->ai,*aj=lu->aj,status; 801677d0b8SKris Buschelman 811677d0b8SKris Buschelman PetscFunctionBegin; 822d4e2982SHong Zhang /* solve Ax = b by umfpack_*_wsolve */ 831677d0b8SKris Buschelman /* ----------------------------------*/ 842d4e2982SHong Zhang 852d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 862d4e2982SHong Zhang ierr = VecConjugate(b); 872d4e2982SHong Zhang #endif 882d4e2982SHong Zhang 891677d0b8SKris Buschelman ierr = VecGetArray(b,&ba); 901677d0b8SKris Buschelman ierr = VecGetArray(x,&xa); 911677d0b8SKris Buschelman 922d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 932d4e2982SHong Zhang status = umfpack_zl_wsolve(UMFPACK_At,ai,aj,(double*)av,NULL,(double*)xa,NULL,(double*)ba,NULL, 942d4e2982SHong Zhang lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 952d4e2982SHong Zhang umfpack_zl_report_info(lu->Control, lu->Info); 961677d0b8SKris Buschelman if (status < 0){ 972d4e2982SHong Zhang umfpack_zl_report_status(lu->Control, status); 982d4e2982SHong Zhang SETERRQ(PETSC_ERR_LIB,"umfpack_zl_wsolve failed"); 991677d0b8SKris Buschelman } 1002d4e2982SHong Zhang #else 1012d4e2982SHong Zhang status = umfpack_dl_wsolve(UMFPACK_At,ai,aj,av,xa,ba, 1022d4e2982SHong Zhang lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W); 1032d4e2982SHong Zhang umfpack_dl_report_info(lu->Control, lu->Info); 1042d4e2982SHong Zhang if (status < 0){ 1052d4e2982SHong Zhang umfpack_dl_report_status(lu->Control, status); 1062d4e2982SHong Zhang SETERRQ(PETSC_ERR_LIB,"umfpack_dl_wsolve failed"); 1072d4e2982SHong Zhang } 1082d4e2982SHong Zhang #endif 1091677d0b8SKris Buschelman 1101677d0b8SKris Buschelman ierr = VecRestoreArray(b,&ba); 1111677d0b8SKris Buschelman ierr = VecRestoreArray(x,&xa); 1122a325a84SHong Zhang 1132d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 1142d4e2982SHong Zhang ierr = VecConjugate(b); 1155b2a9f60SHong Zhang ierr = VecConjugate(x); 1162d4e2982SHong Zhang #endif 1172d4e2982SHong Zhang 1181677d0b8SKris Buschelman PetscFunctionReturn(0); 1191677d0b8SKris Buschelman } 1201677d0b8SKris Buschelman 1211677d0b8SKris Buschelman #undef __FUNCT__ 122f0c56d0fSKris Buschelman #define __FUNCT__ "MatLUFactorNumeric_UMFPACK" 123af281ebdSHong Zhang PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat A,MatFactorInfo *info,Mat *F) 1246849ba73SBarry Smith { 125f0c56d0fSKris Buschelman Mat_UMFPACK *lu=(Mat_UMFPACK*)(*F)->spptr; 1266849ba73SBarry Smith PetscErrorCode ierr; 1272d4e2982SHong Zhang UF_long *ai=lu->ai,*aj=lu->aj,m=A->rmap.n,status; 1281677d0b8SKris Buschelman PetscScalar *av=lu->av; 1291677d0b8SKris Buschelman 1301677d0b8SKris Buschelman PetscFunctionBegin; 1311677d0b8SKris Buschelman /* numeric factorization of A' */ 1321677d0b8SKris Buschelman /* ----------------------------*/ 1332d4e2982SHong Zhang 1342d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 1352a325a84SHong Zhang if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){ 1362d4e2982SHong Zhang umfpack_zl_free_numeric(&lu->Numeric); 1372a325a84SHong Zhang } 1382d4e2982SHong Zhang status = umfpack_zl_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 1392d4e2982SHong Zhang if (status < 0) SETERRQ(PETSC_ERR_LIB,"umfpack_zl_numeric failed"); 1401677d0b8SKris Buschelman /* report numeric factorization of A' when Control[PRL] > 3 */ 1412d4e2982SHong Zhang (void) umfpack_zl_report_numeric(lu->Numeric, lu->Control); 1422d4e2982SHong Zhang #else 1432d4e2982SHong Zhang if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){ 1442d4e2982SHong Zhang umfpack_dl_free_numeric(&lu->Numeric); 1452d4e2982SHong Zhang } 1462d4e2982SHong Zhang status = umfpack_dl_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info); 1472d4e2982SHong Zhang if (status < 0) SETERRQ(PETSC_ERR_LIB,"umfpack_dl_numeric failed"); 1482d4e2982SHong Zhang /* report numeric factorization of A' when Control[PRL] > 3 */ 1492d4e2982SHong Zhang (void) umfpack_dl_report_numeric(lu->Numeric, lu->Control); 1502d4e2982SHong Zhang #endif 1511677d0b8SKris Buschelman 1521677d0b8SKris Buschelman if (lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */ 1531677d0b8SKris Buschelman /* allocate working space to be used by Solve */ 1542d4e2982SHong Zhang ierr = PetscMalloc(m * sizeof(UF_long), &lu->Wi);CHKERRQ(ierr); 1552d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 1562d4e2982SHong Zhang ierr = PetscMalloc(2*5*m * sizeof(double), &lu->W);CHKERRQ(ierr); 1572d4e2982SHong Zhang #else 1581677d0b8SKris Buschelman ierr = PetscMalloc(5*m * sizeof(double), &lu->W);CHKERRQ(ierr); 1592d4e2982SHong Zhang #endif 1601677d0b8SKris Buschelman } 1612d4e2982SHong Zhang 1622a325a84SHong Zhang lu->flg = SAME_NONZERO_PATTERN; 1632a325a84SHong Zhang lu->CleanUpUMFPACK = PETSC_TRUE; 1641677d0b8SKris Buschelman PetscFunctionReturn(0); 1651677d0b8SKris Buschelman } 1661677d0b8SKris Buschelman 1671677d0b8SKris Buschelman /* 1681677d0b8SKris Buschelman Note the r permutation is ignored 1691677d0b8SKris Buschelman */ 1701677d0b8SKris Buschelman #undef __FUNCT__ 171f0c56d0fSKris Buschelman #define __FUNCT__ "MatLUFactorSymbolic_UMFPACK" 1726849ba73SBarry Smith PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 1736849ba73SBarry Smith { 174e1b4c3a1SKris Buschelman Mat B; 1751677d0b8SKris Buschelman Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 176b24902e0SBarry Smith Mat_UMFPACK *lu = (Mat_UMFPACK*)((*F)->spptr); 177dfbe8321SBarry Smith PetscErrorCode ierr; 1782d4e2982SHong Zhang int i,m=A->rmap.n,n=A->cmap.n,*ra,idx; 1792d4e2982SHong Zhang UF_long status; 1801677d0b8SKris Buschelman PetscScalar *av=mat->a; 1811677d0b8SKris Buschelman 1821677d0b8SKris Buschelman PetscFunctionBegin; 1831677d0b8SKris Buschelman if (lu->PetscMatOdering) { 1840f6efc10SHong Zhang ierr = ISGetIndices(r,&ra);CHKERRQ(ierr); 1852d4e2982SHong Zhang ierr = PetscMalloc(m*sizeof(UF_long),&lu->perm_c);CHKERRQ(ierr); 1862d4e2982SHong Zhang /* we cannot simply memcpy on 64 bit archs */ 1872d4e2982SHong Zhang for(i = 0; i < m; i++) lu->perm_c[i] = ra[i]; 1880f6efc10SHong Zhang ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr); 1891677d0b8SKris Buschelman } 1901677d0b8SKris Buschelman 1912d4e2982SHong Zhang if(sizeof(UF_long) != sizeof(int)){ 1922d4e2982SHong Zhang /* we cannot directly use mat->i and mat->j on 64 bit archs */ 1932d4e2982SHong Zhang ierr = PetscMalloc((m+1)*sizeof(UF_long),&lu->ai);CHKERRQ(ierr); 1942d4e2982SHong Zhang ierr = PetscMalloc(mat->nz*sizeof(UF_long),&lu->aj);CHKERRQ(ierr); 1952d4e2982SHong Zhang for(i = 0; i < m + 1; i++) lu->ai[i] = mat->i[i]; 1962d4e2982SHong Zhang for(i = 0; i < mat->nz; i++) lu->aj[i] = mat->j[i]; 1972d4e2982SHong Zhang } 1982d4e2982SHong Zhang else{ 1992d4e2982SHong Zhang lu->ai = (UF_long*)mat->i; 2002d4e2982SHong Zhang lu->aj = (UF_long*)mat->j; 2012d4e2982SHong Zhang } 2022d4e2982SHong Zhang 2031677d0b8SKris Buschelman /* print the control parameters */ 2042d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 2052d4e2982SHong Zhang if(lu->Control[UMFPACK_PRL] > 1) umfpack_zl_report_control(lu->Control); 2062d4e2982SHong Zhang #else 2072d4e2982SHong Zhang if(lu->Control[UMFPACK_PRL] > 1) umfpack_dl_report_control(lu->Control); 2082d4e2982SHong Zhang #endif 2091677d0b8SKris Buschelman 2101677d0b8SKris Buschelman /* symbolic factorization of A' */ 2111677d0b8SKris Buschelman /* ---------------------------------------------------------------------- */ 2122d4e2982SHong Zhang #if defined(PETSC_USE_COMPLEX) 2130f6efc10SHong Zhang if (lu->PetscMatOdering) { /* use Petsc row ordering */ 2142d4e2982SHong Zhang status = umfpack_zl_qsymbolic(n,m,lu->ai,lu->aj,(double*)av,NULL, 2152d4e2982SHong Zhang lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 2160f6efc10SHong Zhang } else { /* use Umfpack col ordering */ 2172d4e2982SHong Zhang status = umfpack_zl_symbolic(n,m,lu->ai,lu->aj,(double*)av,NULL, 2182d4e2982SHong Zhang &lu->Symbolic,lu->Control,lu->Info); 2190f6efc10SHong Zhang } 2201677d0b8SKris Buschelman if (status < 0){ 2212d4e2982SHong Zhang umfpack_zl_report_info(lu->Control, lu->Info); 2222d4e2982SHong Zhang umfpack_zl_report_status(lu->Control, status); 2232d4e2982SHong Zhang SETERRQ(PETSC_ERR_LIB,"umfpack_dl_symbolic failed"); 2241677d0b8SKris Buschelman } 2251677d0b8SKris Buschelman /* report sumbolic factorization of A' when Control[PRL] > 3 */ 2262d4e2982SHong Zhang (void) umfpack_zl_report_symbolic(lu->Symbolic, lu->Control); 2272d4e2982SHong Zhang #else 2282d4e2982SHong Zhang if (lu->PetscMatOdering) { /* use Petsc row ordering */ 2292d4e2982SHong Zhang status = umfpack_dl_qsymbolic(n,m,lu->ai,lu->aj,av, 2302d4e2982SHong Zhang lu->perm_c,&lu->Symbolic,lu->Control,lu->Info); 2312d4e2982SHong Zhang } else { /* use Umfpack col ordering */ 2322d4e2982SHong Zhang status = umfpack_dl_symbolic(n,m,lu->ai,lu->aj,av, 2332d4e2982SHong Zhang &lu->Symbolic,lu->Control,lu->Info); 2342d4e2982SHong Zhang } 2352d4e2982SHong Zhang if (status < 0){ 2362d4e2982SHong Zhang umfpack_dl_report_info(lu->Control, lu->Info); 2372d4e2982SHong Zhang umfpack_dl_report_status(lu->Control, status); 2382d4e2982SHong Zhang SETERRQ(PETSC_ERR_LIB,"umfpack_dl_symbolic failed"); 2392d4e2982SHong Zhang } 2402d4e2982SHong Zhang /* report sumbolic factorization of A' when Control[PRL] > 3 */ 2412d4e2982SHong Zhang (void) umfpack_dl_report_symbolic(lu->Symbolic, lu->Control); 2422d4e2982SHong Zhang #endif 2431677d0b8SKris Buschelman 2441677d0b8SKris Buschelman lu->flg = DIFFERENT_NONZERO_PATTERN; 2451677d0b8SKris Buschelman lu->av = av; 246e1b4c3a1SKris Buschelman lu->CleanUpUMFPACK = PETSC_TRUE; 2471677d0b8SKris Buschelman PetscFunctionReturn(0); 2481677d0b8SKris Buschelman } 2491677d0b8SKris Buschelman 2501677d0b8SKris Buschelman #undef __FUNCT__ 251b24902e0SBarry Smith #define __FUNCT__ "MatGetFactor_seqaij_umfpack" 252*5c9eb25fSBarry Smith PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F) 2536849ba73SBarry Smith { 254b24902e0SBarry Smith Mat B; 255b24902e0SBarry Smith Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data; 256b24902e0SBarry Smith Mat_UMFPACK *lu; 257dfbe8321SBarry Smith PetscErrorCode ierr; 258b24902e0SBarry Smith int i,m=A->rmap.n,n=A->cmap.n,*ra,idx; 259b24902e0SBarry Smith UF_long status; 260b24902e0SBarry Smith 261b24902e0SBarry Smith PetscScalar *av=mat->a; 262b24902e0SBarry Smith const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"}, 263b24902e0SBarry Smith *scale[]={"NONE","SUM","MAX"}; 264b24902e0SBarry Smith PetscTruth flg; 265d844382dSKris Buschelman 2661677d0b8SKris Buschelman PetscFunctionBegin; 267b24902e0SBarry Smith /* Create the factorization matrix F */ 268b24902e0SBarry Smith ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 269b24902e0SBarry Smith ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 270b24902e0SBarry Smith ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 271b24902e0SBarry Smith ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 272b24902e0SBarry Smith ierr = PetscNewLog(B,Mat_UMFPACK,&lu);CHKERRQ(ierr); 273b24902e0SBarry Smith B->spptr = lu; 274b24902e0SBarry Smith B->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK; 275b24902e0SBarry Smith B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK; 276b24902e0SBarry Smith B->ops->solve = MatSolve_UMFPACK; 277b24902e0SBarry Smith B->ops->destroy = MatDestroy_UMFPACK; 278*5c9eb25fSBarry Smith B->factor = MAT_FACTOR_LU; 279b24902e0SBarry Smith B->assembled = PETSC_TRUE; /* required by -ksp_view */ 280*5c9eb25fSBarry Smith B->preallocated = PETSC_TRUE; 281b24902e0SBarry Smith 282b24902e0SBarry Smith /* initializations */ 283b24902e0SBarry Smith /* ------------------------------------------------*/ 284b24902e0SBarry Smith /* get the default control parameters */ 285b24902e0SBarry Smith #if defined(PETSC_USE_COMPLEX) 286b24902e0SBarry Smith umfpack_zl_defaults(lu->Control); 287b24902e0SBarry Smith #else 288b24902e0SBarry Smith umfpack_dl_defaults(lu->Control); 289b24902e0SBarry Smith #endif 290b24902e0SBarry Smith lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */ 291b24902e0SBarry Smith lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */ 292b24902e0SBarry Smith 293b24902e0SBarry Smith ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr); 294b24902e0SBarry Smith /* Control parameters used by reporting routiones */ 295b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr); 296b24902e0SBarry Smith 297b24902e0SBarry Smith /* Control parameters for symbolic factorization */ 298b24902e0SBarry Smith ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);CHKERRQ(ierr); 299b24902e0SBarry Smith if (flg) { 300b24902e0SBarry Smith switch (idx){ 301b24902e0SBarry Smith case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break; 302b24902e0SBarry Smith case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break; 303b24902e0SBarry Smith case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break; 304b24902e0SBarry Smith case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break; 305b24902e0SBarry Smith } 306b24902e0SBarry Smith } 307b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);CHKERRQ(ierr); 308b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);CHKERRQ(ierr); 309b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);CHKERRQ(ierr); 310b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);CHKERRQ(ierr); 311b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_2by2_tolerance","Control[UMFPACK_2BY2_TOLERANCE]","None",lu->Control[UMFPACK_2BY2_TOLERANCE],&lu->Control[UMFPACK_2BY2_TOLERANCE],PETSC_NULL);CHKERRQ(ierr); 312b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr); 313b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr); 314b24902e0SBarry Smith 315b24902e0SBarry Smith /* Control parameters used by numeric factorization */ 316b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr); 317b24902e0SBarry Smith 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],PETSC_NULL);CHKERRQ(ierr); 318b24902e0SBarry Smith ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr); 319b24902e0SBarry Smith if (flg) { 320b24902e0SBarry Smith switch (idx){ 321b24902e0SBarry Smith case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break; 322b24902e0SBarry Smith case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break; 323b24902e0SBarry Smith case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break; 324b24902e0SBarry Smith } 325b24902e0SBarry Smith } 326b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr); 327b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);CHKERRQ(ierr); 328b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);CHKERRQ(ierr); 329b24902e0SBarry Smith 330b24902e0SBarry Smith /* Control parameters used by solve */ 331b24902e0SBarry Smith ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);CHKERRQ(ierr); 332b24902e0SBarry Smith 333b24902e0SBarry Smith /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */ 334b24902e0SBarry Smith ierr = PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);CHKERRQ(ierr); 335b24902e0SBarry Smith PetscOptionsEnd(); 336b24902e0SBarry Smith *F = B; 3371677d0b8SKris Buschelman PetscFunctionReturn(0); 3381677d0b8SKris Buschelman } 3391677d0b8SKris Buschelman 3401677d0b8SKris Buschelman #undef __FUNCT__ 341f0c56d0fSKris Buschelman #define __FUNCT__ "MatFactorInfo_UMFPACK" 3426849ba73SBarry Smith PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer) 3436849ba73SBarry Smith { 344f0c56d0fSKris Buschelman Mat_UMFPACK *lu= (Mat_UMFPACK*)A->spptr; 345dfbe8321SBarry Smith PetscErrorCode ierr; 346f0c56d0fSKris Buschelman 3471677d0b8SKris Buschelman PetscFunctionBegin; 3481677d0b8SKris Buschelman /* check if matrix is UMFPACK type */ 349f0c56d0fSKris Buschelman if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0); 3501677d0b8SKris Buschelman 3511677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr); 3521677d0b8SKris Buschelman /* Control parameters used by reporting routiones */ 3531677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr); 3541677d0b8SKris Buschelman 3551677d0b8SKris Buschelman /* Control parameters used by symbolic factorization */ 3560f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr); 3571677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr); 3581677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr); 3590f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr); 3601677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr); 3610f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);CHKERRQ(ierr); 3620f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr); 3630f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr); 3641677d0b8SKris Buschelman 3651677d0b8SKris Buschelman /* Control parameters used by numeric factorization */ 3661677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr); 3670f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr); 3680f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr); 3691677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr); 3700f6efc10SHong Zhang ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr); 3711677d0b8SKris Buschelman 3721677d0b8SKris Buschelman /* Control parameters used by solve */ 3731677d0b8SKris Buschelman ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr); 3741677d0b8SKris Buschelman 3751677d0b8SKris Buschelman /* mat ordering */ 3761677d0b8SKris Buschelman if(!lu->PetscMatOdering) ierr = PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");CHKERRQ(ierr); 3771677d0b8SKris Buschelman 3781677d0b8SKris Buschelman PetscFunctionReturn(0); 3791677d0b8SKris Buschelman } 3801677d0b8SKris Buschelman 381d844382dSKris Buschelman #undef __FUNCT__ 382f0c56d0fSKris Buschelman #define __FUNCT__ "MatView_UMFPACK" 3836849ba73SBarry Smith PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer) 3846849ba73SBarry Smith { 385dfbe8321SBarry Smith PetscErrorCode ierr; 38632077d6dSBarry Smith PetscTruth iascii; 387d844382dSKris Buschelman PetscViewerFormat format; 388f0c56d0fSKris Buschelman Mat_UMFPACK *lu=(Mat_UMFPACK*)(A->spptr); 389d844382dSKris Buschelman 390d844382dSKris Buschelman PetscFunctionBegin; 391b24902e0SBarry Smith ierr = MatView_SeqAIJ(A,viewer);CHKERRQ(ierr); 392d844382dSKris Buschelman 39332077d6dSBarry Smith ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 39432077d6dSBarry Smith if (iascii) { 395d844382dSKris Buschelman ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 3962f59403fSHong Zhang if (format == PETSC_VIEWER_ASCII_INFO) { 397f0c56d0fSKris Buschelman ierr = MatFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 398d844382dSKris Buschelman } 399d844382dSKris Buschelman } 400d844382dSKris Buschelman PetscFunctionReturn(0); 401d844382dSKris Buschelman } 402d844382dSKris Buschelman 4032f71e704SKris Buschelman /*MC 404*5c9eb25fSBarry Smith MAT_SOLVER_UMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices 4052f71e704SKris Buschelman via the external package UMFPACK. 4062f71e704SKris Buschelman 407*5c9eb25fSBarry Smith config/configure.py --download-umfpack to install PETSc to use UMFPACK 4082f71e704SKris Buschelman 4092f71e704SKris Buschelman Consult UMFPACK documentation for more information about the Control parameters 4102f71e704SKris Buschelman which correspond to the options database keys below. 4112f71e704SKris Buschelman 4122f71e704SKris Buschelman Options Database Keys: 413*5c9eb25fSBarry Smith + -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL] 4142f71e704SKris Buschelman . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL] 4152f71e704SKris Buschelman . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE] 4162f71e704SKris Buschelman . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE] 4172f71e704SKris Buschelman . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT] 4182f71e704SKris Buschelman - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP] 4192f71e704SKris Buschelman 4202f71e704SKris Buschelman Level: beginner 4212f71e704SKris Buschelman 4222f71e704SKris Buschelman .seealso: PCLU 4232f71e704SKris Buschelman M*/ 4242f71e704SKris Buschelman 4252d4e2982SHong Zhang 4262d4e2982SHong Zhang 427