Actual source code: lu.c
petsc-3.7.3 2016-08-01
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*/
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) lu->nonzerosalongdiagonaltol = 1.e-10;
20: else lu->nonzerosalongdiagonaltol = z;
21: return(0);
22: }
26: PetscErrorCode PCFactorSetReuseOrdering_LU(PC pc,PetscBool flag)
27: {
28: PC_LU *lu = (PC_LU*)pc->data;
31: lu->reuseordering = flag;
32: return(0);
33: }
37: PetscErrorCode PCFactorSetReuseFill_LU(PC pc,PetscBool flag)
38: {
39: PC_LU *lu = (PC_LU*)pc->data;
42: lu->reusefill = flag;
43: return(0);
44: }
48: static PetscErrorCode PCSetFromOptions_LU(PetscOptionItems *PetscOptionsObject,PC pc)
49: {
50: PC_LU *lu = (PC_LU*)pc->data;
52: PetscBool flg = PETSC_FALSE;
53: PetscReal tol;
56: PetscOptionsHead(PetscOptionsObject,"LU options");
57: PCSetFromOptions_Factor(PetscOptionsObject,pc);
59: PetscOptionsName("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",&flg);
60: if (flg) {
61: tol = PETSC_DECIDE;
62: PetscOptionsReal("-pc_factor_nonzeros_along_diagonal","Reorder to remove zeros from diagonal","PCFactorReorderForNonzeroDiagonal",lu->nonzerosalongdiagonaltol,&tol,0);
63: PCFactorReorderForNonzeroDiagonal(pc,tol);
64: }
65: PetscOptionsTail();
66: return(0);
67: }
71: static PetscErrorCode PCView_LU(PC pc,PetscViewer viewer)
72: {
73: PC_LU *lu = (PC_LU*)pc->data;
75: PetscBool iascii;
78: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
79: if (iascii) {
80: if (lu->inplace) {
81: PetscViewerASCIIPrintf(viewer," LU: in-place factorization\n");
82: } else {
83: PetscViewerASCIIPrintf(viewer," LU: out-of-place factorization\n");
84: }
86: if (lu->reusefill) {PetscViewerASCIIPrintf(viewer," Reusing fill from past factorization\n");}
87: if (lu->reuseordering) {PetscViewerASCIIPrintf(viewer," Reusing reordering from past factorization\n");}
88: }
89: PCView_Factor(pc,viewer);
90: return(0);
91: }
95: static PetscErrorCode PCSetUp_LU(PC pc)
96: {
98: PC_LU *dir = (PC_LU*)pc->data;
99: const MatSolverPackage stype;
102: if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;
104: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
105: if (dir->inplace) {
106: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
107: ISDestroy(&dir->col);
108: MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
109: if (dir->row) {
110: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
111: PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->col);
112: }
113: MatLUFactor(pc->pmat,dir->row,dir->col,&((PC_Factor*)dir)->info);
114: if (pc->pmat->errortype) { /* Factor() fails */
115: pc->failedreason = (PCFailedReason)pc->pmat->errortype;
116: return(0);
117: }
119: ((PC_Factor*)dir)->fact = pc->pmat;
120: } else {
121: MatInfo info;
122: Mat F;
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((PetscObject)pc,(PetscObject)dir->row);
130: PetscLogObjectParent((PetscObject)pc,(PetscObject)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((PetscObject)pc,(PetscObject)((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((PetscObject)pc,(PetscObject)dir->row);
149: PetscLogObjectParent((PetscObject)pc,(PetscObject)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((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
158: }
159: F = ((PC_Factor*)dir)->fact;
160: if (F->errortype) { /* FactorSymbolic() fails */
161: pc->failedreason = (PCFailedReason)F->errortype;
162: return(0);
163: }
165: MatLUFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
166: if (F->errortype) { /* FactorNumeric() fails */
167: pc->failedreason = (PCFailedReason)F->errortype;
168: }
170: }
172: PCFactorGetMatSolverPackage(pc,&stype);
173: if (!stype) {
174: PCFactorSetMatSolverPackage(pc,((PC_Factor*)dir)->fact->solvertype);
175: }
176: return(0);
177: }
181: static PetscErrorCode PCReset_LU(PC pc)
182: {
183: PC_LU *dir = (PC_LU*)pc->data;
187: if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
188: if (dir->row && dir->col && dir->row != dir->col) {ISDestroy(&dir->row);}
189: ISDestroy(&dir->col);
190: return(0);
191: }
195: static PetscErrorCode PCDestroy_LU(PC pc)
196: {
197: PC_LU *dir = (PC_LU*)pc->data;
201: PCReset_LU(pc);
202: PetscFree(((PC_Factor*)dir)->ordering);
203: PetscFree(((PC_Factor*)dir)->solvertype);
204: PetscFree(pc->data);
205: return(0);
206: }
210: static PetscErrorCode PCApply_LU(PC pc,Vec x,Vec y)
211: {
212: PC_LU *dir = (PC_LU*)pc->data;
216: if (dir->inplace) {
217: MatSolve(pc->pmat,x,y);
218: } else {
219: MatSolve(((PC_Factor*)dir)->fact,x,y);
220: }
221: return(0);
222: }
226: static PetscErrorCode PCApplyTranspose_LU(PC pc,Vec x,Vec y)
227: {
228: PC_LU *dir = (PC_LU*)pc->data;
232: if (dir->inplace) {
233: MatSolveTranspose(pc->pmat,x,y);
234: } else {
235: MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
236: }
237: return(0);
238: }
240: /* -----------------------------------------------------------------------------------*/
244: PetscErrorCode PCFactorSetUseInPlace_LU(PC pc,PetscBool flg)
245: {
246: PC_LU *dir = (PC_LU*)pc->data;
249: dir->inplace = flg;
250: return(0);
251: }
255: PetscErrorCode PCFactorGetUseInPlace_LU(PC pc,PetscBool *flg)
256: {
257: PC_LU *dir = (PC_LU*)pc->data;
260: *flg = dir->inplace;
261: return(0);
262: }
264: /* ------------------------------------------------------------------------ */
266: /*MC
267: PCLU - Uses a direct solver, based on LU factorization, as a preconditioner
269: Options Database Keys:
270: + -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
271: . -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
272: . -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
273: . -pc_factor_fill <fill> - Sets fill amount
274: . -pc_factor_in_place - Activates in-place factorization
275: . -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine
276: . -pc_factor_pivot_in_blocks <true,false> - allow pivoting within the small blocks during factorization (may increase
277: stability of factorization.
278: . -pc_factor_shift_type <shifttype> - Sets shift type or PETSC_DECIDE for the default; use '-help' for a list of available types
279: . -pc_factor_shift_amount <shiftamount> - Sets shift amount or PETSC_DECIDE for the default
280: - -pc_factor_nonzeros_along_diagonal - permutes the rows and columns to try to put nonzero value along the
281: diagonal.
283: Notes: Not all options work for all matrix formats
284: Run with -help to see additional options for particular matrix formats or factorization
285: algorithms
287: Level: beginner
289: Concepts: LU factorization, direct solver
291: Notes: Usually this will compute an "exact" solution in one iteration and does
292: not need a Krylov method (i.e. you can use -ksp_type preonly, or
293: KSPSetType(ksp,KSPPREONLY) for the Krylov method
295: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
296: PCILU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
297: PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
298: PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
299: PCFactorReorderForNonzeroDiagonal()
300: M*/
304: PETSC_EXTERN PetscErrorCode PCCreate_LU(PC pc)
305: {
307: PetscMPIInt size;
308: PC_LU *dir;
311: PetscNewLog(pc,&dir);
313: MatFactorInfoInitialize(&((PC_Factor*)dir)->info);
315: ((PC_Factor*)dir)->fact = NULL;
316: ((PC_Factor*)dir)->factortype = MAT_FACTOR_LU;
317: dir->inplace = PETSC_FALSE;
318: dir->nonzerosalongdiagonal = PETSC_FALSE;
320: ((PC_Factor*)dir)->info.fill = 5.0;
321: ((PC_Factor*)dir)->info.dtcol = 1.e-6; /* default to pivoting; this is only thing PETSc LU supports */
322: ((PC_Factor*)dir)->info.shifttype = (PetscReal)MAT_SHIFT_NONE;
323: ((PC_Factor*)dir)->info.shiftamount = 0.0;
324: ((PC_Factor*)dir)->info.zeropivot = 100.0*PETSC_MACHINE_EPSILON;
325: ((PC_Factor*)dir)->info.pivotinblocks = 1.0;
326: dir->col = 0;
327: dir->row = 0;
329: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
330: if (size == 1) {
331: PetscStrallocpy(MATORDERINGND,(char**)&((PC_Factor*)dir)->ordering);
332: } else {
333: PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
334: }
335: dir->reusefill = PETSC_FALSE;
336: dir->reuseordering = PETSC_FALSE;
337: pc->data = (void*)dir;
339: pc->ops->reset = PCReset_LU;
340: pc->ops->destroy = PCDestroy_LU;
341: pc->ops->apply = PCApply_LU;
342: pc->ops->applytranspose = PCApplyTranspose_LU;
343: pc->ops->setup = PCSetUp_LU;
344: pc->ops->setfromoptions = PCSetFromOptions_LU;
345: pc->ops->view = PCView_LU;
346: pc->ops->applyrichardson = 0;
347: pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;
349: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
350: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
351: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
352: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
353: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
354: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
355: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
356: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_LU);
357: PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_LU);
358: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
359: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_LU);
360: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_LU);
361: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetColumnPivot_C",PCFactorSetColumnPivot_Factor);
362: PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetPivotInBlocks_C",PCFactorSetPivotInBlocks_Factor);
363: PetscObjectComposeFunction((PetscObject)pc,"PCFactorReorderForNonzeroDiagonal_C",PCFactorReorderForNonzeroDiagonal_LU);
364: return(0);
365: }