1 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> /*I "petscpc.h" I*/ 8 9 static PetscErrorCode PCSetUp_QR(PC pc) { 10 PC_QR *dir = (PC_QR *)pc->data; 11 MatSolverType stype; 12 MatFactorError err; 13 const char *prefix; 14 15 PetscFunctionBegin; 16 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 17 PetscCall(MatSetOptionsPrefix(pc->pmat, prefix)); 18 pc->failedreason = PC_NOERROR; 19 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor *)dir)->info.fill = dir->hdr.actualfill; 20 21 PetscCall(MatSetErrorIfFailure(pc->pmat, pc->erroriffailure)); 22 if (dir->hdr.inplace) { 23 MatFactorType ftype; 24 25 PetscCall(MatGetFactorType(pc->pmat, &ftype)); 26 if (ftype == MAT_FACTOR_NONE) { 27 PetscCall(MatQRFactor(pc->pmat, dir->col, &((PC_Factor *)dir)->info)); 28 PetscCall(MatFactorGetError(pc->pmat, &err)); 29 if (err) { /* Factor() fails */ 30 pc->failedreason = (PCFailedReason)err; 31 PetscFunctionReturn(0); 32 } 33 } 34 ((PC_Factor *)dir)->fact = pc->pmat; 35 } else { 36 MatInfo info; 37 38 if (!pc->setupcalled) { 39 if (!((PC_Factor *)dir)->fact) { 40 PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact)); 41 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)((PC_Factor *)dir)->fact)); 42 } 43 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info)); 44 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 45 dir->hdr.actualfill = info.fill_ratio_needed; 46 } else if (pc->flag != SAME_NONZERO_PATTERN) { 47 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info)); 48 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 49 dir->hdr.actualfill = info.fill_ratio_needed; 50 } else { 51 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 52 } 53 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 54 if (err) { /* FactorSymbolic() fails */ 55 pc->failedreason = (PCFailedReason)err; 56 PetscFunctionReturn(0); 57 } 58 59 PetscCall(MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info)); 60 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 61 if (err) { /* FactorNumeric() fails */ 62 pc->failedreason = (PCFailedReason)err; 63 } 64 } 65 66 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 67 if (!stype) { 68 MatSolverType solverpackage; 69 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage)); 70 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 71 } 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode PCReset_QR(PC pc) { 76 PC_QR *dir = (PC_QR *)pc->data; 77 78 PetscFunctionBegin; 79 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 80 PetscCall(ISDestroy(&dir->col)); 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode PCDestroy_QR(PC pc) { 85 PC_QR *dir = (PC_QR *)pc->data; 86 87 PetscFunctionBegin; 88 PetscCall(PCReset_QR(pc)); 89 PetscCall(PetscFree(((PC_Factor *)dir)->ordering)); 90 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype)); 91 PetscCall(PCFactorClearComposedFunctions(pc)); 92 PetscCall(PetscFree(pc->data)); 93 PetscFunctionReturn(0); 94 } 95 96 static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y) { 97 PC_QR *dir = (PC_QR *)pc->data; 98 Mat fact; 99 100 PetscFunctionBegin; 101 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 102 PetscCall(MatSolve(fact, x, y)); 103 PetscFunctionReturn(0); 104 } 105 106 static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y) { 107 PC_QR *dir = (PC_QR *)pc->data; 108 Mat fact; 109 110 PetscFunctionBegin; 111 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 112 PetscCall(MatMatSolve(fact, X, Y)); 113 PetscFunctionReturn(0); 114 } 115 116 static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y) { 117 PC_QR *dir = (PC_QR *)pc->data; 118 Mat fact; 119 120 PetscFunctionBegin; 121 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 122 PetscCall(MatSolveTranspose(fact, x, y)); 123 PetscFunctionReturn(0); 124 } 125 126 /* -----------------------------------------------------------------------------------*/ 127 128 /*MC 129 PCQR - Uses a direct solver, based on QR factorization, as a preconditioner 130 131 Level: beginner 132 133 Notes: 134 Usually this will compute an "exact" solution in one iteration and does 135 not need a Krylov method (i.e. you can use -ksp_type preonly, or 136 KSPSetType(ksp,KSPPREONLY) for the Krylov method 137 138 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 139 `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 140 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`, 141 `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 142 `PCFactorReorderForNonzeroDiagonal()` 143 M*/ 144 145 PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc) { 146 PC_QR *dir; 147 148 PetscFunctionBegin; 149 PetscCall(PetscNewLog(pc, &dir)); 150 pc->data = (void *)dir; 151 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR)); 152 153 dir->col = NULL; 154 pc->ops->reset = PCReset_QR; 155 pc->ops->destroy = PCDestroy_QR; 156 pc->ops->apply = PCApply_QR; 157 pc->ops->matapply = PCMatApply_QR; 158 pc->ops->applytranspose = PCApplyTranspose_QR; 159 pc->ops->setup = PCSetUp_QR; 160 pc->ops->setfromoptions = PCSetFromOptions_Factor; 161 pc->ops->view = PCView_Factor; 162 pc->ops->applyrichardson = NULL; 163 PetscFunctionReturn(0); 164 } 165