Actual source code: qr.c


  2: /*
  3:    Defines a direct QR factorization preconditioner for any Mat implementation
  4:    Note: this need not be considered a preconditioner since it supplies
  5:          a direct solver.
  6: */
  7: #include <../src/ksp/pc/impls/factor/qr/qr.h>

  9: static PetscErrorCode PCSetUp_QR(PC pc)
 10: {
 11:   PC_QR         *dir = (PC_QR *)pc->data;
 12:   MatSolverType  stype;
 13:   MatFactorError err;
 14:   const char    *prefix;

 16:   PCGetOptionsPrefix(pc, &prefix);
 17:   MatSetOptionsPrefix(pc->pmat, prefix);
 18:   pc->failedreason = PC_NOERROR;
 19:   if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill;

 21:   MatSetErrorIfFailure(pc->pmat, pc->erroriffailure);
 22:   if (dir->hdr.inplace) {
 23:     MatFactorType ftype;

 25:     MatGetFactorType(pc->pmat, &ftype);
 26:     if (ftype == MAT_FACTOR_NONE) {
 27:       MatQRFactor(pc->pmat, dir->col, &((PC_Factor *)dir)->info);
 28:       MatFactorGetError(pc->pmat, &err);
 29:       if (err) { /* Factor() fails */
 30:         pc->failedreason = (PCFailedReason)err;
 31:         return 0;
 32:       }
 33:     }
 34:     ((PC_Factor *)dir)->fact = pc->pmat;
 35:   } else {
 36:     MatInfo info;

 38:     if (!pc->setupcalled) {
 39:       if (!((PC_Factor *)dir)->fact) { MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact); }
 40:       MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info);
 41:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
 42:       dir->hdr.actualfill = info.fill_ratio_needed;
 43:     } else if (pc->flag != SAME_NONZERO_PATTERN) {
 44:       MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info);
 45:       MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info);
 46:       dir->hdr.actualfill = info.fill_ratio_needed;
 47:     } else {
 48:       MatFactorGetError(((PC_Factor *)dir)->fact, &err);
 49:     }
 50:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
 51:     if (err) { /* FactorSymbolic() fails */
 52:       pc->failedreason = (PCFailedReason)err;
 53:       return 0;
 54:     }

 56:     MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info);
 57:     MatFactorGetError(((PC_Factor *)dir)->fact, &err);
 58:     if (err) { /* FactorNumeric() fails */
 59:       pc->failedreason = (PCFailedReason)err;
 60:     }
 61:   }

 63:   PCFactorGetMatSolverType(pc, &stype);
 64:   if (!stype) {
 65:     MatSolverType solverpackage;
 66:     MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage);
 67:     PCFactorSetMatSolverType(pc, solverpackage);
 68:   }
 69:   return 0;
 70: }

 72: static PetscErrorCode PCReset_QR(PC pc)
 73: {
 74:   PC_QR *dir = (PC_QR *)pc->data;

 76:   if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) MatDestroy(&((PC_Factor *)dir)->fact);
 77:   ISDestroy(&dir->col);
 78:   return 0;
 79: }

 81: static PetscErrorCode PCDestroy_QR(PC pc)
 82: {
 83:   PC_QR *dir = (PC_QR *)pc->data;

 85:   PCReset_QR(pc);
 86:   PetscFree(((PC_Factor *)dir)->ordering);
 87:   PetscFree(((PC_Factor *)dir)->solvertype);
 88:   PCFactorClearComposedFunctions(pc);
 89:   PetscFree(pc->data);
 90:   return 0;
 91: }

 93: static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y)
 94: {
 95:   PC_QR *dir = (PC_QR *)pc->data;
 96:   Mat    fact;

 98:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
 99:   MatSolve(fact, x, y);
100:   return 0;
101: }

103: static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y)
104: {
105:   PC_QR *dir = (PC_QR *)pc->data;
106:   Mat    fact;

108:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
109:   MatMatSolve(fact, X, Y);
110:   return 0;
111: }

113: static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y)
114: {
115:   PC_QR *dir = (PC_QR *)pc->data;
116:   Mat    fact;

118:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact;
119:   MatSolveTranspose(fact, x, y);
120:   return 0;
121: }

123: /*MC
124:    PCQR - Uses a direct solver, based on QR factorization, as a preconditioner

126:    Level: beginner

128:    Note:
129:    Usually this will compute an "exact" solution in one iteration and does
130:    not need a Krylov method (i.e. you can use -ksp_type preonly, or
131:    `KSPSetType`(ksp,`KSPPREONLY`) for the Krylov method

133: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSVD`,
134:           `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`,
135:           `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`,
136:           `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()`
137:           `PCFactorReorderForNonzeroDiagonal()`
138: M*/

140: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
141: {
142:   PC_QR *dir;

144:   PetscNew(&dir);
145:   pc->data = (void *)dir;
146:   PCFactorInitialize(pc, MAT_FACTOR_QR);

148:   dir->col                 = NULL;
149:   pc->ops->reset           = PCReset_QR;
150:   pc->ops->destroy         = PCDestroy_QR;
151:   pc->ops->apply           = PCApply_QR;
152:   pc->ops->matapply        = PCMatApply_QR;
153:   pc->ops->applytranspose  = PCApplyTranspose_QR;
154:   pc->ops->setup           = PCSetUp_QR;
155:   pc->ops->setfromoptions  = PCSetFromOptions_Factor;
156:   pc->ops->view            = PCView_Factor;
157:   pc->ops->applyrichardson = NULL;
158:   return 0;
159: }