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;
15: lu->nonzerosalongdiagonal = PETSC_TRUE;
16: if (z == PETSC_DECIDE) lu->nonzerosalongdiagonaltol = 1.e-10;
17: else lu->nonzerosalongdiagonaltol = z;
18: return(0);
19: }
21: static PetscErrorCode PCSetFromOptions_LU(PetscOptionItems *PetscOptionsObject,PC pc)
22: {
23: PC_LU *lu = (PC_LU*)pc->data;
25: PetscBool flg = PETSC_FALSE;
26: PetscReal tol;
29: PetscOptionsHead(PetscOptionsObject,"LU options");
30: PCSetFromOptions_Factor(PetscOptionsObject,pc);
32: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
33: if (flg) {
34: tol = PETSC_DECIDE;
35: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,NULL);
36: PCFactorReorderForNonzeroDiagonal(pc,tol);
37: }
38: PetscOptionsTail();
39: return(0);
40: }
42: static PetscErrorCode PCSetUp_LU(PC pc)
43: {
44: PetscErrorCode ierr;
45: PC_LU *dir = (PC_LU*)pc->data;
46: MatSolverType stype;
47: MatFactorError err;
50: pc->failedreason = PC_NOERROR;
51: if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill;
53: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
54: if (dir->hdr.inplace) {
55: MatFactorType ftype;
57: MatGetFactorType(pc->pmat, &ftype);
58: if (ftype == MAT_FACTOR_NONE) {
59: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
60: ISDestroy(&dir->col);
61: /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
62: PCFactorSetDefaultOrdering_Factor(pc);
63: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
64: if (dir->row) {
65: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
66: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
67: }
68: MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
69: MatFactorGetError(pc->pmat,&err);
70: if (err) { /* Factor() fails */
71: pc->failedreason = (PCFailedReason)err;
72: return(0);
73: }
74: }
75: ((PC_Factor*)dir)->fact = pc->pmat;
76: } else {
77: MatInfo info;
79: if (!pc->setupcalled) {
80: PetscBool canuseordering;
81: if (!((PC_Factor*)dir)->fact) {
82: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
83: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
84: }
85: MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
86: if (canuseordering) {
87: PCFactorSetDefaultOrdering_Factor(pc);
88: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
89: if (dir->nonzerosalongdiagonal) {
90: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
91: }
92: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
93: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
94: }
95: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
96: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
97: dir->hdr.actualfill = info.fill_ratio_needed;
98: } else if (pc->flag != SAME_NONZERO_PATTERN) {
99: PetscBool canuseordering;
100: if (!dir->hdr.reuseordering) {
101: MatDestroy(&((PC_Factor*)dir)->fact);
102: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
103: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
104: MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
105: if (canuseordering) {
106: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
107: ISDestroy(&dir->col);
108: PCFactorSetDefaultOrdering_Factor(pc);
109: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
110: if (dir->nonzerosalongdiagonal) {
111: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
112: }
113: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
114: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
115: }
116: }
117: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
118: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
119: dir->hdr.actualfill = info.fill_ratio_needed;
120: } else {
121: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
122: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
123: MatFactorClearError(((PC_Factor*)dir)->fact);
124: pc->failedreason = PC_NOERROR;
125: }
126: }
127: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
128: if (err) { /* FactorSymbolic() fails */
129: pc->failedreason = (PCFailedReason)err;
130: return(0);
131: }
133: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
134: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
135: if (err) { /* FactorNumeric() fails */
136: pc->failedreason = (PCFailedReason)err;
137: }
139: }
141: PCFactorGetMatSolverType(pc,&stype);
142: if (!stype) {
143: MatSolverType solverpackage;
144: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
145: PCFactorSetMatSolverType(pc,solverpackage);
146: }
147: return(0);
148: }
150: static PetscErrorCode PCReset_LU(PC pc)
151: {
152: PC_LU *dir = (PC_LU*)pc->data;
156: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
157: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
158: ISDestroy(&dir->col);
159: return(0);
160: }
162: static PetscErrorCode PCDestroy_LU(PC pc)
163: {
164: PC_LU *dir = (PC_LU*)pc->data;
168: PCReset_LU(pc);
169: PetscFree(((PC_Factor*)dir)->ordering);
170: PetscFree(((PC_Factor*)dir)->solvertype);
171: PetscFree(pc->data);
172: return(0);
173: }
175: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
176: {
177: PC_LU *dir = (PC_LU*)pc->data;
181: if (dir->hdr.inplace) {
182: MatSolve(pc->pmat,x,y);
183: } else {
184: MatSolve(((PC_Factor*)dir)->fact,x,y);
185: }
186: return(0);
187: }
189: static PetscErrorCode PCMatApply_LU(PC pc,Mat X,Mat Y)
190: {
191: PC_LU *dir = (PC_LU*)pc->data;
195: if (dir->hdr.inplace) {
196: MatMatSolve(pc->pmat,X,Y);
197: } else {
198: MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
199: }
200: return(0);
201: }
203: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
204: {
205: PC_LU *dir = (PC_LU*)pc->data;
209: if (dir->hdr.inplace) {
210: MatSolveTranspose(pc->pmat,x,y);
211: } else {
212: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
213: }
214: return(0);
215: }
217: /* -----------------------------------------------------------------------------------*/
219: /*MC
220: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
222: Options Database Keys:
223: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
224: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
225: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
226: . -pc_factor_fill <fill> - Sets fill amount
227: . -pc_factor_in_place - Activates in-place factorization
228: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
229: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
230: stability of factorization.
231: . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
232: . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
233: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
234: diagonal.
236: Notes:
237: Not all options work for all matrix formats
238: Run with -help to see additional options for particular matrix formats or factorization
239: algorithms
241: Level: beginner
243: Notes:
244: Usually this will compute an "exact" solution in one iteration and does
245: not need a Krylov method (i.e. you can use -ksp_type preonly, or
246: KSPSetType(ksp,KSPPREONLY) for the Krylov method
248: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
249: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
250: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
251: PCFactorSetPivotInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
252: PCFactorReorderForNonzeroDiagonal()
253: M*/
255: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
256: {
258: PC_LU *dir;
261: PetscNewLog(pc,&dir);
262: pc->data = (void*)dir;
263: PCFactorInitialize(pc,MAT_FACTOR_LU);
264: dir->nonzerosalongdiagonal = PETSC_FALSE;
266: ((PC_Factor*)dir)->info.fill = 5.0;
267: ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
268: ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
269: dir->col = NULL;
270: dir->row = NULL;
272: pc->ops->reset = PCReset_LU;
273: pc->ops->destroy = PCDestroy_LU;
274: pc->ops->apply = PCApply_LU;
275: pc->ops->matapply = PCMatApply_LU;
276: pc->ops->applytranspose = PCApplyTranspose_LU;
277: pc->ops->setup = PCSetUp_LU;
278: pc->ops->setfromoptions = PCSetFromOptions_LU;
279: pc->ops->view = PCView_Factor;
280: pc->ops->applyrichardson = NULL;
281: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
282: return(0);
283: }