Actual source code: lu.c

petsc-master 2015-09-03
Report Typos and Errors
  2: /*
  3:    Defines a direct factorization preconditioner for any Mat implementation
  4:    Note: this need not be consided a preconditioner since it supplies
  5:          a direct solver.
  6: */

  8: #include <../src/ksp/pc/impls/factor/lu/lu.h>  /*I "petscpc.h" I*/

 12: PetscErrorCode  PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)
 13: {
 14:   PC_LU *lu = (PC_LU*)pc->data;

 17:   lu->nonzerosalongdiagonal = PETSC_TRUE;
 18:   if (z == PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
 19:   else lu->nonzerosalongdiagonaltol = z;
 20:   return(0);
 21: }

 25: PetscErrorCode  PCFactorSetReuseOrdering_LU(PC pc,PetscBool flag)
 26: {
 27:   PC_LU *lu = (PC_LU*)pc->data;

 30:   lu->reuseordering = flag;
 31:   return(0);
 32: }

 36: PetscErrorCode  PCFactorSetReuseFill_LU(PC pc,PetscBool flag)
 37: {
 38:   PC_LU *lu = (PC_LU*)pc->data;

 41:   lu->reusefill = flag;
 42:   return(0);
 43: }

 47: static PetscErrorCode PCSetFromOptions_LU(PetscOptions *PetscOptionsObject,PC pc)
 48: {
 49:   PC_LU          *lu = (PC_LU*)pc->data;
 51:   PetscBool      flg = PETSC_FALSE;
 52:   PetscReal      tol;

 55:   PetscOptionsHead(PetscOptionsObject,"LU options");
 56:   PCSetFromOptions_Factor(PetscOptionsObject,pc);

 58:   PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
 59:   if (flg) {
 60:     tol  = PETSC_DECIDE;
 61:     PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
 62:     PCFactorReorderForNonzeroDiagonal(pc,tol);
 63:   }
 64:   PetscOptionsTail();
 65:   return(0);
 66: }

 70: static PetscErrorCode PCView_LU(PC pc,PetscViewer viewer)
 71: {
 72:   PC_LU          *lu = (PC_LU*)pc->data;
 74:   PetscBool      iascii;

 77:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 78:   if (iascii) {
 79:     if (lu->inplace) {
 80:       PetscViewerASCIIPrintf(viewer,"  LU: in-place factorization\n");
 81:     } else {
 82:       PetscViewerASCIIPrintf(viewer,"  LU: out-of-place factorization\n");
 83:     }

 85:     if (lu->reusefill)    {PetscViewerASCIIPrintf(viewer,"       Reusing fill from past factorization\n");}
 86:     if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer,"       Reusing reordering from past factorization\n");}
 87:   }
 88:   PCView_Factor(pc,viewer);
 89:   return(0);
 90: }

 94: static PetscErrorCode PCSetUp_LU(PC pc)
 95: {
 97:   PC_LU          *dir = (PC_LU*)pc->data;

100:   if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;

102:   if (dir->inplace) {
103:     if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
104:     ISDestroy(&dir->col);
105:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
106:     if (dir->row) {
107:       PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
108:       PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
109:     }
110:     MatSetErrorIfFPE(pc->pmat,pc->erroriffailure);
111:     MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
112:     ((PC_Factor*)dir)->fact = pc->pmat;
113:   } else {
114:     MatInfo info;
115:     if (!pc->setupcalled) {
116:       MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
117:       if (dir->nonzerosalongdiagonal) {
118:         MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
119:       }
120:       if (dir->row) {
121:         PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
122:         PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
123:       }
124:       if (!((PC_Factor*)dir)->fact) {
125:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
126:       }
127:       MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
128:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
129:       dir->actualfill = info.fill_ratio_needed;
130:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
131:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
132:       if (!dir->reuseordering) {
133:         if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
134:         ISDestroy(&dir->col);
135:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
136:         if (dir->nonzerosalongdiagonal) {
137:           MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
138:         }
139:         if (dir->row) {
140:           PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
141:           PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
142:         }
143:       }
144:       MatDestroy(&((PC_Factor*)dir)->fact);
145:       MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
146:       MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
147:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
148:       dir->actualfill = info.fill_ratio_needed;
149:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
150:     }
151:     MatSetErrorIfFPE(pc->pmat,pc->erroriffailure);
152:     MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
153:   }
154:   return(0);
155: }

159: static PetscErrorCode PCReset_LU(PC pc)
160: {
161:   PC_LU          *dir = (PC_LU*)pc->data;

165:   if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
166:   if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
167:   ISDestroy(&dir->col);
168:   return(0);
169: }

173: static PetscErrorCode PCDestroy_LU(PC pc)
174: {
175:   PC_LU          *dir = (PC_LU*)pc->data;

179:   PCReset_LU(pc);
180:   PetscFree(((PC_Factor*)dir)->ordering);
181:   PetscFree(((PC_Factor*)dir)->solvertype);
182:   PetscFree(pc->data);
183:   return(0);
184: }

188: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
189: {
190:   PC_LU          *dir = (PC_LU*)pc->data;

194:   if (dir->inplace) {
195:     MatSolve(pc->pmat,x,y);
196:   } else {
197:     MatSolve(((PC_Factor*)dir)->fact,x,y);
198:   }
199:   return(0);
200: }

204: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
205: {
206:   PC_LU          *dir = (PC_LU*)pc->data;

210:   if (dir->inplace) {
211:     MatSolveTranspose(pc->pmat,x,y);
212:   } else {
213:     MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
214:   }
215:   return(0);
216: }

218: /* -----------------------------------------------------------------------------------*/

222: PetscErrorCode  PCFactorSetUseInPlace_LU(PC pc,PetscBool flg)
223: {
224:   PC_LU *dir = (PC_LU*)pc->data;

227:   dir->inplace = flg;
228:   return(0);
229: }

233: PetscErrorCode  PCFactorGetUseInPlace_LU(PC pc,PetscBool *flg)
234: {
235:   PC_LU *dir = (PC_LU*)pc->data;

238:   *flg = dir->inplace;
239:   return(0);
240: }

242: /* ------------------------------------------------------------------------ */

244: /*MC
245:    PCLU - Uses a direct solver, based on LU factorization, as a preconditioner

247:    Options Database Keys:
248: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
249: .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
250: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
251: .  -pc_factor_fill <fill> - Sets fill amount
252: .  -pc_factor_in_place - Activates in-place factorization
253: .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
254: .  -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
255:                                          stability of factorization.
256: .  -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
257: .  -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
258: -   -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
259:         diagonal.

261:    Notes: Not all options work for all matrix formats
262:           Run with -help to see additional options for particular matrix formats or factorization
263:           algorithms

265:    Level: beginner

267:    Concepts: LU factorization, direct solver

269:    Notes: Usually this will compute an "exact" solution in one iteration and does
270:           not need a Krylov method (i.e. you can use -ksp_type preonly, or
271:           KSPSetType(ksp,KSPPREONLY) for the Krylov method

273: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
274:            PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
275:            PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
276:            PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
277:            PCFactorReorderForNonzeroDiagonal()
278: M*/

282: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
283: {
285:   PetscMPIInt    size;
286:   PC_LU          *dir;

289:   PetscNewLog(pc,&dir);

291:   MatFactorInfoInitialize(&((PC_Factor*)dir)->info);

293:   ((PC_Factor*)dir)->fact       = NULL;
294:   ((PC_Factor*)dir)->factortype = MAT_FACTOR_LU;
295:   dir->inplace                  = PETSC_FALSE;
296:   dir->nonzerosalongdiagonal    = PETSC_FALSE;

298:   ((PC_Factor*)dir)->info.fill          = 5.0;
299:   ((PC_Factor*)dir)->info.dtcol         = 1.e-6;  /* default to pivoting; this is only thing PETSc LU supports */
300:   ((PC_Factor*)dir)->info.shifttype     = (PetscReal)MAT_SHIFT_NONE;
301:   ((PC_Factor*)dir)->info.shiftamount   = 0.0;
302:   ((PC_Factor*)dir)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
303:   ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
304:   dir->col                              = 0;
305:   dir->row                              = 0;

307:   PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype); /* default SolverPackage */
308:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
309:   if (size == 1) {
310:     PetscStrallocpy(MATORDERINGND,(char**)&((PC_Factor*)dir)->ordering);
311:   } else {
312:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
313:   }
314:   dir->reusefill     = PETSC_FALSE;
315:   dir->reuseordering = PETSC_FALSE;
316:   pc->data           = (void*)dir;

318:   pc->ops->reset             = PCReset_LU;
319:   pc->ops->destroy           = PCDestroy_LU;
320:   pc->ops->apply             = PCApply_LU;
321:   pc->ops->applytranspose    = PCApplyTranspose_LU;
322:   pc->ops->setup             = PCSetUp_LU;
323:   pc->ops->setfromoptions    = PCSetFromOptions_LU;
324:   pc->ops->view              = PCView_LU;
325:   pc->ops->applyrichardson   = 0;
326:   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;

328:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
329:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
330:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
331:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
332:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
333:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
334:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
335:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_LU);
336:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_LU);
337:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
338:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_LU);
339:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_LU);
340:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetColumnPivot_C",PCFactorSetColumnPivot_Factor);
341:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);
342:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
343:   return(0);
344: }