Actual source code: superlu.c
petsc-3.7.5 2017-01-01
2: /* --------------------------------------------------------------------
4: This file implements a subclass of the SeqAIJ matrix class that uses
5: the SuperLU sparse solver.
6: */
8: /*
9: Defines the data structure for the base matrix type (SeqAIJ)
10: */
11: #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
13: /*
14: SuperLU include files
15: */
16: EXTERN_C_BEGIN
17: #if defined(PETSC_USE_COMPLEX)
18: #if defined(PETSC_USE_REAL_SINGLE)
19: #include <slu_cdefs.h>
20: #else
21: #include <slu_zdefs.h>
22: #endif
23: #else
24: #if defined(PETSC_USE_REAL_SINGLE)
25: #include <slu_sdefs.h>
26: #else
27: #include <slu_ddefs.h>
28: #endif
29: #endif
30: #include <slu_util.h>
31: EXTERN_C_END
33: /*
34: This is the data we are "ADDING" to the SeqAIJ matrix type to get the SuperLU matrix type
35: */
36: typedef struct {
37: SuperMatrix A,L,U,B,X;
38: superlu_options_t options;
39: PetscInt *perm_c; /* column permutation vector */
40: PetscInt *perm_r; /* row permutations from partial pivoting */
41: PetscInt *etree;
42: PetscReal *R, *C;
43: char equed[1];
44: PetscInt lwork;
45: void *work;
46: PetscReal rpg, rcond;
47: mem_usage_t mem_usage;
48: MatStructure flg;
49: SuperLUStat_t stat;
50: Mat A_dup;
51: PetscScalar *rhs_dup;
52: GlobalLU_t Glu;
54: /* Flag to clean up (non-global) SuperLU objects during Destroy */
55: PetscBool CleanUpSuperLU;
56: } Mat_SuperLU;
58: extern PetscErrorCode MatFactorInfo_SuperLU(Mat,PetscViewer);
59: extern PetscErrorCode MatLUFactorNumeric_SuperLU(Mat,Mat,const MatFactorInfo*);
60: extern PetscErrorCode MatDestroy_SuperLU(Mat);
61: extern PetscErrorCode MatView_SuperLU(Mat,PetscViewer);
62: extern PetscErrorCode MatAssemblyEnd_SuperLU(Mat,MatAssemblyType);
63: extern PetscErrorCode MatSolve_SuperLU(Mat,Vec,Vec);
64: extern PetscErrorCode MatMatSolve_SuperLU(Mat,Mat,Mat);
65: extern PetscErrorCode MatSolveTranspose_SuperLU(Mat,Vec,Vec);
66: extern PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat,Mat,IS,IS,const MatFactorInfo*);
67: extern PetscErrorCode MatDuplicate_SuperLU(Mat, MatDuplicateOption, Mat*);
69: /*
70: Utility function
71: */
74: PetscErrorCode MatFactorInfo_SuperLU(Mat A,PetscViewer viewer)
75: {
76: Mat_SuperLU *lu= (Mat_SuperLU*)A->spptr;
77: PetscErrorCode ierr;
78: superlu_options_t options;
81: /* check if matrix is superlu_dist type */
82: if (A->ops->solve != MatSolve_SuperLU) return(0);
84: options = lu->options;
86: PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
87: PetscViewerASCIIPrintf(viewer," Equil: %s\n",(options.Equil != NO) ? "YES" : "NO");
88: PetscViewerASCIIPrintf(viewer," ColPerm: %D\n",options.ColPerm);
89: PetscViewerASCIIPrintf(viewer," IterRefine: %D\n",options.IterRefine);
90: PetscViewerASCIIPrintf(viewer," SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES" : "NO");
91: PetscViewerASCIIPrintf(viewer," DiagPivotThresh: %g\n",options.DiagPivotThresh);
92: PetscViewerASCIIPrintf(viewer," PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES" : "NO");
93: PetscViewerASCIIPrintf(viewer," ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES" : "NO");
94: PetscViewerASCIIPrintf(viewer," RowPerm: %D\n",options.RowPerm);
95: PetscViewerASCIIPrintf(viewer," ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES" : "NO");
96: PetscViewerASCIIPrintf(viewer," PrintStat: %s\n",(options.PrintStat != NO) ? "YES" : "NO");
97: PetscViewerASCIIPrintf(viewer," lwork: %D\n",lu->lwork);
98: if (A->factortype == MAT_FACTOR_ILU) {
99: PetscViewerASCIIPrintf(viewer," ILU_DropTol: %g\n",options.ILU_DropTol);
100: PetscViewerASCIIPrintf(viewer," ILU_FillTol: %g\n",options.ILU_FillTol);
101: PetscViewerASCIIPrintf(viewer," ILU_FillFactor: %g\n",options.ILU_FillFactor);
102: PetscViewerASCIIPrintf(viewer," ILU_DropRule: %D\n",options.ILU_DropRule);
103: PetscViewerASCIIPrintf(viewer," ILU_Norm: %D\n",options.ILU_Norm);
104: PetscViewerASCIIPrintf(viewer," ILU_MILU: %D\n",options.ILU_MILU);
105: }
106: return(0);
107: }
109: /*
110: These are the methods provided to REPLACE the corresponding methods of the
111: SeqAIJ matrix class. Other methods of SeqAIJ are not replaced
112: */
115: PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info)
116: {
117: Mat_SuperLU *lu = (Mat_SuperLU*)F->spptr;
118: Mat_SeqAIJ *aa;
120: PetscInt sinfo;
121: PetscReal ferr, berr;
122: NCformat *Ustore;
123: SCformat *Lstore;
126: if (lu->flg == SAME_NONZERO_PATTERN) { /* successing numerical factorization */
127: lu->options.Fact = SamePattern;
128: /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
129: Destroy_SuperMatrix_Store(&lu->A);
130: if (lu->options.Equil) {
131: MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN);
132: }
133: if (lu->lwork >= 0) {
134: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
135: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
136: lu->options.Fact = SamePattern;
137: }
138: }
140: /* Create the SuperMatrix for lu->A=A^T:
141: Since SuperLU likes column-oriented matrices,we pass it the transpose,
142: and then solve A^T X = B in MatSolve(). */
143: if (lu->options.Equil) {
144: aa = (Mat_SeqAIJ*)(lu->A_dup)->data;
145: } else {
146: aa = (Mat_SeqAIJ*)(A)->data;
147: }
148: #if defined(PETSC_USE_COMPLEX)
149: #if defined(PETSC_USE_REAL_SINGLE)
150: PetscStackCall("SuperLU:cCreate_CompCol_Matrix",cCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(singlecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_C,SLU_GE));
151: #else
152: PetscStackCall("SuperLU:zCreate_CompCol_Matrix",zCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,(doublecomplex*)aa->a,aa->j,aa->i,SLU_NC,SLU_Z,SLU_GE));
153: #endif
154: #else
155: #if defined(PETSC_USE_REAL_SINGLE)
156: PetscStackCall("SuperLU:sCreate_CompCol_Matrix",sCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_S,SLU_GE));
157: #else
158: PetscStackCall("SuperLU:dCreate_CompCol_Matrix",dCreate_CompCol_Matrix(&lu->A,A->cmap->n,A->rmap->n,aa->nz,aa->a,aa->j,aa->i,SLU_NC,SLU_D,SLU_GE));
159: #endif
160: #endif
162: /* Numerical factorization */
163: lu->B.ncol = 0; /* Indicate not to solve the system */
164: if (F->factortype == MAT_FACTOR_LU) {
165: #if defined(PETSC_USE_COMPLEX)
166: #if defined(PETSC_USE_REAL_SINGLE)
167: PetscStackCall("SuperLU:cgssvx",cgssvx(&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, &ferr, &berr,
169: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
170: #else
171: PetscStackCall("SuperLU:zgssvx",zgssvx(&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, &ferr, &berr,
173: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
174: #endif
175: #else
176: #if defined(PETSC_USE_REAL_SINGLE)
177: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
178: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
179: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
180: #else
181: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
182: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
183: &lu->Glu,&lu->mem_usage, &lu->stat, &sinfo));
184: #endif
185: #endif
186: } else if (F->factortype == MAT_FACTOR_ILU) {
187: /* Compute the incomplete factorization, condition number and pivot growth */
188: #if defined(PETSC_USE_COMPLEX)
189: #if defined(PETSC_USE_REAL_SINGLE)
190: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
191: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
192: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
193: #else
194: PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
195: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
196: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
197: #endif
198: #else
199: #if defined(PETSC_USE_REAL_SINGLE)
200: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
201: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
202: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
203: #else
204: PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
205: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
206: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
207: #endif
208: #endif
209: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
210: if (!sinfo || sinfo == lu->A.ncol+1) {
211: if (lu->options.PivotGrowth) {
212: PetscPrintf(PETSC_COMM_SELF," Recip. pivot growth = %e\n", lu->rpg);
213: }
214: if (lu->options.ConditionNumber) {
215: PetscPrintf(PETSC_COMM_SELF," Recip. condition number = %e\n", lu->rcond);
216: }
217: } else if (sinfo > 0) {
218: if (A->erroriffailure) {
219: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %D",sinfo);
220: } else {
221: if (sinfo <= lu->A.ncol) {
222: if (lu->options.ILU_FillTol == 0.0) {
223: F->errortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
224: }
225: PetscInfo2(F,"Number of zero pivots %D, ILU_FillTol %g\n",sinfo,lu->options.ILU_FillTol);
226: } else if (sinfo == lu->A.ncol + 1) {
227: /*
228: U is nonsingular, but RCOND is less than machine
229: precision, meaning that the matrix is singular to
230: working precision. Nevertheless, the solution and
231: error bounds are computed because there are a number
232: of situations where the computed solution can be more
233: accurate than the value of RCOND would suggest.
234: */
235: PetscInfo1(F,"Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %D",sinfo);
236: } else { /* sinfo > lu->A.ncol + 1 */
237: F->errortype = MAT_FACTOR_OUTMEMORY;
238: PetscInfo1(F,"Number of bytes allocated when memory allocation fails %D\n",sinfo);
239: }
240: }
241: } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", sinfo,-sinfo);
243: if (lu->options.PrintStat) {
244: PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
245: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
246: Lstore = (SCformat*) lu->L.Store;
247: Ustore = (NCformat*) lu->U.Store;
248: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor L = %D\n", Lstore->nnz);
249: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor U = %D\n", Ustore->nnz);
250: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in L+U = %D\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
251: PetscPrintf(PETSC_COMM_SELF," L\\U MB %.3f\ttotal MB needed %.3f\n",
252: lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6);
253: }
255: lu->flg = SAME_NONZERO_PATTERN;
256: F->ops->solve = MatSolve_SuperLU;
257: F->ops->solvetranspose = MatSolveTranspose_SuperLU;
258: F->ops->matsolve = NULL;
259: return(0);
260: }
264: PetscErrorCode MatGetDiagonal_SuperLU(Mat A,Vec v)
265: {
267: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Mat type: SuperLU factor");
268: return(0);
269: }
273: PetscErrorCode MatDestroy_SuperLU(Mat A)
274: {
276: Mat_SuperLU *lu=(Mat_SuperLU*)A->spptr;
279: if (lu && lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
280: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->A));
281: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->B));
282: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->X));
283: PetscStackCall("SuperLU:StatFree",StatFree(&lu->stat));
284: if (lu->lwork >= 0) {
285: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
286: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
287: }
288: }
289: if (lu) {
290: PetscFree(lu->etree);
291: PetscFree(lu->perm_r);
292: PetscFree(lu->perm_c);
293: PetscFree(lu->R);
294: PetscFree(lu->C);
295: PetscFree(lu->rhs_dup);
296: MatDestroy(&lu->A_dup);
297: }
298: PetscFree(A->spptr);
300: /* clear composed functions */
301: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
302: PetscObjectComposeFunction((PetscObject)A,"MatSuperluSetILUDropTol_C",NULL);
304: MatDestroy_SeqAIJ(A);
305: return(0);
306: }
310: PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
311: {
312: PetscErrorCode ierr;
313: PetscBool iascii;
314: PetscViewerFormat format;
317: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
318: if (iascii) {
319: PetscViewerGetFormat(viewer,&format);
320: if (format == PETSC_VIEWER_ASCII_INFO) {
321: MatFactorInfo_SuperLU(A,viewer);
322: }
323: }
324: return(0);
325: }
330: PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x)
331: {
332: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
333: const PetscScalar *barray;
334: PetscScalar *xarray;
335: PetscErrorCode ierr;
336: PetscInt info,i,n;
337: PetscReal ferr,berr;
338: static PetscBool cite = PETSC_FALSE;
341: if (lu->lwork == -1) return(0);
342: PetscCitationsRegister("@article{superlu99,\n author = {James W. Demmel and Stanley C. Eisenstat and\n John R. Gilbert and Xiaoye S. Li and Joseph W. H. Liu},\n title = {A supernodal approach to sparse partial pivoting},\n journal = {SIAM J. Matrix Analysis and Applications},\n year = {1999},\n volume = {20},\n number = {3},\n pages = {720-755}\n}\n",&cite);
344: VecGetLocalSize(x,&n);
345: lu->B.ncol = 1; /* Set the number of right-hand side */
346: if (lu->options.Equil && !lu->rhs_dup) {
347: /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */
348: PetscMalloc1(n,&lu->rhs_dup);
349: }
350: if (lu->options.Equil) {
351: /* Copy b into rsh_dup */
352: VecGetArrayRead(b,&barray);
353: PetscMemcpy(lu->rhs_dup,barray,n*sizeof(PetscScalar));
354: VecRestoreArrayRead(b,&barray);
355: barray = lu->rhs_dup;
356: } else {
357: VecGetArrayRead(b,&barray);
358: }
359: VecGetArray(x,&xarray);
361: #if defined(PETSC_USE_COMPLEX)
362: #if defined(PETSC_USE_REAL_SINGLE)
363: ((DNformat*)lu->B.Store)->nzval = (singlecomplex*)barray;
364: ((DNformat*)lu->X.Store)->nzval = (singlecomplex*)xarray;
365: #else
366: ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
367: ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
368: #endif
369: #else
370: ((DNformat*)lu->B.Store)->nzval = (void*)barray;
371: ((DNformat*)lu->X.Store)->nzval = xarray;
372: #endif
374: lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
375: if (A->factortype == MAT_FACTOR_LU) {
376: #if defined(PETSC_USE_COMPLEX)
377: #if defined(PETSC_USE_REAL_SINGLE)
378: PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
379: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
380: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
381: #else
382: PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
383: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
384: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
385: #endif
386: #else
387: #if defined(PETSC_USE_REAL_SINGLE)
388: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
389: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
390: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
391: #else
392: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
393: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
394: &lu->Glu,&lu->mem_usage, &lu->stat, &info));
395: #endif
396: #endif
397: } else if (A->factortype == MAT_FACTOR_ILU) {
398: #if defined(PETSC_USE_COMPLEX)
399: #if defined(PETSC_USE_REAL_SINGLE)
400: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
401: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
402: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
403: #else
404: PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
405: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
406: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
407: #endif
408: #else
409: #if defined(PETSC_USE_REAL_SINGLE)
410: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
411: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
412: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
413: #else
414: PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
415: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
416: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
417: #endif
418: #endif
419: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
420: if (!lu->options.Equil) {
421: VecRestoreArrayRead(b,&barray);
422: }
423: VecRestoreArray(x,&xarray);
425: if (!info || info == lu->A.ncol+1) {
426: if (lu->options.IterRefine) {
427: PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
428: PetscPrintf(PETSC_COMM_SELF," %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
429: for (i = 0; i < 1; ++i) {
430: PetscPrintf(PETSC_COMM_SELF," %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr);
431: }
432: }
433: } else if (info > 0) {
434: if (lu->lwork == -1) {
435: PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %D bytes\n", info - lu->A.ncol);
436: } else {
437: PetscPrintf(PETSC_COMM_SELF," Warning: gssvx() returns info %D\n",info);
438: }
439: } else if (info < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", info,-info);
441: if (lu->options.PrintStat) {
442: PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
443: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
444: }
445: return(0);
446: }
450: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
451: {
452: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
456: if (A->errortype) {
457: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
458: VecSetInf(x);
459: return(0);
460: }
462: lu->options.Trans = TRANS;
463: MatSolve_SuperLU_Private(A,b,x);
464: return(0);
465: }
469: PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x)
470: {
471: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
475: if (A->errortype) {
476: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
477: VecSetInf(x);
478: return(0);
479: }
481: lu->options.Trans = NOTRANS;
482: MatSolve_SuperLU_Private(A,b,x);
483: return(0);
484: }
488: PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X)
489: {
490: Mat_SuperLU *lu = (Mat_SuperLU*)A->spptr;
491: PetscBool flg;
495: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
496: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
497: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
498: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
499: lu->options.Trans = TRANS;
500: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet");
501: return(0);
502: }
504: /*
505: Note the r permutation is ignored
506: */
509: PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
510: {
511: Mat_SuperLU *lu = (Mat_SuperLU*)(F->spptr);
514: lu->flg = DIFFERENT_NONZERO_PATTERN;
515: lu->CleanUpSuperLU = PETSC_TRUE;
516: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
517: return(0);
518: }
522: static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol)
523: {
524: Mat_SuperLU *lu= (Mat_SuperLU*)F->spptr;
527: lu->options.ILU_DropTol = dtol;
528: return(0);
529: }
533: /*@
534: MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
535: Logically Collective on Mat
537: Input Parameters:
538: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface
539: - dtol - drop tolerance
541: Options Database:
542: . -mat_superlu_ilu_droptol <dtol>
544: Level: beginner
546: References:
547: . SuperLU Users' Guide
549: .seealso: MatGetFactor()
550: @*/
551: PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol)
552: {
558: PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol));
559: return(0);
560: }
564: PetscErrorCode MatFactorGetSolverPackage_seqaij_superlu(Mat A,const MatSolverPackage *type)
565: {
567: *type = MATSOLVERSUPERLU;
568: return(0);
569: }
572: /*MC
573: MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
574: via the external package SuperLU.
576: Use ./configure --download-superlu to have PETSc installed with SuperLU
578: Use -pc_type lu -pc_factor_mat_solver_package superlu to us this direct solver
580: Options Database Keys:
581: + -mat_superlu_equil <FALSE> - Equil (None)
582: . -mat_superlu_colperm <COLAMD> - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD
583: . -mat_superlu_iterrefine <NOREFINE> - (choose one of) NOREFINE SINGLE DOUBLE EXTRA
584: . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None)
585: . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None)
586: . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None)
587: . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None)
588: . -mat_superlu_rowperm <NOROWPERM> - (choose one of) NOROWPERM LargeDiag
589: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
590: . -mat_superlu_printstat <FALSE> - PrintStat (None)
591: . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None)
592: . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None)
593: . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None)
594: . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None)
595: . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None)
596: . -mat_superlu_ilu_norm <0> - ILU_Norm (None)
597: - -mat_superlu_ilu_milu <0> - ILU_MILU (None)
599: Notes: Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves
601: Level: beginner
603: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage
604: M*/
608: static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F)
609: {
610: Mat B;
611: Mat_SuperLU *lu;
613: PetscInt indx,m=A->rmap->n,n=A->cmap->n;
614: PetscBool flg,set;
615: PetscReal real_input;
616: const char *colperm[] ={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
617: const char *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
618: const char *rowperm[] ={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */
621: MatCreate(PetscObjectComm((PetscObject)A),&B);
622: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);
623: MatSetType(B,((PetscObject)A)->type_name);
624: MatSeqAIJSetPreallocation(B,0,NULL);
626: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
627: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
628: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
629: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
631: PetscFree(B->solvertype);
632: PetscStrallocpy(MATSOLVERSUPERLU,&B->solvertype);
634: B->ops->destroy = MatDestroy_SuperLU;
635: B->ops->view = MatView_SuperLU;
636: B->ops->getdiagonal = MatGetDiagonal_SuperLU;
637: B->factortype = ftype;
638: B->assembled = PETSC_TRUE; /* required by -ksp_view */
639: B->preallocated = PETSC_TRUE;
641: PetscNewLog(B,&lu);
643: if (ftype == MAT_FACTOR_LU) {
644: set_default_options(&lu->options);
645: /* Comments from SuperLU_4.0/SRC/dgssvx.c:
646: "Whether or not the system will be equilibrated depends on the
647: scaling of the matrix A, but if equilibration is used, A is
648: overwritten by diag(R)*A*diag(C) and B by diag(R)*B
649: (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
650: We set 'options.Equil = NO' as default because additional space is needed for it.
651: */
652: lu->options.Equil = NO;
653: } else if (ftype == MAT_FACTOR_ILU) {
654: /* Set the default input options of ilu: */
655: PetscStackCall("SuperLU:ilu_set_default_options",ilu_set_default_options(&lu->options));
656: }
657: lu->options.PrintStat = NO;
659: /* Initialize the statistics variables. */
660: PetscStackCall("SuperLU:StatInit",StatInit(&lu->stat));
661: lu->lwork = 0; /* allocate space internally by system malloc */
663: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"SuperLU Options","Mat");
664: PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,NULL);
665: PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
666: if (flg) lu->options.ColPerm = (colperm_t)indx;
667: PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
668: if (flg) lu->options.IterRefine = (IterRefine_t)indx;
669: PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,&set);
670: if (set && flg) lu->options.SymmetricMode = YES;
671: PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&real_input,&flg);
672: if (flg) lu->options.DiagPivotThresh = (double) real_input;
673: PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,&set);
674: if (set && flg) lu->options.PivotGrowth = YES;
675: PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,&set);
676: if (set && flg) lu->options.ConditionNumber = YES;
677: PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg);
678: if (flg) lu->options.RowPerm = (rowperm_t)indx;
679: PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,&set);
680: if (set && flg) lu->options.ReplaceTinyPivot = YES;
681: PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,&set);
682: if (set && flg) lu->options.PrintStat = YES;
683: PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,NULL);
684: if (lu->lwork > 0) {
685: /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/
686: PetscMalloc(lu->lwork,&lu->work);
687: } else if (lu->lwork != 0 && lu->lwork != -1) {
688: PetscPrintf(PETSC_COMM_SELF," Warning: lwork %D is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
689: lu->lwork = 0;
690: }
691: /* ilu options */
692: PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&real_input,&flg);
693: if (flg) lu->options.ILU_DropTol = (double) real_input;
694: PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&real_input,&flg);
695: if (flg) lu->options.ILU_FillTol = (double) real_input;
696: PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&real_input,&flg);
697: if (flg) lu->options.ILU_FillFactor = (double) real_input;
698: PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,NULL);
699: PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg);
700: if (flg) lu->options.ILU_Norm = (norm_t)indx;
701: PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg);
702: if (flg) lu->options.ILU_MILU = (milu_t)indx;
703: PetscOptionsEnd();
704: if (lu->options.Equil == YES) {
705: /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */
706: MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup);
707: }
709: /* Allocate spaces (notice sizes are for the transpose) */
710: PetscMalloc1(m,&lu->etree);
711: PetscMalloc1(n,&lu->perm_r);
712: PetscMalloc1(m,&lu->perm_c);
713: PetscMalloc1(n,&lu->R);
714: PetscMalloc1(m,&lu->C);
716: /* create rhs and solution x without allocate space for .Store */
717: #if defined(PETSC_USE_COMPLEX)
718: #if defined(PETSC_USE_REAL_SINGLE)
719: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
720: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
721: #else
722: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
723: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
724: #endif
725: #else
726: #if defined(PETSC_USE_REAL_SINGLE)
727: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
728: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
729: #else
730: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
731: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
732: #endif
733: #endif
735: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_superlu);
736: PetscObjectComposeFunction((PetscObject)B,"MatSuperluSetILUDropTol_C",MatSuperluSetILUDropTol_SuperLU);
737: B->spptr = lu;
739: *F = B;
740: return(0);
741: }
745: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_SuperLU(void)
746: {
750: MatSolverPackageRegister(MATSOLVERSUPERLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_superlu);
751: MatSolverPackageRegister(MATSOLVERSUPERLU,MATSEQAIJ, MAT_FACTOR_ILU,MatGetFactor_seqaij_superlu);
752: return(0);
753: }