Actual source code: lu.c

petsc-dev 2014-08-28
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(PC pc)
 48: {
 49:   PC_LU          *lu = (PC_LU*)pc->data;
 51:   PetscBool      flg = PETSC_FALSE;
 52:   PetscReal      tol;

 55:   PetscOptionsHead("LU options");
 56:   PCSetFromOptions_Factor(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:     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:     if ((!pc->setupcalled) || pc->setupcalled) {
152:       MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
153:     }
154:   }
155:   return(0);
156: }

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

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

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

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

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

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

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

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

219: /* -----------------------------------------------------------------------------------*/

223: PetscErrorCode  PCFactorSetUseInPlace_LU(PC pc)
224: {
225:   PC_LU *dir = (PC_LU*)pc->data;

228:   dir->inplace = PETSC_TRUE;
229:   return(0);
230: }

232: /* ------------------------------------------------------------------------ */

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

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

251:    Notes: Not all options work for all matrix formats
252:           Run with -help to see additional options for particular matrix formats or factorization
253:           algorithms

255:    Level: beginner

257:    Concepts: LU factorization, direct solver

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

263: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
264:            PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
265:            PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
266:            PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
267:            PCFactorReorderForNonzeroDiagonal()
268: M*/

272: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
273: {
275:   PetscMPIInt    size;
276:   PC_LU          *dir;

279:   PetscNewLog(pc,&dir);

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

283:   ((PC_Factor*)dir)->fact       = NULL;
284:   ((PC_Factor*)dir)->factortype = MAT_FACTOR_LU;
285:   dir->inplace                  = PETSC_FALSE;
286:   dir->nonzerosalongdiagonal    = PETSC_FALSE;

288:   ((PC_Factor*)dir)->info.fill          = 5.0;
289:   ((PC_Factor*)dir)->info.dtcol         = 1.e-6;  /* default to pivoting; this is only thing PETSc LU supports */
290:   ((PC_Factor*)dir)->info.shifttype     = (PetscReal)MAT_SHIFT_NONE;
291:   ((PC_Factor*)dir)->info.shiftamount   = 0.0;
292:   ((PC_Factor*)dir)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
293:   ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
294:   dir->col                              = 0;
295:   dir->row                              = 0;

297:   PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype); /* default SolverPackage */
298:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
299:   if (size == 1) {
300:     PetscStrallocpy(MATORDERINGND,(char**)&((PC_Factor*)dir)->ordering);
301:   } else {
302:     PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
303:   }
304:   dir->reusefill     = PETSC_FALSE;
305:   dir->reuseordering = PETSC_FALSE;
306:   pc->data           = (void*)dir;

308:   pc->ops->reset             = PCReset_LU;
309:   pc->ops->destroy           = PCDestroy_LU;
310:   pc->ops->apply             = PCApply_LU;
311:   pc->ops->applytranspose    = PCApplyTranspose_LU;
312:   pc->ops->setup             = PCSetUp_LU;
313:   pc->ops->setfromoptions    = PCSetFromOptions_LU;
314:   pc->ops->view              = PCView_LU;
315:   pc->ops->applyrichardson   = 0;
316:   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;

318:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
319:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
320:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
321:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
322:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
323:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
324:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
325:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_LU);
326:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
327:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_LU);
328:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_LU);
329:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetColumnPivot_C",PCFactorSetColumnPivot_Factor);
330:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);
331:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
332:   return(0);
333: }