Actual source code: cholesky.c


  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>

  9: typedef struct {
 10:   PC_Factor hdr;
 11:   IS        row,col;                 /* index sets used for reordering */
 12: } PC_Cholesky;

 14: static PetscErrorCode PCSetFromOptions_Cholesky(PetscOptionItems *PetscOptionsObject,PC pc)
 15: {
 16:   PetscOptionsHead(PetscOptionsObject,"Cholesky options");
 17:   PCSetFromOptions_Factor(PetscOptionsObject,pc);
 18:   PetscOptionsTail();
 19:   return 0;
 20: }

 22: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 23: {
 24:   PetscBool              flg;
 25:   PC_Cholesky            *dir = (PC_Cholesky*)pc->data;
 26:   MatSolverType          stype;
 27:   MatFactorError         err;

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

 32:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 33:   if (dir->hdr.inplace) {
 34:     if (dir->row && dir->col && (dir->row != dir->col)) {
 35:       ISDestroy(&dir->row);
 36:     }
 37:     ISDestroy(&dir->col);
 38:     /* should only get reordering if the factor matrix uses it but cannot determine because MatGetFactor() not called */
 39:     PCFactorSetDefaultOrdering_Factor(pc);
 40:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 41:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
 42:       ISDestroy(&dir->col);
 43:     }
 44:     if (dir->row) PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
 45:     MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 46:     MatFactorGetError(pc->pmat,&err);
 47:     if (err) { /* Factor() fails */
 48:       pc->failedreason = (PCFailedReason)err;
 49:       return 0;
 50:     }

 52:     ((PC_Factor*)dir)->fact = pc->pmat;
 53:   } else {
 54:     MatInfo info;

 56:     if (!pc->setupcalled) {
 57:       PetscBool canuseordering;
 58:       if (!((PC_Factor*)dir)->fact) {
 59:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 60:         PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 61:       }
 62:       MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
 63:       if (canuseordering) {
 64:         PCFactorSetDefaultOrdering_Factor(pc);
 65:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 66:         /* check if dir->row == dir->col */
 67:         if (dir->row) {
 68:           ISEqual(dir->row,dir->col,&flg);
 70:         }
 71:         ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 73:         flg  = PETSC_FALSE;
 74:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 75:         if (flg) {
 76:           PetscReal tol = 1.e-10;
 77:           PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 78:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 79:         }
 80:         PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
 81:       }
 82:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 83:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
 84:       dir->hdr.actualfill = info.fill_ratio_needed;
 85:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 86:       if (!dir->hdr.reuseordering) {
 87:         PetscBool canuseordering;
 88:         MatDestroy(&((PC_Factor*)dir)->fact);
 89:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 90:         PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 91:         MatFactorGetCanUseOrdering(((PC_Factor*)dir)->fact,&canuseordering);
 92:         if (canuseordering) {
 93:           ISDestroy(&dir->row);
 94:           PCFactorSetDefaultOrdering_Factor(pc);
 95:           MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 96:           ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */

 98:           flg  = PETSC_FALSE;
 99:           PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
100:           if (flg) {
101:             PetscReal tol = 1.e-10;
102:             PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
103:             MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
104:           }
105:           PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
106:         }
107:       }
108:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
109:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
110:       dir->hdr.actualfill = info.fill_ratio_needed;
111:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
112:     } else {
113:       MatFactorGetError(((PC_Factor*)dir)->fact,&err);
114:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
115:         MatFactorClearError(((PC_Factor*)dir)->fact);
116:         pc->failedreason = PC_NOERROR;
117:       }
118:     }
119:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
120:     if (err) { /* FactorSymbolic() fails */
121:       pc->failedreason = (PCFailedReason)err;
122:       return 0;
123:     }

125:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
126:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
127:     if (err) { /* FactorNumeric() fails */
128:       pc->failedreason = (PCFailedReason)err;
129:     }
130:   }

132:   PCFactorGetMatSolverType(pc,&stype);
133:   if (!stype) {
134:     MatSolverType solverpackage;
135:     MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
136:     PCFactorSetMatSolverType(pc,solverpackage);
137:   }
138:   return 0;
139: }

141: static PetscErrorCode PCReset_Cholesky(PC pc)
142: {
143:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

145:   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) MatDestroy(&((PC_Factor*)dir)->fact);
146:   ISDestroy(&dir->row);
147:   ISDestroy(&dir->col);
148:   return 0;
149: }

151: static PetscErrorCode PCDestroy_Cholesky(PC pc)
152: {
153:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

155:   PCReset_Cholesky(pc);
156:   PetscFree(((PC_Factor*)dir)->ordering);
157:   PetscFree(((PC_Factor*)dir)->solvertype);
158:   PetscFree(pc->data);
159:   return 0;
160: }

162: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
163: {
164:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

166:   if (dir->hdr.inplace) {
167:     MatSolve(pc->pmat,x,y);
168:   } else {
169:     MatSolve(((PC_Factor*)dir)->fact,x,y);
170:   }
171:   return 0;
172: }

174: static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y)
175: {
176:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

178:   if (dir->hdr.inplace) {
179:     MatMatSolve(pc->pmat,X,Y);
180:   } else {
181:     MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
182:   }
183:   return 0;
184: }

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

190:   if (dir->hdr.inplace) {
191:     MatForwardSolve(pc->pmat,x,y);
192:   } else {
193:     MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
194:   }
195:   return 0;
196: }

198: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
199: {
200:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

202:   if (dir->hdr.inplace) {
203:     MatBackwardSolve(pc->pmat,x,y);
204:   } else {
205:     MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
206:   }
207:   return 0;
208: }

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

214:   if (dir->hdr.inplace) {
215:     MatSolveTranspose(pc->pmat,x,y);
216:   } else {
217:     MatSolveTranspose(((PC_Factor*)dir)->fact,x,y);
218:   }
219:   return 0;
220: }

222: /* -----------------------------------------------------------------------------------*/

224: /* -----------------------------------------------------------------------------------*/

226: /*@
227:    PCFactorSetReuseOrdering - When similar matrices are factored, this
228:    causes the ordering computed in the first factor to be used for all
229:    following factors.

231:    Logically Collective on PC

233:    Input Parameters:
234: +  pc - the preconditioner context
235: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

237:    Options Database Key:
238: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

240:    Level: intermediate

242: .seealso: PCFactorSetReuseFill()
243: @*/
244: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
245: {
248:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
249:   return 0;
250: }

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

255:    Options Database Keys:
256: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
257: .  -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
258: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
259: .  -pc_factor_fill <fill> - Sets fill amount
260: .  -pc_factor_in_place - Activates in-place factorization
261: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

263:    Notes:
264:     Not all options work for all matrix formats

266:    Level: beginner

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

273: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
274:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
275:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
276:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

278: M*/

280: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
281: {
282:   PC_Cholesky    *dir;

284:   PetscNewLog(pc,&dir);
285:   pc->data = (void*)dir;
286:   PCFactorInitialize(pc,MAT_FACTOR_CHOLESKY);

288:   ((PC_Factor*)dir)->info.fill  = 5.0;

290:   pc->ops->destroy             = PCDestroy_Cholesky;
291:   pc->ops->reset               = PCReset_Cholesky;
292:   pc->ops->apply               = PCApply_Cholesky;
293:   pc->ops->matapply            = PCMatApply_Cholesky;
294:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
295:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
296:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
297:   pc->ops->setup               = PCSetUp_Cholesky;
298:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
299:   pc->ops->view                = PCView_Factor;
300:   pc->ops->applyrichardson     = NULL;
301:   return 0;
302: }