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