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: {

 19:   PetscOptionsHead(PetscOptionsObject,"Cholesky options");
 20:   PCSetFromOptions_Factor(PetscOptionsObject,pc);
 21:   PetscOptionsTail();
 22:   return(0);
 23: }

 25: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 26: {
 27:   PetscErrorCode         ierr;
 28:   PetscBool              flg;
 29:   PC_Cholesky            *dir = (PC_Cholesky*)pc->data;
 30:   MatSolverType          stype;
 31:   MatFactorError         err;

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

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

 56:     ((PC_Factor*)dir)->fact = pc->pmat;
 57:   } else {
 58:     MatInfo info;

 60:     if (!pc->setupcalled) {
 61:       PetscBool useordering;
 62:       if (!((PC_Factor*)dir)->fact) {
 63:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 64:         PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 65:       }
 66:       MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
 67:       if (useordering) {
 68:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 69:         /* check if dir->row == dir->col */
 70:         ISEqual(dir->row,dir->col,&flg);
 71:         if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
 72:         ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 74:         flg  = PETSC_FALSE;
 75:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 76:         if (flg) {
 77:           PetscReal tol = 1.e-10;
 78:           PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 79:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 80:         }
 81:         PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);
 82:       }
 83:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 84:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
 85:       dir->hdr.actualfill = info.fill_ratio_needed;
 86:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 87:       if (!dir->hdr.reuseordering) {
 88:         PetscBool useordering;
 89:         MatDestroy(&((PC_Factor*)dir)->fact);
 90:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 91:         PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 92:         MatFactorGetUseOrdering(((PC_Factor*)dir)->fact,&useordering);
 93:         if (useordering) {
 94:           ISDestroy(&dir->row);
 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;

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

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

159:   PCReset_Cholesky(pc);
160:   PetscFree(((PC_Factor*)dir)->ordering);
161:   PetscFree(((PC_Factor*)dir)->solvertype);
162:   PetscFree(pc->data);
163:   return(0);
164: }

166: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
167: {
168:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

172:   if (dir->hdr.inplace) {
173:     MatSolve(pc->pmat,x,y);
174:   } else {
175:     MatSolve(((PC_Factor*)dir)->fact,x,y);
176:   }
177:   return(0);
178: }

180: static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y)
181: {
182:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

186:   if (dir->hdr.inplace) {
187:     MatMatSolve(pc->pmat,X,Y);
188:   } else {
189:     MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
190:   }
191:   return(0);
192: }

194: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)
195: {
196:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

208: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
209: {
210:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

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

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

236: /* -----------------------------------------------------------------------------------*/

238: /* -----------------------------------------------------------------------------------*/

240: /*@
241:    PCFactorSetReuseOrdering - When similar matrices are factored, this
242:    causes the ordering computed in the first factor to be used for all
243:    following factors.

245:    Logically Collective on PC

247:    Input Parameters:
248: +  pc - the preconditioner context
249: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

251:    Options Database Key:
252: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

254:    Level: intermediate

256: .seealso: PCFactorSetReuseFill()
257: @*/
258: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
259: {

265:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
266:   return(0);
267: }

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

272:    Options Database Keys:
273: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
274: .  -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
275: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
276: .  -pc_factor_fill <fill> - Sets fill amount
277: .  -pc_factor_in_place - Activates in-place factorization
278: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

280:    Notes:
281:     Not all options work for all matrix formats

283:    Level: beginner

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

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

295: M*/

297: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
298: {
300:   PC_Cholesky    *dir;

303:   PetscNewLog(pc,&dir);
304:   pc->data = (void*)dir;
305:   PCFactorInitialize(pc);

307:   ((PC_Factor*)dir)->factortype         = MAT_FACTOR_CHOLESKY;
308:   ((PC_Factor*)dir)->info.fill          = 5.0;

310:   dir->col = NULL;
311:   dir->row = NULL;

313:   /* MATORDERINGNATURAL_OR_ND allows selecting type based on matrix type sbaij or aij */
314:   PetscStrallocpy(MATORDERINGNATURAL_OR_ND,(char**)&((PC_Factor*)dir)->ordering);

316:   pc->ops->destroy             = PCDestroy_Cholesky;
317:   pc->ops->reset               = PCReset_Cholesky;
318:   pc->ops->apply               = PCApply_Cholesky;
319:   pc->ops->matapply            = PCMatApply_Cholesky;
320:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
321:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
322:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
323:   pc->ops->setup               = PCSetUp_Cholesky;
324:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
325:   pc->ops->view                = PCView_Factor;
326:   pc->ops->applyrichardson     = NULL;
327:   return(0);
328: }