Actual source code: cholesky.c

petsc-dev 2014-04-23
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;

 25:   lu                = (PC_Cholesky*)pc->data;
 26:   lu->reuseordering = flag;
 27:   return(0);
 28: }

 32: static PetscErrorCode  PCFactorSetReuseFill_Cholesky(PC pc,PetscBool flag)
 33: {
 34:   PC_Cholesky *lu;

 37:   lu            = (PC_Cholesky*)pc->data;
 38:   lu->reusefill = flag;
 39:   return(0);
 40: }

 44: static PetscErrorCode PCSetFromOptions_Cholesky(PC pc)
 45: {

 49:   PetscOptionsHead("Cholesky options");
 50:   PCSetFromOptions_Factor(pc);
 51:   PetscOptionsTail();
 52:   return(0);
 53: }

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

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

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


 82: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 83: {
 85:   PetscBool      flg;
 86:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

 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);

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

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

134:         flg  = PETSC_FALSE;
135:         PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
136:         if (flg) {
137:           PetscReal tol = 1.e-10;
138:           PetscOptionsGetReal(((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
139:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
140:         }
141:         if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
142:       }
143:       MatDestroy(&((PC_Factor*)dir)->fact);
144:       MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
145:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
146:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
147:       dir->actualfill = info.fill_ratio_needed;
148:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
149:     }
150:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
151:   }
152:   return(0);
153: }

157: static PetscErrorCode PCReset_Cholesky(PC pc)
158: {
159:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

163:   if (!dir->inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
164:   ISDestroy(&dir->row);
165:   ISDestroy(&dir->col);
166:   return(0);
167: }

171: static PetscErrorCode PCDestroy_Cholesky(PC pc)
172: {
173:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

177:   PCReset_Cholesky(pc);
178:   PetscFree(((PC_Factor*)dir)->ordering);
179:   PetscFree(((PC_Factor*)dir)->solvertype);
180:   PetscFree(pc->data);
181:   return(0);
182: }

186: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
187: {
188:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

192:   if (dir->inplace) {
193:     MatSolve(pc->pmat,x,y);
194:   } else {
195:     MatSolve(((PC_Factor*)dir)->fact,x,y);
196:   }
197:   return(0);
198: }

202: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
203: {
204:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

216: /* -----------------------------------------------------------------------------------*/

220: static PetscErrorCode  PCFactorSetUseInPlace_Cholesky(PC pc)
221: {
222:   PC_Cholesky *dir;

225:   dir          = (PC_Cholesky*)pc->data;
226:   dir->inplace = PETSC_TRUE;
227:   return(0);
228: }

230: /* -----------------------------------------------------------------------------------*/

234: /*@
235:    PCFactorSetReuseOrdering - When similar matrices are factored, this
236:    causes the ordering computed in the first factor to be used for all
237:    following factors.

239:    Logically Collective on PC

241:    Input Parameters:
242: +  pc - the preconditioner context
243: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

245:    Options Database Key:
246: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

248:    Level: intermediate

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

252: .seealso: PCFactorSetReuseFill()
253: @*/
254: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
255: {

261:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
262:   return(0);
263: }


266: /*MC
267:    PCCHOLESKY - Uses a direct solver, based on Cholesky 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

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

279:    Level: beginner

281:    Concepts: Cholesky factorization, direct solver

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

287: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
288:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
289:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
290:            PCFactorSetUseInPlace(), PCFactorSetMatOrderingType()

292: M*/

296: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
297: {
299:   PC_Cholesky    *dir;

302:   PetscNewLog(pc,&dir);

304:   ((PC_Factor*)dir)->fact = 0;
305:   dir->inplace            = PETSC_FALSE;

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

309:   ((PC_Factor*)dir)->factortype         = MAT_FACTOR_CHOLESKY;
310:   ((PC_Factor*)dir)->info.fill          = 5.0;
311:   ((PC_Factor*)dir)->info.shifttype     = (PetscReal) MAT_SHIFT_NONE;
312:   ((PC_Factor*)dir)->info.shiftamount   = 0.0;
313:   ((PC_Factor*)dir)->info.zeropivot     = 100.0*PETSC_MACHINE_EPSILON;
314:   ((PC_Factor*)dir)->info.pivotinblocks = 1.0;

316:   dir->col = 0;
317:   dir->row = 0;

319:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);
320:   PetscStrallocpy(MATSOLVERPETSC,&((PC_Factor*)dir)->solvertype);

322:   dir->reusefill        = PETSC_FALSE;
323:   dir->reuseordering    = PETSC_FALSE;
324:   pc->data              = (void*)dir;

326:   pc->ops->destroy           = PCDestroy_Cholesky;
327:   pc->ops->reset             = PCReset_Cholesky;
328:   pc->ops->apply             = PCApply_Cholesky;
329:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
330:   pc->ops->setup             = PCSetUp_Cholesky;
331:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
332:   pc->ops->view              = PCView_Cholesky;
333:   pc->ops->applyrichardson   = 0;
334:   pc->ops->getfactoredmatrix = PCFactorGetMatrix_Factor;

336:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUpMatSolverPackage_C",PCFactorSetUpMatSolverPackage_Factor);
337:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatSolverPackage_C",PCFactorSetMatSolverPackage_Factor);
338:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorGetMatSolverPackage_C",PCFactorGetMatSolverPackage_Factor);
339:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetZeroPivot_C",PCFactorSetZeroPivot_Factor);
340:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftType_C",PCFactorSetShiftType_Factor);
341:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetShiftAmount_C",PCFactorSetShiftAmount_Factor);
342:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetFill_C",PCFactorSetFill_Factor);
343:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetUseInPlace_C",PCFactorSetUseInPlace_Cholesky);
344:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetMatOrderingType_C",PCFactorSetMatOrderingType_Factor);
345:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseOrdering_C",PCFactorSetReuseOrdering_Cholesky);
346:   PetscObjectComposeFunction((PetscObject)pc,"PCFactorSetReuseFill_C",PCFactorSetReuseFill_Cholesky);
347:   return(0);
348: }