Actual source code: umfpack.c
1: #define PETSCMAT_DLL
3: /*
4: Provides an interface to the UMFPACKv4.3 sparse solver
5: */
7: #include src/mat/impls/aij/seq/aij.h
10: #include "umfpack.h"
13: typedef struct {
14: void *Symbolic, *Numeric;
15: double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W;
16: int *Wi,*ai,*aj,*perm_c;
17: PetscScalar *av;
18: MatStructure flg;
19: PetscTruth PetscMatOdering;
21: /* A few function pointers for inheritance */
22: PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
23: PetscErrorCode (*MatView)(Mat,PetscViewer);
24: PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType);
25: PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
26: PetscErrorCode (*MatDestroy)(Mat);
27:
28: /* Flag to clean up UMFPACK objects during Destroy */
29: PetscTruth CleanUpUMFPACK;
30: } Mat_UMFPACK;
32: EXTERN PetscErrorCode MatDuplicate_UMFPACK(Mat,MatDuplicateOption,Mat*);
37: PetscErrorCode MatConvert_UMFPACK_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat)
38: {
40: Mat B=*newmat;
41: Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr;
44: if (reuse == MAT_INITIAL_MATRIX) {
45: MatDuplicate(A,MAT_COPY_VALUES,&B);
46: }
47: /* Reset the original function pointers */
48: B->ops->duplicate = lu->MatDuplicate;
49: B->ops->view = lu->MatView;
50: B->ops->assemblyend = lu->MatAssemblyEnd;
51: B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
52: B->ops->destroy = lu->MatDestroy;
54: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_umfpack_C","",PETSC_NULL);
55: PetscObjectComposeFunction((PetscObject)B,"MatConvert_umfpack_seqaij_C","",PETSC_NULL);
57: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
58: *newmat = B;
59: return(0);
60: }
65: PetscErrorCode MatDestroy_UMFPACK(Mat A)
66: {
68: Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr;
71: if (lu->CleanUpUMFPACK) {
72: umfpack_di_free_symbolic(&lu->Symbolic);
73: umfpack_di_free_numeric(&lu->Numeric);
74: PetscFree(lu->Wi);
75: PetscFree(lu->W);
76: if (lu->PetscMatOdering) {
77: PetscFree(lu->perm_c);
78: }
79: }
80: MatConvert_UMFPACK_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
81: (*A->ops->destroy)(A);
82: return(0);
83: }
87: PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
88: {
89: Mat_UMFPACK *lu = (Mat_UMFPACK*)A->spptr;
90: PetscScalar *av=lu->av,*ba,*xa;
92: int *ai=lu->ai,*aj=lu->aj,status;
93:
95: /* solve Ax = b by umfpack_di_wsolve */
96: /* ----------------------------------*/
97: VecGetArray(b,&ba);
98: VecGetArray(x,&xa);
100: status = umfpack_di_wsolve(UMFPACK_At,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
101: umfpack_di_report_info(lu->Control, lu->Info);
102: if (status < 0){
103: umfpack_di_report_status(lu->Control, status);
104: SETERRQ(PETSC_ERR_LIB,"umfpack_di_wsolve failed");
105: }
106:
107: VecRestoreArray(b,&ba);
108: VecRestoreArray(x,&xa);
110: return(0);
111: }
115: PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat A,MatFactorInfo *info,Mat *F)
116: {
117: Mat_UMFPACK *lu=(Mat_UMFPACK*)(*F)->spptr;
119: int *ai=lu->ai,*aj=lu->aj,m=A->rmap.n,status;
120: PetscScalar *av=lu->av;
123: /* numeric factorization of A' */
124: /* ----------------------------*/
125: if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric){
126: umfpack_di_free_numeric(&lu->Numeric);
127: }
128: status = umfpack_di_numeric (ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
129: if (status < 0) SETERRQ(PETSC_ERR_LIB,"umfpack_di_numeric failed");
130: /* report numeric factorization of A' when Control[PRL] > 3 */
131: (void) umfpack_di_report_numeric (lu->Numeric, lu->Control);
133: if (lu->flg == DIFFERENT_NONZERO_PATTERN){ /* first numeric factorization */
134: /* allocate working space to be used by Solve */
135: PetscMalloc(m * sizeof(int), &lu->Wi);
136: PetscMalloc(5*m * sizeof(double), &lu->W);
137: }
138: lu->flg = SAME_NONZERO_PATTERN;
139: lu->CleanUpUMFPACK = PETSC_TRUE;
140: return(0);
141: }
143: /*
144: Note the r permutation is ignored
145: */
148: PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
149: {
150: Mat B;
151: Mat_SeqAIJ *mat=(Mat_SeqAIJ*)A->data;
152: Mat_UMFPACK *lu;
154: int m=A->rmap.n,n=A->cmap.n,*ai=mat->i,*aj=mat->j,status,*ra,idx;
155: PetscScalar *av=mat->a;
156: const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC","2BY2"},
157: *scale[]={"NONE","SUM","MAX"};
158: PetscTruth flg;
159:
161: /* Create the factorization matrix F */
162: MatCreate(A->comm,&B);
163: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
164: MatSetType(B,A->type_name);
165: MatSeqAIJSetPreallocation(B,0,PETSC_NULL);
166:
167: B->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
168: B->ops->solve = MatSolve_UMFPACK;
169: B->factor = FACTOR_LU;
170: B->assembled = PETSC_TRUE; /* required by -ksp_view */
172: lu = (Mat_UMFPACK*)(B->spptr);
173:
174: /* initializations */
175: /* ------------------------------------------------*/
176: /* get the default control parameters */
177: umfpack_di_defaults (lu->Control);
178: lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */
179: lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
181: PetscOptionsBegin(A->comm,A->prefix,"UMFPACK Options","Mat");
182: /* Control parameters used by reporting routiones */
183: PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);
185: /* Control parameters for symbolic factorization */
186: PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,4,strategy[0],&idx,&flg);
187: if (flg) {
188: switch (idx){
189: case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
190: case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
191: case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
192: case 3: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_2BY2; break;
193: }
194: }
195: PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);
196: PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);
197: PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);
198: PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);
199: PetscOptionsReal("-mat_umfpack_2by2_tolerance","Control[UMFPACK_2BY2_TOLERANCE]","None",lu->Control[UMFPACK_2BY2_TOLERANCE],&lu->Control[UMFPACK_2BY2_TOLERANCE],PETSC_NULL);
200: PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);
201: PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);
203: /* Control parameters used by numeric factorization */
204: PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);
205: 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);
206: PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);
207: if (flg) {
208: switch (idx){
209: case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
210: case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
211: case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
212: }
213: }
214: PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],PETSC_NULL);
215: 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);
216: PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],PETSC_NULL);
218: /* Control parameters used by solve */
219: PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],PETSC_NULL);
220:
221: /* use Petsc mat ordering (note: size is for the transpose, and PETSc r = Umfpack perm_c) */
222: PetscOptionsHasName(PETSC_NULL,"-pc_factor_mat_ordering_type",&lu->PetscMatOdering);
223: if (lu->PetscMatOdering) {
224: ISGetIndices(r,&ra);
225: PetscMalloc(A->rmap.n*sizeof(int),&lu->perm_c);
226: PetscMemcpy(lu->perm_c,ra,A->rmap.n*sizeof(int));
227: ISRestoreIndices(r,&ra);
228: }
229: PetscOptionsEnd();
231: /* print the control parameters */
232: if( lu->Control[UMFPACK_PRL] > 1 ) umfpack_di_report_control (lu->Control);
234: /* symbolic factorization of A' */
235: /* ---------------------------------------------------------------------- */
236: if (lu->PetscMatOdering) { /* use Petsc row ordering */
237: status = umfpack_di_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
238: } else { /* use Umfpack col ordering */
239: status = umfpack_di_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info);
240: }
241: if (status < 0){
242: umfpack_di_report_info(lu->Control, lu->Info);
243: umfpack_di_report_status(lu->Control, status);
244: SETERRQ(PETSC_ERR_LIB,"umfpack_di_symbolic failed");
245: }
246: /* report sumbolic factorization of A' when Control[PRL] > 3 */
247: (void) umfpack_di_report_symbolic(lu->Symbolic, lu->Control);
249: lu->flg = DIFFERENT_NONZERO_PATTERN;
250: lu->ai = ai;
251: lu->aj = aj;
252: lu->av = av;
253: lu->CleanUpUMFPACK = PETSC_TRUE;
254: *F = B;
255: return(0);
256: }
260: PetscErrorCode MatAssemblyEnd_UMFPACK(Mat A,MatAssemblyType mode)
261: {
263: Mat_UMFPACK *lu=(Mat_UMFPACK*)(A->spptr);
266: (*lu->MatAssemblyEnd)(A,mode);
267: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
268: A->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
269: return(0);
270: }
272: /* used by -ksp_view */
275: PetscErrorCode MatFactorInfo_UMFPACK(Mat A,PetscViewer viewer)
276: {
277: Mat_UMFPACK *lu= (Mat_UMFPACK*)A->spptr;
281: /* check if matrix is UMFPACK type */
282: if (A->ops->solve != MatSolve_UMFPACK) return(0);
284: PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");
285: /* Control parameters used by reporting routiones */
286: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);
288: /* Control parameters used by symbolic factorization */
289: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);
290: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);
291: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);
292: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);
293: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);
294: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_2BY2_TOLERANCE]: %g\n",lu->Control[UMFPACK_2BY2_TOLERANCE]);
295: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);
296: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);
298: /* Control parameters used by numeric factorization */
299: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);
300: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);
301: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);
302: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);
303: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);
305: /* Control parameters used by solve */
306: PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);
308: /* mat ordering */
309: if(!lu->PetscMatOdering) PetscViewerASCIIPrintf(viewer," UMFPACK default matrix ordering is used (not the PETSc matrix ordering) \n");
311: return(0);
312: }
316: PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
317: {
319: PetscTruth iascii;
320: PetscViewerFormat format;
321: Mat_UMFPACK *lu=(Mat_UMFPACK*)(A->spptr);
324: (*lu->MatView)(A,viewer);
326: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
327: if (iascii) {
328: PetscViewerGetFormat(viewer,&format);
329: if (format == PETSC_VIEWER_ASCII_INFO) {
330: MatFactorInfo_UMFPACK(A,viewer);
331: }
332: }
333: return(0);
334: }
339: PetscErrorCode MatConvert_SeqAIJ_UMFPACK(Mat A,MatType type,MatReuse reuse,Mat *newmat)
340: {
342: Mat B=*newmat;
343: Mat_UMFPACK *lu;
346: if (reuse == MAT_INITIAL_MATRIX) {
347: MatDuplicate(A,MAT_COPY_VALUES,&B);
348: }
350: PetscNew(Mat_UMFPACK,&lu);
351: lu->MatDuplicate = A->ops->duplicate;
352: lu->MatView = A->ops->view;
353: lu->MatAssemblyEnd = A->ops->assemblyend;
354: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
355: lu->MatDestroy = A->ops->destroy;
356: lu->CleanUpUMFPACK = PETSC_FALSE;
358: B->spptr = (void*)lu;
359: B->ops->duplicate = MatDuplicate_UMFPACK;
360: B->ops->view = MatView_UMFPACK;
361: B->ops->assemblyend = MatAssemblyEnd_UMFPACK;
362: B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
363: B->ops->destroy = MatDestroy_UMFPACK;
365: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_umfpack_C",
366: "MatConvert_SeqAIJ_UMFPACK",MatConvert_SeqAIJ_UMFPACK);
367: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_umfpack_seqaij_C",
368: "MatConvert_UMFPACK_SeqAIJ",MatConvert_UMFPACK_SeqAIJ);
369: PetscInfo(0,"Using UMFPACK for SeqAIJ LU factorization and solves.\n");
370: PetscObjectChangeTypeName((PetscObject)B,MATUMFPACK);
371: *newmat = B;
372: return(0);
373: }
378: PetscErrorCode MatDuplicate_UMFPACK(Mat A, MatDuplicateOption op, Mat *M)
379: {
381: Mat_UMFPACK *lu=(Mat_UMFPACK*)A->spptr;
384: (*lu->MatDuplicate)(A,op,M);
385: PetscMemcpy((*M)->spptr,lu,sizeof(Mat_UMFPACK));
386: return(0);
387: }
389: /*MC
390: MATUMFPACK - MATUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
391: via the external package UMFPACK.
393: If UMFPACK is installed (see the manual for
394: instructions on how to declare the existence of external packages),
395: a matrix type can be constructed which invokes UMFPACK solvers.
396: After calling MatCreate(...,A), simply call MatSetType(A,UMFPACK).
397: This matrix type is only supported for double precision real.
399: This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is
400: supported for this matrix type. One can also call MatConvert for an inplace conversion to or from
401: the MATSEQAIJ type without data copy.
403: Consult UMFPACK documentation for more information about the Control parameters
404: which correspond to the options database keys below.
406: Options Database Keys:
407: + -mat_type umfpack - sets the matrix type to "umfpack" during a call to MatSetFromOptions()
408: . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
409: . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
410: . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
411: . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
412: . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
413: - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
415: Level: beginner
417: .seealso: PCLU
418: M*/
423: PetscErrorCode MatCreate_UMFPACK(Mat A)
424: {
428: MatSetType(A,MATSEQAIJ);
429: MatConvert_SeqAIJ_UMFPACK(A,MATUMFPACK,MAT_REUSE_MATRIX,&A);
430: return(0);
431: }