Actual source code: cholesky.c

petsc-3.7.3 2016-07-24
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: {
 82:   PetscBool      flg;
 83:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;
 84:   const MatSolverPackage stype;

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

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

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

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

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

160:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
161:     if (F->errortype) { /* FactorNumeric() fails */
162:       pc->failedreason = (PCFailedReason)F->errortype;
163:     }
164:   }

166:   PCFactorGetMatSolverPackage(pc,&stype);
167:   if (!stype) {
168:     PCFactorSetMatSolverPackage(pc,((PC_Factor*)dir)->fact->solvertype);
169:   }
170:   return(0);
171: }

175: static PetscErrorCode PCReset_Cholesky(PC pc)
176: {
177:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

181:   if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
182:   ISDestroy(&dir->row);
183:   ISDestroy(&dir->col);
184:   return(0);
185: }

189: static PetscErrorCode PCDestroy_Cholesky(PC pc)
190: {
191:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

195:   PCReset_Cholesky(pc);
196:   PetscFree(((PC_Factor*)dir)->ordering);
197:   PetscFree(((PC_Factor*)dir)->solvertype);
198:   PetscFree(pc->data);
199:   return(0);
200: }

204: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
205: {
206:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

210:   if (dir->inplace) {
211:     MatSolve(pc->pmat,x,y);
212:   } else {
213:     MatSolve(((PC_Factor*)dir)->fact,x,y);
214:   }
215:   return(0);
216: }

220: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
221: {
222:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

226:   if (dir->inplace) {
227:     MatSolveTranspose(pc->pmat,x,y);
228:   } else {
229:     MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
230:   }
231:   return(0);
232: }

234: /* -----------------------------------------------------------------------------------*/

238: static PetscErrorCode  PCFactorSetUseInPlace_Cholesky(PC pc,PetscBool flg)
239: {
240:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;

243:   dir->inplace = flg;
244:   return(0);
245: }

249: static PetscErrorCode  PCFactorGetUseInPlace_Cholesky(PC pc,PetscBool *flg)
250: {
251:   PC_Cholesky *dir = (PC_Cholesky*)pc->data;

254:   *flg = dir->inplace;
255:   return(0);
256: }

258: /* -----------------------------------------------------------------------------------*/

262: /*@
263:    PCFactorSetReuseOrdering - When similar matrices are factored, this
264:    causes the ordering computed in the first factor to be used for all
265:    following factors.

267:    Logically Collective on PC

269:    Input Parameters:
270: +  pc - the preconditioner context
271: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

273:    Options Database Key:
274: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

276:    Level: intermediate

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

280: .seealso: PCFactorSetReuseFill()
281: @*/
282: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
283: {

289:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
290:   return(0);
291: }

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

296:    Options Database Keys:
297: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
298: .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
299: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
300: .  -pc_factor_fill <fill> - Sets fill amount
301: .  -pc_factor_in_place - Activates in-place factorization
302: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

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

306:    Level: beginner

308:    Concepts: Cholesky factorization, direct solver

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

314: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
315:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
316:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
317:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

319: M*/

323: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
324: {
326:   PC_Cholesky    *dir;

329:   PetscNewLog(pc,&dir);

331:   ((PC_Factor*)dir)->fact = 0;
332:   dir->inplace            = PETSC_FALSE;

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

336:   ((PC_Factor*)dir)->factortype         = MAT_FACTOR_CHOLESKY;
337:   ((PC_Factor*)dir)->info.fill          = 5.0;
338:   ((PC_Factor*)dir)->info.shifttype     = (PetscReal) MAT_SHIFT_NONE;
339:   ((PC_Factor*)dir)->info.shiftamount   = 0.0;
340:   ((PC_Factor*)dir)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
341:   ((PC_Factor*)dir)->info.pivotinblocks = 1.0;

343:   dir->col = 0;
344:   dir->row = 0;

346:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
347:   dir->reusefill        = PETSC_FALSE;
348:   dir->reuseordering    = PETSC_FALSE;
349:   pc->data              = (void*)dir;

351:   pc->ops->destroy           = PCDestroy_Cholesky;
352:   pc->ops->reset             = PCReset_Cholesky;
353:   pc->ops->apply             = PCApply_Cholesky;
354:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
355:   pc->ops->setup             = PCSetUp_Cholesky;
356:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
357:   pc->ops->view              = PCView_Cholesky;
358:   pc->ops->applyrichardson   = 0;
359:   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;

361:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
362:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
363:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
364:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
365:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
366:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
367:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
368:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
369:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetUseInPlace_C",PCFactorGetUseInPlace_Cholesky);
370:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
371:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
372:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
373:   return(0);
374: }