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;

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

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

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

 35:     if (!pc->setupcalled) {
 36:       if (!((PC_Factor*)dir)->fact) {
 37:         MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_QR,&((PC_Factor*)dir)->fact);
 38:         PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact);
 39:       }
 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:   PetscFree(pc->data);
 89:   return 0;
 90: }

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

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

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

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

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

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

122: /* -----------------------------------------------------------------------------------*/

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

127:    Level: beginner

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

134: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
135:            PCILU, PCLU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(),
136:            PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(),
137:            PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount()
138:            PCFactorReorderForNonzeroDiagonal()
139: M*/

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

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

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