Actual source code: superlu.c
2: /*
3: This file implements a subclass of the SeqAIJ matrix class that uses
4: the SuperLU sparse solver.
5: */
7: /*
8: Defines the data structure for the base matrix type (SeqAIJ)
9: */
10: #include <../src/mat/impls/aij/seq/aij.h>
12: /*
13: SuperLU include files
14: */
15: EXTERN_C_BEGIN
16: #if defined(PETSC_USE_COMPLEX)
17: #if defined(PETSC_USE_REAL_SINGLE)
18: #include <slu_cdefs.h>
19: #else
20: #include <slu_zdefs.h>
21: #endif
22: #else
23: #if defined(PETSC_USE_REAL_SINGLE)
24: #include <slu_sdefs.h>
25: #else
26: #include <slu_ddefs.h>
27: #endif
28: #endif
29: #include <slu_util.h>
30: EXTERN_C_END
32: /*
33: This is the data that defines the SuperLU factored matrix type
34: */
35: typedef struct {
36: SuperMatrix A, L, U, B, X;
37: superlu_options_t options;
38: PetscInt *perm_c; /* column permutation vector */
39: PetscInt *perm_r; /* row permutations from partial pivoting */
40: PetscInt *etree;
41: PetscReal *R, *C;
42: char equed[1];
43: PetscInt lwork;
44: void *work;
45: PetscReal rpg, rcond;
46: mem_usage_t mem_usage;
47: MatStructure flg;
48: SuperLUStat_t stat;
49: Mat A_dup;
50: PetscScalar *rhs_dup;
51: GlobalLU_t Glu;
52: PetscBool needconversion;
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 MatView_Info_SuperLU(Mat A, PetscViewer viewer)
62: {
63: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
64: superlu_options_t options;
66: PetscFunctionBegin;
67: options = lu->options;
69: PetscCall(PetscViewerASCIIPrintf(viewer, "SuperLU run parameters:\n"));
70: PetscCall(PetscViewerASCIIPrintf(viewer, " Equil: %s\n", (options.Equil != NO) ? "YES" : "NO"));
71: PetscCall(PetscViewerASCIIPrintf(viewer, " ColPerm: %" PetscInt_FMT "\n", options.ColPerm));
72: PetscCall(PetscViewerASCIIPrintf(viewer, " IterRefine: %" PetscInt_FMT "\n", options.IterRefine));
73: PetscCall(PetscViewerASCIIPrintf(viewer, " SymmetricMode: %s\n", (options.SymmetricMode != NO) ? "YES" : "NO"));
74: PetscCall(PetscViewerASCIIPrintf(viewer, " DiagPivotThresh: %g\n", options.DiagPivotThresh));
75: PetscCall(PetscViewerASCIIPrintf(viewer, " PivotGrowth: %s\n", (options.PivotGrowth != NO) ? "YES" : "NO"));
76: PetscCall(PetscViewerASCIIPrintf(viewer, " ConditionNumber: %s\n", (options.ConditionNumber != NO) ? "YES" : "NO"));
77: PetscCall(PetscViewerASCIIPrintf(viewer, " RowPerm: %" PetscInt_FMT "\n", options.RowPerm));
78: PetscCall(PetscViewerASCIIPrintf(viewer, " ReplaceTinyPivot: %s\n", (options.ReplaceTinyPivot != NO) ? "YES" : "NO"));
79: PetscCall(PetscViewerASCIIPrintf(viewer, " PrintStat: %s\n", (options.PrintStat != NO) ? "YES" : "NO"));
80: PetscCall(PetscViewerASCIIPrintf(viewer, " lwork: %" PetscInt_FMT "\n", lu->lwork));
81: if (A->factortype == MAT_FACTOR_ILU) {
82: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_DropTol: %g\n", options.ILU_DropTol));
83: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_FillTol: %g\n", options.ILU_FillTol));
84: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_FillFactor: %g\n", options.ILU_FillFactor));
85: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_DropRule: %" PetscInt_FMT "\n", options.ILU_DropRule));
86: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_Norm: %" PetscInt_FMT "\n", options.ILU_Norm));
87: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU_MILU: %" PetscInt_FMT "\n", options.ILU_MILU));
88: }
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: PetscErrorCode MatSolve_SuperLU_Private(Mat A, Vec b, Vec x)
93: {
94: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
95: const PetscScalar *barray;
96: PetscScalar *xarray;
97: PetscInt info, i, n;
98: PetscReal ferr, berr;
99: static PetscBool cite = PETSC_FALSE;
101: PetscFunctionBegin;
102: if (lu->lwork == -1) PetscFunctionReturn(PETSC_SUCCESS);
103: PetscCall(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 "
104: "pivoting},\n journal = {SIAM J. Matrix Analysis and Applications},\n year = {1999},\n volume = {20},\n number = {3},\n pages = {720-755}\n}\n",
105: &cite));
107: PetscCall(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: PetscCall(PetscMalloc1(n, &lu->rhs_dup));
112: }
113: if (lu->options.Equil) {
114: /* Copy b into rsh_dup */
115: PetscCall(VecGetArrayRead(b, &barray));
116: PetscCall(PetscArraycpy(lu->rhs_dup, barray, n));
117: PetscCall(VecRestoreArrayRead(b, &barray));
118: barray = lu->rhs_dup;
119: } else {
120: PetscCall(VecGetArrayRead(b, &barray));
121: }
122: PetscCall(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: PetscStackCallExternalVoid("SuperLU:cgssvx", cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
142: #else
143: PetscStackCallExternalVoid("SuperLU:zgssvx", zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
144: #endif
145: #else
146: #if defined(PETSC_USE_REAL_SINGLE)
147: PetscStackCallExternalVoid("SuperLU:sgssvx", sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
148: #else
149: PetscStackCallExternalVoid("SuperLU:dgssvx", dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
150: #endif
151: #endif
152: } else if (A->factortype == MAT_FACTOR_ILU) {
153: #if defined(PETSC_USE_COMPLEX)
154: #if defined(PETSC_USE_REAL_SINGLE)
155: PetscStackCallExternalVoid("SuperLU:cgsisx", cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
156: #else
157: PetscStackCallExternalVoid("SuperLU:zgsisx", zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
158: #endif
159: #else
160: #if defined(PETSC_USE_REAL_SINGLE)
161: PetscStackCallExternalVoid("SuperLU:sgsisx", sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
162: #else
163: PetscStackCallExternalVoid("SuperLU:dgsisx", dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &info));
164: #endif
165: #endif
166: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
167: if (!lu->options.Equil) PetscCall(VecRestoreArrayRead(b, &barray));
168: PetscCall(VecRestoreArray(x, &xarray));
170: if (!info || info == lu->A.ncol + 1) {
171: if (lu->options.IterRefine) {
172: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Iterative Refinement:\n"));
173: PetscCall(PetscPrintf(PETSC_COMM_SELF, " %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR"));
174: for (i = 0; i < 1; ++i) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %8d%8d%16e%16e\n", i + 1, lu->stat.RefineSteps, ferr, berr));
175: }
176: } else if (info > 0) {
177: if (lu->lwork == -1) {
178: PetscCall(PetscPrintf(PETSC_COMM_SELF, " ** Estimated memory: %" PetscInt_FMT " bytes\n", info - lu->A.ncol));
179: } else {
180: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Warning: gssvx() returns info %" PetscInt_FMT "\n", info));
181: }
182: } else PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "info = %" PetscInt_FMT ", the %" PetscInt_FMT "-th argument in gssvx() had an illegal value", info, -info);
184: if (lu->options.PrintStat) {
185: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve__SuperLU():\n"));
186: PetscStackCallExternalVoid("SuperLU:StatPrint", StatPrint(&lu->stat));
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode MatSolve_SuperLU(Mat A, Vec b, Vec x)
192: {
193: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
194: trans_t oldOption;
196: PetscFunctionBegin;
197: if (A->factorerrortype) {
198: PetscCall(PetscInfo(A, "MatSolve is called with singular matrix factor, skip\n"));
199: PetscCall(VecSetInf(x));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: oldOption = lu->options.Trans;
204: lu->options.Trans = TRANS;
205: PetscCall(MatSolve_SuperLU_Private(A, b, x));
206: lu->options.Trans = oldOption;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: PetscErrorCode MatSolveTranspose_SuperLU(Mat A, Vec b, Vec x)
211: {
212: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
213: trans_t oldOption;
215: PetscFunctionBegin;
216: if (A->factorerrortype) {
217: PetscCall(PetscInfo(A, "MatSolve is called with singular matrix factor, skip\n"));
218: PetscCall(VecSetInf(x));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: oldOption = lu->options.Trans;
223: lu->options.Trans = NOTRANS;
224: PetscCall(MatSolve_SuperLU_Private(A, b, x));
225: lu->options.Trans = oldOption;
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F, Mat A, const MatFactorInfo *info)
230: {
231: Mat_SuperLU *lu = (Mat_SuperLU *)F->data;
232: Mat_SeqAIJ *aa;
233: PetscInt sinfo;
234: PetscReal ferr, berr;
235: NCformat *Ustore;
236: SCformat *Lstore;
238: PetscFunctionBegin;
239: if (lu->flg == SAME_NONZERO_PATTERN) { /* successive numerical factorization */
240: lu->options.Fact = SamePattern;
241: /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
242: Destroy_SuperMatrix_Store(&lu->A);
243: if (lu->A_dup) PetscCall(MatCopy_SeqAIJ(A, lu->A_dup, SAME_NONZERO_PATTERN));
245: if (lu->lwork >= 0) {
246: PetscStackCallExternalVoid("SuperLU:Destroy_SuperNode_Matrix", Destroy_SuperNode_Matrix(&lu->L));
247: PetscStackCallExternalVoid("SuperLU:Destroy_CompCol_Matrix", Destroy_CompCol_Matrix(&lu->U));
248: lu->options.Fact = SamePattern;
249: }
250: }
252: /* Create the SuperMatrix for lu->A=A^T:
253: Since SuperLU likes column-oriented matrices,we pass it the transpose,
254: and then solve A^T X = B in MatSolve(). */
255: if (lu->A_dup) {
256: aa = (Mat_SeqAIJ *)(lu->A_dup)->data;
257: } else {
258: aa = (Mat_SeqAIJ *)(A)->data;
259: }
260: #if defined(PETSC_USE_COMPLEX)
261: #if defined(PETSC_USE_REAL_SINGLE)
262: PetscStackCallExternalVoid("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));
263: #else
264: PetscStackCallExternalVoid("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));
265: #endif
266: #else
267: #if defined(PETSC_USE_REAL_SINGLE)
268: PetscStackCallExternalVoid("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));
269: #else
270: PetscStackCallExternalVoid("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));
271: #endif
272: #endif
274: /* Numerical factorization */
275: lu->B.ncol = 0; /* Indicate not to solve the system */
276: if (F->factortype == MAT_FACTOR_LU) {
277: #if defined(PETSC_USE_COMPLEX)
278: #if defined(PETSC_USE_REAL_SINGLE)
279: PetscStackCallExternalVoid("SuperLU:cgssvx", cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
280: #else
281: PetscStackCallExternalVoid("SuperLU:zgssvx", zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
282: #endif
283: #else
284: #if defined(PETSC_USE_REAL_SINGLE)
285: PetscStackCallExternalVoid("SuperLU:sgssvx", sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
286: #else
287: PetscStackCallExternalVoid("SuperLU:dgssvx", dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
288: #endif
289: #endif
290: } else if (F->factortype == MAT_FACTOR_ILU) {
291: /* Compute the incomplete factorization, condition number and pivot growth */
292: #if defined(PETSC_USE_COMPLEX)
293: #if defined(PETSC_USE_REAL_SINGLE)
294: PetscStackCallExternalVoid("SuperLU:cgsisx", cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
295: #else
296: PetscStackCallExternalVoid("SuperLU:zgsisx", zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
297: #endif
298: #else
299: #if defined(PETSC_USE_REAL_SINGLE)
300: PetscStackCallExternalVoid("SuperLU:sgsisx", sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
301: #else
302: PetscStackCallExternalVoid("SuperLU:dgsisx", dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C, &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
303: #endif
304: #endif
305: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
306: if (!sinfo || sinfo == lu->A.ncol + 1) {
307: if (lu->options.PivotGrowth) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Recip. pivot growth = %e\n", lu->rpg));
308: if (lu->options.ConditionNumber) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Recip. condition number = %e\n", lu->rcond));
309: } else if (sinfo > 0) {
310: if (A->erroriffailure) {
311: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot in row %" PetscInt_FMT, sinfo);
312: } else {
313: if (sinfo <= lu->A.ncol) {
314: if (lu->options.ILU_FillTol == 0.0) F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
315: PetscCall(PetscInfo(F, "Number of zero pivots %" PetscInt_FMT ", ILU_FillTol %g\n", sinfo, lu->options.ILU_FillTol));
316: } else if (sinfo == lu->A.ncol + 1) {
317: /*
318: U is nonsingular, but RCOND is less than machine
319: precision, meaning that the matrix is singular to
320: working precision. Nevertheless, the solution and
321: error bounds are computed because there are a number
322: of situations where the computed solution can be more
323: accurate than the value of RCOND would suggest.
324: */
325: PetscCall(PetscInfo(F, "Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %" PetscInt_FMT, sinfo));
326: } else { /* sinfo > lu->A.ncol + 1 */
327: F->factorerrortype = MAT_FACTOR_OUTMEMORY;
328: PetscCall(PetscInfo(F, "Number of bytes allocated when memory allocation fails %" PetscInt_FMT "\n", sinfo));
329: }
330: }
331: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "info = %" PetscInt_FMT ", the %" PetscInt_FMT "-th argument in gssvx() had an illegal value", sinfo, -sinfo);
333: if (lu->options.PrintStat) {
334: PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatLUFactorNumeric_SuperLU():\n"));
335: PetscStackCallExternalVoid("SuperLU:StatPrint", StatPrint(&lu->stat));
336: Lstore = (SCformat *)lu->L.Store;
337: Ustore = (NCformat *)lu->U.Store;
338: PetscCall(PetscPrintf(PETSC_COMM_SELF, " No of nonzeros in factor L = %" PetscInt_FMT "\n", Lstore->nnz));
339: PetscCall(PetscPrintf(PETSC_COMM_SELF, " No of nonzeros in factor U = %" PetscInt_FMT "\n", Ustore->nnz));
340: PetscCall(PetscPrintf(PETSC_COMM_SELF, " No of nonzeros in L+U = %" PetscInt_FMT "\n", Lstore->nnz + Ustore->nnz - lu->A.ncol));
341: PetscCall(PetscPrintf(PETSC_COMM_SELF, " L\\U MB %.3f\ttotal MB needed %.3f\n", lu->mem_usage.for_lu / 1e6, lu->mem_usage.total_needed / 1e6));
342: }
344: lu->flg = SAME_NONZERO_PATTERN;
345: F->ops->solve = MatSolve_SuperLU;
346: F->ops->solvetranspose = MatSolveTranspose_SuperLU;
347: F->ops->matsolve = NULL;
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode MatDestroy_SuperLU(Mat A)
352: {
353: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
355: PetscFunctionBegin;
356: if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
357: PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store", Destroy_SuperMatrix_Store(&lu->A));
358: PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store", Destroy_SuperMatrix_Store(&lu->B));
359: PetscStackCallExternalVoid("SuperLU:Destroy_SuperMatrix_Store", Destroy_SuperMatrix_Store(&lu->X));
360: PetscStackCallExternalVoid("SuperLU:StatFree", StatFree(&lu->stat));
361: if (lu->lwork >= 0) {
362: PetscStackCallExternalVoid("SuperLU:Destroy_SuperNode_Matrix", Destroy_SuperNode_Matrix(&lu->L));
363: PetscStackCallExternalVoid("SuperLU:Destroy_CompCol_Matrix", Destroy_CompCol_Matrix(&lu->U));
364: }
365: }
366: PetscCall(PetscFree(lu->etree));
367: PetscCall(PetscFree(lu->perm_r));
368: PetscCall(PetscFree(lu->perm_c));
369: PetscCall(PetscFree(lu->R));
370: PetscCall(PetscFree(lu->C));
371: PetscCall(PetscFree(lu->rhs_dup));
372: PetscCall(MatDestroy(&lu->A_dup));
373: PetscCall(PetscFree(A->data));
375: /* clear composed functions */
376: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
377: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSuperluSetILUDropTol_C", NULL));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode MatView_SuperLU(Mat A, PetscViewer viewer)
382: {
383: PetscBool iascii;
384: PetscViewerFormat format;
386: PetscFunctionBegin;
387: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
388: if (iascii) {
389: PetscCall(PetscViewerGetFormat(viewer, &format));
390: if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_SuperLU(A, viewer));
391: }
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: PetscErrorCode MatMatSolve_SuperLU(Mat A, Mat B, Mat X)
396: {
397: Mat_SuperLU *lu = (Mat_SuperLU *)A->data;
398: PetscBool flg;
400: PetscFunctionBegin;
401: PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
402: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
403: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
404: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
405: lu->options.Trans = TRANS;
406: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatSolve_SuperLU() is not implemented yet");
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: static PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
411: {
412: Mat_SuperLU *lu = (Mat_SuperLU *)(F->data);
413: PetscInt indx;
414: PetscBool flg, set;
415: PetscReal real_input;
416: const char *colperm[] = {"NATURAL", "MMD_ATA", "MMD_AT_PLUS_A", "COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
417: const char *iterrefine[] = {"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
418: const char *rowperm[] = {"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */
420: PetscFunctionBegin;
421: /* Set options to F */
422: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "SuperLU Options", "Mat");
423: PetscCall(PetscOptionsBool("-mat_superlu_equil", "Equil", "None", (PetscBool)lu->options.Equil, (PetscBool *)&lu->options.Equil, NULL));
424: PetscCall(PetscOptionsEList("-mat_superlu_colperm", "ColPerm", "None", colperm, 4, colperm[3], &indx, &flg));
425: if (flg) lu->options.ColPerm = (colperm_t)indx;
426: PetscCall(PetscOptionsEList("-mat_superlu_iterrefine", "IterRefine", "None", iterrefine, 4, iterrefine[0], &indx, &flg));
427: if (flg) lu->options.IterRefine = (IterRefine_t)indx;
428: PetscCall(PetscOptionsBool("-mat_superlu_symmetricmode", "SymmetricMode", "None", (PetscBool)lu->options.SymmetricMode, &flg, &set));
429: if (set && flg) lu->options.SymmetricMode = YES;
430: PetscCall(PetscOptionsReal("-mat_superlu_diagpivotthresh", "DiagPivotThresh", "None", lu->options.DiagPivotThresh, &real_input, &flg));
431: if (flg) lu->options.DiagPivotThresh = (double)real_input;
432: PetscCall(PetscOptionsBool("-mat_superlu_pivotgrowth", "PivotGrowth", "None", (PetscBool)lu->options.PivotGrowth, &flg, &set));
433: if (set && flg) lu->options.PivotGrowth = YES;
434: PetscCall(PetscOptionsBool("-mat_superlu_conditionnumber", "ConditionNumber", "None", (PetscBool)lu->options.ConditionNumber, &flg, &set));
435: if (set && flg) lu->options.ConditionNumber = YES;
436: PetscCall(PetscOptionsEList("-mat_superlu_rowperm", "rowperm", "None", rowperm, 2, rowperm[lu->options.RowPerm], &indx, &flg));
437: if (flg) lu->options.RowPerm = (rowperm_t)indx;
438: PetscCall(PetscOptionsBool("-mat_superlu_replacetinypivot", "ReplaceTinyPivot", "None", (PetscBool)lu->options.ReplaceTinyPivot, &flg, &set));
439: if (set && flg) lu->options.ReplaceTinyPivot = YES;
440: PetscCall(PetscOptionsBool("-mat_superlu_printstat", "PrintStat", "None", (PetscBool)lu->options.PrintStat, &flg, &set));
441: if (set && flg) lu->options.PrintStat = YES;
442: PetscCall(PetscOptionsInt("-mat_superlu_lwork", "size of work array in bytes used by factorization", "None", lu->lwork, &lu->lwork, NULL));
443: if (lu->lwork > 0) {
444: /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/
445: PetscCall(PetscMalloc(lu->lwork, &lu->work));
446: } else if (lu->lwork != 0 && lu->lwork != -1) {
447: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Warning: lwork %" PetscInt_FMT " is not supported by SUPERLU. The default lwork=0 is used.\n", lu->lwork));
448: lu->lwork = 0;
449: }
450: /* ilu options */
451: PetscCall(PetscOptionsReal("-mat_superlu_ilu_droptol", "ILU_DropTol", "None", lu->options.ILU_DropTol, &real_input, &flg));
452: if (flg) lu->options.ILU_DropTol = (double)real_input;
453: PetscCall(PetscOptionsReal("-mat_superlu_ilu_filltol", "ILU_FillTol", "None", lu->options.ILU_FillTol, &real_input, &flg));
454: if (flg) lu->options.ILU_FillTol = (double)real_input;
455: PetscCall(PetscOptionsReal("-mat_superlu_ilu_fillfactor", "ILU_FillFactor", "None", lu->options.ILU_FillFactor, &real_input, &flg));
456: if (flg) lu->options.ILU_FillFactor = (double)real_input;
457: PetscCall(PetscOptionsInt("-mat_superlu_ilu_droprull", "ILU_DropRule", "None", lu->options.ILU_DropRule, &lu->options.ILU_DropRule, NULL));
458: PetscCall(PetscOptionsInt("-mat_superlu_ilu_norm", "ILU_Norm", "None", lu->options.ILU_Norm, &indx, &flg));
459: if (flg) lu->options.ILU_Norm = (norm_t)indx;
460: PetscCall(PetscOptionsInt("-mat_superlu_ilu_milu", "ILU_MILU", "None", lu->options.ILU_MILU, &indx, &flg));
461: if (flg) lu->options.ILU_MILU = (milu_t)indx;
462: PetscOptionsEnd();
464: lu->flg = DIFFERENT_NONZERO_PATTERN;
465: lu->CleanUpSuperLU = PETSC_TRUE;
466: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;
468: /* if we are here, the nonzero pattern has changed unless the user explicitly called MatLUFactorSymbolic */
469: PetscCall(MatDestroy(&lu->A_dup));
470: if (lu->needconversion) PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &lu->A_dup));
471: if (lu->options.Equil == YES && !lu->A_dup) { /* superlu overwrites input matrix and rhs when Equil is used, thus create A_dup to keep user's A unchanged */
472: PetscCall(MatDuplicate_SeqAIJ(A, MAT_COPY_VALUES, &lu->A_dup));
473: }
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F, PetscReal dtol)
478: {
479: Mat_SuperLU *lu = (Mat_SuperLU *)F->data;
481: PetscFunctionBegin;
482: lu->options.ILU_DropTol = dtol;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /*@
487: MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
489: Logically Collective
491: Input Parameters:
492: + F - the factored matrix obtained by calling `MatGetFactor()`
493: - dtol - drop tolerance
495: Options Database Key:
496: . -mat_superlu_ilu_droptol <dtol> - the drop tolerance
498: Level: beginner
500: References:
501: . * - SuperLU Users' Guide
503: .seealso: [](ch_matrices), `Mat`, `MatGetFactor()`, `MATSOLVERSUPERLU`
504: @*/
505: PetscErrorCode MatSuperluSetILUDropTol(Mat F, PetscReal dtol)
506: {
507: PetscFunctionBegin;
510: PetscTryMethod(F, "MatSuperluSetILUDropTol_C", (Mat, PetscReal), (F, dtol));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: PetscErrorCode MatFactorGetSolverType_seqaij_superlu(Mat A, MatSolverType *type)
515: {
516: PetscFunctionBegin;
517: *type = MATSOLVERSUPERLU;
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*MC
522: MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
523: via the external package SuperLU.
525: Use `./configure --download-superlu` to have PETSc installed with SuperLU
527: Use `-pc_type lu` `-pc_factor_mat_solver_type superlu` to use this direct solver
529: Options Database Keys:
530: + -mat_superlu_equil <FALSE> - Equil (None)
531: . -mat_superlu_colperm <COLAMD> - (choose one of) `NATURAL`, `MMD_ATA MMD_AT_PLUS_A`, `COLAMD`
532: . -mat_superlu_iterrefine <NOREFINE> - (choose one of) `NOREFINE`, `SINGLE`, `DOUBLE`, `EXTRA`
533: . -mat_superlu_symmetricmode: <FALSE> - SymmetricMode (None)
534: . -mat_superlu_diagpivotthresh <1> - DiagPivotThresh (None)
535: . -mat_superlu_pivotgrowth <FALSE> - PivotGrowth (None)
536: . -mat_superlu_conditionnumber <FALSE> - ConditionNumber (None)
537: . -mat_superlu_rowperm <NOROWPERM> - (choose one of) `NOROWPERM`, `LargeDiag`
538: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
539: . -mat_superlu_printstat <FALSE> - PrintStat (None)
540: . -mat_superlu_lwork <0> - size of work array in bytes used by factorization (None)
541: . -mat_superlu_ilu_droptol <0> - ILU_DropTol (None)
542: . -mat_superlu_ilu_filltol <0> - ILU_FillTol (None)
543: . -mat_superlu_ilu_fillfactor <0> - ILU_FillFactor (None)
544: . -mat_superlu_ilu_droprull <0> - ILU_DropRule (None)
545: . -mat_superlu_ilu_norm <0> - ILU_Norm (None)
546: - -mat_superlu_ilu_milu <0> - ILU_MILU (None)
548: Level: beginner
550: Notes:
551: Do not confuse this with `MATSOLVERSUPERLU_DIST` which is for parallel sparse solves
553: Cannot use ordering provided by PETSc, provides its own.
555: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCILU`, `MATSOLVERSUPERLU_DIST`, `MATSOLVERMUMPS`, `PCFactorSetMatSolverType()`, `MatSolverType`
556: M*/
558: static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A, MatFactorType ftype, Mat *F)
559: {
560: Mat B;
561: Mat_SuperLU *lu;
562: PetscInt m = A->rmap->n, n = A->cmap->n;
564: PetscFunctionBegin;
565: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
566: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, PETSC_DETERMINE, PETSC_DETERMINE));
567: PetscCall(PetscStrallocpy("superlu", &((PetscObject)B)->type_name));
568: PetscCall(MatSetUp(B));
569: B->trivialsymbolic = PETSC_TRUE;
570: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
571: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU;
572: B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
573: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Factor type not supported");
575: PetscCall(PetscFree(B->solvertype));
576: PetscCall(PetscStrallocpy(MATSOLVERSUPERLU, &B->solvertype));
578: B->ops->getinfo = MatGetInfo_External;
579: B->ops->destroy = MatDestroy_SuperLU;
580: B->ops->view = MatView_SuperLU;
581: B->factortype = ftype;
582: B->assembled = PETSC_TRUE; /* required by -ksp_view */
583: B->preallocated = PETSC_TRUE;
585: PetscCall(PetscNew(&lu));
587: if (ftype == MAT_FACTOR_LU) {
588: set_default_options(&lu->options);
589: /* Comments from SuperLU_4.0/SRC/dgssvx.c:
590: "Whether or not the system will be equilibrated depends on the
591: scaling of the matrix A, but if equilibration is used, A is
592: overwritten by diag(R)*A*diag(C) and B by diag(R)*B
593: (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
594: We set 'options.Equil = NO' as default because additional space is needed for it.
595: */
596: lu->options.Equil = NO;
597: } else if (ftype == MAT_FACTOR_ILU) {
598: /* Set the default input options of ilu: */
599: PetscStackCallExternalVoid("SuperLU:ilu_set_default_options", ilu_set_default_options(&lu->options));
600: }
601: lu->options.PrintStat = NO;
603: /* Initialize the statistics variables. */
604: PetscStackCallExternalVoid("SuperLU:StatInit", StatInit(&lu->stat));
605: lu->lwork = 0; /* allocate space internally by system malloc */
607: /* Allocate spaces (notice sizes are for the transpose) */
608: PetscCall(PetscMalloc1(m, &lu->etree));
609: PetscCall(PetscMalloc1(n, &lu->perm_r));
610: PetscCall(PetscMalloc1(m, &lu->perm_c));
611: PetscCall(PetscMalloc1(n, &lu->R));
612: PetscCall(PetscMalloc1(m, &lu->C));
614: /* create rhs and solution x without allocate space for .Store */
615: #if defined(PETSC_USE_COMPLEX)
616: #if defined(PETSC_USE_REAL_SINGLE)
617: PetscStackCallExternalVoid("SuperLU:cCreate_Dense_Matrix(", cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
618: PetscStackCallExternalVoid("SuperLU:cCreate_Dense_Matrix(", cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
619: #else
620: PetscStackCallExternalVoid("SuperLU:zCreate_Dense_Matrix", zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
621: PetscStackCallExternalVoid("SuperLU:zCreate_Dense_Matrix", zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
622: #endif
623: #else
624: #if defined(PETSC_USE_REAL_SINGLE)
625: PetscStackCallExternalVoid("SuperLU:sCreate_Dense_Matrix", sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
626: PetscStackCallExternalVoid("SuperLU:sCreate_Dense_Matrix", sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
627: #else
628: PetscStackCallExternalVoid("SuperLU:dCreate_Dense_Matrix", dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
629: PetscStackCallExternalVoid("SuperLU:dCreate_Dense_Matrix", dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
630: #endif
631: #endif
633: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_superlu));
634: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSuperluSetILUDropTol_C", MatSuperluSetILUDropTol_SuperLU));
635: B->data = lu;
637: *F = B;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: static PetscErrorCode MatGetFactor_seqsell_superlu(Mat A, MatFactorType ftype, Mat *F)
642: {
643: Mat_SuperLU *lu;
645: PetscFunctionBegin;
646: PetscCall(MatGetFactor_seqaij_superlu(A, ftype, F));
647: lu = (Mat_SuperLU *)((*F)->data);
648: lu->needconversion = PETSC_TRUE;
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuperLU(void)
653: {
654: PetscFunctionBegin;
655: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_superlu));
656: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU, MATSEQAIJ, MAT_FACTOR_ILU, MatGetFactor_seqaij_superlu));
657: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU, MATSEQSELL, MAT_FACTOR_LU, MatGetFactor_seqsell_superlu));
658: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU, MATSEQSELL, MAT_FACTOR_ILU, MatGetFactor_seqsell_superlu));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }