Actual source code: klu.c

      2: /*
      3:    Provides an interface to the KLUv1.2 sparse solver

      5:    When build with PETSC_USE_64BIT_INDICES this will use SuiteSparse_long as the
      6:    integer type in KLU, otherwise it will use int. This means
      7:    all integers in this file are simply declared as PetscInt. Also it means
      8:    that KLU SuiteSparse_long version MUST be built with 64 bit integers when used.

     10: */
     11: #include <../src/mat/impls/aij/seq/aij.h>

     13: #if defined(PETSC_USE_64BIT_INDICES)
     14: #define klu_K_defaults                   klu_l_defaults
     15: #define klu_K_analyze(a,b,c,d)           klu_l_analyze((SuiteSparse_long)a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d)
     16: #define klu_K_analyze_given(a,b,c,d,e,f) klu_l_analyze_given((SuiteSparse_long)a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,(SuiteSparse_long*)e,f)
     17: #define klu_K_free_symbolic              klu_l_free_symbolic
     18: #define klu_K_free_numeric               klu_l_free_numeric
     19: #define klu_K_common                     klu_l_common
     20: #define klu_K_symbolic                   klu_l_symbolic
     21: #define klu_K_numeric                    klu_l_numeric
     22: #if defined(PETSC_USE_COMPLEX)
     23: #define klu_K_factor(a,b,c,d,e)       klu_zl_factor((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e);
     24: #define klu_K_solve                   klu_zl_solve
     25: #define klu_K_tsolve                  klu_zl_tsolve
     26: #define klu_K_refactor                klu_zl_refactor
     27: #define klu_K_sort                    klu_zl_sort
     28: #define klu_K_flops                   klu_zl_flops
     29: #define klu_K_rgrowth                 klu_zl_rgrowth
     30: #define klu_K_condest                 klu_zl_condest
     31: #define klu_K_rcond                   klu_zl_rcond
     32: #define klu_K_scale                   klu_zl_scale
     33: #else
     34: #define klu_K_factor(a,b,c,d,e)       klu_l_factor((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e);
     35: #define klu_K_solve                   klu_l_solve
     36: #define klu_K_tsolve                  klu_l_tsolve
     37: #define klu_K_refactor                klu_l_refactor
     38: #define klu_K_sort                    klu_l_sort
     39: #define klu_K_flops                   klu_l_flops
     40: #define klu_K_rgrowth                 klu_l_rgrowth
     41: #define klu_K_condest                 klu_l_condest
     42: #define klu_K_rcond                   klu_l_rcond
     43: #define klu_K_scale                   klu_l_scale
     44: #endif
     45: #else
     46: #define klu_K_defaults                klu_defaults
     47: #define klu_K_analyze                 klu_analyze
     48: #define klu_K_analyze_given           klu_analyze_given
     49: #define klu_K_free_symbolic           klu_free_symbolic
     50: #define klu_K_free_numeric            klu_free_numeric
     51: #define klu_K_common                  klu_common
     52: #define klu_K_symbolic                klu_symbolic
     53: #define klu_K_numeric                 klu_numeric
     54: #if defined(PETSC_USE_COMPLEX)
     55: #define klu_K_factor                  klu_z_factor
     56: #define klu_K_solve                   klu_z_solve
     57: #define klu_K_tsolve                  klu_z_tsolve
     58: #define klu_K_refactor                klu_z_refactor
     59: #define klu_K_sort                    klu_z_sort
     60: #define klu_K_flops                   klu_z_flops
     61: #define klu_K_rgrowth                 klu_z_rgrowth
     62: #define klu_K_condest                 klu_z_condest
     63: #define klu_K_rcond                   klu_z_rcond
     64: #define klu_K_scale                   klu_z_scale
     65: #else
     66: #define klu_K_factor                  klu_factor
     67: #define klu_K_solve                   klu_solve
     68: #define klu_K_tsolve                  klu_tsolve
     69: #define klu_K_refactor                klu_refactor
     70: #define klu_K_sort                    klu_sort
     71: #define klu_K_flops                   klu_flops
     72: #define klu_K_rgrowth                 klu_rgrowth
     73: #define klu_K_condest                 klu_condest
     74: #define klu_K_rcond                   klu_rcond
     75: #define klu_K_scale                   klu_scale
     76: #endif
     77: #endif

     79: EXTERN_C_BEGIN
     80: #include <klu.h>
     81: EXTERN_C_END

     83: static const char *KluOrderingTypes[] = {"AMD","COLAMD","PETSC"};
     84: static const char *scale[] ={"NONE","SUM","MAX"};

     86: typedef struct {
     87:   klu_K_common   Common;
     88:   klu_K_symbolic *Symbolic;
     89:   klu_K_numeric  *Numeric;
     90:   PetscInt       *perm_c,*perm_r;
     91:   MatStructure   flg;
     92:   PetscBool      PetscMatOrdering;
     93:   PetscBool      CleanUpKLU;
     94: } Mat_KLU;

     96: static PetscErrorCode MatDestroy_KLU(Mat A)
     97: {
     99:   Mat_KLU    *lu=(Mat_KLU*)A->data;

    102:   if (lu->CleanUpKLU) {
    103:     klu_K_free_symbolic(&lu->Symbolic,&lu->Common);
    104:     klu_K_free_numeric(&lu->Numeric,&lu->Common);
    105:     PetscFree2(lu->perm_r,lu->perm_c);
    106:   }
    107:   PetscFree(A->data);
    108:   return(0);
    109: }

    111: static PetscErrorCode MatSolveTranspose_KLU(Mat A,Vec b,Vec x)
    112: {
    113:   Mat_KLU       *lu = (Mat_KLU*)A->data;
    114:   PetscScalar    *xa;
    116:   PetscInt       status;

    119:   /* KLU uses a column major format, solve Ax = b by klu_*_solve */
    120:   /* ----------------------------------*/
    121:   VecCopy(b,x); /* klu_solve stores the solution in rhs */
    122:   VecGetArray(x,&xa);
    123:   status = klu_K_solve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,&lu->Common);
    124:   if (status != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed");
    125:   VecRestoreArray(x,&xa);
    126:   return(0);
    127: }

    129: static PetscErrorCode MatSolve_KLU(Mat A,Vec b,Vec x)
    130: {
    131:   Mat_KLU       *lu = (Mat_KLU*)A->data;
    132:   PetscScalar    *xa;
    134:   PetscInt       status;

    137:   /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */
    138:   /* ----------------------------------*/
    139:   VecCopy(b,x); /* klu_solve stores the solution in rhs */
    140:   VecGetArray(x,&xa);
    141: #if defined(PETSC_USE_COMPLEX)
    142:   PetscInt conj_solve=1;
    143:   status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,conj_solve,&lu->Common); /* conjugate solve */
    144: #else
    145:   status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,xa,&lu->Common);
    146: #endif
    147:   if (status != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed");
    148:   VecRestoreArray(x,&xa);
    149:   return(0);
    150: }

    152: static PetscErrorCode MatLUFactorNumeric_KLU(Mat F,Mat A,const MatFactorInfo *info)
    153: {
    154:   Mat_KLU        *lu = (Mat_KLU*)(F)->data;
    155:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
    156:   PetscInt       *ai = a->i,*aj=a->j;
    157:   PetscScalar    *av = a->a;

    160:   /* numeric factorization of A' */
    161:   /* ----------------------------*/

    163:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) {
    164:     klu_K_free_numeric(&lu->Numeric,&lu->Common);
    165:   }
    166:   lu->Numeric = klu_K_factor(ai,aj,(PetscReal*)av,lu->Symbolic,&lu->Common);
    167:   if (!lu->Numeric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Numeric factorization failed");

    169:   lu->flg                = SAME_NONZERO_PATTERN;
    170:   lu->CleanUpKLU         = PETSC_TRUE;
    171:   F->ops->solve          = MatSolve_KLU;
    172:   F->ops->solvetranspose = MatSolveTranspose_KLU;
    173:   return(0);
    174: }

    176: static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
    177: {
    178:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
    179:   Mat_KLU       *lu = (Mat_KLU*)(F->data);
    181:   PetscInt       i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
    182:   const PetscInt *ra,*ca;

    185:   if (lu->PetscMatOrdering) {
    186:     ISGetIndices(r,&ra);
    187:     ISGetIndices(c,&ca);
    188:     PetscMalloc2(m,&lu->perm_r,n,&lu->perm_c);
    189:     /* we cannot simply memcpy on 64 bit archs */
    190:     for (i = 0; i < m; i++) lu->perm_r[i] = ra[i];
    191:     for (i = 0; i < n; i++) lu->perm_c[i] = ca[i];
    192:     ISRestoreIndices(r,&ra);
    193:     ISRestoreIndices(c,&ca);
    194:   }

    196:   /* symbolic factorization of A' */
    197:   /* ---------------------------------------------------------------------- */
    198:   if (lu->PetscMatOrdering) { /* use Petsc ordering */
    199:     lu->Symbolic = klu_K_analyze_given(n,ai,aj,lu->perm_c,lu->perm_r,&lu->Common);
    200:   } else { /* use klu internal ordering */
    201:     lu->Symbolic = klu_K_analyze(n,ai,aj,&lu->Common);
    202:   }
    203:   if (!lu->Symbolic) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Symbolic Factorization failed");

    205:   lu->flg                   = DIFFERENT_NONZERO_PATTERN;
    206:   lu->CleanUpKLU            = PETSC_TRUE;
    207:   (F)->ops->lufactornumeric = MatLUFactorNumeric_KLU;
    208:   return(0);
    209: }

    211: static PetscErrorCode MatView_Info_KLU(Mat A,PetscViewer viewer)
    212: {
    213:   Mat_KLU       *lu= (Mat_KLU*)A->data;
    214:   klu_K_numeric *Numeric=(klu_K_numeric*)lu->Numeric;

    218:   /* check if matrix is KLU type */
    219:   if (A->ops->solve != MatSolve_KLU) return(0);

    221:   PetscViewerASCIIPrintf(viewer,"KLU stats:\n");
    222:   PetscViewerASCIIPrintf(viewer,"  Number of diagonal blocks: %d\n",Numeric->nblocks);
    223:   PetscViewerASCIIPrintf(viewer,"  Total nonzeros=%d\n",Numeric->lnz+Numeric->unz);

    225:   PetscViewerASCIIPrintf(viewer,"KLU runtime parameters:\n");

    227:   /* Control parameters used by numeric factorization */
    228:   PetscViewerASCIIPrintf(viewer,"  Partial pivoting tolerance: %g\n",lu->Common.tol);
    229:   /* BTF preordering */
    230:   PetscViewerASCIIPrintf(viewer,"  BTF preordering enabled: %d\n",lu->Common.btf);
    231:   /* mat ordering */
    232:   if (!lu->PetscMatOrdering) {
    233:     PetscViewerASCIIPrintf(viewer,"  Ordering: %s (not using the PETSc ordering)\n",KluOrderingTypes[(int)lu->Common.ordering]);
    234:   } else {
    235:     PetscViewerASCIIPrintf(viewer,"  Using PETSc ordering\n");
    236:   }
    237:   /* matrix row scaling */
    238:   PetscViewerASCIIPrintf(viewer, "  Matrix row scaling: %s\n",scale[(int)lu->Common.scale]);
    239:   return(0);
    240: }

    242: static PetscErrorCode MatView_KLU(Mat A,PetscViewer viewer)
    243: {
    244:   PetscErrorCode    ierr;
    245:   PetscBool         iascii;
    246:   PetscViewerFormat format;

    249:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    250:   if (iascii) {
    251:     PetscViewerGetFormat(viewer,&format);
    252:     if (format == PETSC_VIEWER_ASCII_INFO) {
    253:       MatView_Info_KLU(A,viewer);
    254:     }
    255:   }
    256:   return(0);
    257: }

    259: PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A,MatSolverType *type)
    260: {
    262:   *type = MATSOLVERKLU;
    263:   return(0);
    264: }

    267: /*MC
    268:   MATSOLVERKLU = "klu" - A matrix type providing direct solvers (LU) for sequential matrices
    269:   via the external package KLU.

    271:   ./configure --download-suitesparse to install PETSc to use KLU

    273:   Use -pc_type lu -pc_factor_mat_solver_type klu to use this direct solver

    275:   Consult KLU documentation for more information on the options database keys below.

    277:   Options Database Keys:
    278: + -mat_klu_pivot_tol <0.001>                  - Partial pivoting tolerance
    279: . -mat_klu_use_btf <1>                        - Use BTF preordering
    280: . -mat_klu_ordering <AMD>                     - KLU reordering scheme to reduce fill-in (choose one of) AMD COLAMD PETSC
    281: - -mat_klu_row_scale <NONE>                   - Matrix row scaling (choose one of) NONE SUM MAX

    283:    Note: KLU is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

    285:    Level: beginner

    287: .seealso: PCLU, MATSOLVERUMFPACK, MATSOLVERCHOLMOD, PCFactorSetMatSolverType(), MatSolverType
    288: M*/

    290: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A,MatFactorType ftype,Mat *F)
    291: {
    292:   Mat            B;
    293:   Mat_KLU       *lu;
    295:   PetscInt       m=A->rmap->n,n=A->cmap->n,idx,status;
    296:   PetscBool      flg;

    299:   /* Create the factorization matrix F */
    300:   MatCreate(PetscObjectComm((PetscObject)A),&B);
    301:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
    302:   PetscStrallocpy("klu",&((PetscObject)B)->type_name);
    303:   MatSetUp(B);

    305:   PetscNewLog(B,&lu);

    307:   B->data                  = lu;
    308:   B->ops->getinfo          = MatGetInfo_External;
    309:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
    310:   B->ops->destroy          = MatDestroy_KLU;
    311:   B->ops->view             = MatView_KLU;

    313:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_klu);

    315:   B->factortype   = MAT_FACTOR_LU;
    316:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */
    317:   B->preallocated = PETSC_TRUE;

    319:   PetscFree(B->solvertype);
    320:   PetscStrallocpy(MATSOLVERKLU,&B->solvertype);
    321:   B->useordering = PETSC_TRUE;

    323:   /* initializations */
    324:   /* ------------------------------------------------*/
    325:   /* get the default control parameters */
    326:   status = klu_K_defaults(&lu->Common);
    327:   if (status <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Initialization failed");

    329:   lu->Common.scale = 0; /* No row scaling */

    331:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"KLU Options","Mat");
    332:   /* Partial pivoting tolerance */
    333:   PetscOptionsReal("-mat_klu_pivot_tol","Partial pivoting tolerance","None",lu->Common.tol,&lu->Common.tol,NULL);
    334:   /* BTF pre-ordering */
    335:   PetscOptionsInt("-mat_klu_use_btf","Enable BTF preordering","None",(PetscInt)lu->Common.btf,(PetscInt*)&lu->Common.btf,NULL);
    336:   /* Matrix reordering */
    337:   PetscOptionsEList("-mat_klu_ordering","Internal ordering method","None",KluOrderingTypes,sizeof(KluOrderingTypes)/sizeof(KluOrderingTypes[0]),KluOrderingTypes[0],&idx,&flg);
    338:   if (flg) {
    339:     if ((int)idx == 2) lu->PetscMatOrdering = PETSC_TRUE;   /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Klu perm_c) */
    340:     else lu->Common.ordering = (int)idx;
    341:   }
    342:   /* Matrix row scaling */
    343:   PetscOptionsEList("-mat_klu_row_scale","Matrix row scaling","None",scale,3,scale[0],&idx,&flg);
    344:   PetscOptionsEnd();
    345:   *F = B;
    346:   return(0);
    347: }