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