Actual source code: lu.c


  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>

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

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

 20: static PetscErrorCode PCSetFromOptions_LU(PC pc, PetscOptionItems *PetscOptionsObject)
 21: {
 22:   PC_LU    *lu  = (PC_LU *)pc->data;
 23:   PetscBool flg = PETSC_FALSE;
 24:   PetscReal tol;

 26:   PetscOptionsHeadBegin(PetscOptionsObject, "LU options");
 27:   PCSetFromOptions_Factor(pc, PetscOptionsObject);

 29:   PetscOptionsName("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", &flg);
 30:   if (flg) {
 31:     tol = PETSC_DECIDE;
 32:     PetscOptionsReal("-pc_factor_nonzeros_along_diagonal", "Reorder to remove zeros from diagonal", "PCFactorReorderForNonzeroDiagonal", lu->nonzerosalongdiagonaltol, &tol, NULL);
 33:     PCFactorReorderForNonzeroDiagonal(pc, tol);
 34:   }
 35:   PetscOptionsHeadEnd();
 36:   return 0;
 37: }

 39: static PetscErrorCode PCSetUp_LU(PC pc)
 40: {
 41:   PC_LU         *dir = (PC_LU *)pc->data;
 42:   MatSolverType  stype;
 43:   MatFactorError err;
 44:   const char    *prefix;

 46:   pc->failedreason = PC_NOERROR;
 47:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

 49:   PCGetOptionsPrefix(pc, &prefix);
 50:   MatSetOptionsPrefixFactor(pc->pmat, prefix);

 52:   MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
 53:   if (dir->hdr.inplace) {
 54:     MatFactorType ftype;

 56:     MatGetFactorType(pc->pmat, &ftype);
 57:     if (ftype == MAT_FACTOR_NONE) {
 58:       if (dir->row && dir->col && dir->row != dir->col) ISDestroy(&dir->row);
 59:       ISDestroy(&dir->col);
 60:       /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
 61:       PCFactorSetDefaultOrdering_Factor(pc);
 62:       MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 63:       if (dir->row) { }
 64:       MatLUFactor(pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info);
 65:       MatFactorGetError(pc->pmat, &err);
 66:       if (err) { /* Factor() fails */
 67:         pc->failedreason = (PCFailedReason)err;
 68:         return 0;
 69:       }
 70:     }
 71:     ((PC_Factor *)dir)->fact = pc->pmat;
 72:   } else {
 73:     MatInfo info;

 75:     if (!pc->setupcalled) {
 76:       PetscBool canuseordering;
 77:       if (!((PC_Factor *)dir)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_LU, &((PC_Factor *)dir)->fact); }
 78:       MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
 79:       if (canuseordering) {
 80:         PCFactorSetDefaultOrdering_Factor(pc);
 81:         MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 82:         if (dir->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col);
 83:       }
 84:       MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info);
 85:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
 86:       dir->hdr.actualfill = info.fill_ratio_needed;
 87:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 88:       PetscBool canuseordering;
 89:       if (!dir->hdr.reuseordering) {
 90:         MatDestroy(&((PC_Factor *)dir)->fact);
 91:         MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_LU, &((PC_Factor *)dir)->fact);
 92:         MatFactorGetCanUseOrdering(((PC_Factor *)dir)->fact, &canuseordering);
 93:         if (canuseordering) {
 94:           if (dir->row && dir->col && dir->row != dir->col) ISDestroy(&dir->row);
 95:           ISDestroy(&dir->col);
 96:           PCFactorSetDefaultOrdering_Factor(pc);
 97:           MatGetOrdering(pc->pmat, ((PC_Factor *)dir)->ordering, &dir->row, &dir->col);
 98:           if (dir->nonzerosalongdiagonal) MatReorderForNonzeroDiagonal(pc->pmat, dir->nonzerosalongdiagonaltol, dir->row, dir->col);
 99:         }
100:       }
101:       MatLUFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->row, dir->col, &((PC_Factor *)dir)->info);
102:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
103:       dir->hdr.actualfill = info.fill_ratio_needed;
104:     } else {
105:       MatFactorGetError(((PC_Factor *)dir)->fact, &err);
106:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
107:         MatFactorClearError(((PC_Factor *)dir)->fact);
108:         pc->failedreason = PC_NOERROR;
109:       }
110:     }
111:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
112:     if (err) { /* FactorSymbolic() fails */
113:       pc->failedreason = (PCFailedReason)err;
114:       return 0;
115:     }

117:     MatLUFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info);
118:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
119:     if (err) { /* FactorNumeric() fails */
120:       pc->failedreason = (PCFailedReason)err;
121:     }
122:   }

124:   PCFactorGetMatSolverType(pc, &stype);
125:   if (!stype) {
126:     MatSolverType solverpackage;
127:     MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage);
128:     PCFactorSetMatSolverType(pc, solverpackage);
129:   }
130:   return 0;
131: }

133: static PetscErrorCode PCReset_LU(PC pc)
134: {
135:   PC_LU *dir = (PC_LU *)pc->data;

137:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) MatDestroy(&((PC_Factor *)dir)->fact);
138:   if (dir->row && dir->col && dir->row != dir->col) ISDestroy(&dir->row);
139:   ISDestroy(&dir->col);
140:   return 0;
141: }

143: static PetscErrorCode PCDestroy_LU(PC pc)
144: {
145:   PC_LU *dir = (PC_LU *)pc->data;

147:   PCReset_LU(pc);
148:   PetscFree(((PC_Factor *)dir)->ordering);
149:   PetscFree(((PC_Factor *)dir)->solvertype);
150:   PCFactorClearComposedFunctions(pc);
151:   PetscFree(pc->data);
152:   return 0;
153: }

155: static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y)
156: {
157:   PC_LU *dir = (PC_LU *)pc->data;

159:   if (dir->hdr.inplace) {
160:     MatSolve(pc->pmat, x, y);
161:   } else {
162:     MatSolve(((PC_Factor *)dir)->fact, x, y);
163:   }
164:   return 0;
165: }

167: static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y)
168: {
169:   PC_LU *dir = (PC_LU *)pc->data;

171:   if (dir->hdr.inplace) {
172:     MatMatSolve(pc->pmat, X, Y);
173:   } else {
174:     MatMatSolve(((PC_Factor *)dir)->fact, X, Y);
175:   }
176:   return 0;
177: }

179: static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y)
180: {
181:   PC_LU *dir = (PC_LU *)pc->data;

183:   if (dir->hdr.inplace) {
184:     MatSolveTranspose(pc->pmat, x, y);
185:   } else {
186:     MatSolveTranspose(((PC_Factor *)dir)->fact, x, y);
187:   }
188:   return 0;
189: }

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

194:    Options Database Keys:
195: +  -pc_factor_reuse_ordering - Activate `PCFactorSetReuseOrdering()`
196: .  -pc_factor_mat_solver_type - Actives `PCFactorSetMatSolverType()` to choose the direct solver, like superlu
197: .  -pc_factor_reuse_fill - Activates `PCFactorSetReuseFill()`
198: .  -pc_factor_fill <fill> - Sets fill amount
199: .  -pc_factor_in_place - Activates in-place factorization
200: .  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
201: .  -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
202:                                          stability of factorization.
203: .  -pc_factor_shift_type <shifttype> - Sets shift type or -1 for the default; use '-help' for a list of available types
204: .  -pc_factor_shift_amount <shiftamount> - Sets shift amount or -1 for the default
205: .  -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the diagonal.
206: .  -pc_factor_mat_solver_type <packagename> - use an external package for the solve, see `MatSolverType` for possibilities
207: -  -mat_solvertype_optionname - options for a specific solver package, for example -mat_mumps_cntl_1

209:    Level: beginner

211:    Notes:
212:    Not all options work for all matrix formats

214:    Run with -help to see additional options for particular matrix formats or factorization algorithms

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

220: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MatSolverType`, `MatGetFactor()`, `PCQR`, `PCSVD`,
221:           `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
222:           `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
223:           `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
224:           `PCFactorReorderForNonzeroDiagonal()`
225: M*/

227: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
228: {
229:   PC_LU *dir;

231:   PetscNew(&dir);
232:   pc->data = (void *)dir;
233:   PCFactorInitialize(pc, MAT_FACTOR_LU);
234:   dir->nonzerosalongdiagonal = PETSC_FALSE;

236:   ((PC_Factor *)dir)->info.fill      = 5.0;
237:   ((PC_Factor *)dir)->info.dtcol     = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
238:   ((PC_Factor *)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
239:   dir->col                           = NULL;
240:   dir->row                           = NULL;

242:   pc->ops->reset           = PCReset_LU;
243:   pc->ops->destroy         = PCDestroy_LU;
244:   pc->ops->apply           = PCApply_LU;
245:   pc->ops->matapply        = PCMatApply_LU;
246:   pc->ops->applytranspose  = PCApplyTranspose_LU;
247:   pc->ops->setup           = PCSetUp_LU;
248:   pc->ops->setfromoptions  = PCSetFromOptions_LU;
249:   pc->ops->view            = PCView_Factor;
250:   pc->ops->applyrichardson = NULL;
251:   PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_LU);
252:   return 0;
253: }