1
2 /*
3 Provides an interface to the UMFPACK sparse solver available through SuiteSparse version 4.2.1
4
5 When build with PETSC_USE_64BIT_INDICES this will use Suitesparse_long as the
6 integer type in UMFPACK, otherwise it will use int. This means
7 all integers in this file as simply declared as PetscInt. Also it means
8 that one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only]
9
10 */
11 #include <../src/mat/impls/aij/seq/aij.h>
12
13 #if defined(PETSC_USE_64BIT_INDICES)
14 #if defined(PETSC_USE_COMPLEX)
15 #define umfpack_UMF_free_symbolic umfpack_zl_free_symbolic
16 #define umfpack_UMF_free_numeric umfpack_zl_free_numeric
17 /* the type casts are needed because PetscInt is long long while SuiteSparse_long is long and compilers warn even when they are identical */
18 #define umfpack_UMF_wsolve(a,b,c,d,e,f,g,h,i,j,k,l,m,n) umfpack_zl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k,l,(SuiteSparse_long*)m,n)
19 #define umfpack_UMF_numeric(a,b,c,d,e,f,g,h) umfpack_zl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g,h)
20 #define umfpack_UMF_report_numeric umfpack_zl_report_numeric
21 #define umfpack_UMF_report_control umfpack_zl_report_control
22 #define umfpack_UMF_report_status umfpack_zl_report_status
23 #define umfpack_UMF_report_info umfpack_zl_report_info
24 #define umfpack_UMF_report_symbolic umfpack_zl_report_symbolic
25 #define umfpack_UMF_qsymbolic(a,b,c,d,e,f,g,h,i,j) umfpack_zl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,(SuiteSparse_long*)g,h,i,j)
26 #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h,i) umfpack_zl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h,i)
27 #define umfpack_UMF_defaults umfpack_zl_defaults
28
29 #else
30 #define umfpack_UMF_free_symbolic umfpack_dl_free_symbolic
31 #define umfpack_UMF_free_numeric umfpack_dl_free_numeric
32 #define umfpack_UMF_wsolve(a,b,c,d,e,f,g,h,i,j,k) umfpack_dl_wsolve(a,(SuiteSparse_long*)b,(SuiteSparse_long*)c,d,e,f,g,h,i,(SuiteSparse_long*)j,k)
33 #define umfpack_UMF_numeric(a,b,c,d,e,f,g) umfpack_dl_numeric((SuiteSparse_long*)a,(SuiteSparse_long*)b,c,d,e,f,g)
34 #define umfpack_UMF_report_numeric umfpack_dl_report_numeric
35 #define umfpack_UMF_report_control umfpack_dl_report_control
36 #define umfpack_UMF_report_status umfpack_dl_report_status
37 #define umfpack_UMF_report_info umfpack_dl_report_info
38 #define umfpack_UMF_report_symbolic umfpack_dl_report_symbolic
39 #define umfpack_UMF_qsymbolic(a,b,c,d,e,f,g,h,i) umfpack_dl_qsymbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,(SuiteSparse_long*)f,g,h,i)
40 #define umfpack_UMF_symbolic(a,b,c,d,e,f,g,h) umfpack_dl_symbolic(a,b,(SuiteSparse_long*)c,(SuiteSparse_long*)d,e,f,g,h)
41 #define umfpack_UMF_defaults umfpack_dl_defaults
42 #endif
43
44 #else
45 #if defined(PETSC_USE_COMPLEX)
46 #define umfpack_UMF_free_symbolic umfpack_zi_free_symbolic
47 #define umfpack_UMF_free_numeric umfpack_zi_free_numeric
48 #define umfpack_UMF_wsolve umfpack_zi_wsolve
49 #define umfpack_UMF_numeric umfpack_zi_numeric
50 #define umfpack_UMF_report_numeric umfpack_zi_report_numeric
51 #define umfpack_UMF_report_control umfpack_zi_report_control
52 #define umfpack_UMF_report_status umfpack_zi_report_status
53 #define umfpack_UMF_report_info umfpack_zi_report_info
54 #define umfpack_UMF_report_symbolic umfpack_zi_report_symbolic
55 #define umfpack_UMF_qsymbolic umfpack_zi_qsymbolic
56 #define umfpack_UMF_symbolic umfpack_zi_symbolic
57 #define umfpack_UMF_defaults umfpack_zi_defaults
58
59 #else
60 #define umfpack_UMF_free_symbolic umfpack_di_free_symbolic
61 #define umfpack_UMF_free_numeric umfpack_di_free_numeric
62 #define umfpack_UMF_wsolve umfpack_di_wsolve
63 #define umfpack_UMF_numeric umfpack_di_numeric
64 #define umfpack_UMF_report_numeric umfpack_di_report_numeric
65 #define umfpack_UMF_report_control umfpack_di_report_control
66 #define umfpack_UMF_report_status umfpack_di_report_status
67 #define umfpack_UMF_report_info umfpack_di_report_info
68 #define umfpack_UMF_report_symbolic umfpack_di_report_symbolic
69 #define umfpack_UMF_qsymbolic umfpack_di_qsymbolic
70 #define umfpack_UMF_symbolic umfpack_di_symbolic
71 #define umfpack_UMF_defaults umfpack_di_defaults
72 #endif
73 #endif
74
75 EXTERN_C_BEGIN
76 #include <umfpack.h>
77 EXTERN_C_END
78
79 static const char *const UmfpackOrderingTypes[] = {"CHOLMOD","AMD","GIVEN","METIS","BEST","NONE","USER","UmfpackOrderingTypes","UMFPACK_ORDERING_",0};
80
81 typedef struct {
82 void *Symbolic, *Numeric;
83 double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL],*W;
84 PetscInt *Wi,*perm_c;
85 Mat A; /* Matrix used for factorization */
86 MatStructure flg;
87
88 /* Flag to clean up UMFPACK objects during Destroy */
89 PetscBool CleanUpUMFPACK;
90 } Mat_UMFPACK;
91
MatDestroy_UMFPACK(Mat A)92 static PetscErrorCode MatDestroy_UMFPACK(Mat A)
93 {
94 PetscErrorCode ierr;
95 Mat_UMFPACK *lu=(Mat_UMFPACK*)A->data;
96
97 PetscFunctionBegin;
98 if (lu->CleanUpUMFPACK) {
99 umfpack_UMF_free_symbolic(&lu->Symbolic);
100 umfpack_UMF_free_numeric(&lu->Numeric);
101 ierr = PetscFree(lu->Wi);CHKERRQ(ierr);
102 ierr = PetscFree(lu->W);CHKERRQ(ierr);
103 ierr = PetscFree(lu->perm_c);CHKERRQ(ierr);
104 }
105 ierr = MatDestroy(&lu->A);CHKERRQ(ierr);
106 ierr = PetscFree(A->data);CHKERRQ(ierr);
107 PetscFunctionReturn(0);
108 }
109
MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag)110 static PetscErrorCode MatSolve_UMFPACK_Private(Mat A,Vec b,Vec x,int uflag)
111 {
112 Mat_UMFPACK *lu = (Mat_UMFPACK*)A->data;
113 Mat_SeqAIJ *a = (Mat_SeqAIJ*)lu->A->data;
114 PetscScalar *av = a->a,*xa;
115 const PetscScalar *ba;
116 PetscErrorCode ierr;
117 PetscInt *ai = a->i,*aj = a->j,status;
118 static PetscBool cite = PETSC_FALSE;
119
120 PetscFunctionBegin;
121 if (!A->rmap->n) PetscFunctionReturn(0);
122 ierr = PetscCitationsRegister("@article{davis2004algorithm,\n title={Algorithm 832: {UMFPACK} V4.3---An Unsymmetric-Pattern Multifrontal Method},\n author={Davis, Timothy A},\n journal={ACM Transactions on Mathematical Software (TOMS)},\n volume={30},\n number={2},\n pages={196--199},\n year={2004},\n publisher={ACM}\n}\n",&cite);CHKERRQ(ierr);
123 /* solve Ax = b by umfpack_*_wsolve */
124 /* ----------------------------------*/
125
126 if (!lu->Wi) { /* first time, allocate working space for wsolve */
127 ierr = PetscMalloc1(A->rmap->n,&lu->Wi);CHKERRQ(ierr);
128 ierr = PetscMalloc1(5*A->rmap->n,&lu->W);CHKERRQ(ierr);
129 }
130
131 ierr = VecGetArrayRead(b,&ba);
132 ierr = VecGetArray(x,&xa);
133 #if defined(PETSC_USE_COMPLEX)
134 status = umfpack_UMF_wsolve(uflag,ai,aj,(PetscReal*)av,NULL,(PetscReal*)xa,NULL,(PetscReal*)ba,NULL,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
135 #else
136 status = umfpack_UMF_wsolve(uflag,ai,aj,av,xa,ba,lu->Numeric,lu->Control,lu->Info,lu->Wi,lu->W);
137 #endif
138 umfpack_UMF_report_info(lu->Control, lu->Info);
139 if (status < 0) {
140 umfpack_UMF_report_status(lu->Control, status);
141 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_wsolve failed");
142 }
143
144 ierr = VecRestoreArrayRead(b,&ba);CHKERRQ(ierr);
145 ierr = VecRestoreArray(x,&xa);CHKERRQ(ierr);
146 PetscFunctionReturn(0);
147 }
148
MatSolve_UMFPACK(Mat A,Vec b,Vec x)149 static PetscErrorCode MatSolve_UMFPACK(Mat A,Vec b,Vec x)
150 {
151 PetscErrorCode ierr;
152
153 PetscFunctionBegin;
154 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
155 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_Aat);CHKERRQ(ierr);
156 PetscFunctionReturn(0);
157 }
158
MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)159 static PetscErrorCode MatSolveTranspose_UMFPACK(Mat A,Vec b,Vec x)
160 {
161 PetscErrorCode ierr;
162
163 PetscFunctionBegin;
164 /* We gave UMFPACK the algebraic transpose (because it assumes column alignment) */
165 ierr = MatSolve_UMFPACK_Private(A,b,x,UMFPACK_A);CHKERRQ(ierr);
166 PetscFunctionReturn(0);
167 }
168
MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo * info)169 static PetscErrorCode MatLUFactorNumeric_UMFPACK(Mat F,Mat A,const MatFactorInfo *info)
170 {
171 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F)->data;
172 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
173 PetscInt *ai = a->i,*aj=a->j,status;
174 PetscScalar *av = a->a;
175 PetscErrorCode ierr;
176
177 PetscFunctionBegin;
178 if (!A->rmap->n) PetscFunctionReturn(0);
179 /* numeric factorization of A' */
180 /* ----------------------------*/
181
182 if (lu->flg == SAME_NONZERO_PATTERN && lu->Numeric) {
183 umfpack_UMF_free_numeric(&lu->Numeric);
184 }
185 #if defined(PETSC_USE_COMPLEX)
186 status = umfpack_UMF_numeric(ai,aj,(double*)av,NULL,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
187 #else
188 status = umfpack_UMF_numeric(ai,aj,av,lu->Symbolic,&lu->Numeric,lu->Control,lu->Info);
189 #endif
190 if (status < 0) {
191 umfpack_UMF_report_status(lu->Control, status);
192 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_numeric failed");
193 }
194 /* report numeric factorization of A' when Control[PRL] > 3 */
195 (void) umfpack_UMF_report_numeric(lu->Numeric, lu->Control);
196
197 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
198 ierr = MatDestroy(&lu->A);CHKERRQ(ierr);
199
200 lu->A = A;
201 lu->flg = SAME_NONZERO_PATTERN;
202 lu->CleanUpUMFPACK = PETSC_TRUE;
203 F->ops->solve = MatSolve_UMFPACK;
204 F->ops->solvetranspose = MatSolveTranspose_UMFPACK;
205 PetscFunctionReturn(0);
206 }
207
208 /*
209 Note the r permutation is ignored
210 */
MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)211 static PetscErrorCode MatLUFactorSymbolic_UMFPACK(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
212 {
213 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
214 Mat_UMFPACK *lu = (Mat_UMFPACK*)(F->data);
215 PetscErrorCode ierr;
216 PetscInt i,*ai = a->i,*aj = a->j,m=A->rmap->n,n=A->cmap->n;
217 #if !defined(PETSC_USE_COMPLEX)
218 PetscScalar *av = a->a;
219 #endif
220 const PetscInt *ra;
221 PetscInt status;
222
223 PetscFunctionBegin;
224 (F)->ops->lufactornumeric = MatLUFactorNumeric_UMFPACK;
225 if (!n) PetscFunctionReturn(0);
226 if (r) {
227 ierr = ISGetIndices(r,&ra);CHKERRQ(ierr);
228 ierr = PetscMalloc1(m,&lu->perm_c);CHKERRQ(ierr);
229 /* we cannot simply memcpy on 64 bit archs */
230 for (i = 0; i < m; i++) lu->perm_c[i] = ra[i];
231 ierr = ISRestoreIndices(r,&ra);CHKERRQ(ierr);
232 }
233
234 /* print the control parameters */
235 if (lu->Control[UMFPACK_PRL] > 1) umfpack_UMF_report_control(lu->Control);
236
237 /* symbolic factorization of A' */
238 /* ---------------------------------------------------------------------- */
239 if (r) { /* use Petsc row ordering */
240 #if !defined(PETSC_USE_COMPLEX)
241 status = umfpack_UMF_qsymbolic(n,m,ai,aj,av,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
242 #else
243 status = umfpack_UMF_qsymbolic(n,m,ai,aj,NULL,NULL,lu->perm_c,&lu->Symbolic,lu->Control,lu->Info);
244 #endif
245 } else { /* use Umfpack col ordering */
246 #if !defined(PETSC_USE_COMPLEX)
247 status = umfpack_UMF_symbolic(n,m,ai,aj,av,&lu->Symbolic,lu->Control,lu->Info);
248 #else
249 status = umfpack_UMF_symbolic(n,m,ai,aj,NULL,NULL,&lu->Symbolic,lu->Control,lu->Info);
250 #endif
251 }
252 if (status < 0) {
253 umfpack_UMF_report_info(lu->Control, lu->Info);
254 umfpack_UMF_report_status(lu->Control, status);
255 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"umfpack_UMF_symbolic failed");
256 }
257 /* report sumbolic factorization of A' when Control[PRL] > 3 */
258 (void) umfpack_UMF_report_symbolic(lu->Symbolic, lu->Control);
259
260 lu->flg = DIFFERENT_NONZERO_PATTERN;
261 lu->CleanUpUMFPACK = PETSC_TRUE;
262 PetscFunctionReturn(0);
263 }
264
MatView_Info_UMFPACK(Mat A,PetscViewer viewer)265 static PetscErrorCode MatView_Info_UMFPACK(Mat A,PetscViewer viewer)
266 {
267 Mat_UMFPACK *lu= (Mat_UMFPACK*)A->data;
268 PetscErrorCode ierr;
269
270 PetscFunctionBegin;
271 /* check if matrix is UMFPACK type */
272 if (A->ops->solve != MatSolve_UMFPACK) PetscFunctionReturn(0);
273
274 ierr = PetscViewerASCIIPrintf(viewer,"UMFPACK run parameters:\n");CHKERRQ(ierr);
275 /* Control parameters used by reporting routiones */
276 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PRL]: %g\n",lu->Control[UMFPACK_PRL]);CHKERRQ(ierr);
277
278 /* Control parameters used by symbolic factorization */
279 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_STRATEGY]: %g\n",lu->Control[UMFPACK_STRATEGY]);CHKERRQ(ierr);
280 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_COL]: %g\n",lu->Control[UMFPACK_DENSE_COL]);CHKERRQ(ierr);
281 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DENSE_ROW]: %g\n",lu->Control[UMFPACK_DENSE_ROW]);CHKERRQ(ierr);
282 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AMD_DENSE]: %g\n",lu->Control[UMFPACK_AMD_DENSE]);CHKERRQ(ierr);
283 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_BLOCK_SIZE]: %g\n",lu->Control[UMFPACK_BLOCK_SIZE]);CHKERRQ(ierr);
284 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_FIXQ]: %g\n",lu->Control[UMFPACK_FIXQ]);CHKERRQ(ierr);
285 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_AGGRESSIVE]: %g\n",lu->Control[UMFPACK_AGGRESSIVE]);CHKERRQ(ierr);
286
287 /* Control parameters used by numeric factorization */
288 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_PIVOT_TOLERANCE]);CHKERRQ(ierr);
289 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SYM_PIVOT_TOLERANCE]: %g\n",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE]);CHKERRQ(ierr);
290 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_SCALE]: %g\n",lu->Control[UMFPACK_SCALE]);CHKERRQ(ierr);
291 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ALLOC_INIT]: %g\n",lu->Control[UMFPACK_ALLOC_INIT]);CHKERRQ(ierr);
292 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_DROPTOL]: %g\n",lu->Control[UMFPACK_DROPTOL]);CHKERRQ(ierr);
293
294 /* Control parameters used by solve */
295 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_IRSTEP]: %g\n",lu->Control[UMFPACK_IRSTEP]);CHKERRQ(ierr);
296
297 /* mat ordering */
298 if (!lu->perm_c) {
299 ierr = PetscViewerASCIIPrintf(viewer," Control[UMFPACK_ORDERING]: AMD (not using the PETSc ordering)\n",UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]]);CHKERRQ(ierr);
300 }
301 PetscFunctionReturn(0);
302 }
303
MatView_UMFPACK(Mat A,PetscViewer viewer)304 static PetscErrorCode MatView_UMFPACK(Mat A,PetscViewer viewer)
305 {
306 PetscErrorCode ierr;
307 PetscBool iascii;
308 PetscViewerFormat format;
309
310 PetscFunctionBegin;
311 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
312 if (iascii) {
313 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
314 if (format == PETSC_VIEWER_ASCII_INFO) {
315 ierr = MatView_Info_UMFPACK(A,viewer);CHKERRQ(ierr);
316 }
317 }
318 PetscFunctionReturn(0);
319 }
320
MatFactorGetSolverType_seqaij_umfpack(Mat A,MatSolverType * type)321 PetscErrorCode MatFactorGetSolverType_seqaij_umfpack(Mat A,MatSolverType *type)
322 {
323 PetscFunctionBegin;
324 *type = MATSOLVERUMFPACK;
325 PetscFunctionReturn(0);
326 }
327
328
329 /*MC
330 MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
331 via the external package UMFPACK.
332
333 Use ./configure --download-suitesparse to install PETSc to use UMFPACK
334
335 Use -pc_type lu -pc_factor_mat_solver_type umfpack to use this direct solver
336
337 Consult UMFPACK documentation for more information about the Control parameters
338 which correspond to the options database keys below.
339
340 Options Database Keys:
341 + -mat_umfpack_ordering - CHOLMOD, AMD, GIVEN, METIS, BEST, NONE
342 . -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
343 . -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
344 . -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
345 . -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW]
346 . -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE]
347 . -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
348 . -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE]
349 . -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ]
350 . -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE]
351 . -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
352 . -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
353 . -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX
354 . -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
355 . -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL]
356 - -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
357
358 Level: beginner
359
360 Note: UMFPACK is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html
361
362 .seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
363 M*/
364
MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat * F)365 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
366 {
367 Mat B;
368 Mat_UMFPACK *lu;
369 PetscErrorCode ierr;
370 PetscInt m=A->rmap->n,n=A->cmap->n,idx;
371
372 const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC"};
373 const char *scale[] ={"NONE","SUM","MAX"};
374 PetscBool flg;
375
376 PetscFunctionBegin;
377 /* Create the factorization matrix F */
378 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr);
379 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
380 ierr = PetscStrallocpy("umfpack",&((PetscObject)B)->type_name);CHKERRQ(ierr);
381 ierr = MatSetUp(B);CHKERRQ(ierr);
382
383 ierr = PetscNewLog(B,&lu);CHKERRQ(ierr);
384
385 B->data = lu;
386 B->ops->getinfo = MatGetInfo_External;
387 B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
388 B->ops->destroy = MatDestroy_UMFPACK;
389 B->ops->view = MatView_UMFPACK;
390 B->ops->matsolve = NULL;
391
392 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_umfpack);CHKERRQ(ierr);
393
394 B->factortype = MAT_FACTOR_LU;
395 B->assembled = PETSC_TRUE; /* required by -ksp_view */
396 B->preallocated = PETSC_TRUE;
397
398 ierr = PetscFree(B->solvertype);CHKERRQ(ierr);
399 ierr = PetscStrallocpy(MATSOLVERUMFPACK,&B->solvertype);CHKERRQ(ierr);
400 B->useordering = PETSC_TRUE;
401
402 /* initializations */
403 /* ------------------------------------------------*/
404 /* get the default control parameters */
405 umfpack_UMF_defaults(lu->Control);
406 lu->perm_c = NULL; /* use defaul UMFPACK col permutation */
407 lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
408
409 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr);
410 /* Control parameters used by reporting routiones */
411 ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],NULL);CHKERRQ(ierr);
412
413 /* Control parameters for symbolic factorization */
414 ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg);CHKERRQ(ierr);
415 if (flg) {
416 switch (idx) {
417 case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
418 case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
419 case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
420 }
421 }
422 ierr = PetscOptionsEList("-mat_umfpack_ordering","Internal ordering method","None",UmfpackOrderingTypes,sizeof(UmfpackOrderingTypes)/sizeof(UmfpackOrderingTypes[0]),UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]],&idx,&flg);CHKERRQ(ierr);
423 if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
424 ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],NULL);CHKERRQ(ierr);
425 ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],NULL);CHKERRQ(ierr);
426 ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],NULL);CHKERRQ(ierr);
427 ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],NULL);CHKERRQ(ierr);
428 ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],NULL);CHKERRQ(ierr);
429 ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],NULL);CHKERRQ(ierr);
430
431 /* Control parameters used by numeric factorization */
432 ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],NULL);CHKERRQ(ierr);
433 ierr = PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],NULL);CHKERRQ(ierr);
434 ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr);
435 if (flg) {
436 switch (idx) {
437 case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
438 case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
439 case 2: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_MAX; break;
440 }
441 }
442 ierr = PetscOptionsReal("-mat_umfpack_alloc_init","Control[UMFPACK_ALLOC_INIT]","None",lu->Control[UMFPACK_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);CHKERRQ(ierr);
443 ierr = PetscOptionsReal("-mat_umfpack_front_alloc_init","Control[UMFPACK_FRONT_ALLOC_INIT]","None",lu->Control[UMFPACK_FRONT_ALLOC_INIT],&lu->Control[UMFPACK_ALLOC_INIT],NULL);CHKERRQ(ierr);
444 ierr = PetscOptionsReal("-mat_umfpack_droptol","Control[UMFPACK_DROPTOL]","None",lu->Control[UMFPACK_DROPTOL],&lu->Control[UMFPACK_DROPTOL],NULL);CHKERRQ(ierr);
445
446 /* Control parameters used by solve */
447 ierr = PetscOptionsReal("-mat_umfpack_irstep","Control[UMFPACK_IRSTEP]","None",lu->Control[UMFPACK_IRSTEP],&lu->Control[UMFPACK_IRSTEP],NULL);CHKERRQ(ierr);
448 PetscOptionsEnd();
449 *F = B;
450 PetscFunctionReturn(0);
451 }
452
453 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*);
454 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*);
455 PETSC_INTERN PetscErrorCode MatGetFactor_seqaij_klu(Mat,MatFactorType,Mat*);
456
MatSolverTypeRegister_SuiteSparse(void)457 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuiteSparse(void)
458 {
459 PetscErrorCode ierr;
460
461 PetscFunctionBegin;
462 ierr = MatSolverTypeRegister(MATSOLVERUMFPACK,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_umfpack);CHKERRQ(ierr);
463 ierr = MatSolverTypeRegister(MATSOLVERCHOLMOD,MATSEQAIJ, MAT_FACTOR_CHOLESKY,MatGetFactor_seqaij_cholmod);CHKERRQ(ierr);
464 ierr = MatSolverTypeRegister(MATSOLVERCHOLMOD,MATSEQSBAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr);
465 ierr = MatSolverTypeRegister(MATSOLVERKLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_klu);CHKERRQ(ierr);
466 PetscFunctionReturn(0);
467 }
468