Actual source code: cholesky.c

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

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

 59:     if (!pc->setupcalled) {
 60:       MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 61:       /* check if dir->row == dir->col */
 62:       ISEqual(dir->row,dir->col,&flg);
 63:       if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
 64:       ISDestroy(&dir->col); /* only pass one ordering into CholeskyFactor */

 66:       flg  = PETSC_FALSE;
 67:       PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 68:       if (flg) {
 69:         PetscReal tol = 1.e-10;
 70:         PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 71:         MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 72:       }
 73:       if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 74:       if (!((PC_Factor*)dir)->fact) {
 75:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 76:       }
 77:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 78:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
 79:       dir->hdr.actualfill = info.fill_ratio_needed;
 80:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 81:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 82:       if (!dir->hdr.reuseordering) {
 83:         ISDestroy(&dir->row);
 84:         MatGetOrdering(pc->pmat,((PC_Factor*)dir)->ordering,&dir->row,&dir->col);
 85:         ISDestroy(&dir->col); /* only use dir->row ordering in CholeskyFactor */

 87:         flg  = PETSC_FALSE;
 88:         PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&flg,NULL);
 89:         if (flg) {
 90:           PetscReal tol = 1.e-10;
 91:           PetscOptionsGetReal(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_factor_nonzeros_along_diagonal",&tol,NULL);
 92:           MatReorderForNonzeroDiagonal(pc->pmat,tol,dir->row,dir->row);
 93:         }
 94:         if (dir->row) {PetscLogObjectParent((PetscObject)pc,(PetscObject)dir->row);}
 95:       }
 96:       MatDestroy(&((PC_Factor*)dir)->fact);
 97:       MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_CHOLESKY,&((PC_Factor*)dir)->fact);
 98:       MatCholeskyFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->row,&((PC_Factor*)dir)->info);
 99:       MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info);
100:       dir->hdr.actualfill = info.fill_ratio_needed;
101:       PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
102:     } else {
103:       MatFactorGetError(((PC_Factor*)dir)->fact,&err);
104:       if (err == MAT_FACTOR_NUMERIC_ZEROPIVOT) {
105:         MatFactorClearError(((PC_Factor*)dir)->fact);
106:         pc->failedreason = PC_NOERROR;
107:       }
108:     }
109:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
110:     if (err) { /* FactorSymbolic() fails */
111:       pc->failedreason = (PCFailedReason)err;
112:       return(0);
113:     }

115:     MatCholeskyFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info);
116:     MatFactorGetError(((PC_Factor*)dir)->fact,&err);
117:     if (err) { /* FactorNumeric() fails */
118:       pc->failedreason = (PCFailedReason)err;
119:     }
120:   }

122:   PCFactorGetMatSolverType(pc,&stype);
123:   if (!stype) {
124:     MatSolverType solverpackage;
125:     MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage);
126:     PCFactorSetMatSolverType(pc,solverpackage);
127:   }
128:   return(0);
129: }

131: static PetscErrorCode PCReset_Cholesky(PC pc)
132: {
133:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

137:   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
138:   ISDestroy(&dir->row);
139:   ISDestroy(&dir->col);
140:   return(0);
141: }

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

149:   PCReset_Cholesky(pc);
150:   PetscFree(((PC_Factor*)dir)->ordering);
151:   PetscFree(((PC_Factor*)dir)->solvertype);
152:   PetscFree(pc->data);
153:   return(0);
154: }

156: static PetscErrorCode PCApply_Cholesky(PC pc,Vec x,Vec y)
157: {
158:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

162:   if (dir->hdr.inplace) {
163:     MatSolve(pc->pmat,x,y);
164:   } else {
165:     MatSolve(((PC_Factor*)dir)->fact,x,y);
166:   }
167:   return(0);
168: }

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

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

184: static PetscErrorCode PCApplySymmetricRight_Cholesky(PC pc,Vec x,Vec y)
185: {
186:   PC_Cholesky    *dir = (PC_Cholesky*)pc->data;

190:   if (dir->hdr.inplace) {
191:     MatBackwardSolve(pc->pmat,x,y);
192:   } else {
193:     MatBackwardSolve(((PC_Factor*)dir)->fact,x,y);
194:   }
195:   return(0);
196: }

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

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

212: /* -----------------------------------------------------------------------------------*/

214: /* -----------------------------------------------------------------------------------*/

216: /*@
217:    PCFactorSetReuseOrdering - When similar matrices are factored, this
218:    causes the ordering computed in the first factor to be used for all
219:    following factors.

221:    Logically Collective on PC

223:    Input Parameters:
224: +  pc - the preconditioner context
225: -  flag - PETSC_TRUE to reuse else PETSC_FALSE

227:    Options Database Key:
228: .  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()

230:    Level: intermediate

232: .seealso: PCFactorSetReuseFill()
233: @*/
234: PetscErrorCode  PCFactorSetReuseOrdering(PC pc,PetscBool flag)
235: {

241:   PetscTryMethod(pc,"PCFactorSetReuseOrdering_C",(PC,PetscBool),(pc,flag));
242:   return(0);
243: }

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

248:    Options Database Keys:
249: +  -pc_factor_reuse_ordering - Activate PCFactorSetReuseOrdering()
250: .  -pc_factor_mat_solver_type - Actives PCFactorSetMatSolverType() to choose the direct solver, like superlu
251: .  -pc_factor_reuse_fill - Activates PCFactorSetReuseFill()
252: .  -pc_factor_fill <fill> - Sets fill amount
253: .  -pc_factor_in_place - Activates in-place factorization
254: -  -pc_factor_mat_ordering_type <nd,rcm,...> - Sets ordering routine

256:    Notes:
257:     Not all options work for all matrix formats

259:    Level: beginner

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

266: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
267:            PCILU, PCLU, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
268:            PCFactorSetFill(), PCFactorSetShiftNonzero(), PCFactorSetShiftType(), PCFactorSetShiftAmount()
269:            PCFactorSetUseInPlace(), PCFactorGetUseInPlace(), PCFactorSetMatOrderingType()

271: M*/

273: PETSC_EXTERN PetscErrorCode PCCreate_Cholesky(PC pc)
274: {
276:   PC_Cholesky    *dir;

279:   PetscNewLog(pc,&dir);
280:   pc->data = (void*)dir;
281:   PCFactorInitialize(pc);

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

286:   dir->col = 0;
287:   dir->row = 0;

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

291:   pc->ops->destroy             = PCDestroy_Cholesky;
292:   pc->ops->reset               = PCReset_Cholesky;
293:   pc->ops->apply               = PCApply_Cholesky;
294:   pc->ops->applysymmetricleft  = PCApplySymmetricLeft_Cholesky;
295:   pc->ops->applysymmetricright = PCApplySymmetricRight_Cholesky;
296:   pc->ops->applytranspose      = PCApplyTranspose_Cholesky;
297:   pc->ops->setup               = PCSetUp_Cholesky;
298:   pc->ops->setfromoptions      = PCSetFromOptions_Cholesky;
299:   pc->ops->view                = PCView_Factor;
300:   pc->ops->applyrichardson     = 0;
301:   return(0);
302: }