Actual source code: cholesky.c

petsc-3.10.4 2019-02-26
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 PCView_Cholesky(PC pc,PetscViewer viewer)
 26: {
 28:   PetscBool      iascii;

 31:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 32:   PCView_Factor(pc,viewer);
 33:   return(0);
 34: }

 36: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 37: {
 38:   PetscErrorCode         ierr;
 39:   PetscBool              flg;
 40:   PC_Cholesky            *dir = (PC_Cholesky*)pc->data;
 41:   MatSolverType          stype;
 42:   MatFactorError         err;

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

 48:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 49:   if (dir->hdr.inplace) {
 50:     if (dir->row && dir->col && (dir->row != dir->col)) {
 51:       ISDestroy(&dir->row);
 52:     }
 53:     ISDestroy(&dir->col);
 54:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 55:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
 56:       ISDestroy(&dir->col);
 57:     }
 58:     if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 59:     MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 60:     MatFactorGetError(pc->pmat,&err);
 61:     if (err) { /* Factor() fails */
 62:       pc->failedreason = (PCFailedReason)err;
 63:       return(0);
 64:     }

 66:     ((PC_Factor*)dir)->fact = pc->pmat;
 67:   } else {
 68:     MatInfo info;

 70:     if (!pc->setupcalled) {
 71:       MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 72:       /* check if dir->row == dir->col */
 73:       ISEqual(dir->row,dir->col,&flg);
 74:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
 75:       ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 77:       flg  = PETSC_FALSE;
 78:       PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 79:       if (flg) {
 80:         PetscReal tol = 1.e-10;
 81:         PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 82:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 83:       }
 84:       if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 85:       if (!((PC_Factor*)dir)->fact) {
 86:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 87:       }
 88:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 89:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
 90:       dir->hdr.actualfill = info.fill_ratio_needed;
 91:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 92:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 93:       if (!dir->hdr.reuseordering) {
 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:         if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
106:       }
107:       MatDestroy(&((PC_Factor*)dir)->fact);
108:       MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
109:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
110:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
111:       dir->hdr.actualfill = info.fill_ratio_needed;
112:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
113:     } else {
114:       MatFactorGetError(((PC_Factor*)dir)->fact,&err);
115:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
116:         MatFactorClearError(((PC_Factor*)dir)->fact);
117:         pc->failedreason = PC_NOERROR;
118:       }
119:     }
120:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
121:     if (err) { /* FactorSymbolic() fails */
122:       pc->failedreason = (PCFailedReason)err;
123:       return(0);
124:     }

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

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

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

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

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

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

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

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

181: static PetscErrorCode PCApplySymmetricLeft_Cholesky(PC pc,Vec x,Vec y)
182: {
183:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

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

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

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

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

223: /* -----------------------------------------------------------------------------------*/

225: /* -----------------------------------------------------------------------------------*/

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

232:    Logically Collective on PC

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

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

241:    Level: intermediate

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

245: .seealso: PCFactorSetReuseFill()
246: @*/
247: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
248: {

254:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
255:   return(0);
256: }

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

261:    Options Database Keys:
262: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
263: .  -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
264: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
265: .  -pc_factor_fill <fill> - Sets fill amount
266: .  -pc_factor_in_place - Activates in-place factorization
267: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

269:    Notes:
270:     Not all options work for all matrix formats

272:    Level: beginner

274:    Concepts: Cholesky factorization, direct solver

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

281: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
282:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
283:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
284:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

286: M*/

288: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
289: {
291:   PC_Cholesky    *dir;

294:   PetscNewLog(pc,&dir);
295:   pc->data = (void*)dir;
296:   PCFactorInitialize(pc);

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

301:   dir->col = 0;
302:   dir->row = 0;

304:   PetscStrallocpy(MATORDERINGNATURAL,(char**)&((PC_Factor*)dir)->ordering);

306:   pc->ops->destroy             = PCDestroy_Cholesky;
307:   pc->ops->reset               = PCReset_Cholesky;
308:   pc->ops->apply               = PCApply_Cholesky;
309:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
310:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
311:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
312:   pc->ops->setup               = PCSetUp_Cholesky;
313:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
314:   pc->ops->view                = PCView_Cholesky;
315:   pc->ops->applyrichardson     = 0;
316:   return(0);
317: }