Actual source code: cholesky.c

petsc-main 2021-04-20
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>

  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:     PCFactorSetDefaultOrdering_Factor(pc);
 45:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 46:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
 47:       ISDestroy(&dir->col);
 48:     }
 49:     if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 50:     MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 51:     MatFactorGetError(pc->pmat,&err);
 52:     if (err) { /* Factor() fails */
 53:       pc->failedreason = (PCFailedReason)err;
 54:       return(0);
 55:     }

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

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

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

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

130:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
131:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
132:     if (err) { /* FactorNumeric() fails */
133:       pc->failedreason = (PCFailedReason)err;
134:     }
135:   }

137:   PCFactorGetMatSolverType(pc,&stype);
138:   if (!stype) {
139:     MatSolverType solverpackage;
140:     MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
141:     PCFactorSetMatSolverType(pc,solverpackage);
142:   }
143:   return(0);
144: }

146: static PetscErrorCode PCReset_Cholesky(PC pc)
147: {
148:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

152:   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
153:   ISDestroy(&dir->row);
154:   ISDestroy(&dir->col);
155:   return(0);
156: }

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

164:   PCReset_Cholesky(pc);
165:   PetscFree(((PC_Factor*)dir)->ordering);
166:   PetscFree(((PC_Factor*)dir)->solvertype);
167:   PetscFree(pc->data);
168:   return(0);
169: }

171: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
172: {
173:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

177:   if (dir->hdr.inplace) {
178:     MatSolve(pc->pmat,x,y);
179:   } else {
180:     MatSolve(((PC_Factor*)dir)->fact,x,y);
181:   }
182:   return(0);
183: }

185: static PetscErrorCode PCMatApply_Cholesky(PC pc,Mat X,Mat Y)
186: {
187:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

191:   if (dir->hdr.inplace) {
192:     MatMatSolve(pc->pmat,X,Y);
193:   } else {
194:     MatMatSolve(((PC_Factor*)dir)->fact,X,Y);
195:   }
196:   return(0);
197: }

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

205:   if (dir->hdr.inplace) {
206:     MatForwardSolve(pc->pmat,x,y);
207:   } else {
208:     MatForwardSolve(((PC_Factor*)dir)->fact,x,y);
209:   }
210:   return(0);
211: }

213: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
214: {
215:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

219:   if (dir->hdr.inplace) {
220:     MatBackwardSolve(pc->pmat,x,y);
221:   } else {
222:     MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
223:   }
224:   return(0);
225: }

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

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

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

243: /* -----------------------------------------------------------------------------------*/

245: /*@
246:    PCFactorSetReuseOrdering - When similar matrices are factored, this
247:    causes the ordering computed in the first factor to be used for all
248:    following factors.

250:    Logically Collective on PC

252:    Input Parameters:
253: +  pc - the preconditioner context
254: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

256:    Options Database Key:
257: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

259:    Level: intermediate

261: .seealso: PCFactorSetReuseFill()
262: @*/
263: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
264: {

270:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
271:   return(0);
272: }

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

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

285:    Notes:
286:     Not all options work for all matrix formats

288:    Level: beginner

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

295: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
296:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
297:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
298:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

300: M*/

302: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
303: {
305:   PC_Cholesky    *dir;

308:   PetscNewLog(pc,&dir);
309:   pc->data = (void*)dir;
310:   PCFactorInitialize(pc,MAT_FACTOR_CHOLESKY);

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

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