Actual source code: superlu.c

petsc-master 2018-05-20
Report Typos and Errors

  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;
 53:   PetscBool         needconversion;

 55:   /* Flag to clean up (non-global) SuperLU objects during Destroy */
 56:   PetscBool CleanUpSuperLU;
 57: } Mat_SuperLU;

 59: /*
 60:     Utility function
 61: */
 62: static PetscErrorCode MatView_Info_SuperLU(Mat A,PetscViewer viewer)
 63: {
 64:   Mat_SuperLU       *lu= (Mat_SuperLU*)A->data;
 65:   PetscErrorCode    ierr;
 66:   superlu_options_t options;

 69:   options = lu->options;

 71:   PetscViewerASCIIPrintf(viewer,"SuperLU run parameters:\n");
 72:   PetscViewerASCIIPrintf(viewer,"  Equil: %s\n",(options.Equil != NO) ? "YES" : "NO");
 73:   PetscViewerASCIIPrintf(viewer,"  ColPerm: %D\n",options.ColPerm);
 74:   PetscViewerASCIIPrintf(viewer,"  IterRefine: %D\n",options.IterRefine);
 75:   PetscViewerASCIIPrintf(viewer,"  SymmetricMode: %s\n",(options.SymmetricMode != NO) ? "YES" : "NO");
 76:   PetscViewerASCIIPrintf(viewer,"  DiagPivotThresh: %g\n",options.DiagPivotThresh);
 77:   PetscViewerASCIIPrintf(viewer,"  PivotGrowth: %s\n",(options.PivotGrowth != NO) ? "YES" : "NO");
 78:   PetscViewerASCIIPrintf(viewer,"  ConditionNumber: %s\n",(options.ConditionNumber != NO) ? "YES" : "NO");
 79:   PetscViewerASCIIPrintf(viewer,"  RowPerm: %D\n",options.RowPerm);
 80:   PetscViewerASCIIPrintf(viewer,"  ReplaceTinyPivot: %s\n",(options.ReplaceTinyPivot != NO) ? "YES" : "NO");
 81:   PetscViewerASCIIPrintf(viewer,"  PrintStat: %s\n",(options.PrintStat != NO) ? "YES" : "NO");
 82:   PetscViewerASCIIPrintf(viewer,"  lwork: %D\n",lu->lwork);
 83:   if (A->factortype == MAT_FACTOR_ILU) {
 84:     PetscViewerASCIIPrintf(viewer,"  ILU_DropTol: %g\n",options.ILU_DropTol);
 85:     PetscViewerASCIIPrintf(viewer,"  ILU_FillTol: %g\n",options.ILU_FillTol);
 86:     PetscViewerASCIIPrintf(viewer,"  ILU_FillFactor: %g\n",options.ILU_FillFactor);
 87:     PetscViewerASCIIPrintf(viewer,"  ILU_DropRule: %D\n",options.ILU_DropRule);
 88:     PetscViewerASCIIPrintf(viewer,"  ILU_Norm: %D\n",options.ILU_Norm);
 89:     PetscViewerASCIIPrintf(viewer,"  ILU_MILU: %D\n",options.ILU_MILU);
 90:   }
 91:   return(0);
 92: }

 94: PetscErrorCode MatSolve_SuperLU_Private(Mat A,Vec b,Vec x)
 95: {
 96:   Mat_SuperLU       *lu = (Mat_SuperLU*)A->data;
 97:   const PetscScalar *barray;
 98:   PetscScalar       *xarray;
 99:   PetscErrorCode    ierr;
100:   PetscInt          info,i,n;
101:   PetscReal         ferr,berr;
102:   static PetscBool  cite = PETSC_FALSE;

105:   if (lu->lwork == -1) return(0);
106:   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);

108:   VecGetLocalSize(x,&n);
109:   lu->B.ncol = 1;   /* Set the number of right-hand side */
110:   if (lu->options.Equil && !lu->rhs_dup) {
111:     /* superlu overwrites b when Equil is used, thus create rhs_dup to keep user's b unchanged */
112:     PetscMalloc1(n,&lu->rhs_dup);
113:   }
114:   if (lu->options.Equil) {
115:     /* Copy b into rsh_dup */
116:     VecGetArrayRead(b,&barray);
117:     PetscMemcpy(lu->rhs_dup,barray,n*sizeof(PetscScalar));
118:     VecRestoreArrayRead(b,&barray);
119:     barray = lu->rhs_dup;
120:   } else {
121:     VecGetArrayRead(b,&barray);
122:   }
123:   VecGetArray(x,&xarray);

125: #if defined(PETSC_USE_COMPLEX)
126: #if defined(PETSC_USE_REAL_SINGLE)
127:   ((DNformat*)lu->B.Store)->nzval = (singlecomplex*)barray;
128:   ((DNformat*)lu->X.Store)->nzval = (singlecomplex*)xarray;
129: #else
130:   ((DNformat*)lu->B.Store)->nzval = (doublecomplex*)barray;
131:   ((DNformat*)lu->X.Store)->nzval = (doublecomplex*)xarray;
132: #endif
133: #else
134:   ((DNformat*)lu->B.Store)->nzval = (void*)barray;
135:   ((DNformat*)lu->X.Store)->nzval = xarray;
136: #endif

138:   lu->options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
139:   if (A->factortype == MAT_FACTOR_LU) {
140: #if defined(PETSC_USE_COMPLEX)
141: #if defined(PETSC_USE_REAL_SINGLE)
142:     PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
143:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
144:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
145: #else
146:     PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
147:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
148:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
149: #endif
150: #else
151: #if defined(PETSC_USE_REAL_SINGLE)
152:     PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
153:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
154:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
155: #else
156:     PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
157:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
158:                                      &lu->Glu,&lu->mem_usage, &lu->stat, &info));
159: #endif
160: #endif
161:   } else if (A->factortype == MAT_FACTOR_ILU) {
162: #if defined(PETSC_USE_COMPLEX)
163: #if defined(PETSC_USE_REAL_SINGLE)
164:     PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
165:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
166:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
167: #else
168:     PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
169:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
170:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
171: #endif
172: #else
173: #if defined(PETSC_USE_REAL_SINGLE)
174:     PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
175:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
176:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
177: #else
178:     PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
179:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
180:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &info));
181: #endif
182: #endif
183:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
184:   if (!lu->options.Equil) {
185:     VecRestoreArrayRead(b,&barray);
186:   }
187:   VecRestoreArray(x,&xarray);

189:   if (!info || info == lu->A.ncol+1) {
190:     if (lu->options.IterRefine) {
191:       PetscPrintf(PETSC_COMM_SELF,"Iterative Refinement:\n");
192:       PetscPrintf(PETSC_COMM_SELF,"  %8s%8s%16s%16s\n", "rhs", "Steps", "FERR", "BERR");
193:       for (i = 0; i < 1; ++i) {
194:         PetscPrintf(PETSC_COMM_SELF,"  %8d%8d%16e%16e\n", i+1, lu->stat.RefineSteps, ferr, berr);
195:       }
196:     }
197:   } else if (info > 0) {
198:     if (lu->lwork == -1) {
199:       PetscPrintf(PETSC_COMM_SELF,"  ** Estimated memory: %D bytes\n", info - lu->A.ncol);
200:     } else {
201:       PetscPrintf(PETSC_COMM_SELF,"  Warning: gssvx() returns info %D\n",info);
202:     }
203:   } 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);

205:   if (lu->options.PrintStat) {
206:     PetscPrintf(PETSC_COMM_SELF,"MatSolve__SuperLU():\n");
207:     PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
208:   }
209:   return(0);
210: }

212: PetscErrorCode MatSolve_SuperLU(Mat A,Vec b,Vec x)
213: {
214:   Mat_SuperLU    *lu = (Mat_SuperLU*)A->data;

218:   if (A->factorerrortype) {
219:     PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
220:     VecSetInf(x);
221:     return(0);
222:   }

224:   lu->options.Trans = TRANS;
225:   MatSolve_SuperLU_Private(A,b,x);
226:   return(0);
227: }

229: PetscErrorCode MatSolveTranspose_SuperLU(Mat A,Vec b,Vec x)
230: {
231:   Mat_SuperLU    *lu = (Mat_SuperLU*)A->data;

235:   if (A->factorerrortype) {
236:     PetscInfo(A,"MatSolve is called with singular matrix factor, skip\n");
237:     VecSetInf(x);
238:     return(0);
239:   }

241:   lu->options.Trans = NOTRANS;
242:   MatSolve_SuperLU_Private(A,b,x);
243:   return(0);
244: }

246: static PetscErrorCode MatLUFactorNumeric_SuperLU(Mat F,Mat A,const MatFactorInfo *info)
247: {
248:   Mat_SuperLU    *lu = (Mat_SuperLU*)F->data;
249:   Mat_SeqAIJ     *aa;
251:   PetscInt       sinfo;
252:   PetscReal      ferr, berr;
253:   NCformat       *Ustore;
254:   SCformat       *Lstore;

257:   if (lu->flg == SAME_NONZERO_PATTERN) { /* successing numerical factorization */
258:     lu->options.Fact = SamePattern;
259:     /* Ref: ~SuperLU_3.0/EXAMPLE/dlinsolx2.c */
260:     Destroy_SuperMatrix_Store(&lu->A);
261:     if (lu->A_dup) {
262:       MatCopy_SeqAIJ(A,lu->A_dup,SAME_NONZERO_PATTERN);
263:     }

265:     if (lu->lwork >= 0) {
266:       PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
267:       PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
268:       lu->options.Fact = SamePattern;
269:     }
270:   }

272:   /* Create the SuperMatrix for lu->A=A^T:
273:        Since SuperLU likes column-oriented matrices,we pass it the transpose,
274:        and then solve A^T X = B in MatSolve(). */
275:   if (lu->A_dup) {
276:     aa = (Mat_SeqAIJ*)(lu->A_dup)->data;
277:   } else {
278:     aa = (Mat_SeqAIJ*)(A)->data;
279:   }
280: #if defined(PETSC_USE_COMPLEX)
281: #if defined(PETSC_USE_REAL_SINGLE)
282:   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));
283: #else
284:   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));
285: #endif
286: #else
287: #if defined(PETSC_USE_REAL_SINGLE)
288:   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));
289: #else
290:   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));
291: #endif
292: #endif

294:   /* Numerical factorization */
295:   lu->B.ncol = 0;  /* Indicate not to solve the system */
296:   if (F->factortype == MAT_FACTOR_LU) {
297: #if defined(PETSC_USE_COMPLEX)
298: #if defined(PETSC_USE_REAL_SINGLE)
299:     PetscStackCall("SuperLU:cgssvx",cgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
300:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
301:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
302: #else
303:     PetscStackCall("SuperLU:zgssvx",zgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
304:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
305:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
306: #endif
307: #else
308: #if defined(PETSC_USE_REAL_SINGLE)
309:     PetscStackCall("SuperLU:sgssvx",sgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
310:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
311:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
312: #else
313:     PetscStackCall("SuperLU:dgssvx",dgssvx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
314:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond, &ferr, &berr,
315:                                      &lu->Glu,&lu->mem_usage, &lu->stat, &sinfo));
316: #endif
317: #endif
318:   } else if (F->factortype == MAT_FACTOR_ILU) {
319:     /* Compute the incomplete factorization, condition number and pivot growth */
320: #if defined(PETSC_USE_COMPLEX)
321: #if defined(PETSC_USE_REAL_SINGLE)
322:     PetscStackCall("SuperLU:cgsisx",cgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
323:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
324:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
325: #else
326:     PetscStackCall("SuperLU:zgsisx",zgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r,lu->etree, lu->equed, lu->R, lu->C,
327:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
328:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
329: #endif
330: #else
331: #if defined(PETSC_USE_REAL_SINGLE)
332:     PetscStackCall("SuperLU:sgsisx",sgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
333:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
334:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
335: #else
336:     PetscStackCall("SuperLU:dgsisx",dgsisx(&lu->options, &lu->A, lu->perm_c, lu->perm_r, lu->etree, lu->equed, lu->R, lu->C,
337:                                      &lu->L, &lu->U, lu->work, lu->lwork, &lu->B, &lu->X, &lu->rpg, &lu->rcond,
338:                                      &lu->Glu, &lu->mem_usage, &lu->stat, &sinfo));
339: #endif
340: #endif
341:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");
342:   if (!sinfo || sinfo == lu->A.ncol+1) {
343:     if (lu->options.PivotGrowth) {
344:       PetscPrintf(PETSC_COMM_SELF,"  Recip. pivot growth = %e\n", lu->rpg);
345:     }
346:     if (lu->options.ConditionNumber) {
347:       PetscPrintf(PETSC_COMM_SELF,"  Recip. condition number = %e\n", lu->rcond);
348:     }
349:   } else if (sinfo > 0) {
350:     if (A->erroriffailure) {
351:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %D",sinfo);
352:     } else {
353:       if (sinfo <= lu->A.ncol) {
354:         if (lu->options.ILU_FillTol == 0.0) {
355:           F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
356:         }
357:         PetscInfo2(F,"Number of zero pivots %D, ILU_FillTol %g\n",sinfo,lu->options.ILU_FillTol);
358:       } else if (sinfo == lu->A.ncol + 1) {
359:         /*
360:          U is nonsingular, but RCOND is less than machine
361:                        precision, meaning that the matrix is singular to
362:                        working precision. Nevertheless, the solution and
363:                        error bounds are computed because there are a number
364:                        of situations where the computed solution can be more
365:                        accurate than the value of RCOND would suggest.
366:          */
367:         PetscInfo1(F,"Matrix factor U is nonsingular, but is singular to working precision. The solution is computed. info %D",sinfo);
368:       } else { /* sinfo > lu->A.ncol + 1 */
369:         F->factorerrortype = MAT_FACTOR_OUTMEMORY;
370:         PetscInfo1(F,"Number of bytes allocated when memory allocation fails %D\n",sinfo);
371:       }
372:     }
373:   } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, the %D-th argument in gssvx() had an illegal value", sinfo,-sinfo);

375:   if (lu->options.PrintStat) {
376:     PetscPrintf(PETSC_COMM_SELF,"MatLUFactorNumeric_SuperLU():\n");
377:     PetscStackCall("SuperLU:StatPrint",StatPrint(&lu->stat));
378:     Lstore = (SCformat*) lu->L.Store;
379:     Ustore = (NCformat*) lu->U.Store;
380:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in factor L = %D\n", Lstore->nnz);
381:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in factor U = %D\n", Ustore->nnz);
382:     PetscPrintf(PETSC_COMM_SELF,"  No of nonzeros in L+U = %D\n", Lstore->nnz + Ustore->nnz - lu->A.ncol);
383:     PetscPrintf(PETSC_COMM_SELF,"  L\\U MB %.3f\ttotal MB needed %.3f\n",
384:                          lu->mem_usage.for_lu/1e6, lu->mem_usage.total_needed/1e6);
385:   }

387:   lu->flg                = SAME_NONZERO_PATTERN;
388:   F->ops->solve          = MatSolve_SuperLU;
389:   F->ops->solvetranspose = MatSolveTranspose_SuperLU;
390:   F->ops->matsolve       = NULL;
391:   return(0);
392: }

394: static PetscErrorCode MatDestroy_SuperLU(Mat A)
395: {
397:   Mat_SuperLU    *lu=(Mat_SuperLU*)A->data;

400:   if (lu->CleanUpSuperLU) { /* Free the SuperLU datastructures */
401:     PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->A));
402:     PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->B));
403:     PetscStackCall("SuperLU:Destroy_SuperMatrix_Store",Destroy_SuperMatrix_Store(&lu->X));
404:     PetscStackCall("SuperLU:StatFree",StatFree(&lu->stat));
405:     if (lu->lwork >= 0) {
406:       PetscStackCall("SuperLU:Destroy_SuperNode_Matrix",Destroy_SuperNode_Matrix(&lu->L));
407:       PetscStackCall("SuperLU:Destroy_CompCol_Matrix",Destroy_CompCol_Matrix(&lu->U));
408:     }
409:   }
410:   PetscFree(lu->etree);
411:   PetscFree(lu->perm_r);
412:   PetscFree(lu->perm_c);
413:   PetscFree(lu->R);
414:   PetscFree(lu->C);
415:   PetscFree(lu->rhs_dup);
416:   MatDestroy(&lu->A_dup);
417:   PetscFree(A->data);

419:   /* clear composed functions */
420:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
421:   PetscObjectComposeFunction((PetscObject)A,"MatSuperluSetILUDropTol_C",NULL);
422:   return(0);
423: }

425: static PetscErrorCode MatView_SuperLU(Mat A,PetscViewer viewer)
426: {
427:   PetscErrorCode    ierr;
428:   PetscBool         iascii;
429:   PetscViewerFormat format;

432:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
433:   if (iascii) {
434:     PetscViewerGetFormat(viewer,&format);
435:     if (format == PETSC_VIEWER_ASCII_INFO) {
436:       MatView_Info_SuperLU(A,viewer);
437:     }
438:   }
439:   return(0);
440: }

442: PetscErrorCode MatMatSolve_SuperLU(Mat A,Mat B,Mat X)
443: {
444:   Mat_SuperLU    *lu = (Mat_SuperLU*)A->data;
445:   PetscBool      flg;

449:   PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
450:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
451:   PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
452:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
453:   lu->options.Trans = TRANS;
454:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatSolve_SuperLU() is not implemented yet");
455:   return(0);
456: }

458: /*
459:    Note the r permutation is ignored
460: */
461: static PetscErrorCode MatLUFactorSymbolic_SuperLU(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
462: {
463:   Mat_SuperLU    *lu = (Mat_SuperLU*)(F->data);

467:   lu->flg                 = DIFFERENT_NONZERO_PATTERN;
468:   lu->CleanUpSuperLU      = PETSC_TRUE;
469:   F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU;

471:   /* if we are here, the nonzero pattern has changed unless the user explicitly called MatLUFactorSymbolic */
472:   MatDestroy(&lu->A_dup);
473:   if (lu->needconversion) {
474:     MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&lu->A_dup);
475:   }
476:   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 */
477:     MatDuplicate_SeqAIJ(A,MAT_COPY_VALUES,&lu->A_dup);
478:   }
479:   return(0);
480: }

482: static PetscErrorCode MatSuperluSetILUDropTol_SuperLU(Mat F,PetscReal dtol)
483: {
484:   Mat_SuperLU *lu= (Mat_SuperLU*)F->data;

487:   lu->options.ILU_DropTol = dtol;
488:   return(0);
489: }

491: /*@
492:   MatSuperluSetILUDropTol - Set SuperLU ILU drop tolerance
493:    Logically Collective on Mat

495:    Input Parameters:
496: +  F - the factored matrix obtained by calling MatGetFactor() from PETSc-SuperLU interface
497: -  dtol - drop tolerance

499:   Options Database:
500: .   -mat_superlu_ilu_droptol <dtol>

502:    Level: beginner

504:    References:
505: .      SuperLU Users' Guide

507: .seealso: MatGetFactor()
508: @*/
509: PetscErrorCode MatSuperluSetILUDropTol(Mat F,PetscReal dtol)
510: {

516:   PetscTryMethod(F,"MatSuperluSetILUDropTol_C",(Mat,PetscReal),(F,dtol));
517:   return(0);
518: }

520: PetscErrorCode MatFactorGetSolverType_seqaij_superlu(Mat A,MatSolverType *type)
521: {
523:   *type = MATSOLVERSUPERLU;
524:   return(0);
525: }

527: /*MC
528:   MATSOLVERSUPERLU = "superlu" - A solver package providing solvers LU and ILU for sequential matrices
529:   via the external package SuperLU.

531:   Use ./configure --download-superlu to have PETSc installed with SuperLU

533:   Use -pc_type lu -pc_factor_mat_solver_type superlu to use this direct solver

535:   Options Database Keys:
536: + -mat_superlu_equil <FALSE>            - Equil (None)
537: . -mat_superlu_colperm <COLAMD>         - (choose one of) NATURAL MMD_ATA MMD_AT_PLUS_A COLAMD
538: . -mat_superlu_iterrefine <NOREFINE>    - (choose one of) NOREFINE SINGLE DOUBLE EXTRA
539: . -mat_superlu_symmetricmode: <FALSE>   - SymmetricMode (None)
540: . -mat_superlu_diagpivotthresh <1>      - DiagPivotThresh (None)
541: . -mat_superlu_pivotgrowth <FALSE>      - PivotGrowth (None)
542: . -mat_superlu_conditionnumber <FALSE>  - ConditionNumber (None)
543: . -mat_superlu_rowperm <NOROWPERM>      - (choose one of) NOROWPERM LargeDiag
544: . -mat_superlu_replacetinypivot <FALSE> - ReplaceTinyPivot (None)
545: . -mat_superlu_printstat <FALSE>        - PrintStat (None)
546: . -mat_superlu_lwork <0>                - size of work array in bytes used by factorization (None)
547: . -mat_superlu_ilu_droptol <0>          - ILU_DropTol (None)
548: . -mat_superlu_ilu_filltol <0>          - ILU_FillTol (None)
549: . -mat_superlu_ilu_fillfactor <0>       - ILU_FillFactor (None)
550: . -mat_superlu_ilu_droprull <0>         - ILU_DropRule (None)
551: . -mat_superlu_ilu_norm <0>             - ILU_Norm (None)
552: - -mat_superlu_ilu_milu <0>             - ILU_MILU (None)

554:    Notes:
555:     Do not confuse this with MATSOLVERSUPERLU_DIST which is for parallel sparse solves

557:    Level: beginner

559: .seealso: PCLU, PCILU, MATSOLVERSUPERLU_DIST, MATSOLVERMUMPS, PCFactorSetMatSolverType(), MatSolverType
560: M*/

562: static PetscErrorCode MatGetFactor_seqaij_superlu(Mat A,MatFactorType ftype,Mat *F)
563: {
564:   Mat            B;
565:   Mat_SuperLU    *lu;
567:   PetscInt       indx,m=A->rmap->n,n=A->cmap->n;
568:   PetscBool      flg,set;
569:   PetscReal      real_input;
570:   const char     *colperm[]   ={"NATURAL","MMD_ATA","MMD_AT_PLUS_A","COLAMD"}; /* MY_PERMC - not supported by the petsc interface yet */
571:   const char     *iterrefine[]={"NOREFINE", "SINGLE", "DOUBLE", "EXTRA"};
572:   const char     *rowperm[]   ={"NOROWPERM", "LargeDiag"}; /* MY_PERMC - not supported by the petsc interface yet */

575:   MatCreate(PetscObjectComm((PetscObject)A),&B);
576:   MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DETERMINE,PETSC_DETERMINE);
577:   PetscStrallocpy("superlu",&((PetscObject)B)->type_name);
578:   MatSetUp(B);
579:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
580:     B->ops->lufactorsymbolic  = MatLUFactorSymbolic_SuperLU;
581:     B->ops->ilufactorsymbolic = MatLUFactorSymbolic_SuperLU;
582:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported");

584:   PetscFree(B->solvertype);
585:   PetscStrallocpy(MATSOLVERSUPERLU,&B->solvertype);

587:   B->ops->getinfo     = MatGetInfo_External;
588:   B->ops->destroy     = MatDestroy_SuperLU;
589:   B->ops->view        = MatView_SuperLU;
590:   B->factortype       = ftype;
591:   B->assembled        = PETSC_TRUE;           /* required by -ksp_view */
592:   B->preallocated     = PETSC_TRUE;

594:   PetscNewLog(B,&lu);

596:   if (ftype == MAT_FACTOR_LU) {
597:     set_default_options(&lu->options);
598:     /* Comments from SuperLU_4.0/SRC/dgssvx.c:
599:       "Whether or not the system will be equilibrated depends on the
600:        scaling of the matrix A, but if equilibration is used, A is
601:        overwritten by diag(R)*A*diag(C) and B by diag(R)*B
602:        (if options->Trans=NOTRANS) or diag(C)*B (if options->Trans = TRANS or CONJ)."
603:      We set 'options.Equil = NO' as default because additional space is needed for it.
604:     */
605:     lu->options.Equil = NO;
606:   } else if (ftype == MAT_FACTOR_ILU) {
607:     /* Set the default input options of ilu: */
608:     PetscStackCall("SuperLU:ilu_set_default_options",ilu_set_default_options(&lu->options));
609:   }
610:   lu->options.PrintStat = NO;

612:   /* Initialize the statistics variables. */
613:   PetscStackCall("SuperLU:StatInit",StatInit(&lu->stat));
614:   lu->lwork = 0;   /* allocate space internally by system malloc */

616:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"SuperLU Options","Mat");
617:   PetscOptionsBool("-mat_superlu_equil","Equil","None",(PetscBool)lu->options.Equil,(PetscBool*)&lu->options.Equil,NULL);
618:   PetscOptionsEList("-mat_superlu_colperm","ColPerm","None",colperm,4,colperm[3],&indx,&flg);
619:   if (flg) lu->options.ColPerm = (colperm_t)indx;
620:   PetscOptionsEList("-mat_superlu_iterrefine","IterRefine","None",iterrefine,4,iterrefine[0],&indx,&flg);
621:   if (flg) lu->options.IterRefine = (IterRefine_t)indx;
622:   PetscOptionsBool("-mat_superlu_symmetricmode","SymmetricMode","None",(PetscBool)lu->options.SymmetricMode,&flg,&set);
623:   if (set && flg) lu->options.SymmetricMode = YES;
624:   PetscOptionsReal("-mat_superlu_diagpivotthresh","DiagPivotThresh","None",lu->options.DiagPivotThresh,&real_input,&flg);
625:   if (flg) lu->options.DiagPivotThresh = (double) real_input;
626:   PetscOptionsBool("-mat_superlu_pivotgrowth","PivotGrowth","None",(PetscBool)lu->options.PivotGrowth,&flg,&set);
627:   if (set && flg) lu->options.PivotGrowth = YES;
628:   PetscOptionsBool("-mat_superlu_conditionnumber","ConditionNumber","None",(PetscBool)lu->options.ConditionNumber,&flg,&set);
629:   if (set && flg) lu->options.ConditionNumber = YES;
630:   PetscOptionsEList("-mat_superlu_rowperm","rowperm","None",rowperm,2,rowperm[lu->options.RowPerm],&indx,&flg);
631:   if (flg) lu->options.RowPerm = (rowperm_t)indx;
632:   PetscOptionsBool("-mat_superlu_replacetinypivot","ReplaceTinyPivot","None",(PetscBool)lu->options.ReplaceTinyPivot,&flg,&set);
633:   if (set && flg) lu->options.ReplaceTinyPivot = YES;
634:   PetscOptionsBool("-mat_superlu_printstat","PrintStat","None",(PetscBool)lu->options.PrintStat,&flg,&set);
635:   if (set && flg) lu->options.PrintStat = YES;
636:   PetscOptionsInt("-mat_superlu_lwork","size of work array in bytes used by factorization","None",lu->lwork,&lu->lwork,NULL);
637:   if (lu->lwork > 0) {
638:     /* lwork is in bytes, hence PetscMalloc() is used here, not PetscMalloc1()*/
639:     PetscMalloc(lu->lwork,&lu->work);
640:   } else if (lu->lwork != 0 && lu->lwork != -1) {
641:     PetscPrintf(PETSC_COMM_SELF,"   Warning: lwork %D is not supported by SUPERLU. The default lwork=0 is used.\n",lu->lwork);
642:     lu->lwork = 0;
643:   }
644:   /* ilu options */
645:   PetscOptionsReal("-mat_superlu_ilu_droptol","ILU_DropTol","None",lu->options.ILU_DropTol,&real_input,&flg);
646:   if (flg) lu->options.ILU_DropTol = (double) real_input;
647:   PetscOptionsReal("-mat_superlu_ilu_filltol","ILU_FillTol","None",lu->options.ILU_FillTol,&real_input,&flg);
648:   if (flg) lu->options.ILU_FillTol = (double) real_input;
649:   PetscOptionsReal("-mat_superlu_ilu_fillfactor","ILU_FillFactor","None",lu->options.ILU_FillFactor,&real_input,&flg);
650:   if (flg) lu->options.ILU_FillFactor = (double) real_input;
651:   PetscOptionsInt("-mat_superlu_ilu_droprull","ILU_DropRule","None",lu->options.ILU_DropRule,&lu->options.ILU_DropRule,NULL);
652:   PetscOptionsInt("-mat_superlu_ilu_norm","ILU_Norm","None",lu->options.ILU_Norm,&indx,&flg);
653:   if (flg) lu->options.ILU_Norm = (norm_t)indx;
654:   PetscOptionsInt("-mat_superlu_ilu_milu","ILU_MILU","None",lu->options.ILU_MILU,&indx,&flg);
655:   if (flg) lu->options.ILU_MILU = (milu_t)indx;
656:   PetscOptionsEnd();

658:   /* Allocate spaces (notice sizes are for the transpose) */
659:   PetscMalloc1(m,&lu->etree);
660:   PetscMalloc1(n,&lu->perm_r);
661:   PetscMalloc1(m,&lu->perm_c);
662:   PetscMalloc1(n,&lu->R);
663:   PetscMalloc1(m,&lu->C);

665:   /* create rhs and solution x without allocate space for .Store */
666: #if defined(PETSC_USE_COMPLEX)
667: #if defined(PETSC_USE_REAL_SINGLE)
668:   PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
669:   PetscStackCall("SuperLU:cCreate_Dense_Matrix(",cCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_C, SLU_GE));
670: #else
671:   PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
672:   PetscStackCall("SuperLU:zCreate_Dense_Matrix",zCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_Z, SLU_GE));
673: #endif
674: #else
675: #if defined(PETSC_USE_REAL_SINGLE)
676:   PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
677:   PetscStackCall("SuperLU:sCreate_Dense_Matrix",sCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_S, SLU_GE));
678: #else
679:   PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->B, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
680:   PetscStackCall("SuperLU:dCreate_Dense_Matrix",dCreate_Dense_Matrix(&lu->X, m, 1, NULL, m, SLU_DN, SLU_D, SLU_GE));
681: #endif
682: #endif

684:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqaij_superlu);
685:   PetscObjectComposeFunction((PetscObject)B,"MatSuperluSetILUDropTol_C",MatSuperluSetILUDropTol_SuperLU);
686:   B->data  = lu;

688:   *F       = B;
689:   return(0);
690: }

692: static PetscErrorCode MatGetFactor_seqsell_superlu(Mat A,MatFactorType ftype,Mat *F)
693: {
694:   Mat_SuperLU    *lu;

698:   MatGetFactor_seqaij_superlu(A,ftype,F);
699:   lu   = (Mat_SuperLU*)((*F)->data);
700:   lu->needconversion = PETSC_TRUE;
701:   return(0);
702: }

704: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuperLU(void)
705: {

709:   MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_seqaij_superlu);
710:   MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQAIJ,MAT_FACTOR_ILU,MatGetFactor_seqaij_superlu);
711:   MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_LU,MatGetFactor_seqsell_superlu);
712:   MatSolverTypeRegister(MATSOLVERSUPERLU,MATSEQSELL,MAT_FACTOR_ILU,MatGetFactor_seqsell_superlu);
713:   return(0);
714: }