Actual source code: klu.c
petsc-3.7.7 2017-09-25
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
80: EXTERN_C_BEGIN
81: #include <klu.h>
82: EXTERN_C_END
84: static const char *KluOrderingTypes[] = {"AMD","COLAMD","PETSC"};
85: static const char *scale[] ={"NONE","SUM","MAX"};
87: typedef struct {
88: klu_K_common Common;
89: klu_K_symbolic *Symbolic;
90: klu_K_numeric *Numeric;
91: PetscInt *perm_c,*perm_r;
92: MatStructure flg;
93: PetscBool PetscMatOrdering;
95: /* Flag to clean up KLU objects during Destroy */
96: PetscBool CleanUpKLU;
97: } Mat_KLU;
101: static PetscErrorCode MatDestroy_KLU(Mat A)
102: {
104: Mat_KLU *lu=(Mat_KLU*)A->spptr;
107: if (lu && lu->CleanUpKLU) {
108: klu_K_free_symbolic(&lu->Symbolic,&lu->Common);
109: klu_K_free_numeric(&lu->Numeric,&lu->Common);
110: PetscFree2(lu->perm_r,lu->perm_c);
111: }
112: PetscFree(A->spptr);
113: MatDestroy_SeqAIJ(A);
114: return(0);
115: }
119: static PetscErrorCode MatSolveTranspose_KLU(Mat A,Vec b,Vec x)
120: {
121: Mat_KLU *lu = (Mat_KLU*)A->spptr;
122: PetscScalar *xa;
124: PetscInt status;
127: /* KLU uses a column major format, solve Ax = b by klu_*_solve */
128: /* ----------------------------------*/
129: VecCopy(b,x); /* klu_solve stores the solution in rhs */
130: VecGetArray(x,&xa);
131: status = klu_K_solve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,&lu->Common);
132: if (status != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed");
133: VecRestoreArray(x,&xa);
134: return(0);
135: }
139: static PetscErrorCode MatSolve_KLU(Mat A,Vec b,Vec x)
140: {
141: Mat_KLU *lu = (Mat_KLU*)A->spptr;
142: PetscScalar *xa;
144: PetscInt status;
147: /* KLU uses a column major format, solve A^Tx = b by klu_*_tsolve */
148: /* ----------------------------------*/
149: VecCopy(b,x); /* klu_solve stores the solution in rhs */
150: VecGetArray(x,&xa);
151: #if defined(PETSC_USE_COMPLEX)
152: PetscInt conj_solve=1;
153: status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,(PetscReal*)xa,conj_solve,&lu->Common); /* conjugate solve */
154: #else
155: status = klu_K_tsolve(lu->Symbolic,lu->Numeric,A->rmap->n,1,xa,&lu->Common);
156: #endif
157: if (status != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Solve failed");
158: VecRestoreArray(x,&xa);
159: return(0);
160: }
164: static PetscErrorCode MatLUFactorNumeric_KLU(Mat F,Mat A,const MatFactorInfo *info)
165: {
166: Mat_KLU *lu = (Mat_KLU*)(F)->spptr;
167: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
168: PetscInt *ai = a->i,*aj=a->j;
169: PetscScalar *av = a->a;
172: /* numeric factorization of A' */
173: /* ----------------------------*/
175: if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) {
176: klu_K_free_numeric(&lu->Numeric,&lu->Common);
177: }
178: lu->Numeric = klu_K_factor(ai,aj,(PetscReal*)av,lu->Symbolic,&lu->Common);
179: if(!lu->Numeric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Numeric factorization failed");
181: lu->flg = SAME_NONZERO_PATTERN;
182: lu->CleanUpKLU = PETSC_TRUE;
183: F->ops->solve = MatSolve_KLU;
184: F->ops->solvetranspose = MatSolveTranspose_KLU;
185: return(0);
186: }
190: static PetscErrorCode MatLUFactorSymbolic_KLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
191: {
192: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
193: Mat_KLU *lu = (Mat_KLU*)(F->spptr);
195: PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
196: const PetscInt *ra,*ca;
199: if (lu->PetscMatOrdering) {
200: ISGetIndices(r,&ra);
201: ISGetIndices(c,&ca);
202: PetscMalloc2(m,&lu->perm_r,n,&lu->perm_c);
203: /* we cannot simply memcpy on 64 bit archs */
204: for (i = 0; i < m; i++) lu->perm_r[i] = ra[i];
205: for (i = 0; i < n; i++) lu->perm_c[i] = ca[i];
206: ISRestoreIndices(r,&ra);
207: ISRestoreIndices(c,&ca);
208: }
210: /* symbolic factorization of A' */
211: /* ---------------------------------------------------------------------- */
212: if (lu->PetscMatOrdering) { /* use Petsc ordering */
213: lu->Symbolic = klu_K_analyze_given(n,ai,aj,lu->perm_c,lu->perm_r,&lu->Common);
214: } else { /* use klu internal ordering */
215: lu->Symbolic = klu_K_analyze(n,ai,aj,&lu->Common);
216: }
217: if (!lu->Symbolic) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Symbolic Factorization failed");
219: lu->flg = DIFFERENT_NONZERO_PATTERN;
220: lu->CleanUpKLU = PETSC_TRUE;
221: (F)->ops->lufactornumeric = MatLUFactorNumeric_KLU;
222: return(0);
223: }
227: static PetscErrorCode MatFactorInfo_KLU(Mat A,PetscViewer viewer)
228: {
229: Mat_KLU *lu= (Mat_KLU*)A->spptr;
230: klu_K_numeric *Numeric=(klu_K_numeric*)lu->Numeric;
234: /* check if matrix is KLU type */
235: if (A->ops->solve != MatSolve_KLU) return(0);
237: PetscViewerASCIIPrintf(viewer,"KLU stats:\n");
238: PetscViewerASCIIPrintf(viewer," Number of diagonal blocks: %d\n",Numeric->nblocks);
239: PetscViewerASCIIPrintf(viewer," Total nonzeros=%d\n",Numeric->lnz+Numeric->unz);
241: PetscViewerASCIIPrintf(viewer,"KLU runtime parameters:\n");
243: /* Control parameters used by numeric factorization */
244: PetscViewerASCIIPrintf(viewer," Partial pivoting tolerance: %g\n",lu->Common.tol);
245: /* BTF preordering */
246: PetscViewerASCIIPrintf(viewer," BTF preordering enabled: %d\n",lu->Common.btf);
247: /* mat ordering */
248: if (!lu->PetscMatOrdering) {
249: PetscViewerASCIIPrintf(viewer," Ordering: %s (not using the PETSc ordering)\n",KluOrderingTypes[(int)lu->Common.ordering]);
250: } else {
251: PetscViewerASCIIPrintf(viewer," Using PETSc ordering\n");
252: }
253: /* matrix row scaling */
254: PetscViewerASCIIPrintf(viewer, " Matrix row scaling: %s\n",scale[(int)lu->Common.scale]);
255: return(0);
256: }
260: static PetscErrorCode MatView_KLU(Mat A,PetscViewer viewer)
261: {
262: PetscErrorCode ierr;
263: PetscBool iascii;
264: PetscViewerFormat format;
267: MatView_SeqAIJ(A,viewer);
269: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
270: if (iascii) {
271: PetscViewerGetFormat(viewer,&format);
272: if (format == PETSC_VIEWER_ASCII_INFO) {
273: MatFactorInfo_KLU(A,viewer);
274: }
275: }
276: return(0);
277: }
281: PetscErrorCode MatFactorGetSolverPackage_seqaij_klu(Mat A,const MatSolverPackage *type)
282: {
284: *type = MATSOLVERKLU;
285: return(0);
286: }
289: /*MC
290: MATSOLVERKLU = "klu" - A matrix type providing direct solvers (LU) for sequential matrices
291: via the external package KLU.
293: ./configure --download-suitesparse to install PETSc to use KLU
295: Use -pc_type lu -pc_factor_mat_solver_package klu to us this direct solver
297: Consult KLU documentation for more information on the options database keys below.
299: Options Database Keys:
300: + -mat_klu_pivot_tol <0.001> - Partial pivoting tolerance
301: . -mat_klu_use_btf <1> - Use BTF preordering
302: . -mat_klu_ordering <AMD> - KLU reordering scheme to reduce fill-in (choose one of) AMD COLAMD PETSC
303: - -mat_klu_row_scale <NONE> - Matrix row scaling (choose one of) NONE SUM MAX
305: Note: KLU is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html
307: Level: beginner
309: .seealso: PCLU, MATSOLVERUMFPACK, MATSOLVERCHOLMOD, PCFactorSetMatSolverPackage(), MatSolverPackage
310: M*/
314: PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat A,MatFactorType ftype,Mat *F)
315: {
316: Mat B;
317: Mat_KLU *lu;
319: PetscInt m=A->rmap->n,n=A->cmap->n,idx,status;
320: PetscBool flg;
323: /* Create the factorization matrix F */
324: MatCreate(PetscObjectComm((PetscObject)A),&B);
325: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
326: MatSetType(B,((PetscObject)A)->type_name);
327: MatSeqAIJSetPreallocation(B,0,NULL);
328: PetscNewLog(B,&lu);
330: B->spptr = lu;
331: B->ops->lufactorsymbolic = MatLUFactorSymbolic_KLU;
332: B->ops->destroy = MatDestroy_KLU;
333: B->ops->view = MatView_KLU;
335: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_klu);
337: B->factortype = MAT_FACTOR_LU;
338: B->assembled = PETSC_TRUE; /* required by -ksp_view */
339: B->preallocated = PETSC_TRUE;
341: PetscFree(B->solvertype);
342: PetscStrallocpy(MATSOLVERKLU,&B->solvertype);
344: /* initializations */
345: /* ------------------------------------------------*/
346: /* get the default control parameters */
347: status = klu_K_defaults(&lu->Common);
348: if(status <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"KLU Initialization failed");
350: lu->Common.scale = 0; /* No row scaling */
352: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"KLU Options","Mat");
353: /* Partial pivoting tolerance */
354: PetscOptionsReal("-mat_klu_pivot_tol","Partial pivoting tolerance","None",lu->Common.tol,&lu->Common.tol,NULL);
355: /* BTF pre-ordering */
356: PetscOptionsInt("-mat_klu_use_btf","Enable BTF preordering","None",(PetscInt)lu->Common.btf,(PetscInt*)&lu->Common.btf,NULL);
357: /* Matrix reordering */
358: PetscOptionsEList("-mat_klu_ordering","Internal ordering method","None",KluOrderingTypes,sizeof(KluOrderingTypes)/sizeof(KluOrderingTypes[0]),KluOrderingTypes[0],&idx,&flg);
359: if (flg) {
360: if ((int)idx == 2) lu->PetscMatOrdering = PETSC_TRUE; /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Klu perm_c) */
361: else lu->Common.ordering = (int)idx;
362: }
363: /* Matrix row scaling */
364: PetscOptionsEList("-mat_klu_row_scale","Matrix row scaling","None",scale,3,scale[0],&idx,&flg);
365: PetscOptionsEnd();
366: *F = B;
367: return(0);
368: }