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(PetscFree(pc->data)); 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode PCApply_QR(PC pc,Vec x,Vec y) 99 { 100 PC_QR *dir = (PC_QR*)pc->data; 101 Mat fact; 102 103 PetscFunctionBegin; 104 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact; 105 PetscCall(MatSolve(fact,x,y)); 106 PetscFunctionReturn(0); 107 } 108 109 static PetscErrorCode PCMatApply_QR(PC pc,Mat X,Mat Y) 110 { 111 PC_QR *dir = (PC_QR*)pc->data; 112 Mat fact; 113 114 PetscFunctionBegin; 115 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact; 116 PetscCall(MatMatSolve(fact,X,Y)); 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode PCApplyTranspose_QR(PC pc,Vec x,Vec y) 121 { 122 PC_QR *dir = (PC_QR*)pc->data; 123 Mat fact; 124 125 PetscFunctionBegin; 126 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact; 127 PetscCall(MatSolveTranspose(fact,x,y)); 128 PetscFunctionReturn(0); 129 } 130 131 /* -----------------------------------------------------------------------------------*/ 132 133 /*MC 134 PCQR - Uses a direct solver, based on QR factorization, as a preconditioner 135 136 Level: beginner 137 138 Notes: 139 Usually this will compute an "exact" solution in one iteration and does 140 not need a Krylov method (i.e. you can use -ksp_type preonly, or 141 KSPSetType(ksp,KSPPREONLY) for the Krylov method 142 143 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 144 `PCILU`, `PCLU`, `PCCHOLESKY`, `PCICC`, `PCFactorSetReuseOrdering()`, `PCFactorSetReuseFill()`, `PCFactorGetMatrix()`, 145 `PCFactorSetFill()`, `PCFactorSetUseInPlace()`, `PCFactorSetMatOrderingType()`, `PCFactorSetColumnPivot()`, 146 `PCFactorSetPivotingInBlocks()`, `PCFactorSetShiftType()`, `PCFactorSetShiftAmount()` 147 `PCFactorReorderForNonzeroDiagonal()` 148 M*/ 149 150 PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc) 151 { 152 PC_QR *dir; 153 154 PetscFunctionBegin; 155 PetscCall(PetscNewLog(pc,&dir)); 156 pc->data = (void*)dir; 157 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR)); 158 159 dir->col = NULL; 160 pc->ops->reset = PCReset_QR; 161 pc->ops->destroy = PCDestroy_QR; 162 pc->ops->apply = PCApply_QR; 163 pc->ops->matapply = PCMatApply_QR; 164 pc->ops->applytranspose = PCApplyTranspose_QR; 165 pc->ops->setup = PCSetUp_QR; 166 pc->ops->setfromoptions = PCSetFromOptions_Factor; 167 pc->ops->view = PCView_Factor; 168 pc->ops->applyrichardson = NULL; 169 PetscFunctionReturn(0); 170 } 171