Actual source code: cholesky.c

petsc-master 2017-03-28
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: {
 27:   PC_Cholesky    *chol = (PC_Cholesky*)pc->data;
 29:   PetscBool      iascii;

 32:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 33:   if (iascii) {
 34:     if (chol->hdr.inplace) {
 35:       PetscViewerASCIIPrintf(viewer,"  Cholesky: in-place factorization\n");
 36:     } else {
 37:       PetscViewerASCIIPrintf(viewer,"  Cholesky: out-of-place factorization\n");
 38:     }

 40:     if (chol->hdr.reusefill)    {PetscViewerASCIIPrintf(viewer,"  Reusing fill from past factorization\n");}
 41:     if (chol->hdr.reuseordering) {PetscViewerASCIIPrintf(viewer,"  Reusing reordering from past factorization\n");}
 42:   }
 43:   PCView_Factor(pc,viewer);
 44:   return(0);
 45: }

 47: static PetscErrorCode PCSetUp_Cholesky(PC pc)
 48: {
 49:   PetscErrorCode         ierr;
 50:   PetscBool              flg;
 51:   PC_Cholesky            *dir = (PC_Cholesky*)pc->data;
 52:   const MatSolverPackage stype;
 53:   MatFactorError         err;

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

 59:   MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
 60:   if (dir->hdr.inplace) {
 61:     if (dir->row && dir->col && (dir->row != dir->col)) {
 62:       ISDestroy(&dir->row);
 63:     }
 64:     ISDestroy(&dir->col);
 65:     MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 66:     if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
 67:       ISDestroy(&dir->col);
 68:     }
 69:     if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 70:     MatCholeskyFactor(pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 71:     MatFactorGetError(pc->pmat,&err);
 72:     if (err) { /* Factor() fails */
 73:       pc->failedreason = (PCFailedReason)err;
 74:       return(0);
 75:     }

 77:     ((PC_Factor*)dir)->fact = pc->pmat;
 78:   } else {
 79:     MatInfo info;

 81:     if (!pc->setupcalled) {
 82:       MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 83:       /* check if dir->row == dir->col */
 84:       ISEqual(dir->row,dir->col,&flg);
 85:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
 86:       ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 88:       flg  = PETSC_FALSE;
 89:       PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 90:       if (flg) {
 91:         PetscReal tol = 1.e-10;
 92:         PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 93:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 94:       }
 95:       if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 96:       if (!((PC_Factor*)dir)->fact) {
 97:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 98:       }
 99:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
100:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
101:       dir->hdr.actualfill = info.fill_ratio_needed;
102:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
103:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
104:       if (!dir->hdr.reuseordering) {
105:         ISDestroy(&dir->row);
106:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
107:         ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */

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

137:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
138:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
139:     if (err) { /* FactorNumeric() fails */
140:       pc->failedreason = (PCFailedReason)err;
141:     }
142:   }

144:   PCFactorGetMatSolverPackage(pc,&stype);
145:   if (!stype) {
146:     const MatSolverPackage solverpackage;
147:     MatFactorGetSolverPackage(((PC_Factor*)dir)->fact,&solverpackage);
148:     PCFactorSetMatSolverPackage(pc,solverpackage);
149:   }
150:   return(0);
151: }

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

159:   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
160:   ISDestroy(&dir->row);
161:   ISDestroy(&dir->col);
162:   return(0);
163: }

165: static PetscErrorCode PCDestroy_Cholesky(PC pc)
166: {
167:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

171:   PCReset_Cholesky(pc);
172:   PetscFree(((PC_Factor*)dir)->ordering);
173:   PetscFree(((PC_Factor*)dir)->solvertype);
174:   PetscFree(pc->data);
175:   return(0);
176: }

178: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
179: {
180:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

184:   if (dir->hdr.inplace) {
185:     MatSolve(pc->pmat,x,y);
186:   } else {
187:     MatSolve(((PC_Factor*)dir)->fact,x,y);
188:   }
189:   return(0);
190: }

192: static PetscErrorCode PCApplyTranspose_Cholesky(PC pc,Vec x,Vec y)
193: {
194:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

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

206: /* -----------------------------------------------------------------------------------*/

208: /* -----------------------------------------------------------------------------------*/

210: /*@
211:    PCFactorSetReuseOrdering - When similar matrices are factored, this
212:    causes the ordering computed in the first factor to be used for all
213:    following factors.

215:    Logically Collective on PC

217:    Input Parameters:
218: +  pc - the preconditioner context
219: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

221:    Options Database Key:
222: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

224:    Level: intermediate

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

228: .seealso: PCFactorSetReuseFill()
229: @*/
230: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
231: {

237:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
238:   return(0);
239: }

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

244:    Options Database Keys:
245: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
246: .  -pc_factor_mat_solver_package - Actives PCFactorSetMatSolverPackage() to choose the direct solver, like superlu
247: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
248: .  -pc_factor_fill <fill> - Sets fill amount
249: .  -pc_factor_in_place - Activates in-place factorization
250: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

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

254:    Level: beginner

256:    Concepts: Cholesky factorization, direct solver

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

262: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
263:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
264:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
265:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

267: M*/

269: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
270: {
272:   PC_Cholesky    *dir;

275:   PetscNewLog(pc,&dir);
276:   pc->data = (void*)dir;
277:   PCFactorInitialize(pc);

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

282:   dir->col = 0;
283:   dir->row = 0;

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

287:   pc->ops->destroy           = PCDestroy_Cholesky;
288:   pc->ops->reset             = PCReset_Cholesky;
289:   pc->ops->apply             = PCApply_Cholesky;
290:   pc->ops->applytranspose    = PCApplyTranspose_Cholesky;
291:   pc->ops->setup             = PCSetUp_Cholesky;
292:   pc->ops->setfromoptions    = PCSetFromOptions_Cholesky;
293:   pc->ops->view              = PCView_Cholesky;
294:   pc->ops->applyrichardson   = 0;
295:   return(0);
296: }