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