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"};
 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: {
 98:   Mat_KLU    *lu=(Mat_KLU*)A->data;

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

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

115:   /* KLU uses a column major format, solve Ax = b by klu_*_solve */
116:   /* ----------------------------------*/
117:   VecCopy(b,x); /* klu_solve stores the solution in rhs */
118:   VecGetArray(x,&xa);
119:   status = klu_K_solve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,&lu->Common);
121:   VecRestoreArray(x,&xa);
122:   return 0;
123: }

125: static PetscErrorCode MatSolve_KLU(Mat A,Vec b,Vec x)
126: {
127:   Mat_KLU     *lu = (Mat_KLU*)A->data;
128:   PetscScalar *xa;
129:   PetscInt     status;

131:   /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */
132:   /* ----------------------------------*/
133:   VecCopy(b,x); /* klu_solve stores the solution in rhs */
134:   VecGetArray(x,&xa);
135: #if defined(PETSC_USE_COMPLEX)
136:   PetscInt conj_solve=1;
137:   status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,conj_solve,&lu->Common); /* conjugate solve */
138: #else
139:   status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,xa,&lu->Common);
140: #endif
142:   VecRestoreArray(x,&xa);
143:   return 0;
144: }

146: static PetscErrorCode MatLUFactorNumeric_KLU(Mat F,Mat A,const MatFactorInfo *info)
147: {
148:   Mat_KLU        *lu = (Mat_KLU*)(F)->data;
149:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
150:   PetscInt       *ai = a->i,*aj=a->j;
151:   PetscScalar    *av = a->a;

153:   /* numeric factorization of A' */
154:   /* ----------------------------*/

156:   if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) {
157:     klu_K_free_numeric(&lu->Numeric,&lu->Common);
158:   }
159:   lu->Numeric = klu_K_factor(ai,aj,(PetscReal*)av,lu->Symbolic,&lu->Common);

162:   lu->flg                = SAME_NONZERO_PATTERN;
163:   lu->CleanUpKLU         = PETSC_TRUE;
164:   F->ops->solve          = MatSolve_KLU;
165:   F->ops->solvetranspose = MatSolveTranspose_KLU;
166:   return 0;
167: }

169: static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
170: {
171:   Mat_SeqAIJ     *a  = (Mat_SeqAIJ*)A->data;
172:   Mat_KLU        *lu = (Mat_KLU*)(F->data);
173:   PetscInt       i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
174:   const PetscInt *ra,*ca;

176:   if (lu->PetscMatOrdering) {
177:     ISGetIndices(r,&ra);
178:     ISGetIndices(c,&ca);
179:     PetscMalloc2(m,&lu->perm_r,n,&lu->perm_c);
180:     /* we cannot simply memcpy on 64 bit archs */
181:     for (i = 0; i < m; i++) lu->perm_r[i] = ra[i];
182:     for (i = 0; i < n; i++) lu->perm_c[i] = ca[i];
183:     ISRestoreIndices(r,&ra);
184:     ISRestoreIndices(c,&ca);
185:   }

187:   /* symbolic factorization of A' */
188:   /* ---------------------------------------------------------------------- */
189:   if (r) {
190:     lu->PetscMatOrdering = PETSC_TRUE;
191:     lu->Symbolic = klu_K_analyze_given(n,ai,aj,lu->perm_c,lu->perm_r,&lu->Common);
192:   } else { /* use klu internal ordering */
193:     lu->Symbolic = klu_K_analyze(n,ai,aj,&lu->Common);
194:   }

197:   lu->flg                   = DIFFERENT_NONZERO_PATTERN;
198:   lu->CleanUpKLU            = PETSC_TRUE;
199:   (F)->ops->lufactornumeric = MatLUFactorNumeric_KLU;
200:   return 0;
201: }

203: static PetscErrorCode MatView_Info_KLU(Mat A,PetscViewer viewer)
204: {
205:   Mat_KLU       *lu= (Mat_KLU*)A->data;
206:   klu_K_numeric *Numeric=(klu_K_numeric*)lu->Numeric;

208:   PetscViewerASCIIPrintf(viewer,"KLU stats:\n");
209:   PetscViewerASCIIPrintf(viewer,"  Number of diagonal blocks: %" PetscInt_FMT "\n",(PetscInt)(Numeric->nblocks));
210:   PetscViewerASCIIPrintf(viewer,"  Total nonzeros=%" PetscInt_FMT "\n",(PetscInt)(Numeric->lnz+Numeric->unz));
211:   PetscViewerASCIIPrintf(viewer,"KLU runtime parameters:\n");
212:   /* Control parameters used by numeric factorization */
213:   PetscViewerASCIIPrintf(viewer,"  Partial pivoting tolerance: %g\n",lu->Common.tol);
214:   /* BTF preordering */
215:   PetscViewerASCIIPrintf(viewer,"  BTF preordering enabled: %" PetscInt_FMT "\n",(PetscInt)(lu->Common.btf));
216:   /* mat ordering */
217:   if (!lu->PetscMatOrdering) {
218:     PetscViewerASCIIPrintf(viewer,"  Ordering: %s (not using the PETSc ordering)\n",KluOrderingTypes[(int)lu->Common.ordering]);
219:   }
220:   /* matrix row scaling */
221:   PetscViewerASCIIPrintf(viewer, "  Matrix row scaling: %s\n",scale[(int)lu->Common.scale]);
222:   return 0;
223: }

225: static PetscErrorCode MatView_KLU(Mat A,PetscViewer viewer)
226: {
227:   PetscBool         iascii;
228:   PetscViewerFormat format;

230:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
231:   if (iascii) {
232:     PetscViewerGetFormat(viewer,&format);
233:     if (format == PETSC_VIEWER_ASCII_INFO) {
234:       MatView_Info_KLU(A,viewer);
235:     }
236:   }
237:   return 0;
238: }

240: PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A,MatSolverType *type)
241: {
242:   *type = MATSOLVERKLU;
243:   return 0;
244: }

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

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

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

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

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

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

264:    Level: beginner

266: .seealso: PCLU, MATSOLVERUMFPACK, MATSOLVERCHOLMOD, PCFactorSetMatSolverType(), MatSolverType
267: M*/

269: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A,MatFactorType ftype,Mat *F)
270: {
271:   Mat            B;
272:   Mat_KLU       *lu;
274:   PetscInt       m=A->rmap->n,n=A->cmap->n,idx = 0,status;
275:   PetscBool      flg;

277:   /* Create the factorization matrix F */
278:   MatCreate(PetscObjectComm((PetscObject)A),&B);
279:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
280:   PetscStrallocpy("klu",&((PetscObject)B)->type_name);
281:   MatSetUp(B);

283:   PetscNewLog(B,&lu);

285:   B->data                  = lu;
286:   B->ops->getinfo          = MatGetInfo_External;
287:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
288:   B->ops->destroy          = MatDestroy_KLU;
289:   B->ops->view             = MatView_KLU;

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

293:   B->factortype   = MAT_FACTOR_LU;
294:   B->assembled    = PETSC_TRUE;           /* required by -ksp_view */
295:   B->preallocated = PETSC_TRUE;

297:   PetscFree(B->solvertype);
298:   PetscStrallocpy(MATSOLVERKLU,&B->solvertype);
299:   B->canuseordering = PETSC_TRUE;
300:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU]);

302:   /* initializations */
303:   /* ------------------------------------------------*/
304:   /* get the default control parameters */
305:   status = klu_K_defaults(&lu->Common);

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

310:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"KLU Options","Mat");
311:   /* Partial pivoting tolerance */
312:   PetscOptionsReal("-mat_klu_pivot_tol","Partial pivoting tolerance","None",lu->Common.tol,&lu->Common.tol,NULL);
313:   /* BTF pre-ordering */
314:   PetscOptionsInt("-mat_klu_use_btf","Enable BTF preordering","None",(PetscInt)lu->Common.btf,(PetscInt*)&lu->Common.btf,NULL);
315:   /* Matrix reordering */
316:   PetscOptionsEList("-mat_klu_ordering","Internal ordering method","None",KluOrderingTypes,sizeof(KluOrderingTypes)/sizeof(KluOrderingTypes[0]),KluOrderingTypes[0],&idx,&flg);
317:   lu->Common.ordering = (int)idx;
318:   /* Matrix row scaling */
319:   PetscOptionsEList("-mat_klu_row_scale","Matrix row scaling","None",scale,3,scale[0],&idx,&flg);
320:   PetscOptionsEnd();
321:   *F = B;
322:   return 0;
323: }