Actual source code: superlu_dist.c

petsc-master 2019-12-10
Report Typos and Errors
  1: /*
  2:         Provides an interface to the SuperLU_DIST sparse solver
  3: */

  5: #include <../src/mat/impls/aij/seq/aij.h>
  6: #include <../src/mat/impls/aij/mpi/mpiaij.h>

  8: EXTERN_C_BEGIN
  9: #if defined(PETSC_USE_COMPLEX)
 10: #include <superlu_zdefs.h>
 11: #else
 12: #include <superlu_ddefs.h>
 13: #endif
 14: EXTERN_C_END

 16: typedef struct {
 17:   int_t                  nprow,npcol,*row,*col;
 18:   gridinfo_t             grid;
 19:   superlu_dist_options_t options;
 20:   SuperMatrix            A_sup;
 21:   ScalePermstruct_t      ScalePermstruct;
 22:   LUstruct_t             LUstruct;
 23:   int                    StatPrint;
 24:   SOLVEstruct_t          SOLVEstruct;
 25:   fact_t                 FactPattern;
 26:   MPI_Comm               comm_superlu;
 27: #if defined(PETSC_USE_COMPLEX)
 28:   doublecomplex          *val;
 29: #else
 30:   double                 *val;
 31: #endif
 32:   PetscBool              matsolve_iscalled,matmatsolve_iscalled;
 33:   PetscBool              CleanUpSuperLU_Dist;  /* Flag to clean up (non-global) SuperLU objects during Destroy */
 34: } Mat_SuperLU_DIST;


 37: PetscErrorCode MatSuperluDistGetDiagU_SuperLU_DIST(Mat F,PetscScalar *diagU)
 38: {
 39:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)F->data;

 42: #if defined(PETSC_USE_COMPLEX)
 43:   PetscStackCall("SuperLU_DIST:pzGetDiagU",pzGetDiagU(F->rmap->N,&lu->LUstruct,&lu->grid,(doublecomplex*)diagU));
 44: #else
 45:   PetscStackCall("SuperLU_DIST:pdGetDiagU",pdGetDiagU(F->rmap->N,&lu->LUstruct,&lu->grid,diagU));
 46: #endif
 47:   return(0);
 48: }

 50: PetscErrorCode MatSuperluDistGetDiagU(Mat F,PetscScalar *diagU)
 51: {

 56:   PetscTryMethod(F,"MatSuperluDistGetDiagU_C",(Mat,PetscScalar*),(F,diagU));
 57:   return(0);
 58: }

 60: static PetscErrorCode MatDestroy_SuperLU_DIST(Mat A)
 61: {
 62:   PetscErrorCode   ierr;
 63:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)A->data;

 66:   if (lu->CleanUpSuperLU_Dist) {
 67:     /* Deallocate SuperLU_DIST storage */
 68:     PetscStackCall("SuperLU_DIST:Destroy_CompRowLoc_Matrix_dist",Destroy_CompRowLoc_Matrix_dist(&lu->A_sup));
 69:     if (lu->options.SolveInitialized) {
 70: #if defined(PETSC_USE_COMPLEX)
 71:       PetscStackCall("SuperLU_DIST:zSolveFinalize",zSolveFinalize(&lu->options, &lu->SOLVEstruct));
 72: #else
 73:       PetscStackCall("SuperLU_DIST:dSolveFinalize",dSolveFinalize(&lu->options, &lu->SOLVEstruct));
 74: #endif
 75:     }
 76:     PetscStackCall("SuperLU_DIST:Destroy_LU",Destroy_LU(A->cmap->N, &lu->grid, &lu->LUstruct));
 77:     PetscStackCall("SuperLU_DIST:ScalePermstructFree",ScalePermstructFree(&lu->ScalePermstruct));
 78:     PetscStackCall("SuperLU_DIST:LUstructFree",LUstructFree(&lu->LUstruct));

 80:     /* Release the SuperLU_DIST process grid. */
 81:     PetscStackCall("SuperLU_DIST:superlu_gridexit",superlu_gridexit(&lu->grid));
 82:     MPI_Comm_free(&(lu->comm_superlu));
 83:   }
 84:   PetscFree(A->data);
 85:   /* clear composed functions */
 86:   PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverType_C",NULL);
 87:   PetscObjectComposeFunction((PetscObject)A,"MatSuperluDistGetDiagU_C",NULL);

 89:   return(0);
 90: }

 92: static PetscErrorCode MatSolve_SuperLU_DIST(Mat A,Vec b_mpi,Vec x)
 93: {
 94:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)A->data;
 95:   PetscErrorCode   ierr;
 96:   PetscInt         m=A->rmap->n;
 97:   SuperLUStat_t    stat;
 98:   double           berr[1];
 99:   PetscScalar      *bptr=NULL;
100:   int              info; /* SuperLU_Dist info code is ALWAYS an int, even with long long indices */
101:   static PetscBool cite = PETSC_FALSE;

104:   if (lu->options.Fact != FACTORED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"SuperLU_DIST options.Fact must equal FACTORED");
105:   PetscCitationsRegister("@article{lidemmel03,\n  author = {Xiaoye S. Li and James W. Demmel},\n  title = {{SuperLU_DIST}: A Scalable Distributed-Memory Sparse Direct\n           Solver for Unsymmetric Linear Systems},\n  journal = {ACM Trans. Mathematical Software},\n  volume = {29},\n  number = {2},\n  pages = {110-140},\n  year = 2003\n}\n",&cite);

107:   if (lu->options.SolveInitialized && !lu->matsolve_iscalled) {
108:     /* see comments in MatMatSolve() */
109: #if defined(PETSC_USE_COMPLEX)
110:     PetscStackCall("SuperLU_DIST:zSolveFinalize",zSolveFinalize(&lu->options, &lu->SOLVEstruct));
111: #else
112:     PetscStackCall("SuperLU_DIST:dSolveFinalize",dSolveFinalize(&lu->options, &lu->SOLVEstruct));
113: #endif
114:     lu->options.SolveInitialized = NO;
115:   }
116:   VecCopy(b_mpi,x);
117:   VecGetArray(x,&bptr);

119:   PetscStackCall("SuperLU_DIST:PStatInit",PStatInit(&stat));        /* Initialize the statistics variables. */
120: #if defined(PETSC_USE_COMPLEX)
121:   PetscStackCall("SuperLU_DIST:pzgssvx",pzgssvx(&lu->options,&lu->A_sup,&lu->ScalePermstruct,(doublecomplex*)bptr,m,1,&lu->grid,&lu->LUstruct,&lu->SOLVEstruct,berr,&stat,&info));
122: #else
123:   PetscStackCall("SuperLU_DIST:pdgssvx",pdgssvx(&lu->options,&lu->A_sup,&lu->ScalePermstruct,bptr,m,1,&lu->grid,&lu->LUstruct,&lu->SOLVEstruct,berr,&stat,&info));
124: #endif
125:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"pdgssvx fails, info: %d\n",info);

127:   if (lu->options.PrintStat) PetscStackCall("SuperLU_DIST:PStatPrint",PStatPrint(&lu->options, &stat, &lu->grid));  /* Print the statistics. */
128:   PetscStackCall("SuperLU_DIST:PStatFree",PStatFree(&stat));

130:   VecRestoreArray(x,&bptr);
131:   lu->matsolve_iscalled    = PETSC_TRUE;
132:   lu->matmatsolve_iscalled = PETSC_FALSE;
133:   return(0);
134: }

136: static PetscErrorCode MatMatSolve_SuperLU_DIST(Mat A,Mat B_mpi,Mat X)
137: {
138:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)A->data;
139:   PetscErrorCode   ierr;
140:   PetscInt         m=A->rmap->n,nrhs;
141:   SuperLUStat_t    stat;
142:   double           berr[1];
143:   PetscScalar      *bptr;
144:   int              info; /* SuperLU_Dist info code is ALWAYS an int, even with long long indices */
145:   PetscBool        flg;

148:   if (lu->options.Fact != FACTORED) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"SuperLU_DIST options.Fact must equal FACTORED");
149:   PetscObjectTypeCompareAny((PetscObject)B_mpi,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
150:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix");
151:   if (X != B_mpi) {
152:     PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);
153:     if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix");
154:   }

156:   if (lu->options.SolveInitialized && !lu->matmatsolve_iscalled) {
157:     /* communication pattern of SOLVEstruct is unlikely created for matmatsolve,
158:        thus destroy it and create a new SOLVEstruct.
159:        Otherwise it may result in memory corruption or incorrect solution
160:        See src/mat/examples/tests/ex125.c */
161: #if defined(PETSC_USE_COMPLEX)
162:     PetscStackCall("SuperLU_DIST:zSolveFinalize",zSolveFinalize(&lu->options, &lu->SOLVEstruct));
163: #else
164:     PetscStackCall("SuperLU_DIST:dSolveFinalize",dSolveFinalize(&lu->options, &lu->SOLVEstruct));
165: #endif
166:     lu->options.SolveInitialized = NO;
167:   }
168:   if (X != B_mpi) {
169:     MatCopy(B_mpi,X,SAME_NONZERO_PATTERN);
170:   }

172:   MatGetSize(B_mpi,NULL,&nrhs);

174:   PetscStackCall("SuperLU_DIST:PStatInit",PStatInit(&stat));        /* Initialize the statistics variables. */
175:   MatDenseGetArray(X,&bptr);

177: #if defined(PETSC_USE_COMPLEX)
178:   PetscStackCall("SuperLU_DIST:pzgssvx",pzgssvx(&lu->options,&lu->A_sup,&lu->ScalePermstruct,(doublecomplex*)bptr,m,nrhs,&lu->grid, &lu->LUstruct,&lu->SOLVEstruct,berr,&stat,&info));
179: #else
180:   PetscStackCall("SuperLU_DIST:pdgssvx",pdgssvx(&lu->options,&lu->A_sup,&lu->ScalePermstruct,bptr,m,nrhs,&lu->grid,&lu->LUstruct,&lu->SOLVEstruct,berr,&stat,&info));
181: #endif

183:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"pdgssvx fails, info: %d\n",info);
184:   MatDenseRestoreArray(X,&bptr);

186:   if (lu->options.PrintStat) PetscStackCall("SuperLU_DIST:PStatPrint",PStatPrint(&lu->options, &stat, &lu->grid));  /* Print the statistics. */
187:   PetscStackCall("SuperLU_DIST:PStatFree",PStatFree(&stat));
188:   lu->matsolve_iscalled    = PETSC_FALSE;
189:   lu->matmatsolve_iscalled = PETSC_TRUE;
190:   return(0);
191: }

193: /*
194:   input:
195:    F:        numeric Cholesky factor
196:   output:
197:    nneg:     total number of negative pivots
198:    nzero:    total number of zero pivots
199:    npos:     (global dimension of F) - nneg - nzero
200: */
201: static PetscErrorCode MatGetInertia_SuperLU_DIST(Mat F,PetscInt *nneg,PetscInt *nzero,PetscInt *npos)
202: {
203:   PetscErrorCode   ierr;
204:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)F->data;
205:   PetscScalar      *diagU=NULL;
206:   PetscInt         M,i,neg=0,zero=0,pos=0;
207:   PetscReal        r;

210:   if (!F->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix factor F is not assembled");
211:   if (lu->options.RowPerm != NOROWPERM) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must set NOROWPERM");
212:   MatGetSize(F,&M,NULL);
213:   PetscMalloc1(M,&diagU);
214:   MatSuperluDistGetDiagU(F,diagU);
215:   for (i=0; i<M; i++) {
216: #if defined(PETSC_USE_COMPLEX)
217:     r = PetscImaginaryPart(diagU[i])/10.0;
218:     if (r< -PETSC_MACHINE_EPSILON || r>PETSC_MACHINE_EPSILON) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"diagU[%d]=%g + i %g is non-real",i,PetscRealPart(diagU[i]),r*10.0);
219:     r = PetscRealPart(diagU[i]);
220: #else
221:     r = diagU[i];
222: #endif
223:     if (r > 0) {
224:       pos++;
225:     } else if (r < 0) {
226:       neg++;
227:     } else zero++;
228:   }

230:   PetscFree(diagU);
231:   if (nneg)  *nneg  = neg;
232:   if (nzero) *nzero = zero;
233:   if (npos)  *npos  = pos;
234:   return(0);
235: }

237: static PetscErrorCode MatLUFactorNumeric_SuperLU_DIST(Mat F,Mat A,const MatFactorInfo *info)
238: {
239:   Mat_SuperLU_DIST  *lu = (Mat_SuperLU_DIST*)F->data;
240:   Mat               Aloc;
241:   const PetscScalar *av;
242:   const PetscInt    *ai=NULL,*aj=NULL;
243:   PetscInt          nz,dummy;
244:   int               sinfo;   /* SuperLU_Dist info flag is always an int even with long long indices */
245:   SuperLUStat_t     stat;
246:   double            *berr=0;
247:   PetscBool         ismpiaij,isseqaij,flg;
248:   PetscErrorCode    ierr;

251:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
252:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
253:   if (ismpiaij) {
254:     MatMPIAIJGetLocalMat(A,MAT_INITIAL_MATRIX,&Aloc);
255:   } else if (isseqaij) {
256:     PetscObjectReference((PetscObject)A);
257:     Aloc = A;
258:   } else SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for type %s",((PetscObject)A)->type_name);

260:   MatGetRowIJ(Aloc,0,PETSC_FALSE,PETSC_FALSE,&dummy,&ai,&aj,&flg);
261:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GetRowIJ failed");
262:   MatSeqAIJGetArrayRead(Aloc,&av);
263:   nz   = ai[Aloc->rmap->n];

265:   /* Allocations for A_sup */
266:   if (lu->options.Fact == DOFACT) { /* first numeric factorization */
267: #if defined(PETSC_USE_COMPLEX)
268:     PetscStackCall("SuperLU_DIST:zallocateA_dist",zallocateA_dist(Aloc->rmap->n, nz, &lu->val, &lu->col, &lu->row));
269: #else
270:     PetscStackCall("SuperLU_DIST:dallocateA_dist",dallocateA_dist(Aloc->rmap->n, nz, &lu->val, &lu->col, &lu->row));
271: #endif
272:   } else { /* successive numeric factorization, sparsity pattern and perm_c are reused. */
273:     if (lu->FactPattern == SamePattern_SameRowPerm) {
274:       lu->options.Fact = SamePattern_SameRowPerm; /* matrix has similar numerical values */
275:     } else if (lu->FactPattern == SamePattern) {
276:       PetscStackCall("SuperLU_DIST:Destroy_LU",Destroy_LU(A->rmap->N, &lu->grid, &lu->LUstruct)); /* Deallocate L and U matrices. */
277:       lu->options.Fact = SamePattern;
278:     } else if (lu->FactPattern == DOFACT) {
279:       PetscStackCall("SuperLU_DIST:Destroy_CompRowLoc_Matrix_dist",Destroy_CompRowLoc_Matrix_dist(&lu->A_sup));
280:       PetscStackCall("SuperLU_DIST:Destroy_LU",Destroy_LU(A->rmap->N, &lu->grid, &lu->LUstruct));
281:       lu->options.Fact = DOFACT;

283: #if defined(PETSC_USE_COMPLEX)
284:       PetscStackCall("SuperLU_DIST:zallocateA_dist",zallocateA_dist(Aloc->rmap->n, nz, &lu->val, &lu->col, &lu->row));
285: #else
286:       PetscStackCall("SuperLU_DIST:dallocateA_dist",dallocateA_dist(Aloc->rmap->n, nz, &lu->val, &lu->col, &lu->row));
287: #endif
288:     } else {
289:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"options.Fact must be one of SamePattern SamePattern_SameRowPerm DOFACT");
290:     }
291:   }

293:   /* Copy AIJ matrix to superlu_dist matrix */
294:   PetscArraycpy(lu->row,ai,Aloc->rmap->n+1);
295:   PetscArraycpy(lu->col,aj,nz);
296:   PetscArraycpy(lu->val,av,nz);
297:   MatRestoreRowIJ(Aloc,0,PETSC_FALSE,PETSC_FALSE,&dummy,&ai,&aj,&flg);
298:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"RestoreRowIJ failed");
299:   MatSeqAIJRestoreArrayRead(Aloc,&av);
300:   MatDestroy(&Aloc);

302:   /* Create and setup A_sup */
303:   if (lu->options.Fact == DOFACT) {
304: #if defined(PETSC_USE_COMPLEX)
305:     PetscStackCall("SuperLU_DIST:zCreate_CompRowLoc_Matrix_dist",zCreate_CompRowLoc_Matrix_dist(&lu->A_sup, A->rmap->N, A->cmap->N, nz, A->rmap->n, A->rmap->rstart, lu->val, lu->col, lu->row, SLU_NR_loc, SLU_Z, SLU_GE));
306: #else
307:     PetscStackCall("SuperLU_DIST:dCreate_CompRowLoc_Matrix_dist",dCreate_CompRowLoc_Matrix_dist(&lu->A_sup, A->rmap->N, A->cmap->N, nz, A->rmap->n, A->rmap->rstart, lu->val, lu->col, lu->row, SLU_NR_loc, SLU_D, SLU_GE));
308: #endif
309:   }

311:   /* Factor the matrix. */
312:   PetscStackCall("SuperLU_DIST:PStatInit",PStatInit(&stat));   /* Initialize the statistics variables. */
313: #if defined(PETSC_USE_COMPLEX)
314:   PetscStackCall("SuperLU_DIST:pzgssvx",pzgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, A->rmap->n, 0, &lu->grid, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &sinfo));
315: #else
316:   PetscStackCall("SuperLU_DIST:pdgssvx",pdgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, A->rmap->n, 0, &lu->grid, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &sinfo));
317: #endif

319:   if (sinfo > 0) {
320:     if (A->erroriffailure) {
321:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot in row %D",sinfo);
322:     } else {
323:       if (sinfo <= lu->A_sup.ncol) {
324:         F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
325:         PetscInfo1(F,"U(i,i) is exactly zero, i= %D\n",sinfo);
326:       } else if (sinfo > lu->A_sup.ncol) {
327:         /*
328:          number of bytes allocated when memory allocation
329:          failure occurred, plus A->ncol.
330:          */
331:         F->factorerrortype = MAT_FACTOR_OUTMEMORY;
332:         PetscInfo1(F,"Number of bytes allocated when memory allocation fails %D\n",sinfo);
333:       }
334:     }
335:   } else if (sinfo < 0) {
336:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB, "info = %D, argument in p*gssvx() had an illegal value", sinfo);
337:   }

339:   if (lu->options.PrintStat) {
340:     PetscStackCall("SuperLU_DIST:PStatPrint",PStatPrint(&lu->options, &stat, &lu->grid));  /* Print the statistics. */
341:   }
342:   PetscStackCall("SuperLU_DIST:PStatFree",PStatFree(&stat));
343:   F->assembled     = PETSC_TRUE;
344:   F->preallocated  = PETSC_TRUE;
345:   lu->options.Fact = FACTORED; /* The factored form of A is supplied. Local option used by this func. only */
346:   return(0);
347: }

349: /* Note the Petsc r and c permutations are ignored */
350: static PetscErrorCode MatLUFactorSymbolic_SuperLU_DIST(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
351: {
352:   Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST*)F->data;
353:   PetscInt         M   = A->rmap->N,N=A->cmap->N;

356:   /* Initialize the SuperLU process grid. */
357:   PetscStackCall("SuperLU_DIST:superlu_gridinit",superlu_gridinit(lu->comm_superlu, lu->nprow, lu->npcol, &lu->grid));

359:   /* Initialize ScalePermstruct and LUstruct. */
360:   PetscStackCall("SuperLU_DIST:ScalePermstructInit",ScalePermstructInit(M, N, &lu->ScalePermstruct));
361:   PetscStackCall("SuperLU_DIST:LUstructInit",LUstructInit(N, &lu->LUstruct));
362:   F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU_DIST;
363:   F->ops->solve           = MatSolve_SuperLU_DIST;
364:   F->ops->matsolve        = MatMatSolve_SuperLU_DIST;
365:   F->ops->getinertia      = NULL;

367:   if (A->symmetric || A->hermitian) F->ops->getinertia = MatGetInertia_SuperLU_DIST;
368:   lu->CleanUpSuperLU_Dist = PETSC_TRUE;
369:   return(0);
370: }

372: static PetscErrorCode MatCholeskyFactorSymbolic_SuperLU_DIST(Mat F,Mat A,IS r,const MatFactorInfo *info)
373: {

377:   MatLUFactorSymbolic_SuperLU_DIST(F,A,r,r,info);
378:   F->ops->choleskyfactornumeric = MatLUFactorNumeric_SuperLU_DIST;
379:   return(0);
380: }

382: static PetscErrorCode MatFactorGetSolverType_aij_superlu_dist(Mat A,MatSolverType *type)
383: {
385:   *type = MATSOLVERSUPERLU_DIST;
386:   return(0);
387: }

389: static PetscErrorCode MatView_Info_SuperLU_DIST(Mat A,PetscViewer viewer)
390: {
391:   Mat_SuperLU_DIST       *lu=(Mat_SuperLU_DIST*)A->data;
392:   superlu_dist_options_t options;
393:   PetscErrorCode         ierr;

396:   /* check if matrix is superlu_dist type */
397:   if (A->ops->solve != MatSolve_SuperLU_DIST) return(0);

399:   options = lu->options;
400:   PetscViewerASCIIPrintf(viewer,"SuperLU_DIST run parameters:\n");
401:   PetscViewerASCIIPrintf(viewer,"  Process grid nprow %D x npcol %D \n",lu->nprow,lu->npcol);
402:   PetscViewerASCIIPrintf(viewer,"  Equilibrate matrix %s \n",PetscBools[options.Equil != NO]);
403:   PetscViewerASCIIPrintf(viewer,"  Replace tiny pivots %s \n",PetscBools[options.ReplaceTinyPivot != NO]);
404:   PetscViewerASCIIPrintf(viewer,"  Use iterative refinement %s \n",PetscBools[options.IterRefine == SLU_DOUBLE]);
405:   PetscViewerASCIIPrintf(viewer,"  Processors in row %d col partition %d \n",lu->nprow,lu->npcol);

407:   switch (options.RowPerm) {
408:   case NOROWPERM:
409:     PetscViewerASCIIPrintf(viewer,"  Row permutation NOROWPERM\n");
410:     break;
411:   case LargeDiag_MC64:
412:     PetscViewerASCIIPrintf(viewer,"  Row permutation LargeDiag_MC64\n");
413:     break;
414:   case LargeDiag_AWPM:
415:     PetscViewerASCIIPrintf(viewer,"  Row permutation LargeDiag_AWPM\n");
416:     break;
417:   case MY_PERMR:
418:     PetscViewerASCIIPrintf(viewer,"  Row permutation MY_PERMR\n");
419:     break;
420:   default:
421:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown column permutation");
422:   }

424:   switch (options.ColPerm) {
425:   case NATURAL:
426:     PetscViewerASCIIPrintf(viewer,"  Column permutation NATURAL\n");
427:     break;
428:   case MMD_AT_PLUS_A:
429:     PetscViewerASCIIPrintf(viewer,"  Column permutation MMD_AT_PLUS_A\n");
430:     break;
431:   case MMD_ATA:
432:     PetscViewerASCIIPrintf(viewer,"  Column permutation MMD_ATA\n");
433:     break;
434:   /*  Even though this is called METIS, the SuperLU_DIST code sets this by default if PARMETIS is defined, not METIS */
435:   case METIS_AT_PLUS_A:
436:     PetscViewerASCIIPrintf(viewer,"  Column permutation METIS_AT_PLUS_A\n");
437:     break;
438:   case PARMETIS:
439:     PetscViewerASCIIPrintf(viewer,"  Column permutation PARMETIS\n");
440:     break;
441:   default:
442:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown column permutation");
443:   }

445:   PetscViewerASCIIPrintf(viewer,"  Parallel symbolic factorization %s \n",PetscBools[options.ParSymbFact != NO]);

447:   if (lu->FactPattern == SamePattern) {
448:     PetscViewerASCIIPrintf(viewer,"  Repeated factorization SamePattern\n");
449:   } else if (lu->FactPattern == SamePattern_SameRowPerm) {
450:     PetscViewerASCIIPrintf(viewer,"  Repeated factorization SamePattern_SameRowPerm\n");
451:   } else if (lu->FactPattern == DOFACT) {
452:     PetscViewerASCIIPrintf(viewer,"  Repeated factorization DOFACT\n");
453:   } else {
454:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown factorization pattern");
455:   }
456:   return(0);
457: }

459: static PetscErrorCode MatView_SuperLU_DIST(Mat A,PetscViewer viewer)
460: {
461:   PetscErrorCode    ierr;
462:   PetscBool         iascii;
463:   PetscViewerFormat format;

466:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
467:   if (iascii) {
468:     PetscViewerGetFormat(viewer,&format);
469:     if (format == PETSC_VIEWER_ASCII_INFO) {
470:       MatView_Info_SuperLU_DIST(A,viewer);
471:     }
472:   }
473:   return(0);
474: }

476: static PetscErrorCode MatGetFactor_aij_superlu_dist(Mat A,MatFactorType ftype,Mat *F)
477: {
478:   Mat                    B;
479:   Mat_SuperLU_DIST       *lu;
480:   PetscErrorCode         ierr;
481:   PetscInt               M=A->rmap->N,N=A->cmap->N,indx;
482:   PetscMPIInt            size;
483:   superlu_dist_options_t options;
484:   PetscBool              flg;
485:   const char             *colperm[]     = {"NATURAL","MMD_AT_PLUS_A","MMD_ATA","METIS_AT_PLUS_A","PARMETIS"};
486:   const char             *rowperm[]     = {"NOROWPERM","LargeDiag_MC64","LargeDiag_AWPM","MY_PERMR"};
487:   const char             *factPattern[] = {"SamePattern","SamePattern_SameRowPerm","DOFACT"};
488:   PetscBool              set;

491:   /* Create the factorization matrix */
492:   MatCreate(PetscObjectComm((PetscObject)A),&B);
493:   MatSetSizes(B,A->rmap->n,A->cmap->n,M,N);
494:   PetscStrallocpy("superlu_dist",&((PetscObject)B)->type_name);
495:   MatSetUp(B);
496:   B->ops->getinfo = MatGetInfo_External;
497:   B->ops->view    = MatView_SuperLU_DIST;
498:   B->ops->destroy = MatDestroy_SuperLU_DIST;

500:   /* Set the default input options:
501:      options.Fact              = DOFACT;
502:      options.Equil             = YES;
503:      options.ParSymbFact       = NO;
504:      options.ColPerm           = METIS_AT_PLUS_A;
505:      options.RowPerm           = LargeDiag_MC64;
506:      options.ReplaceTinyPivot  = YES;
507:      options.IterRefine        = DOUBLE;
508:      options.Trans             = NOTRANS;
509:      options.SolveInitialized  = NO; -hold the communication pattern used MatSolve() and MatMatSolve()
510:      options.RefineInitialized = NO;
511:      options.PrintStat         = YES;
512:      options.SymPattern        = NO;
513:   */
514:   set_default_options_dist(&options);

516:   if (ftype == MAT_FACTOR_LU) {
517:     B->factortype = MAT_FACTOR_LU;
518:     B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU_DIST;
519:   } else {
520:     B->factortype = MAT_FACTOR_CHOLESKY;
521:     B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SuperLU_DIST;
522:     options.SymPattern = YES;
523:   }

525:   /* set solvertype */
526:   PetscFree(B->solvertype);
527:   PetscStrallocpy(MATSOLVERSUPERLU_DIST,&B->solvertype);

529:   PetscNewLog(B,&lu);
530:   B->data = lu;

532:   MPI_Comm_dup(PetscObjectComm((PetscObject)A),&(lu->comm_superlu));
533:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
534:   /* Default num of process columns and rows */
535:   lu->nprow = (int_t) (0.5 + PetscSqrtReal((PetscReal)size));
536:   if (!lu->nprow) lu->nprow = 1;
537:   while (lu->nprow > 0) {
538:     lu->npcol = (int_t) (size/lu->nprow);
539:     if (size == lu->nprow * lu->npcol) break;
540:     lu->nprow--;
541:   }

543:   PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"SuperLU_Dist Options","Mat");
544:   PetscOptionsInt("-mat_superlu_dist_r","Number rows in processor partition","None",lu->nprow,(PetscInt*)&lu->nprow,NULL);
545:   PetscOptionsInt("-mat_superlu_dist_c","Number columns in processor partition","None",lu->npcol,(PetscInt*)&lu->npcol,NULL);
546:   if (size != lu->nprow * lu->npcol) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number of processes %d must equal to nprow %d * npcol %d",size,lu->nprow,lu->npcol);

548:   PetscOptionsBool("-mat_superlu_dist_equil","Equilibrate matrix","None",options.Equil ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
549:   if (set && !flg) options.Equil = NO;

551:   PetscOptionsEList("-mat_superlu_dist_rowperm","Row permutation","None",rowperm,4,rowperm[1],&indx,&flg);
552:   if (flg) {
553:     switch (indx) {
554:     case 0:
555:       options.RowPerm = NOROWPERM;
556:       break;
557:     case 1:
558:       options.RowPerm = LargeDiag_MC64;
559:       break;
560:     case 2:
561:       options.RowPerm = LargeDiag_AWPM;
562:       break;
563:     case 3:
564:       options.RowPerm = MY_PERMR;
565:       break;
566:     default:
567:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown row permutation");
568:     }
569:   }

571:   PetscOptionsEList("-mat_superlu_dist_colperm","Column permutation","None",colperm,5,colperm[3],&indx,&flg);
572:   if (flg) {
573:     switch (indx) {
574:     case 0:
575:       options.ColPerm = NATURAL;
576:       break;
577:     case 1:
578:       options.ColPerm = MMD_AT_PLUS_A;
579:       break;
580:     case 2:
581:       options.ColPerm = MMD_ATA;
582:       break;
583:     case 3:
584:       options.ColPerm = METIS_AT_PLUS_A;
585:       break;
586:     case 4:
587:       options.ColPerm = PARMETIS;   /* only works for np>1 */
588:       break;
589:     default:
590:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown column permutation");
591:     }
592:   }

594:   options.ReplaceTinyPivot = NO;
595:   PetscOptionsBool("-mat_superlu_dist_replacetinypivot","Replace tiny pivots","None",options.ReplaceTinyPivot ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
596:   if (set && flg) options.ReplaceTinyPivot = YES;

598:   options.ParSymbFact = NO;
599:   PetscOptionsBool("-mat_superlu_dist_parsymbfact","Parallel symbolic factorization","None",PETSC_FALSE,&flg,&set);
600:   if (set && flg && size>1) {
601: #if defined(PETSC_HAVE_PARMETIS)
602:     options.ParSymbFact = YES;
603:     options.ColPerm     = PARMETIS;   /* in v2.2, PARMETIS is forced for ParSymbFact regardless of user ordering setting */
604: #else
605:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"parsymbfact needs PARMETIS");
606: #endif
607:   }

609:   lu->FactPattern = SamePattern;
610:   PetscOptionsEList("-mat_superlu_dist_fact","Sparsity pattern for repeated matrix factorization","None",factPattern,3,factPattern[0],&indx,&flg);
611:   if (flg) {
612:     switch (indx) {
613:     case 0:
614:       lu->FactPattern = SamePattern;
615:       break;
616:     case 1:
617:       lu->FactPattern = SamePattern_SameRowPerm;
618:       break;
619:     case 2:
620:       lu->FactPattern = DOFACT;
621:       break;
622:     }
623:   }

625:   options.IterRefine = NOREFINE;
626:   PetscOptionsBool("-mat_superlu_dist_iterrefine","Use iterative refinement","None",options.IterRefine == NOREFINE ? PETSC_FALSE : PETSC_TRUE ,&flg,&set);
627:   if (set) {
628:     if (flg) options.IterRefine = SLU_DOUBLE;
629:     else options.IterRefine = NOREFINE;
630:   }

632:   if (PetscLogPrintInfo) options.PrintStat = YES;
633:   else options.PrintStat = NO;
634:   PetscOptionsBool("-mat_superlu_dist_statprint","Print factorization information","None",(PetscBool)options.PrintStat,(PetscBool*)&options.PrintStat,NULL);
635:   PetscOptionsEnd();

637:   lu->options              = options;
638:   lu->options.Fact         = DOFACT;
639:   lu->matsolve_iscalled    = PETSC_FALSE;
640:   lu->matmatsolve_iscalled = PETSC_FALSE;

642:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_aij_superlu_dist);
643:   PetscObjectComposeFunction((PetscObject)B,"MatSuperluDistGetDiagU_C",MatSuperluDistGetDiagU_SuperLU_DIST);

645:   *F = B;
646:   return(0);
647: }

649: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_SuperLU_DIST(void)
650: {
653:   MatSolverTypeRegister(MATSOLVERSUPERLU_DIST,MATMPIAIJ,MAT_FACTOR_LU,MatGetFactor_aij_superlu_dist);
654:   MatSolverTypeRegister(MATSOLVERSUPERLU_DIST,MATSEQAIJ,MAT_FACTOR_LU,MatGetFactor_aij_superlu_dist);
655:   MatSolverTypeRegister(MATSOLVERSUPERLU_DIST,MATMPIAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_superlu_dist);
656:   MatSolverTypeRegister(MATSOLVERSUPERLU_DIST,MATSEQAIJ,MAT_FACTOR_CHOLESKY,MatGetFactor_aij_superlu_dist);
657:   return(0);
658: }

660: /*MC
661:   MATSOLVERSUPERLU_DIST - Parallel direct solver package for LU factorization

663:   Use ./configure --download-superlu_dist --download-parmetis --download-metis --download-ptscotch  to have PETSc installed with SuperLU_DIST

665:   Use -pc_type lu -pc_factor_mat_solver_type superlu_dist to use this direct solver

667:    Works with AIJ matrices

669:   Options Database Keys:
670: + -mat_superlu_dist_r <n> - number of rows in processor partition
671: . -mat_superlu_dist_c <n> - number of columns in processor partition
672: . -mat_superlu_dist_equil - equilibrate the matrix
673: . -mat_superlu_dist_rowperm <NOROWPERM,LargeDiag_MC64,LargeDiag_AWPM,MY_PERMR> - row permutation
674: . -mat_superlu_dist_colperm <MMD_AT_PLUS_A,MMD_ATA,NATURAL> - column permutation
675: . -mat_superlu_dist_replacetinypivot - replace tiny pivots
676: . -mat_superlu_dist_fact <SamePattern> - (choose one of) SamePattern SamePattern_SameRowPerm DOFACT
677: . -mat_superlu_dist_iterrefine - use iterative refinement
678: - -mat_superlu_dist_statprint - print factorization information

680:    Level: beginner

682: .seealso: PCLU

684: .seealso: PCFactorSetMatSolverType(), MatSolverType

686: M*/