xref: /petsc/src/mat/impls/aij/seq/umfpack/umfpack.c (revision 5c9eb25f753630cfd361293e05e7358a00a954ac)
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