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:   PetscErrorCode         ierr;
 12:   PC_QR                  *dir = (PC_QR*)pc->data;
 13:   MatSolverType          stype;
 14:   MatFactorError         err;

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

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

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

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

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

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

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

 80:   if (!dir->hdr.inplace && ((PC_Factor*)dir)->fact) {MatDestroy(&((PC_Factor*)dir)->fact);}
 81:   ISDestroy(&dir->col);
 82:   return(0);
 83: }

 85: static PetscErrorCode PCDestroy_QR(PC pc)
 86: {
 87:   PC_QR          *dir = (PC_QR*)pc->data;

 91:   PCReset_QR(pc);
 92:   PetscFree(((PC_Factor*)dir)->ordering);
 93:   PetscFree(((PC_Factor*)dir)->solvertype);
 94:   PetscFree(pc->data);
 95:   return(0);
 96: }

 98: static PetscErrorCode PCApply_QR(PC pc,Vec x,Vec y)
 99: {
100:   PC_QR          *dir = (PC_QR*)pc->data;
101:   Mat            fact;

105:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
106:   MatSolve(fact,x,y);
107:   return(0);
108: }

110: static PetscErrorCode PCMatApply_QR(PC pc,Mat X,Mat Y)
111: {
112:   PC_QR          *dir = (PC_QR*)pc->data;
113:   Mat            fact;

117:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
118:   MatMatSolve(fact,X,Y);
119:   return(0);
120: }

122: static PetscErrorCode PCApplyTranspose_QR(PC pc,Vec x,Vec y)
123: {
124:   PC_QR          *dir = (PC_QR*)pc->data;
125:   Mat            fact;

129:   fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact;
130:   MatSolveTranspose(fact,x,y);
131:   return(0);
132: }

134: /* -----------------------------------------------------------------------------------*/

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

139:    Level: beginner

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

146: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
147:            PCILU, PCLU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
148:            PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
149:            PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
150:            PCFactorReorderForNonzeroDiagonal()
151: M*/

153: PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc)
154: {
156:   PC_QR          *dir;

159:   PetscNewLog(pc,&dir);
160:   pc->data = (void*)dir;
161:   PCFactorInitialize(pc, MAT_FACTOR_QR);

163:   dir->col                   = NULL;
164:   pc->ops->reset             = PCReset_QR;
165:   pc->ops->destroy           = PCDestroy_QR;
166:   pc->ops->apply             = PCApply_QR;
167:   pc->ops->matapply          = PCMatApply_QR;
168:   pc->ops->applytranspose    = PCApplyTranspose_QR;
169:   pc->ops->setup             = PCSetUp_QR;
170:   pc->ops->setfromoptions    = PCSetFromOptions_Factor;
171:   pc->ops->view              = PCView_Factor;
172:   pc->ops->applyrichardson   = NULL;
173:   return(0);
174: }