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: }