Actual source code: lu.c
petsc-3.3-p7 2013-05-11
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> /*I "petscpc.h" I*/
10: EXTERN_C_BEGIN
13: PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)
14: {
15: PC_LU *lu = (PC_LU*)pc->data;
18: lu->nonzerosalongdiagonal = PETSC_TRUE;
19: if (z == PETSC_DECIDE) {
20: lu->nonzerosalongdiagonaltol = 1.e-10;
21: } else {
22: lu->nonzerosalongdiagonaltol = z;
23: }
24: return(0);
25: }
26: EXTERN_C_END
28: EXTERN_C_BEGIN
31: PetscErrorCode PCFactorSetReuseOrdering_LU(PC pc,PetscBool flag)
32: {
33: PC_LU *lu = (PC_LU*)pc->data;
36: lu->reuseordering = flag;
37: return(0);
38: }
39: EXTERN_C_END
41: EXTERN_C_BEGIN
44: PetscErrorCode PCFactorSetReuseFill_LU(PC pc,PetscBool flag)
45: {
46: PC_LU *lu = (PC_LU*)pc->data;
49: lu->reusefill = flag;
50: return(0);
51: }
52: EXTERN_C_END
56: static PetscErrorCode PCSetFromOptions_LU(PC pc)
57: {
58: PC_LU *lu = (PC_LU*)pc->data;
59: PetscErrorCode ierr;
60: PetscBool flg = PETSC_FALSE;
61: PetscReal tol;
64: PetscOptionsHead("LU options");
65: PCSetFromOptions_Factor(pc);
67: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
68: if (flg) {
69: tol = PETSC_DECIDE;
70: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
71: PCFactorReorderForNonzeroDiagonal(pc,tol);
72: }
73: PetscOptionsTail();
74: return(0);
75: }
79: static PetscErrorCode PCView_LU(PC pc,PetscViewer viewer)
80: {
81: PC_LU *lu = (PC_LU*)pc->data;
83: PetscBool iascii;
86: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
87: if (iascii) {
88: if (lu->inplace) {
89: PetscViewerASCIIPrintf(viewer," LU: in-place factorization\n");
90: } else {
91: PetscViewerASCIIPrintf(viewer," LU: out-of-place factorization\n");
92: }
93:
94: if (lu->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
95: if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
96: }
97: PCView_Factor(pc,viewer);
98: return(0);
99: }
103: static PetscErrorCode PCSetUp_LU(PC pc)
104: {
106: PC_LU *dir = (PC_LU*)pc->data;
109: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
111: if (dir->inplace) {
112: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
113: ISDestroy(&dir->col);
114: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
115: if (dir->row) {
116: PetscLogObjectParent(pc,dir->row);
117: PetscLogObjectParent(pc,dir->col);
118: }
119: MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
120: ((PC_Factor*)dir)->fact = pc->pmat;
121: } else {
122: MatInfo info;
123: if (!pc->setupcalled) {
124: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
125: if (dir->nonzerosalongdiagonal) {
126: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
127: }
128: if (dir->row) {
129: PetscLogObjectParent(pc,dir->row);
130: PetscLogObjectParent(pc,dir->col);
131: }
132: if (!((PC_Factor*)dir)->fact){
133: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
134: }
135: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
136: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
137: dir->actualfill = info.fill_ratio_needed;
138: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
139: } else if (pc->flag != SAME_NONZERO_PATTERN) {
140: if (!dir->reuseordering) {
141: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
142: ISDestroy(&dir->col);
143: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
144: if (dir->nonzerosalongdiagonal) {
145: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
146: }
147: if (dir->row) {
148: PetscLogObjectParent(pc,dir->row);
149: PetscLogObjectParent(pc,dir->col);
150: }
151: }
152: MatDestroy(&((PC_Factor*)dir)->fact);
153: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
154: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
155: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
156: dir->actualfill = info.fill_ratio_needed;
157: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
158: }
159: if((!pc->setupcalled) || (pc->setupcalled && !pc->reuse)) {
160: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
161: }
162: }
163: return(0);
164: }
168: static PetscErrorCode PCReset_LU(PC pc)
169: {
170: PC_LU *dir = (PC_LU*)pc->data;
174: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
175: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
176: ISDestroy(&dir->col);
177: return(0);
178: }
182: static PetscErrorCode PCDestroy_LU(PC pc)
183: {
184: PC_LU *dir = (PC_LU*)pc->data;
188: PCReset_LU(pc);
189: PetscFree(((PC_Factor*)dir)->ordering);
190: PetscFree(((PC_Factor*)dir)->solvertype);
191: PetscFree(pc->data);
192: return(0);
193: }
197: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
198: {
199: PC_LU *dir = (PC_LU*)pc->data;
203: if (dir->inplace) {MatSolve(pc->pmat,x,y);}
204: else {MatSolve(((PC_Factor*)dir)->fact,x,y);}
205: return(0);
206: }
210: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
211: {
212: PC_LU *dir = (PC_LU*)pc->data;
216: if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
217: else {MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);}
218: return(0);
219: }
221: /* -----------------------------------------------------------------------------------*/
223: EXTERN_C_BEGIN
226: PetscErrorCode PCFactorSetUseInPlace_LU(PC pc)
227: {
228: PC_LU *dir = (PC_LU*)pc->data;
231: dir->inplace = PETSC_TRUE;
232: return(0);
233: }
234: EXTERN_C_END
236: /* ------------------------------------------------------------------------ */
238: /*MC
239: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
241: Options Database Keys:
242: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
243: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles
244: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
245: . -pc_factor_fill <fill> - Sets fill amount
246: . -pc_factor_in_place - Activates in-place factorization
247: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
248: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
249: stability of factorization.
250: . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
251: . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
252: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
253: diagonal.
255: Notes: Not all options work for all matrix formats
256: Run with -help to see additional options for particular matrix formats or factorization
257: algorithms
259: Level: beginner
261: Concepts: LU factorization, direct solver
263: Notes: Usually this will compute an "exact" solution in one iteration and does
264: not need a Krylov method (i.e. you can use -ksp_type preonly, or
265: KSPSetType(ksp,KSPPREONLY) for the Krylov method
267: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
268: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
269: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
270: PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
271: PCFactorReorderForNonzeroDiagonal()
272: M*/
274: EXTERN_C_BEGIN
277: PetscErrorCode PCCreate_LU(PC pc)
278: {
280: PetscMPIInt size;
281: PC_LU *dir;
284: PetscNewLog(pc,PC_LU,&dir);
286: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
287: ((PC_Factor*)dir)->fact = PETSC_NULL;
288: ((PC_Factor*)dir)->factortype = MAT_FACTOR_LU;
289: dir->inplace = PETSC_FALSE;
290: dir->nonzerosalongdiagonal = PETSC_FALSE;
292: ((PC_Factor*)dir)->info.fill = 5.0;
293: ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
294: ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
295: ((PC_Factor*)dir)->info.shiftamount = 0.0;
296: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
297: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
298: dir->col = 0;
299: dir->row = 0;
301: PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype); /* default SolverPackage */
302: MPI_Comm_size(((PetscObject)pc)->comm,&size);
303: if (size == 1) {
304: PetscStrallocpy(MATORDERINGND,&((PC_Factor*)dir)->ordering);
305: } else {
306: PetscStrallocpy(MATORDERINGNATURAL,&((PC_Factor*)dir)->ordering);
307: }
308: dir->reusefill = PETSC_FALSE;
309: dir->reuseordering = PETSC_FALSE;
310: pc->data = (void*)dir;
312: pc->ops->reset = PCReset_LU;
313: pc->ops->destroy = PCDestroy_LU;
314: pc->ops->apply = PCApply_LU;
315: pc->ops->applytranspose = PCApplyTranspose_LU;
316: pc->ops->setup = PCSetUp_LU;
317: pc->ops->setfromoptions = PCSetFromOptions_LU;
318: pc->ops->view = PCView_LU;
319: pc->ops->applyrichardson = 0;
320: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
322: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C","PCFactorSetUpMatSolverPackage_Factor",
323: PCFactorSetUpMatSolverPackage_Factor);
324: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
325: PCFactorGetMatSolverPackage_Factor);
326: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
327: PCFactorSetMatSolverPackage_Factor);
328: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
329: PCFactorSetZeroPivot_Factor);
330: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftType_C","PCFactorSetShiftType_Factor",
331: PCFactorSetShiftType_Factor);
332: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftAmount_C","PCFactorSetShiftAmount_Factor",
333: PCFactorSetShiftAmount_Factor);
334: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
335: PCFactorSetFill_Factor);
336: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_LU",
337: PCFactorSetUseInPlace_LU);
338: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
339: PCFactorSetMatOrderingType_Factor);
340: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_LU",
341: PCFactorSetReuseOrdering_LU);
342: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_LU",
343: PCFactorSetReuseFill_LU);
344: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetColumnPivot_C","PCFactorSetColumnPivot_Factor",
345: PCFactorSetColumnPivot_Factor);
346: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",
347: PCFactorSetPivotInBlocks_Factor);
348: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_LU",
349: PCFactorReorderForNonzeroDiagonal_LU);
350: return(0);
351: }
352: EXTERN_C_END