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>
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,NULL);
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: /* This should only get the ordering if needed, but since MatGetFactor() is not called we can't know if it is needed */
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 useordering;
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: MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
86: if (useordering) {
87: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
88: if (dir->nonzerosalongdiagonal) {
89: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
90: }
91: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
92: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
93: }
94: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
95: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
96: dir->hdr.actualfill = info.fill_ratio_needed;
97: } else if (pc->flag != SAME_NONZERO_PATTERN) {
98: PetscBool useordering;
99: if (!dir->hdr.reuseordering) {
100: MatDestroy(&((PC_Factor*)dir)->fact);
101: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
102: PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
103: MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
104: if (useordering) {
105: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
106: ISDestroy(&dir->col);
107: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
108: if (dir->nonzerosalongdiagonal) {
109: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
110: }
111: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
112: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
113: }
114: }
115: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
116: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
117: dir->hdr.actualfill = info.fill_ratio_needed;
118: } else {
119: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
120: if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
121: MatFactorClearError(((PC_Factor*)dir)->fact);
122: pc->failedreason = PC_NOERROR;
123: }
124: }
125: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
126: if (err) { /* FactorSymbolic() fails */
127: pc->failedreason = (PCFailedReason)err;
128: return(0);
129: }
131: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
132: MatFactorGetError(((PC_Factor*)dir)->fact,&err);
133: if (err) { /* FactorNumeric() fails */
134: pc->failedreason = (PCFailedReason)err;
135: }
137: }
139: PCFactorGetMatSolverType(pc,&stype);
140: if (!stype) {
141: MatSolverType solverpackage;
142: MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
143: PCFactorSetMatSolverType(pc,solverpackage);
144: }
145: return(0);
146: }
148: static PetscErrorCode PCReset_LU(PC pc)
149: {
150: PC_LU *dir = (PC_LU*)pc->data;
154: if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
155: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
156: ISDestroy(&dir->col);
157: return(0);
158: }
160: static PetscErrorCode PCDestroy_LU(PC pc)
161: {
162: PC_LU *dir = (PC_LU*)pc->data;
166: PCReset_LU(pc);
167: PetscFree(((PC_Factor*)dir)->ordering);
168: PetscFree(((PC_Factor*)dir)->solvertype);
169: PetscFree(pc->data);
170: return(0);
171: }
173: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
174: {
175: PC_LU *dir = (PC_LU*)pc->data;
179: if (dir->hdr.inplace) {
180: MatSolve(pc->pmat,x,y);
181: } else {
182: MatSolve(((PC_Factor*)dir)->fact,x,y);
183: }
184: return(0);
185: }
187: static PetscErrorCode PCMatApply_LU(PC pc,Mat X,Mat Y)
188: {
189: PC_LU *dir = (PC_LU*)pc->data;
193: if (dir->hdr.inplace) {
194: MatMatSolve(pc->pmat,X,Y);
195: } else {
196: MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
197: }
198: return(0);
199: }
201: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
202: {
203: PC_LU *dir = (PC_LU*)pc->data;
207: if (dir->hdr.inplace) {
208: MatSolveTranspose(pc->pmat,x,y);
209: } else {
210: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
211: }
212: return(0);
213: }
215: /* -----------------------------------------------------------------------------------*/
217: /*MC
218: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
220: Options Database Keys:
221: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
222: . -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
223: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
224: . -pc_factor_fill <fill> - Sets fill amount
225: . -pc_factor_in_place - Activates in-place factorization
226: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
227: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
228: stability of factorization.
229: . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
230: . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
231: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
232: diagonal.
234: Notes:
235: Not all options work for all matrix formats
236: Run with -help to see additional options for particular matrix formats or factorization
237: algorithms
239: Level: beginner
241: Notes:
242: Usually this will compute an "exact" solution in one iteration and does
243: not need a Krylov method (i.e. you can use -ksp_type preonly, or
244: KSPSetType(ksp,KSPPREONLY) for the Krylov method
246: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
247: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
248: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
249: PCFactorSetPivotInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
250: PCFactorReorderForNonzeroDiagonal()
251: M*/
253: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
254: {
256: PetscMPIInt size;
257: PC_LU *dir;
260: PetscNewLog(pc,&dir);
261: pc->data = (void*)dir;
262: PCFactorInitialize(pc);
263: ((PC_Factor*)dir)->factortype = 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: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
273: if (size == 1) {
274: PetscStrallocpy(MATORDERINGND,(char**)&((PC_Factor*)dir)->ordering);
275: } else {
276: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
277: }
279: pc->ops->reset = PCReset_LU;
280: pc->ops->destroy = PCDestroy_LU;
281: pc->ops->apply = PCApply_LU;
282: pc->ops->matapply = PCMatApply_LU;
283: pc->ops->applytranspose = PCApplyTranspose_LU;
284: pc->ops->setup = PCSetUp_LU;
285: pc->ops->setfromoptions = PCSetFromOptions_LU;
286: pc->ops->view = PCView_Factor;
287: pc->ops->applyrichardson = NULL;
288: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
289: return(0);
290: }