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:   PetscFunctionBegin;
 15:   lu->nonzerosalongdiagonal = PETSC_TRUE;
 16:   if (z == (PetscReal)PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
 17:   else lu->nonzerosalongdiagonaltol = z;
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

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

 27:   PetscFunctionBegin;
 28:   PetscOptionsHeadBegin(PetscOptionsObject, "LU options");
 29:   PetscCall(PCSetFromOptions_Factor(pc, PetscOptionsObject));

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

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

 48:   PetscFunctionBegin;
 49:   pc->failedreason = PC_NOERROR;
 50:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

 52:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
 53:   PetscCall(MatSetOptionsPrefixFactor(pc->pmat, prefix));

 55:   PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure));
 56:   if (dir->hdr.inplace) {
 57:     MatFactorType ftype;

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

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

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

127:   PetscCall(PCFactorGetMatSolverType(pc, &stype));
128:   if (!stype) {
129:     MatSolverType solverpackage;
130:     PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage));
131:     PetscCall(PCFactorSetMatSolverType(pc, solverpackage));
132:   }
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: static PetscErrorCode PCReset_LU(PC pc)
137: {
138:   PC_LU *dir = (PC_LU *)pc->data;

140:   PetscFunctionBegin;
141:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact));
142:   if (dir->row && dir->col && dir->row != dir->col) PetscCall(ISDestroy(&dir->row));
143:   PetscCall(ISDestroy(&dir->col));
144:   PetscFunctionReturn(PETSC_SUCCESS);
145: }

147: static PetscErrorCode PCDestroy_LU(PC pc)
148: {
149:   PC_LU *dir = (PC_LU *)pc->data;

151:   PetscFunctionBegin;
152:   PetscCall(PCReset_LU(pc));
153:   PetscCall(PetscFree(((PC_Factor *)dir)->ordering));
154:   PetscCall(PetscFree(((PC_Factor *)dir)->solvertype));
155:   PetscCall(PCFactorClearComposedFunctions(pc));
156:   PetscCall(PetscFree(pc->data));
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: static PetscErrorCode PCApply_LU(PC pc, Vec x, Vec y)
161: {
162:   PC_LU *dir = (PC_LU *)pc->data;

164:   PetscFunctionBegin;
165:   if (dir->hdr.inplace) {
166:     PetscCall(MatSolve(pc->pmat, x, y));
167:   } else {
168:     PetscCall(MatSolve(((PC_Factor *)dir)->fact, x, y));
169:   }
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: static PetscErrorCode PCMatApply_LU(PC pc, Mat X, Mat Y)
174: {
175:   PC_LU *dir = (PC_LU *)pc->data;

177:   PetscFunctionBegin;
178:   if (dir->hdr.inplace) {
179:     PetscCall(MatMatSolve(pc->pmat, X, Y));
180:   } else {
181:     PetscCall(MatMatSolve(((PC_Factor *)dir)->fact, X, Y));
182:   }
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

186: static PetscErrorCode PCApplyTranspose_LU(PC pc, Vec x, Vec y)
187: {
188:   PC_LU *dir = (PC_LU *)pc->data;

190:   PetscFunctionBegin;
191:   if (dir->hdr.inplace) {
192:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
193:   } else {
194:     PetscCall(MatSolveTranspose(((PC_Factor *)dir)->fact, x, y));
195:   }
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

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

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

217:    Level: beginner

219:    Notes:
220:    Not all options work for all matrix formats

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

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

228: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MatSolverType`, `MatGetFactor()`, `PCQR`, `PCSVD`,
229:           `PCILU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
230:           `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
231:           `PCFactorSetPivotInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
232:           `PCFactorReorderForNonzeroDiagonal()`
233: M*/

235: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
236: {
237:   PC_LU *dir;

239:   PetscFunctionBegin;
240:   PetscCall(PetscNew(&dir));
241:   pc->data = (void *)dir;
242:   PetscCall(PCFactorInitialize(pc, MAT_FACTOR_LU));
243:   dir->nonzerosalongdiagonal = PETSC_FALSE;

245:   ((PC_Factor *)dir)->info.fill      = 5.0;
246:   ((PC_Factor *)dir)->info.dtcol     = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
247:   ((PC_Factor *)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
248:   dir->col                           = NULL;
249:   dir->row                           = NULL;

251:   pc->ops->reset           = PCReset_LU;
252:   pc->ops->destroy         = PCDestroy_LU;
253:   pc->ops->apply           = PCApply_LU;
254:   pc->ops->matapply        = PCMatApply_LU;
255:   pc->ops->applytranspose  = PCApplyTranspose_LU;
256:   pc->ops->setup           = PCSetUp_LU;
257:   pc->ops->setfromoptions  = PCSetFromOptions_LU;
258:   pc->ops->view            = PCView_Factor;
259:   pc->ops->applyrichardson = NULL;
260:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFactorReorderForNonzeroDiagonal_C", PCFactorReorderForNonzeroDiagonal_LU));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }