Actual source code: lu.c
1: #define PETSCKSP_DLL
3: /*
4: Defines a direct factorization preconditioner for any Mat implementation
5: Note: this need not be consided a preconditioner since it supplies
6: a direct solver.
7: */
9: #include ../src/ksp/pc/impls/factor/lu/lu.h
14: PetscErrorCode PCFactorReorderForNonzeroDiagonal_LU(PC pc,PetscReal z)
15: {
16: PC_LU *lu = (PC_LU*)pc->data;
19: lu->nonzerosalongdiagonal = PETSC_TRUE;
20: if (z == PETSC_DECIDE) {
21: lu->nonzerosalongdiagonaltol = 1.e-10;
22: } else {
23: lu->nonzerosalongdiagonaltol = z;
24: }
25: return(0);
26: }
32: PetscErrorCode PCFactorSetReuseOrdering_LU(PC pc,PetscTruth flag)
33: {
34: PC_LU *lu = (PC_LU*)pc->data;
37: lu->reuseordering = flag;
38: return(0);
39: }
45: PetscErrorCode PCFactorSetReuseFill_LU(PC pc,PetscTruth flag)
46: {
47: PC_LU *lu = (PC_LU*)pc->data;
50: lu->reusefill = flag;
51: return(0);
52: }
57: static PetscErrorCode PCSetFromOptions_LU(PC pc)
58: {
59: PC_LU *lu = (PC_LU*)pc->data;
60: PetscErrorCode ierr;
61: PetscTruth flg,set;
62: char tname[256], solvertype[64];
63: PetscFList ordlist;
64: PetscReal tol;
67: if (!MatOrderingRegisterAllCalled) {MatOrderingRegisterAll(PETSC_NULL);}
68: PetscOptionsHead("LU options");
69: PetscOptionsName("-pc_factor_in_place","Form LU in the same memory as the matrix","PCFactorSetUseInPlace",&flg);
70: if (flg) {
71: PCFactorSetUseInPlace(pc);
72: }
73: PetscOptionsReal("-pc_factor_fill","Expected non-zeros in LU/non-zeros in matrix","PCFactorSetFill",((PC_Factor*)lu)->info.fill,&((PC_Factor*)lu)->info.fill,0);
75: PetscOptionsName("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",&flg);
76: if (flg) {
77: PCFactorSetShiftNonzero(pc,(PetscReal) PETSC_DECIDE);
78: }
79: PetscOptionsReal("-pc_factor_shift_nonzero","Shift added to diagonal","PCFactorSetShiftNonzero",((PC_Factor*)lu)->info.shiftnz,&((PC_Factor*)lu)->info.shiftnz,0);
80: PetscOptionsName("-pc_factor_shift_positive_definite","Manteuffel shift applied to diagonal","PCFactorSetShiftPd",&flg);
81: if (flg) {
82: PCFactorSetShiftPd(pc,PETSC_TRUE);
83: }
84: PetscOptionsReal("-pc_factor_zeropivot","Pivot is considered zero if less than","PCFactorSetZeroPivot",((PC_Factor*)lu)->info.zeropivot,&((PC_Factor*)lu)->info.zeropivot,0);
86: PetscOptionsName("-pc_factor_reuse_fill","Use fill from previous factorization","PCFactorSetReuseFill",&flg);
87: if (flg) {
88: PCFactorSetReuseFill(pc,PETSC_TRUE);
89: }
90: PetscOptionsName("-pc_factor_reuse_ordering","Reuse ordering from previous factorization","PCFactorSetReuseOrdering",&flg);
91: if (flg) {
92: PCFactorSetReuseOrdering(pc,PETSC_TRUE);
93: }
95: MatGetOrderingList(&ordlist);
96: PetscOptionsList("-pc_factor_mat_ordering_type","Reordering to reduce nonzeros in LU","PCFactorSetMatOrderingType",ordlist,((PC_Factor*)lu)->ordering,tname,256,&flg);
97: if (flg) {
98: PCFactorSetMatOrderingType(pc,tname);
99: }
101: /* maybe should have MatGetSolverTypes(Mat,&list) like the ordering list */
102: PetscOptionsString("-pc_factor_mat_solver_package","Specific LU solver to use","MatGetFactor",((PC_Factor*)lu)->solvertype,solvertype,64,&flg);
103: if (flg) {
104: PCFactorSetMatSolverPackage(pc,solvertype);
105: }
107: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
108: if (flg) {
109: tol = PETSC_DECIDE;
110: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
111: PCFactorReorderForNonzeroDiagonal(pc,tol);
112: }
114: PetscOptionsReal("-pc_factor_pivoting","Pivoting tolerance (used only for some factorization)","PCFactorSetPivoting",((PC_Factor*)lu)->info.dtcol,&((PC_Factor*)lu)->info.dtcol,&flg);
116: flg = ((PC_Factor*)lu)->info.pivotinblocks ? PETSC_TRUE : PETSC_FALSE;
117: PetscOptionsTruth("-pc_factor_pivot_in_blocks","Pivot inside matrix blocks for BAIJ and SBAIJ","PCFactorSetPivotInBlocks",flg,&flg,&set);
118: if (set) {
119: PCFactorSetPivotInBlocks(pc,flg);
120: }
121: PetscOptionsTail();
122: return(0);
123: }
127: static PetscErrorCode PCView_LU(PC pc,PetscViewer viewer)
128: {
129: PC_LU *lu = (PC_LU*)pc->data;
131: PetscTruth iascii,isstring;
134: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
135: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
136: if (iascii) {
138: if (lu->inplace) {PetscViewerASCIIPrintf(viewer," LU: in-place factorization\n");}
139: else {PetscViewerASCIIPrintf(viewer," LU: out-of-place factorization\n");}
140: PetscViewerASCIIPrintf(viewer," matrix ordering: %s\n",((PC_Factor*)lu)->ordering);
141: PetscViewerASCIIPrintf(viewer," LU: tolerance for zero pivot %G\n",((PC_Factor*)lu)->info.zeropivot);
142: if (((PC_Factor*)lu)->info.shiftpd) {PetscViewerASCIIPrintf(viewer," LU: using Manteuffel shift\n");}
143: if (((PC_Factor*)lu)->fact) {
144: PetscViewerASCIIPrintf(viewer," LU: factor fill ratio needed %G\n",lu->actualfill);
145: PetscViewerASCIIPrintf(viewer," Factored matrix follows\n");
146: PetscViewerASCIIPushTab(viewer);
147: PetscViewerASCIIPushTab(viewer);
148: PetscViewerASCIIPushTab(viewer);
149: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
150: MatView(((PC_Factor*)lu)->fact,viewer);
151: PetscViewerPopFormat(viewer);
152: PetscViewerASCIIPopTab(viewer);
153: PetscViewerASCIIPopTab(viewer);
154: PetscViewerASCIIPopTab(viewer);
155: }
156: if (lu->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
157: if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
158: } else if (isstring) {
159: PetscViewerStringSPrintf(viewer," order=%s",((PC_Factor*)lu)->ordering);
160: } else {
161: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCLU",((PetscObject)viewer)->type_name);
162: }
163: return(0);
164: }
168: static PetscErrorCode PCSetUp_LU(PC pc)
169: {
171: PC_LU *dir = (PC_LU*)pc->data;
174: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
176: if (dir->inplace) {
177: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
178: if (dir->col) {ISDestroy(dir->col);}
179: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
180: if (dir->row) {
181: PetscLogObjectParent(pc,dir->row);
182: PetscLogObjectParent(pc,dir->col);
183: }
184: MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
185: ((PC_Factor*)dir)->fact = pc->pmat;
186: } else {
187: MatInfo info;
188: if (!pc->setupcalled) {
189: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
190: if (dir->nonzerosalongdiagonal) {
191: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
192: }
193: if (dir->row) {
194: PetscLogObjectParent(pc,dir->row);
195: PetscLogObjectParent(pc,dir->col);
196: }
197: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
198: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
199: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
200: dir->actualfill = info.fill_ratio_needed;
201: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
202: } else if (pc->flag != SAME_NONZERO_PATTERN) {
203: if (!dir->reuseordering) {
204: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
205: if (dir->col) {ISDestroy(dir->col);}
206: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
207: if (dir->nonzerosalongdiagonal) {
208: MatReorderForNonzeroDiagonal(pc->pmat,dir->nonzerosalongdiagonaltol,dir->row,dir->col);
209: }
210: if (dir->row) {
211: PetscLogObjectParent(pc,dir->row);
212: PetscLogObjectParent(pc,dir->col);
213: }
214: }
215: MatDestroy(((PC_Factor*)dir)->fact);
216: MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_LU,&((PC_Factor*)dir)->fact);
217: MatLUFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
218: MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
219: dir->actualfill = info.fill_ratio_needed;
220: PetscLogObjectParent(pc,((PC_Factor*)dir)->fact);
221: }
222: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
223: }
224: return(0);
225: }
229: static PetscErrorCode PCDestroy_LU(PC pc)
230: {
231: PC_LU *dir = (PC_LU*)pc->data;
235: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(((PC_Factor*)dir)->fact);}
236: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(dir->row);}
237: if (dir->col) {ISDestroy(dir->col);}
238: PetscStrfree(((PC_Factor*)dir)->ordering);
239: PetscStrfree(((PC_Factor*)dir)->solvertype);
240: PetscFree(dir);
241: return(0);
242: }
246: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
247: {
248: PC_LU *dir = (PC_LU*)pc->data;
252: if (dir->inplace) {MatSolve(pc->pmat,x,y);}
253: else {MatSolve(((PC_Factor*)dir)->fact,x,y);}
254: return(0);
255: }
259: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
260: {
261: PC_LU *dir = (PC_LU*)pc->data;
265: if (dir->inplace) {MatSolveTranspose(pc->pmat,x,y);}
266: else {MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);}
267: return(0);
268: }
270: /* -----------------------------------------------------------------------------------*/
275: PetscErrorCode PCFactorSetUseInPlace_LU(PC pc)
276: {
277: PC_LU *dir = (PC_LU*)pc->data;
280: dir->inplace = PETSC_TRUE;
281: return(0);
282: }
285: /* ------------------------------------------------------------------------ */
287: /*MC
288: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
290: Options Database Keys:
291: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
292: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like spooles
293: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
294: . -pc_factor_fill <fill> - Sets fill amount
295: . -pc_factor_in_place - Activates in-place factorization
296: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
297: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
298: stability of factorization.
299: . -pc_factor_shift_nonzero <shift> - Sets shift amount or PETSC_DECIDE for the default
300: . -pc_factor_shift_positive_definite [PETSC_TRUE/PETSC_FALSE] - Activate/Deactivate PCFactorSetShiftPd(); the value
301: is optional with PETSC_TRUE being the default
302: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
303: diagonal.
305: Notes: Not all options work for all matrix formats
306: Run with -help to see additional options for particular matrix formats or factorization
307: algorithms
309: Level: beginner
311: Concepts: LU factorization, direct solver
313: Notes: Usually this will compute an "exact" solution in one iteration and does
314: not need a Krylov method (i.e. you can use -ksp_type preonly, or
315: KSPSetType(ksp,KSPPREONLY) for the Krylov method
317: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
318: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
319: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetPivoting(),
320: PCFactorSetPivotingInBlocks(),PCFactorSetShiftNonzero(),PCFactorSetShiftPd(), PCFactorReorderForNonzeroDiagonal()
321: M*/
326: PetscErrorCode PCCreate_LU(PC pc)
327: {
329: PetscMPIInt size;
330: PC_LU *dir;
333: PetscNewLog(pc,PC_LU,&dir);
335: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
336: ((PC_Factor*)dir)->fact = 0;
337: dir->inplace = PETSC_FALSE;
338: dir->nonzerosalongdiagonal = PETSC_FALSE;
340: ((PC_Factor*)dir)->info.fill = 5.0;
341: ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
342: ((PC_Factor*)dir)->info.shiftnz = 0.0;
343: ((PC_Factor*)dir)->info.zeropivot = 1.e-12;
344: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
345: ((PC_Factor*)dir)->info.shiftpd = 0.0; /* false */
346: dir->col = 0;
347: dir->row = 0;
349: PetscStrallocpy(MAT_SOLVER_PETSC,&((PC_Factor*)dir)->solvertype);
350: MPI_Comm_size(((PetscObject)pc)->comm,&size);
351: if (size == 1) {
352: PetscStrallocpy(MATORDERING_ND,&((PC_Factor*)dir)->ordering);
353: } else {
354: PetscStrallocpy(MATORDERING_NATURAL,&((PC_Factor*)dir)->ordering);
355: }
356: dir->reusefill = PETSC_FALSE;
357: dir->reuseordering = PETSC_FALSE;
358: pc->data = (void*)dir;
360: pc->ops->destroy = PCDestroy_LU;
361: pc->ops->apply = PCApply_LU;
362: pc->ops->applytranspose = PCApplyTranspose_LU;
363: pc->ops->setup = PCSetUp_LU;
364: pc->ops->setfromoptions = PCSetFromOptions_LU;
365: pc->ops->view = PCView_LU;
366: pc->ops->applyrichardson = 0;
367: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
369: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorGetMatSolverPackage_C","PCFactorGetMatSolverPackage_Factor",
370: PCFactorGetMatSolverPackage_Factor);
371: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatSolverPackage_C","PCFactorSetMatSolverPackage_Factor",
372: PCFactorSetMatSolverPackage_Factor);
373: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetZeroPivot_C","PCFactorSetZeroPivot_Factor",
374: PCFactorSetZeroPivot_Factor);
375: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftNonzero_C","PCFactorSetShiftNonzero_Factor",
376: PCFactorSetShiftNonzero_Factor);
377: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetShiftPd_C","PCFactorSetShiftPd_Factor",
378: PCFactorSetShiftPd_Factor);
380: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetFill_C","PCFactorSetFill_Factor",
381: PCFactorSetFill_Factor);
382: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetUseInPlace_C","PCFactorSetUseInPlace_LU",
383: PCFactorSetUseInPlace_LU);
384: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetMatOrderingType_C","PCFactorSetMatOrderingType_Factor",
385: PCFactorSetMatOrderingType_Factor);
386: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseOrdering_C","PCFactorSetReuseOrdering_LU",
387: PCFactorSetReuseOrdering_LU);
388: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetReuseFill_C","PCFactorSetReuseFill_LU",
389: PCFactorSetReuseFill_LU);
390: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivoting_C","PCFactorSetPivoting_Factor",
391: PCFactorSetPivoting_Factor);
392: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorSetPivotInBlocks_C","PCFactorSetPivotInBlocks_Factor",
393: PCFactorSetPivotInBlocks_Factor);
394: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C","PCFactorReorderForNonzeroDiagonal_LU",
395: PCFactorReorderForNonzeroDiagonal_LU);
396: return(0);
397: }