Actual source code: superlu.c
petsc-3.8.4 2018-03-24
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>
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 that defines the SuperLU factored 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: /*
59: Utility function
60: */
61: static PetscErrorCode MatFactorInfo_SuperLU(Mat A,PetscViewer viewer)
62: {
63: Mat_SuperLU *lu= (Mat_SuperLU*)A->data;
64: PetscErrorCode ierr;
65: superlu_options_t options;
68: options = lu->options;
70: PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
71: PetscViewerASCIIPrintf(viewer," Equil: %s\n",(options.Equil != NO) ? "YES" : "NO");
72: PetscViewerASCIIPrintf(viewer," ColPerm: %D\n",options.ColPerm);
73: PetscViewerASCIIPrintf(viewer," IterRefine: %D\n",options.IterRefine);
74: PetscViewerASCIIPrintf(viewer," SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES" : "NO");
75: PetscViewerASCIIPrintf(viewer," DiagPivotThresh: %g\n",options.DiagPivotThresh);
76: PetscViewerASCIIPrintf(viewer," PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES" : "NO");
77: PetscViewerASCIIPrintf(viewer," ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES" : "NO");
78: PetscViewerASCIIPrintf(viewer," RowPerm: %D\n",options.RowPerm);
79: PetscViewerASCIIPrintf(viewer," ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES" : "NO");
80: PetscViewerASCIIPrintf(viewer," PrintStat: %s\n",(options.PrintStat != NO) ? "YES" : "NO");
81: PetscViewerASCIIPrintf(viewer," lwork: %D\n",lu->lwork);
82: if (A->factortype == MAT_FACTOR_ILU) {
83: PetscViewerASCIIPrintf(viewer," ILU_DropTol: %g\n",options.ILU_DropTol);
84: PetscViewerASCIIPrintf(viewer," ILU_FillTol: %g\n",options.ILU_FillTol);
85: PetscViewerASCIIPrintf(viewer," ILU_FillFactor: %g\n",options.ILU_FillFactor);
86: PetscViewerASCIIPrintf(viewer," ILU_DropRule: %D\n",options.ILU_DropRule);
87: PetscViewerASCIIPrintf(viewer," ILU_Norm: %D\n",options.ILU_Norm);
88: PetscViewerASCIIPrintf(viewer," ILU_MILU: %D\n",options.ILU_MILU);
89: }
90: return(0);
91: }
93: PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x)
94: {
95: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
96: const PetscScalar *barray;
97: PetscScalar *xarray;
98: PetscErrorCode ierr;
99: PetscInt info,i,n;
100: PetscReal ferr,berr;
101: static PetscBool cite = PETSC_FALSE;
104: if (lu->lwork == -1) return(0);
105: 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);
107: VecGetLocalSize(x,&n);
108: lu->B.ncol = 1; /* Set the number of right-hand side */
109: if (lu->options.Equil && !lu->rhs_dup) {
110: /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */
111: PetscMalloc1(n,&lu->rhs_dup);
112: }
113: if (lu->options.Equil) {
114: /* Copy b into rsh_dup */
115: VecGetArrayRead(b,&barray);
116: PetscMemcpy(lu->rhs_dup,barray,n*sizeof(PetscScalar));
117: VecRestoreArrayRead(b,&barray);
118: barray = lu->rhs_dup;
119: } else {
120: VecGetArrayRead(b,&barray);
121: }
122: VecGetArray(x,&xarray);
124: #if defined(PETSC_USE_COMPLEX)
125: #if defined(PETSC_USE_REAL_SINGLE)
126: ((DNformat*)lu->B.Store)->nzval = (singlecomplex*)barray;
127: ((DNformat*)lu->X.Store)->nzval = (singlecomplex*)xarray;
128: #else
129: ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
130: ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
131: #endif
132: #else
133: ((DNformat*)lu->B.Store)->nzval = (void*)barray;
134: ((DNformat*)lu->X.Store)->nzval = xarray;
135: #endif
137: lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
138: if (A->factortype == MAT_FACTOR_LU) {
139: #if defined(PETSC_USE_COMPLEX)
140: #if defined(PETSC_USE_REAL_SINGLE)
141: PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
142: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
143: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
144: #else
145: PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
146: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
147: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
148: #endif
149: #else
150: #if defined(PETSC_USE_REAL_SINGLE)
151: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
152: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
153: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
154: #else
155: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
156: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
157: &lu->Glu,&lu->mem_usage, &lu->stat, &info));
158: #endif
159: #endif
160: } else if (A->factortype == MAT_FACTOR_ILU) {
161: #if defined(PETSC_USE_COMPLEX)
162: #if defined(PETSC_USE_REAL_SINGLE)
163: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
164: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
165: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
166: #else
167: PetscStackCall("SuperLU:zgsisx",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->Glu, &lu->mem_usage, &lu->stat, &info));
170: #endif
171: #else
172: #if defined(PETSC_USE_REAL_SINGLE)
173: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
174: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
175: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
176: #else
177: PetscStackCall("SuperLU:dgsisx",dgsisx(&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,
179: &lu->Glu, &lu->mem_usage, &lu->stat, &info));
180: #endif
181: #endif
182: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
183: if (!lu->options.Equil) {
184: VecRestoreArrayRead(b,&barray);
185: }
186: VecRestoreArray(x,&xarray);
188: if (!info || info == lu->A.ncol+1) {
189: if (lu->options.IterRefine) {
190: PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
191: PetscPrintf(PETSC_COMM_SELF," %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
192: for (i = 0; i < 1; ++i) {
193: PetscPrintf(PETSC_COMM_SELF," %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr);
194: }
195: }
196: } else if (info > 0) {
197: if (lu->lwork == -1) {
198: PetscPrintf(PETSC_COMM_SELF," ** Estimated memory: %D bytes\n", info - lu->A.ncol);
199: } else {
200: PetscPrintf(PETSC_COMM_SELF," Warning: gssvx() returns info %D\n",info);
201: }
202: } 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);
204: if (lu->options.PrintStat) {
205: PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
206: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
207: }
208: return(0);
209: }
211: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
212: {
213: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
217: if (A->factorerrortype) {
218: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
219: VecSetInf(x);
220: return(0);
221: }
223: lu->options.Trans = TRANS;
224: MatSolve_SuperLU_Private(A,b,x);
225: return(0);
226: }
228: PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x)
229: {
230: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
234: if (A->factorerrortype) {
235: PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
236: VecSetInf(x);
237: return(0);
238: }
240: lu->options.Trans = NOTRANS;
241: MatSolve_SuperLU_Private(A,b,x);
242: return(0);
243: }
245: static PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info)
246: {
247: Mat_SuperLU *lu = (Mat_SuperLU*)F->data;
248: Mat_SeqAIJ *aa;
250: PetscInt sinfo;
251: PetscReal ferr, berr;
252: NCformat *Ustore;
253: SCformat *Lstore;
256: if (lu->flg == SAME_NONZERO_PATTERN) { /* successing numerical factorization */
257: lu->options.Fact = SamePattern;
258: /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
259: Destroy_SuperMatrix_Store(&lu->A);
260: if (lu->options.Equil) {
261: MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN);
262: }
263: if (lu->lwork >= 0) {
264: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
265: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
266: lu->options.Fact = SamePattern;
267: }
268: }
270: /* Create the SuperMatrix for lu->A=A^T:
271: Since SuperLU likes column-oriented matrices,we pass it the transpose,
272: and then solve A^T X = B in MatSolve(). */
273: if (lu->options.Equil) {
274: aa = (Mat_SeqAIJ*)(lu->A_dup)->data;
275: } else {
276: aa = (Mat_SeqAIJ*)(A)->data;
277: }
278: #if defined(PETSC_USE_COMPLEX)
279: #if defined(PETSC_USE_REAL_SINGLE)
280: 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));
281: #else
282: 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));
283: #endif
284: #else
285: #if defined(PETSC_USE_REAL_SINGLE)
286: 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));
287: #else
288: 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));
289: #endif
290: #endif
292: /* Numerical factorization */
293: lu->B.ncol = 0; /* Indicate not to solve the system */
294: if (F->factortype == MAT_FACTOR_LU) {
295: #if defined(PETSC_USE_COMPLEX)
296: #if defined(PETSC_USE_REAL_SINGLE)
297: PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
298: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
299: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
300: #else
301: PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
302: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
303: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
304: #endif
305: #else
306: #if defined(PETSC_USE_REAL_SINGLE)
307: PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
308: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
309: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
310: #else
311: PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
312: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
313: &lu->Glu,&lu->mem_usage, &lu->stat, &sinfo));
314: #endif
315: #endif
316: } else if (F->factortype == MAT_FACTOR_ILU) {
317: /* Compute the incomplete factorization, condition number and pivot growth */
318: #if defined(PETSC_USE_COMPLEX)
319: #if defined(PETSC_USE_REAL_SINGLE)
320: PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
321: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
322: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
323: #else
324: PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
325: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
326: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
327: #endif
328: #else
329: #if defined(PETSC_USE_REAL_SINGLE)
330: PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
331: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
332: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
333: #else
334: PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
335: &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
336: &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
337: #endif
338: #endif
339: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
340: if (!sinfo || sinfo == lu->A.ncol+1) {
341: if (lu->options.PivotGrowth) {
342: PetscPrintf(PETSC_COMM_SELF," Recip. pivot growth = %e\n", lu->rpg);
343: }
344: if (lu->options.ConditionNumber) {
345: PetscPrintf(PETSC_COMM_SELF," Recip. condition number = %e\n", lu->rcond);
346: }
347: } else if (sinfo > 0) {
348: if (A->erroriffailure) {
349: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %D",sinfo);
350: } else {
351: if (sinfo <= lu->A.ncol) {
352: if (lu->options.ILU_FillTol == 0.0) {
353: F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
354: }
355: PetscInfo2(F,"Number of zero pivots %D, ILU_FillTol %g\n",sinfo,lu->options.ILU_FillTol);
356: } else if (sinfo == lu->A.ncol + 1) {
357: /*
358: U is nonsingular, but RCOND is less than machine
359: precision, meaning that the matrix is singular to
360: working precision. Nevertheless, the solution and
361: error bounds are computed because there are a number
362: of situations where the computed solution can be more
363: accurate than the value of RCOND would suggest.
364: */
365: PetscInfo1(F,"Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %D",sinfo);
366: } else { /* sinfo > lu->A.ncol + 1 */
367: F->factorerrortype = MAT_FACTOR_OUTMEMORY;
368: PetscInfo1(F,"Number of bytes allocated when memory allocation fails %D\n",sinfo);
369: }
370: }
371: } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", sinfo,-sinfo);
373: if (lu->options.PrintStat) {
374: PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
375: PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
376: Lstore = (SCformat*) lu->L.Store;
377: Ustore = (NCformat*) lu->U.Store;
378: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor L = %D\n", Lstore->nnz);
379: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in factor U = %D\n", Ustore->nnz);
380: PetscPrintf(PETSC_COMM_SELF," No of nonzeros in L+U = %D\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
381: PetscPrintf(PETSC_COMM_SELF," L\\U MB %.3f\ttotal MB needed %.3f\n",
382: lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6);
383: }
385: lu->flg = SAME_NONZERO_PATTERN;
386: F->ops->solve = MatSolve_SuperLU;
387: F->ops->solvetranspose = MatSolveTranspose_SuperLU;
388: F->ops->matsolve = NULL;
389: return(0);
390: }
392: static PetscErrorCode MatDestroy_SuperLU(Mat A)
393: {
395: Mat_SuperLU *lu=(Mat_SuperLU*)A->data;
398: if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
399: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->A));
400: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->B));
401: PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->X));
402: PetscStackCall("SuperLU:StatFree",StatFree(&lu->stat));
403: if (lu->lwork >= 0) {
404: PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
405: PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
406: }
407: }
408: PetscFree(lu->etree);
409: PetscFree(lu->perm_r);
410: PetscFree(lu->perm_c);
411: PetscFree(lu->R);
412: PetscFree(lu->C);
413: PetscFree(lu->rhs_dup);
414: MatDestroy(&lu->A_dup);
415: PetscFree(A->data);
417: /* clear composed functions */
418: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
419: PetscObjectComposeFunction((PetscObject)A,"MatSuperluSetILUDropTol_C",NULL);
420: return(0);
421: }
423: static PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
424: {
425: PetscErrorCode ierr;
426: PetscBool iascii;
427: PetscViewerFormat format;
430: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
431: if (iascii) {
432: PetscViewerGetFormat(viewer,&format);
433: if (format == PETSC_VIEWER_ASCII_INFO) {
434: MatFactorInfo_SuperLU(A,viewer);
435: }
436: }
437: return(0);
438: }
440: PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X)
441: {
442: Mat_SuperLU *lu = (Mat_SuperLU*)A->data;
443: PetscBool flg;
447: PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
448: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
449: PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
450: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
451: lu->options.Trans = TRANS;
452: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet");
453: return(0);
454: }
456: /*
457: Note the r permutation is ignored
458: */
459: static PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
460: {
461: Mat_SuperLU *lu = (Mat_SuperLU*)(F->data);
464: lu->flg = DIFFERENT_NONZERO_PATTERN;
465: lu->CleanUpSuperLU = PETSC_TRUE;
466: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
467: return(0);
468: }
470: static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol)
471: {
472: Mat_SuperLU *lu= (Mat_SuperLU*)F->data;
475: lu->options.ILU_DropTol = dtol;
476: return(0);
477: }
479: /*@
480: MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
481: Logically Collective on Mat
483: Input Parameters:
484: + F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface
485: - dtol - drop tolerance
487: Options Database:
488: . -mat_superlu_ilu_droptol <dtol>
490: Level: beginner
492: References:
493: . SuperLU Users' Guide
495: .seealso: MatGetFactor()
496: @*/
497: PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol)
498: {
504: PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol));
505: return(0);
506: }
508: PetscErrorCode MatFactorGetSolverPackage_seqaij_superlu(Mat A,const MatSolverPackage *type)
509: {
511: *type = MATSOLVERSUPERLU;
512: return(0);
513: }
515: /*MC
516: MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
517: via the external package SuperLU.
519: Use ./configure --download-superlu to have PETSc installed with SuperLU
521: Use -pc_type lu -pc_factor_mat_solver_package superlu to use this direct solver
523: Options Database Keys:
524: + -mat_superlu_equil <FALSE> - Equil (None)
525: . -mat_superlu_colperm <COLAMD> - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD
526: . -mat_superlu_iterrefine <NOREFINE> - (choose one of) NOREFINE SINGLE DOUBLE EXTRA
527: . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None)
528: . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None)
529: . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None)
530: . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None)
531: . -mat_superlu_rowperm <NOROWPERM> - (choose one of) NOROWPERM LargeDiag
532: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
533: . -mat_superlu_printstat <FALSE> - PrintStat (None)
534: . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None)
535: . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None)
536: . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None)
537: . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None)
538: . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None)
539: . -mat_superlu_ilu_norm <0> - ILU_Norm (None)
540: - -mat_superlu_ilu_milu <0> - ILU_MILU (None)
542: Notes: Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves
544: Level: beginner
546: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverPackage(), MatSolverPackage
547: M*/
549: static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F)
550: {
551: Mat B;
552: Mat_SuperLU *lu;
554: PetscInt indx,m=A->rmap->n,n=A->cmap->n;
555: PetscBool flg,set;
556: PetscReal real_input;
557: const char *colperm[] ={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
558: const char *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
559: const char *rowperm[] ={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */
562: MatCreate(PetscObjectComm((PetscObject)A),&B);
563: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);
564: PetscStrallocpy("superlu",&((PetscObject)B)->type_name);
565: MatSetUp(B);
566: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
567: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
568: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
569: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
571: PetscFree(B->solvertype);
572: PetscStrallocpy(MATSOLVERSUPERLU,&B->solvertype);
574: B->ops->getinfo = MatGetInfo_External;
575: B->ops->destroy = MatDestroy_SuperLU;
576: B->ops->view = MatView_SuperLU;
577: B->factortype = ftype;
578: B->assembled = PETSC_TRUE; /* required by -ksp_view */
579: B->preallocated = PETSC_TRUE;
581: PetscNewLog(B,&lu);
583: if (ftype == MAT_FACTOR_LU) {
584: set_default_options(&lu->options);
585: /* Comments from SuperLU_4.0/SRC/dgssvx.c:
586: "Whether or not the system will be equilibrated depends on the
587: scaling of the matrix A, but if equilibration is used, A is
588: overwritten by diag(R)*A*diag(C) and B by diag(R)*B
589: (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
590: We set 'options.Equil = NO' as default because additional space is needed for it.
591: */
592: lu->options.Equil = NO;
593: } else if (ftype == MAT_FACTOR_ILU) {
594: /* Set the default input options of ilu: */
595: PetscStackCall("SuperLU:ilu_set_default_options",ilu_set_default_options(&lu->options));
596: }
597: lu->options.PrintStat = NO;
599: /* Initialize the statistics variables. */
600: PetscStackCall("SuperLU:StatInit",StatInit(&lu->stat));
601: lu->lwork = 0; /* allocate space internally by system malloc */
603: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"SuperLU Options","Mat");
604: PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,NULL);
605: PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
606: if (flg) lu->options.ColPerm = (colperm_t)indx;
607: PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
608: if (flg) lu->options.IterRefine = (IterRefine_t)indx;
609: PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,&set);
610: if (set && flg) lu->options.SymmetricMode = YES;
611: PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&real_input,&flg);
612: if (flg) lu->options.DiagPivotThresh = (double) real_input;
613: PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,&set);
614: if (set && flg) lu->options.PivotGrowth = YES;
615: PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,&set);
616: if (set && flg) lu->options.ConditionNumber = YES;
617: PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg);
618: if (flg) lu->options.RowPerm = (rowperm_t)indx;
619: PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,&set);
620: if (set && flg) lu->options.ReplaceTinyPivot = YES;
621: PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,&set);
622: if (set && flg) lu->options.PrintStat = YES;
623: PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,NULL);
624: if (lu->lwork > 0) {
625: /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/
626: PetscMalloc(lu->lwork,&lu->work);
627: } else if (lu->lwork != 0 && lu->lwork != -1) {
628: PetscPrintf(PETSC_COMM_SELF," Warning: lwork %D is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
629: lu->lwork = 0;
630: }
631: /* ilu options */
632: PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&real_input,&flg);
633: if (flg) lu->options.ILU_DropTol = (double) real_input;
634: PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&real_input,&flg);
635: if (flg) lu->options.ILU_FillTol = (double) real_input;
636: PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&real_input,&flg);
637: if (flg) lu->options.ILU_FillFactor = (double) real_input;
638: PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,NULL);
639: PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg);
640: if (flg) lu->options.ILU_Norm = (norm_t)indx;
641: PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg);
642: if (flg) lu->options.ILU_MILU = (milu_t)indx;
643: PetscOptionsEnd();
644: if (lu->options.Equil == YES) {
645: /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */
646: MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup);
647: }
649: /* Allocate spaces (notice sizes are for the transpose) */
650: PetscMalloc1(m,&lu->etree);
651: PetscMalloc1(n,&lu->perm_r);
652: PetscMalloc1(m,&lu->perm_c);
653: PetscMalloc1(n,&lu->R);
654: PetscMalloc1(m,&lu->C);
656: /* create rhs and solution x without allocate space for .Store */
657: #if defined(PETSC_USE_COMPLEX)
658: #if defined(PETSC_USE_REAL_SINGLE)
659: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
660: PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
661: #else
662: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
663: PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
664: #endif
665: #else
666: #if defined(PETSC_USE_REAL_SINGLE)
667: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
668: PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
669: #else
670: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
671: PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
672: #endif
673: #endif
675: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqaij_superlu);
676: PetscObjectComposeFunction((PetscObject)B,"MatSuperluSetILUDropTol_C",MatSuperluSetILUDropTol_SuperLU);
677: B->data = lu;
679: *F = B;
680: return(0);
681: }
683: PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_SuperLU(void)
684: {
688: MatSolverPackageRegister(MATSOLVERSUPERLU,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_superlu);
689: MatSolverPackageRegister(MATSOLVERSUPERLU,MATSEQAIJ, MAT_FACTOR_ILU,MatGetFactor_seqaij_superlu);
690: return(0);
691: }