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:   PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL);
106:   PetscFree(A->data);
107:   return 0;
108: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

222: static PetscErrorCode MatView_KLU(Mat A, PetscViewer viewer)
223: {
224:   PetscBool         iascii;
225:   PetscViewerFormat format;

227:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
228:   if (iascii) {
229:     PetscViewerGetFormat(viewer, &format);
230:     if (format == PETSC_VIEWER_ASCII_INFO) MatView_Info_KLU(A, viewer);
231:   }
232:   return 0;
233: }

235: PetscErrorCode MatFactorGetSolverType_seqaij_klu(Mat A, MatSolverType *type)
236: {
237:   *type = MATSOLVERKLU;
238:   return 0;
239: }

241: /*MC
242:   MATSOLVERKLU = "klu" - A matrix type providing direct solvers, LU, for sequential matrices
243:   via the external package KLU.

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

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

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

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

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

259:    Level: beginner

261: .seealso: `PCLU`, `MATSOLVERUMFPACK`, `MATSOLVERCHOLMOD`, `PCFactorSetMatSolverType()`, `MatSolverType`
262: M*/

264: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A, MatFactorType ftype, Mat *F)
265: {
266:   Mat       B;
267:   Mat_KLU  *lu;
268:   PetscInt  m = A->rmap->n, n = A->cmap->n, idx = 0, status;
269:   PetscBool flg;

271:   /* Create the factorization matrix F */
272:   MatCreate(PetscObjectComm((PetscObject)A), &B);
273:   MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n);
274:   PetscStrallocpy("klu", &((PetscObject)B)->type_name);
275:   MatSetUp(B);

277:   PetscNew(&lu);

279:   B->data                  = lu;
280:   B->ops->getinfo          = MatGetInfo_External;
281:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
282:   B->ops->destroy          = MatDestroy_KLU;
283:   B->ops->view             = MatView_KLU;

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

287:   B->factortype   = MAT_FACTOR_LU;
288:   B->assembled    = PETSC_TRUE; /* required by -ksp_view */
289:   B->preallocated = PETSC_TRUE;

291:   PetscFree(B->solvertype);
292:   PetscStrallocpy(MATSOLVERKLU, &B->solvertype);
293:   B->canuseordering = PETSC_TRUE;
294:   PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]);

296:   /* initializations */
297:   /* ------------------------------------------------*/
298:   /* get the default control parameters */
299:   status = klu_K_defaults(&lu->Common);

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

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