Actual source code: cholesky.c

petsc-master 2017-07-18
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:   const MatSolverPackage 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:   PCFactorGetMatSolverPackage(pc,&stype);
134:   if (!stype) {
135:     const MatSolverPackage solverpackage;
136:     MatFactorGetSolverPackage(((PC_Factor*)dir)->fact,&solverpackage);
137:     PCFactorSetMatSolverPackage(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 PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
182: {
183:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

195: /* -----------------------------------------------------------------------------------*/

197: /* -----------------------------------------------------------------------------------*/

199: /*@
200:    PCFactorSetReuseOrdering - When similar matrices are factored, this
201:    causes the ordering computed in the first factor to be used for all
202:    following factors.

204:    Logically Collective on PC

206:    Input Parameters:
207: +  pc - the preconditioner context
208: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

210:    Options Database Key:
211: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

213:    Level: intermediate

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

217: .seealso: PCFactorSetReuseFill()
218: @*/
219: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
220: {

226:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
227:   return(0);
228: }

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

233:    Options Database Keys:
234: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
235: .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
236: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
237: .  -pc_factor_fill <fill> - Sets fill amount
238: .  -pc_factor_in_place - Activates in-place factorization
239: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

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

243:    Level: beginner

245:    Concepts: Cholesky factorization, direct solver

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

251: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
252:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
253:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
254:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

256: M*/

258: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
259: {
261:   PC_Cholesky    *dir;

264:   PetscNewLog(pc,&dir);
265:   pc->data = (void*)dir;
266:   PCFactorInitialize(pc);

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

271:   dir->col = 0;
272:   dir->row = 0;

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

276:   pc->ops->destroy           = PCDestroy_Cholesky;
277:   pc->ops->reset             = PCReset_Cholesky;
278:   pc->ops->apply             = PCApply_Cholesky;
279:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
280:   pc->ops->setup             = PCSetUp_Cholesky;
281:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
282:   pc->ops->view              = PCView_Cholesky;
283:   pc->ops->applyrichardson   = 0;
284:   return(0);
285: }