Actual source code: superlu.c
petsc-3.3-p7 2013-05-11
2: /* --------------------------------------------------------------------
4: This file implements a subclass of the SeqAIJ matrix class that uses
5: the SuperLU sparse solver. You can use this as a starting point for
6: implementing your own subclass of a PETSc matrix class.
8: This demonstrates a way to make an implementation inheritence of a PETSc
9: matrix type. This means constructing a new matrix type (SuperLU) by changing some
10: of the methods of the previous type (SeqAIJ), adding additional data, and possibly
11: additional method. (See any book on object oriented programming).
12: */
14: /*
15: Defines the data structure for the base matrix type (SeqAIJ)
16: */
17: #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
19: /*
20: SuperLU include files
21: */
22: EXTERN_C_BEGIN
23: #if defined(PETSC_USE_COMPLEX)
24: #include <slu_zdefs.h>
25: #else
26: #include <slu_ddefs.h>
27: #endif
28: #include <slu_util.h>
29: EXTERN_C_END
31: /*
32: This is the data we are "ADDING" to the SeqAIJ matrix type to get the SuperLU matrix type
33: */
34: typedef struct {
35: SuperMatrix A,L,U,B,X;
36: superlu_options_t options;
37: PetscInt *perm_c; /* column permutation vector */
38: PetscInt *perm_r; /* row permutations from partial pivoting */
39: PetscInt *etree;
40: PetscReal *R, *C;
41: char equed[1];
42: PetscInt lwork;
43: void *work;
44: PetscReal rpg, rcond;
45: mem_usage_t mem_usage;
46: MatStructure flg;
47: SuperLUStat_t stat;
48: Mat A_dup;
49: PetscScalar *rhs_dup;
51: /* Flag to clean up (non-global) SuperLU objects during Destroy */
52: PetscBool CleanUpSuperLU;
53: } Mat_SuperLU;
55: extern PetscErrorCode MatFactorInfo_SuperLU(Mat,PetscViewer);
56: extern PetscErrorCode MatLUFactorNumeric_SuperLU(Mat,Mat,const MatFactorInfo *);
57: extern PetscErrorCode MatDestroy_SuperLU(Mat);
58: extern PetscErrorCode MatView_SuperLU(Mat,PetscViewer);
59: extern PetscErrorCode MatAssemblyEnd_SuperLU(Mat,MatAssemblyType);
60: extern PetscErrorCode MatSolve_SuperLU(Mat,Vec,Vec);
61: extern PetscErrorCode MatMatSolve_SuperLU(Mat,Mat,Mat);
62: extern PetscErrorCode MatSolveTranspose_SuperLU(Mat,Vec,Vec);
63: extern PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat,Mat,IS,IS,const MatFactorInfo*);
64: extern PetscErrorCode MatDuplicate_SuperLU(Mat, MatDuplicateOption, Mat *);
66: /*
67: Utility function
68: */
71: PetscErrorCode MatFactorInfo_SuperLU(Mat A,PetscViewer viewer)
72: {
73: Mat_SuperLU *lu= (Mat_SuperLU*)A->spptr;
74: PetscErrorCode ierr;
75: superlu_options_t options;
78: /* check if matrix is superlu_dist type */
79: if (A->ops->solve != MatSolve_SuperLU) return(0);
81: options = lu->options;
82: PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
83: PetscViewerASCIIPrintf(viewer," Equil: %s\n",(options.Equil != NO) ? "YES": "NO");
84: PetscViewerASCIIPrintf(viewer," ColPerm: %D\n",options.ColPerm);
85: PetscViewerASCIIPrintf(viewer," IterRefine: %D\n",options.IterRefine);
86: PetscViewerASCIIPrintf(viewer," SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES": "NO");
87: PetscViewerASCIIPrintf(viewer," DiagPivotThresh: %g\n",options.DiagPivotThresh);
88: PetscViewerASCIIPrintf(viewer," PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES": "NO");
89: PetscViewerASCIIPrintf(viewer," ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES": "NO");
90: PetscViewerASCIIPrintf(viewer," RowPerm: %D\n",options.RowPerm);
91: PetscViewerASCIIPrintf(viewer," ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES": "NO");
92: PetscViewerASCIIPrintf(viewer," PrintStat: %s\n",(options.PrintStat != NO) ? "YES": "NO");
93: PetscViewerASCIIPrintf(viewer," lwork: %D\n",lu->lwork);
94: if (A->factortype == MAT_FACTOR_ILU){
95: PetscViewerASCIIPrintf(viewer," ILU_DropTol: %g\n",options.ILU_DropTol);
96: PetscViewerASCIIPrintf(viewer," ILU_FillTol: %g\n",options.ILU_FillTol);
97: PetscViewerASCIIPrintf(viewer," ILU_FillFactor: %g\n",options.ILU_FillFactor);
98: PetscViewerASCIIPrintf(viewer," ILU_DropRule: %D\n",options.ILU_DropRule);
99: PetscViewerASCIIPrintf(viewer," ILU_Norm: %D\n",options.ILU_Norm);
100: PetscViewerASCIIPrintf(viewer," ILU_MILU: %D\n",options.ILU_MILU);
101: }
102: return(0);
103: }
105: /*
106: These are the methods provided to REPLACE the corresponding methods of the
107: SeqAIJ matrix class. Other methods of SeqAIJ are not replaced
108: */
111: PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info)
112: {
113: Mat_SuperLU *lu = (Mat_SuperLU*)F->spptr;
114: Mat_SeqAIJ *aa;
116: PetscInt sinfo;
117: PetscReal ferr, berr;
118: NCformat *Ustore;
119: SCformat *Lstore;
120:
122: if (lu->flg == SAME_NONZERO_PATTERN){ /* successing numerical factorization */
123: lu->options.Fact = SamePattern;
124: /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
125: Destroy_SuperMatrix_Store(&lu->A);
126: if (lu->options.Equil){
127: MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN);
128: }
129: if ( lu->lwork >= 0 ) {
130: Destroy_SuperNode_Matrix(&lu->L);
131: Destroy_CompCol_Matrix(&lu->U);
132: lu->options.Fact = SamePattern;
133: }
134: }
136: /* Create the SuperMatrix for lu->A=A^T:
137: Since SuperLU likes column-oriented matrices,we pass it the transpose,
138: and then solve A^T X = B in MatSolve(). */
139: if (lu->options.Equil){
140: aa = (Mat_SeqAIJ*)(lu->A_dup)->data;
141: } else {
142: aa = (Mat_SeqAIJ*)(A)->data;
143: }
144: #if defined(PETSC_USE_COMPLEX)
145: zCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,
146: SLU_NC,SLU_Z,SLU_GE);
147: #else
148: dCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,
149: SLU_NC,SLU_D,SLU_GE);
150: #endif
152: /* Numerical factorization */
153: lu->B.ncol = 0; /* Indicate not to solve the system */
154: if (F->factortype == MAT_FACTOR_LU){
155: #if defined(PETSC_USE_COMPLEX)
156: zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
157: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
158: &lu->mem_usage, &lu->stat, &sinfo);
159: #else
160: dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
161: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
162: &lu->mem_usage, &lu->stat, &sinfo);
163: #endif
164: } else if (F->factortype == MAT_FACTOR_ILU){
165: /* Compute the incomplete factorization, condition number and pivot growth */
166: #if defined(PETSC_USE_COMPLEX)
167: zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
168: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
169: &lu->mem_usage, &lu->stat, &sinfo);
170: #else
171: dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
172: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
173: &lu->mem_usage, &lu->stat, &sinfo);
174: #endif
175: } else {
176: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
177: }
178: if ( !sinfo || sinfo == lu->A.ncol+1 ) {
179: if ( lu->options.PivotGrowth )
180: PetscPrintf(PETSC_COMM_SELF," Recip. pivot growth = %e\n", lu->rpg);
181: if ( lu->options.ConditionNumber )
182: PetscPrintf(PETSC_COMM_SELF," Recip. condition number = %e\n", lu->rcond);
183: } else if ( sinfo > 0 ){
184: if ( lu->lwork == -1 ) {
185: PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %D bytes\n", sinfo - lu->A.ncol);
186: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %D",sinfo);
187: } else { /* sinfo < 0 */
188: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", sinfo,-sinfo);
189: }
191: if ( lu->options.PrintStat ) {
192: PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
193: StatPrint(&lu->stat);
194: Lstore = (SCformat *) lu->L.Store;
195: Ustore = (NCformat *) lu->U.Store;
196: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor L = %D\n", Lstore->nnz);
197: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor U = %D\n", Ustore->nnz);
198: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in L+U = %D\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
199: PetscPrintf(PETSC_COMM_SELF," L\\U MB %.3f\ttotal MB needed %.3f\n",
200: lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6);
201: }
203: lu->flg = SAME_NONZERO_PATTERN;
204: F->ops->solve = MatSolve_SuperLU;
205: F->ops->solvetranspose = MatSolveTranspose_SuperLU;
206: F->ops->matsolve = MatMatSolve_SuperLU;
207: return(0);
208: }
212: PetscErrorCode MatDestroy_SuperLU(Mat A)
213: {
215: Mat_SuperLU *lu=(Mat_SuperLU*)A->spptr;
218: if (lu && lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
219: Destroy_SuperMatrix_Store(&lu->A);
220: Destroy_SuperMatrix_Store(&lu->B);
221: Destroy_SuperMatrix_Store(&lu->X);
222: StatFree(&lu->stat);
223: if (lu->lwork >= 0) {
224: Destroy_SuperNode_Matrix(&lu->L);
225: Destroy_CompCol_Matrix(&lu->U);
226: }
227: }
228: if (lu) {
229: PetscFree(lu->etree);
230: PetscFree(lu->perm_r);
231: PetscFree(lu->perm_c);
232: PetscFree(lu->R);
233: PetscFree(lu->C);
234: PetscFree(lu->rhs_dup);
235: MatDestroy(&lu->A_dup);
236: }
237: PetscFree(A->spptr);
239: /* clear composed functions */
240: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatFactorGetSolverPackage_C","",PETSC_NULL);
241: PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSuperluSetILUDropTol_C","",PETSC_NULL);
243: MatDestroy_SeqAIJ(A);
244: return(0);
245: }
249: PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
250: {
251: PetscErrorCode ierr;
252: PetscBool iascii;
253: PetscViewerFormat format;
256: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
257: if (iascii) {
258: PetscViewerGetFormat(viewer,&format);
259: if (format == PETSC_VIEWER_ASCII_INFO) {
260: MatFactorInfo_SuperLU(A,viewer);
261: }
262: }
263: return(0);
264: }
269: PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x)
270: {
271: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
272: PetscScalar *barray,*xarray;
274: PetscInt info,i,n=x->map->n;
275: PetscReal ferr,berr;
276:
278: if ( lu->lwork == -1 ) {
279: return(0);
280: }
282: lu->B.ncol = 1; /* Set the number of right-hand side */
283: if (lu->options.Equil && !lu->rhs_dup){
284: /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */
285: PetscMalloc(n*sizeof(PetscScalar),&lu->rhs_dup);
286: }
287: if (lu->options.Equil){
288: /* Copy b into rsh_dup */
289: VecGetArray(b,&barray);
290: PetscMemcpy(lu->rhs_dup,barray,n*sizeof(PetscScalar));
291: VecRestoreArray(b,&barray);
292: barray = lu->rhs_dup;
293: } else {
294: VecGetArray(b,&barray);
295: }
296: VecGetArray(x,&xarray);
298: #if defined(PETSC_USE_COMPLEX)
299: ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
300: ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
301: #else
302: ((DNformat*)lu->B.Store)->nzval = barray;
303: ((DNformat*)lu->X.Store)->nzval = xarray;
304: #endif
306: lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
307: if (A->factortype == MAT_FACTOR_LU){
308: #if defined(PETSC_USE_COMPLEX)
309: zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
310: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
311: &lu->mem_usage, &lu->stat, &info);
312: #else
313: dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
314: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
315: &lu->mem_usage, &lu->stat, &info);
316: #endif
317: } else if (A->factortype == MAT_FACTOR_ILU){
318: #if defined(PETSC_USE_COMPLEX)
319: zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
320: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
321: &lu->mem_usage, &lu->stat, &info);
322: #else
323: dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
324: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
325: &lu->mem_usage, &lu->stat, &info);
326: #endif
327: } else {
328: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
329: }
330: if (!lu->options.Equil){
331: VecRestoreArray(b,&barray);
332: }
333: VecRestoreArray(x,&xarray);
335: if ( !info || info == lu->A.ncol+1 ) {
336: if ( lu->options.IterRefine ) {
337: PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
338: PetscPrintf(PETSC_COMM_SELF," %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
339: for (i = 0; i < 1; ++i)
340: PetscPrintf(PETSC_COMM_SELF," %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr);
341: }
342: } else if ( info > 0 ){
343: if ( lu->lwork == -1 ) {
344: PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %D bytes\n", info - lu->A.ncol);
345: } else {
346: PetscPrintf(PETSC_COMM_SELF," Warning: gssvx() returns info %D\n",info);
347: }
348: } else if (info < 0){
349: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", info,-info);
350: }
352: if ( lu->options.PrintStat ) {
353: PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
354: StatPrint(&lu->stat);
355: }
356: return(0);
357: }
361: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
362: {
363: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
367: lu->options.Trans = TRANS;
368: MatSolve_SuperLU_Private(A,b,x);
369: return(0);
370: }
374: PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x)
375: {
376: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
380: lu->options.Trans = NOTRANS;
381: MatSolve_SuperLU_Private(A,b,x);
382: return(0);
383: }
387: PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X)
388: {
389: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
390: PetscBool flg;
394: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
395: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
396: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,PETSC_NULL);
397: if (!flg) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); lu->options.Trans = TRANS;
398: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet");
399: return(0);
400: }
402: /*
403: Note the r permutation is ignored
404: */
407: PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
408: {
409: Mat_SuperLU *lu = (Mat_SuperLU*)(F->spptr);
410:
412: lu->flg = DIFFERENT_NONZERO_PATTERN;
413: lu->CleanUpSuperLU = PETSC_TRUE;
414: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
415: return(0);
416: }
418: EXTERN_C_BEGIN
421: PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol)
422: {
423: Mat_SuperLU *lu= (Mat_SuperLU*)F->spptr;
426: lu->options.ILU_DropTol = dtol;
427: return(0);
428: }
429: EXTERN_C_END
433: /*@
434: MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
435: Logically Collective on Mat
437: Input Parameters:
438: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface
439: - dtol - drop tolerance
441: Options Database:
442: . -mat_superlu_ilu_droptol <dtol>
444: Level: beginner
446: References: SuperLU Users' Guide
448: .seealso: MatGetFactor()
449: @*/
450: PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol)
451: {
457: PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol));
458: return(0);
459: }
461: EXTERN_C_BEGIN
464: PetscErrorCode MatFactorGetSolverPackage_seqaij_superlu(Mat A,const MatSolverPackage *type)
465: {
467: *type = MATSOLVERSUPERLU;
468: return(0);
469: }
470: EXTERN_C_END
471:
473: /*MC
474: MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
475: via the external package SuperLU.
477: Use ./configure --download-superlu to have PETSc installed with SuperLU
479: Options Database Keys:
480: + -mat_superlu_equil <FALSE> - Equil (None)
481: . -mat_superlu_colperm <COLAMD> - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD
482: . -mat_superlu_iterrefine <NOREFINE> - (choose one of) NOREFINE SINGLE DOUBLE EXTRA
483: . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None)
484: . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None)
485: . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None)
486: . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None)
487: . -mat_superlu_rowperm <NOROWPERM> - (choose one of) NOROWPERM LargeDiag
488: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
489: . -mat_superlu_printstat <FALSE> - PrintStat (None)
490: . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None)
491: . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None)
492: . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None)
493: . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None)
494: . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None)
495: . -mat_superlu_ilu_norm <0> - ILU_Norm (None)
496: - -mat_superlu_ilu_milu <0> - ILU_MILU (None)
498: Notes: Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves
500: Level: beginner
502: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, MATSOLVERSPOOLES, PCFactorSetMatSolverPackage(), MatSolverPackage
503: M*/
505: EXTERN_C_BEGIN
508: PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F)
509: {
510: Mat B;
511: Mat_SuperLU *lu;
513: PetscInt indx,m=A->rmap->n,n=A->cmap->n;
514: PetscBool flg;
515: const char *colperm[]={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
516: const char *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
517: const char *rowperm[]={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */
520: MatCreate(((PetscObject)A)->comm,&B);
521: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);
522: MatSetType(B,((PetscObject)A)->type_name);
523: MatSeqAIJSetPreallocation(B,0,PETSC_NULL);
525: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU){
526: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
527: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
528: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
530: B->ops->destroy = MatDestroy_SuperLU;
531: B->ops->view = MatView_SuperLU;
532: B->factortype = ftype;
533: B->assembled = PETSC_TRUE; /* required by -ksp_view */
534: B->preallocated = PETSC_TRUE;
535:
536: PetscNewLog(B,Mat_SuperLU,&lu);
537:
538: if (ftype == MAT_FACTOR_LU){
539: set_default_options(&lu->options);
540: /* Comments from SuperLU_4.0/SRC/dgssvx.c:
541: "Whether or not the system will be equilibrated depends on the
542: scaling of the matrix A, but if equilibration is used, A is
543: overwritten by diag(R)*A*diag(C) and B by diag(R)*B
544: (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
545: We set 'options.Equil = NO' as default because additional space is needed for it.
546: */
547: lu->options.Equil = NO;
548: } else if (ftype == MAT_FACTOR_ILU){
549: /* Set the default input options of ilu: */
550: ilu_set_default_options(&lu->options);
551: }
552: lu->options.PrintStat = NO;
553:
554: /* Initialize the statistics variables. */
555: StatInit(&lu->stat);
556: lu->lwork = 0; /* allocate space internally by system malloc */
558: PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"SuperLU Options","Mat");
559: PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,0);
560: PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
561: if (flg) {lu->options.ColPerm = (colperm_t)indx;}
562: PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
563: if (flg) { lu->options.IterRefine = (IterRefine_t)indx;}
564: PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,0);
565: if (flg) lu->options.SymmetricMode = YES;
566: PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&lu->options.DiagPivotThresh,PETSC_NULL);
567: PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,0);
568: if (flg) lu->options.PivotGrowth = YES;
569: PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,0);
570: if (flg) lu->options.ConditionNumber = YES;
571: PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg);
572: if (flg) {lu->options.RowPerm = (rowperm_t)indx;}
573: PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,0);
574: if (flg) lu->options.ReplaceTinyPivot = YES;
575: PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,0);
576: if (flg) lu->options.PrintStat = YES;
577: PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,PETSC_NULL);
578: if (lu->lwork > 0 ){
579: PetscMalloc(lu->lwork,&lu->work);
580: } else if (lu->lwork != 0 && lu->lwork != -1){
581: PetscPrintf(PETSC_COMM_SELF," Warning: lwork %D is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
582: lu->lwork = 0;
583: }
584: /* ilu options */
585: PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&lu->options.ILU_DropTol,PETSC_NULL);
586: PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&lu->options.ILU_FillTol,PETSC_NULL);
587: PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&lu->options.ILU_FillFactor,PETSC_NULL);
588: PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,PETSC_NULL);
589: PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg);
590: if (flg){
591: lu->options.ILU_Norm = (norm_t)indx;
592: }
593: PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg);
594: if (flg){
595: lu->options.ILU_MILU = (milu_t)indx;
596: }
597: PetscOptionsEnd();
598: if (lu->options.Equil == YES) {
599: /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */
600: MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup);
601: }
603: /* Allocate spaces (notice sizes are for the transpose) */
604: PetscMalloc(m*sizeof(PetscInt),&lu->etree);
605: PetscMalloc(n*sizeof(PetscInt),&lu->perm_r);
606: PetscMalloc(m*sizeof(PetscInt),&lu->perm_c);
607: PetscMalloc(n*sizeof(PetscScalar),&lu->R);
608: PetscMalloc(m*sizeof(PetscScalar),&lu->C);
609:
610: /* create rhs and solution x without allocate space for .Store */
611: #if defined(PETSC_USE_COMPLEX)
612: zCreate_Dense_Matrix(&lu->B, m, 1, PETSC_NULL, m, SLU_DN, SLU_Z, SLU_GE);
613: zCreate_Dense_Matrix(&lu->X, m, 1, PETSC_NULL, m, SLU_DN, SLU_Z, SLU_GE);
614: #else
615: dCreate_Dense_Matrix(&lu->B, m, 1, PETSC_NULL, m, SLU_DN, SLU_D, SLU_GE);
616: dCreate_Dense_Matrix(&lu->X, m, 1, PETSC_NULL, m, SLU_DN, SLU_D, SLU_GE);
617: #endif
619: #ifdef SUPERLU2
620: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatCreateNull","MatCreateNull_SuperLU",(void(*)(void))MatCreateNull_SuperLU);
621: #endif
622: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_superlu",MatFactorGetSolverPackage_seqaij_superlu);
623: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSuperluSetILUDropTol_C","MatSuperluSetILUDropTol_SuperLU",MatSuperluSetILUDropTol_SuperLU);
624: B->spptr = lu;
625: *F = B;
626: return(0);
627: }
628: EXTERN_C_END