Actual source code: cholesky.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  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: */
  7: #include <../src/ksp/pc/impls/factor/factor.h>         /*I "petscpc.h" I*/

  9: typedef struct {
 10:   PC_Factor hdr;
 11:   PetscReal actualfill;              /* actual fill in factor */
 12:   PetscBool inplace;                 /* flag indicating in-place factorization */
 13:   IS        row,col;                 /* index sets used for reordering */
 14:   PetscBool reuseordering;           /* reuses previous reordering computed */
 15:   PetscBool reusefill;               /* reuse fill from previous Cholesky */
 16: } PC_Cholesky;

 20: static PetscErrorCode  PCFactorSetReuseOrdering_Cholesky(PC pc,PetscBool flag)
 21: {
 22:   PC_Cholesky *lu = (PC_Cholesky*)pc->data;

 25:   lu->reuseordering = flag;
 26:   return(0);
 27: }

 31: static PetscErrorCode  PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag)
 32: {
 33:   PC_Cholesky *lu = (PC_Cholesky*)pc->data;

 36:   lu->reusefill = flag;
 37:   return(0);
 38: }

 42: static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptionItems *PetscOptionsObject,PC pc)
 43: {

 47:   PetscOptionsHead(PetscOptionsObject,"Cholesky options");
 48:   PCSetFromOptions_Factor(PetscOptionsObject,pc);
 49:   PetscOptionsTail();
 50:   return(0);
 51: }

 55: static PetscErrorCode PCView_Cholesky(PC pc,PetscViewer viewer)
 56: {
 57:   PC_Cholesky    *chol = (PC_Cholesky*)pc->data;
 59:   PetscBool      iascii;

 62:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 63:   if (iascii) {
 64:     if (chol->inplace) {
 65:       PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorization\n");
 66:     } else {
 67:       PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorization\n");
 68:     }

 70:     if (chol->reusefill)    {PetscViewerASCIIPrintf(viewer,"  Reusing fill from past factorization\n");}
 71:     if (chol->reuseordering) {PetscViewerASCIIPrintf(viewer,"  Reusing reordering from past factorization\n");}
 72:   }
 73:   PCView_Factor(pc,viewer);
 74:   return(0);
 75: }

 79: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 80: {
 81:   PetscErrorCode         ierr;
 82:   PetscBool              flg;
 83:   PC_Cholesky            *dir = (PC_Cholesky*)pc->data;
 84:   const MatSolverPackage stype;

 87:   pc->failedreason = PC_NOERROR;
 88:   if (dir->reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->actualfill;

 90:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 91:   if (dir->inplace) {
 92:     if (dir->row && dir->col && (dir->row != dir->col)) {
 93:       ISDestroy(&dir->row);
 94:     }
 95:     ISDestroy(&dir->col);
 96:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 97:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
 98:       ISDestroy(&dir->col);
 99:     }
100:     if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
101:     MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
102:     if (pc->pmat->errortype) { /* Factor() fails */
103:       pc->failedreason = (PCFailedReason)pc->pmat->errortype;
104:       return(0);
105:     }

107:     ((PC_Factor*)dir)->fact = pc->pmat;
108:   } else {
109:     MatInfo info;
110:     Mat     F;
111:     if (!pc->setupcalled) {
112:       MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
113:       /* check if dir->row == dir->col */
114:       ISEqual(dir->row,dir->col,&flg);
115:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
116:       ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

118:       flg  = PETSC_FALSE;
119:       PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
120:       if (flg) {
121:         PetscReal tol = 1.e-10;
122:         PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
123:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
124:       }
125:       if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
126:       if (!((PC_Factor*)dir)->fact) {
127:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
128:       }
129:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
130:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
131:       dir->actualfill = info.fill_ratio_needed;
132:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
133:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
134:       if (!dir->reuseordering) {
135:         ISDestroy(&dir->row);
136:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
137:         ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */

139:         flg  = PETSC_FALSE;
140:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
141:         if (flg) {
142:           PetscReal tol = 1.e-10;
143:           PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
144:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
145:         }
146:         if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
147:       }
148:       MatDestroy(&((PC_Factor*)dir)->fact);
149:       MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
150:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
151:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
152:       dir->actualfill = info.fill_ratio_needed;
153:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
154:     } else {
155:       F = ((PC_Factor*)dir)->fact;
156:       if ((PCFailedReason)F->errortype == PC_FACTOR_NUMERIC_ZEROPIVOT) {
157:         F->errortype     = MAT_FACTOR_NOERROR;
158:         pc->failedreason = (PCFailedReason)F->errortype;
159:       }
160:     }
161:     F = ((PC_Factor*)dir)->fact;
162:     if (F->errortype) { /* FactorSymbolic() fails */
163:       pc->failedreason = (PCFailedReason)F->errortype;
164:       return(0);
165:     }

167:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
168:     if (F->errortype) { /* FactorNumeric() fails */
169:       pc->failedreason = (PCFailedReason)F->errortype;
170:     }
171:   }

173:   PCFactorGetMatSolverPackage(pc,&stype);
174:   if (!stype) {
175:     PCFactorSetMatSolverPackage(pc,((PC_Factor*)dir)->fact->solvertype);
176:   }
177:   return(0);
178: }

182: static PetscErrorCode PCReset_Cholesky(PC pc)
183: {
184:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

188:   if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
189:   ISDestroy(&dir->row);
190:   ISDestroy(&dir->col);
191:   return(0);
192: }

196: static PetscErrorCode PCDestroy_Cholesky(PC pc)
197: {
198:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

202:   PCReset_Cholesky(pc);
203:   PetscFree(((PC_Factor*)dir)->ordering);
204:   PetscFree(((PC_Factor*)dir)->solvertype);
205:   PetscFree(pc->data);
206:   return(0);
207: }

211: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
212: {
213:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

217:   if (dir->inplace) {
218:     MatSolve(pc->pmat,x,y);
219:   } else {
220:     MatSolve(((PC_Factor*)dir)->fact,x,y);
221:   }
222:   return(0);
223: }

227: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
228: {
229:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

233:   if (dir->inplace) {
234:     MatSolveTranspose(pc->pmat,x,y);
235:   } else {
236:     MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
237:   }
238:   return(0);
239: }

241: /* -----------------------------------------------------------------------------------*/

245: static PetscErrorCode  PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg)
246: {
247:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;

250:   dir->inplace = flg;
251:   return(0);
252: }

256: static PetscErrorCode  PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg)
257: {
258:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;

261:   *flg = dir->inplace;
262:   return(0);
263: }

265: /* -----------------------------------------------------------------------------------*/

269: /*@
270:    PCFactorSetReuseOrdering - When similar matrices are factored, this
271:    causes the ordering computed in the first factor to be used for all
272:    following factors.

274:    Logically Collective on PC

276:    Input Parameters:
277: +  pc - the preconditioner context
278: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

280:    Options Database Key:
281: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

283:    Level: intermediate

285: .keywords: PC, levels, reordering, factorization, incomplete, LU

287: .seealso: PCFactorSetReuseFill()
288: @*/
289: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
290: {

296:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
297:   return(0);
298: }

300: /*MC
301:    PCCHOLESKY - Uses a direct solver, based on Cholesky factorization, as a preconditioner

303:    Options Database Keys:
304: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
305: .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
306: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
307: .  -pc_factor_fill <fill> - Sets fill amount
308: .  -pc_factor_in_place - Activates in-place factorization
309: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

311:    Notes: Not all options work for all matrix formats

313:    Level: beginner

315:    Concepts: Cholesky factorization, direct solver

317:    Notes: Usually this will compute an "exact" solution in one iteration and does
318:           not need a Krylov method (i.e. you can use -ksp_type preonly, or
319:           KSPSetType(ksp,KSPPREONLY) for the Krylov method

321: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
322:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
323:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
324:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

326: M*/

330: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
331: {
333:   PC_Cholesky    *dir;

336:   PetscNewLog(pc,&dir);

338:   ((PC_Factor*)dir)->fact = 0;
339:   dir->inplace            = PETSC_FALSE;

341:   MatFactorInfoInitialize(&((PC_Factor*)dir)->info);

343:   ((PC_Factor*)dir)->factortype         = MAT_FACTOR_CHOLESKY;
344:   ((PC_Factor*)dir)->info.fill          = 5.0;
345:   ((PC_Factor*)dir)->info.shifttype     = (PetscReal) MAT_SHIFT_NONE;
346:   ((PC_Factor*)dir)->info.shiftamount   = 0.0;
347:   ((PC_Factor*)dir)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
348:   ((PC_Factor*)dir)->info.pivotinblocks = 1.0;

350:   dir->col = 0;
351:   dir->row = 0;

353:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
354:   dir->reusefill        = PETSC_FALSE;
355:   dir->reuseordering    = PETSC_FALSE;
356:   pc->data              = (void*)dir;

358:   pc->ops->destroy           = PCDestroy_Cholesky;
359:   pc->ops->reset             = PCReset_Cholesky;
360:   pc->ops->apply             = PCApply_Cholesky;
361:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
362:   pc->ops->setup             = PCSetUp_Cholesky;
363:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
364:   pc->ops->view              = PCView_Cholesky;
365:   pc->ops->applyrichardson   = 0;
366:   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;

368:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
369:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetZeroPivot_C",PCFactorGetZeroPivot_Factor);
370:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
371:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftType_C",PCFactorGetShiftType_Factor);
372:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
373:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetShiftAmount_C",PCFactorGetShiftAmount_Factor);
374:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
375:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
376:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
377:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
378:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
379:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
380:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);
381:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
382:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
383:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
384:   return(0);
385: }