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) { PetscCall(MatGetFactor(pc->pmat, ((PC_Factor *)dir)->solvertype, MAT_FACTOR_QR, &((PC_Factor *)dir)->fact)); } 40 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info)); 41 PetscCall(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 PetscCall(MatQRFactorSymbolic(((PC_Factor *)dir)->fact, pc->pmat, dir->col, &((PC_Factor *)dir)->info)); 45 PetscCall(MatGetInfo(((PC_Factor *)dir)->fact, MAT_LOCAL, &info)); 46 dir->hdr.actualfill = info.fill_ratio_needed; 47 } else { 48 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 49 } 50 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 51 if (err) { /* FactorSymbolic() fails */ 52 pc->failedreason = (PCFailedReason)err; 53 PetscFunctionReturn(0); 54 } 55 56 PetscCall(MatQRFactorNumeric(((PC_Factor *)dir)->fact, pc->pmat, &((PC_Factor *)dir)->info)); 57 PetscCall(MatFactorGetError(((PC_Factor *)dir)->fact, &err)); 58 if (err) { /* FactorNumeric() fails */ 59 pc->failedreason = (PCFailedReason)err; 60 } 61 } 62 63 PetscCall(PCFactorGetMatSolverType(pc, &stype)); 64 if (!stype) { 65 MatSolverType solverpackage; 66 PetscCall(MatFactorGetSolverType(((PC_Factor *)dir)->fact, &solverpackage)); 67 PetscCall(PCFactorSetMatSolverType(pc, solverpackage)); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode PCReset_QR(PC pc) { 73 PC_QR *dir = (PC_QR *)pc->data; 74 75 PetscFunctionBegin; 76 if (!dir->hdr.inplace && ((PC_Factor *)dir)->fact) PetscCall(MatDestroy(&((PC_Factor *)dir)->fact)); 77 PetscCall(ISDestroy(&dir->col)); 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode PCDestroy_QR(PC pc) { 82 PC_QR *dir = (PC_QR *)pc->data; 83 84 PetscFunctionBegin; 85 PetscCall(PCReset_QR(pc)); 86 PetscCall(PetscFree(((PC_Factor *)dir)->ordering)); 87 PetscCall(PetscFree(((PC_Factor *)dir)->solvertype)); 88 PetscCall(PCFactorClearComposedFunctions(pc)); 89 PetscCall(PetscFree(pc->data)); 90 PetscFunctionReturn(0); 91 } 92 93 static PetscErrorCode PCApply_QR(PC pc, Vec x, Vec y) { 94 PC_QR *dir = (PC_QR *)pc->data; 95 Mat fact; 96 97 PetscFunctionBegin; 98 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 99 PetscCall(MatSolve(fact, x, y)); 100 PetscFunctionReturn(0); 101 } 102 103 static PetscErrorCode PCMatApply_QR(PC pc, Mat X, Mat Y) { 104 PC_QR *dir = (PC_QR *)pc->data; 105 Mat fact; 106 107 PetscFunctionBegin; 108 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 109 PetscCall(MatMatSolve(fact, X, Y)); 110 PetscFunctionReturn(0); 111 } 112 113 static PetscErrorCode PCApplyTranspose_QR(PC pc, Vec x, Vec y) { 114 PC_QR *dir = (PC_QR *)pc->data; 115 Mat fact; 116 117 PetscFunctionBegin; 118 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor *)dir)->fact; 119 PetscCall(MatSolveTranspose(fact, x, y)); 120 PetscFunctionReturn(0); 121 } 122 123 /*MC 124 PCQR - Uses a direct solver, based on QR factorization, as a preconditioner 125 126 Level: beginner 127 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 132 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*/ 139 140 PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc) { 141 PC_QR *dir; 142 143 PetscFunctionBegin; 144 PetscCall(PetscNew(&dir)); 145 pc->data = (void *)dir; 146 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR)); 147 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 PetscFunctionReturn(0); 159 } 160