Actual source code: lu.c
petsc-3.12.5 2020-03-29
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>
11: PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)
12: {
13: PC_LU *lu = (PC_LU*)pc->data;
16: lu->nonzerosalongdiagonal = PETSC_TRUE;
17: if (z == PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
18: else lu->nonzerosalongdiagonaltol = z;
19: return(0);
20: }
22: static PetscErrorCode PCSetFromOptions_LU(PetscOptionItems *PetscOptionsObject,PC pc)
23: {
24: PC_LU *lu = (PC_LU*)pc->data;
26: PetscBool flg = PETSC_FALSE;
27: PetscReal tol;
30: PetscOptionsHead(PetscOptionsObject,"LU options");
31: PCSetFromOptions_Factor(PetscOptionsObject,pc);
33: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
34: if (flg) {
35: tol = PETSC_DECIDE;
36: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
37: PCFactorReorderForNonzeroDiagonal(pc,tol);
38: }
39: PetscOptionsTail();
40: return(0);
41: }
43: static PetscErrorCode PCSetUp_LU(PC pc)
44: {
45: PetscErrorCode ierr;
46: PC_LU *dir = (PC_LU*)pc->data;
47: MatSolverType stype;
48: MatFactorError err;
51: pc->failedreason = PC_NOERROR;
52: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
54: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
55: if (dir->hdr.inplace) {
56: MatFactorType ftype;
58: MatGetFactorType(pc->pmat, &ftype);
59: if (ftype == MAT_FACTOR_NONE) {
60: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
61: ISDestroy(&dir->col);
62: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
63: if (dir->row) {
64: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
65: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
66: }
67: MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
68: MatFactorGetError(pc->pmat,&err);
69: if (err) { /* Factor() fails */
70: pc->failedreason = (PCFailedReason)err;
71: return(0);
72: }
73: }
74: ((PC_Factor*)dir)->fact = pc->pmat;
75: } else {
76: MatInfo info;
78: if (!pc->setupcalled) {
79: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
80: if (dir->nonzerosalongdiagonal) {
81: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
82: }
83: if (dir->row) {
84: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
85: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
86: }
87: if (!((PC_Factor*)dir)->fact) {
88: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
89: }
90: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
91: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
92: dir->hdr.actualfill = info.fill_ratio_needed;
93: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
94: } else if (pc->flag != SAME_NONZERO_PATTERN) {
95: if (!dir->hdr.reuseordering) {
96: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
97: ISDestroy(&dir->col);
98: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
99: if (dir->nonzerosalongdiagonal) {
100: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
101: }
102: if (dir->row) {
103: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
104: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
105: }
106: }
107: MatDestroy(&((PC_Factor*)dir)->fact);
108: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
109: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
110: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
111: dir->hdr.actualfill = info.fill_ratio_needed;
112: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
113: } else {
114: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
115: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
116: MatFactorClearError(((PC_Factor*)dir)->fact);
117: pc->failedreason = PC_NOERROR;
118: }
119: }
120: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
121: if (err) { /* FactorSymbolic() fails */
122: pc->failedreason = (PCFailedReason)err;
123: return(0);
124: }
126: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
127: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
128: if (err) { /* FactorNumeric() fails */
129: pc->failedreason = (PCFailedReason)err;
130: }
132: }
134: PCFactorGetMatSolverType(pc,&stype);
135: if (!stype) {
136: MatSolverType solverpackage;
137: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
138: PCFactorSetMatSolverType(pc,solverpackage);
139: }
140: return(0);
141: }
143: static PetscErrorCode PCReset_LU(PC pc)
144: {
145: PC_LU *dir = (PC_LU*)pc->data;
149: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
150: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
151: ISDestroy(&dir->col);
152: return(0);
153: }
155: static PetscErrorCode PCDestroy_LU(PC pc)
156: {
157: PC_LU *dir = (PC_LU*)pc->data;
161: PCReset_LU(pc);
162: PetscFree(((PC_Factor*)dir)->ordering);
163: PetscFree(((PC_Factor*)dir)->solvertype);
164: PetscFree(pc->data);
165: return(0);
166: }
168: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
169: {
170: PC_LU *dir = (PC_LU*)pc->data;
174: if (dir->hdr.inplace) {
175: MatSolve(pc->pmat,x,y);
176: } else {
177: MatSolve(((PC_Factor*)dir)->fact,x,y);
178: }
179: return(0);
180: }
182: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
183: {
184: PC_LU *dir = (PC_LU*)pc->data;
188: if (dir->hdr.inplace) {
189: MatSolveTranspose(pc->pmat,x,y);
190: } else {
191: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
192: }
193: return(0);
194: }
196: /* -----------------------------------------------------------------------------------*/
198: /*MC
199: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
201: Options Database Keys:
202: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
203: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
204: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
205: . -pc_factor_fill <fill> - Sets fill amount
206: . -pc_factor_in_place - Activates in-place factorization
207: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
208: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
209: stability of factorization.
210: . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
211: . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
212: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
213: diagonal.
215: Notes:
216: Not all options work for all matrix formats
217: Run with -help to see additional options for particular matrix formats or factorization
218: algorithms
220: Level: beginner
222: Notes:
223: Usually this will compute an "exact" solution in one iteration and does
224: not need a Krylov method (i.e. you can use -ksp_type preonly, or
225: KSPSetType(ksp,KSPPREONLY) for the Krylov method
227: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
228: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
229: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
230: PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
231: PCFactorReorderForNonzeroDiagonal()
232: M*/
234: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
235: {
237: PetscMPIInt size;
238: PC_LU *dir;
241: PetscNewLog(pc,&dir);
242: pc->data = (void*)dir;
243: PCFactorInitialize(pc);
244: ((PC_Factor*)dir)->factortype = MAT_FACTOR_LU;
245: dir->nonzerosalongdiagonal = PETSC_FALSE;
247: ((PC_Factor*)dir)->info.fill = 5.0;
248: ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
249: ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
250: dir->col = 0;
251: dir->row = 0;
253: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
254: if (size == 1) {
255: PetscStrallocpy(MATORDERINGND,(char**)&((PC_Factor*)dir)->ordering);
256: } else {
257: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
258: }
260: pc->ops->reset = PCReset_LU;
261: pc->ops->destroy = PCDestroy_LU;
262: pc->ops->apply = PCApply_LU;
263: pc->ops->applytranspose = PCApplyTranspose_LU;
264: pc->ops->setup = PCSetUp_LU;
265: pc->ops->setfromoptions = PCSetFromOptions_LU;
266: pc->ops->view = PCView_Factor;
267: pc->ops->applyrichardson = 0;
268: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
269: return(0);
270: }