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 15 PetscFunctionBegin; 16 pc->failedreason = PC_NOERROR; 17 if (dir->hdr.reusefill && pc->setupcalled) ((PC_Factor*)dir)->info.fill = dir->hdr.actualfill; 18 19 PetscCall(MatSetErrorIfFailure(pc->pmat,pc->erroriffailure)); 20 if (dir->hdr.inplace) { 21 MatFactorType ftype; 22 23 PetscCall(MatGetFactorType(pc->pmat, &ftype)); 24 if (ftype == MAT_FACTOR_NONE) { 25 PetscCall(MatQRFactor(pc->pmat,dir->col,&((PC_Factor*)dir)->info)); 26 PetscCall(MatFactorGetError(pc->pmat,&err)); 27 if (err) { /* Factor() fails */ 28 pc->failedreason = (PCFailedReason)err; 29 PetscFunctionReturn(0); 30 } 31 } 32 ((PC_Factor*)dir)->fact = pc->pmat; 33 } else { 34 MatInfo info; 35 36 if (!pc->setupcalled) { 37 if (!((PC_Factor*)dir)->fact) { 38 PetscCall(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_QR,&((PC_Factor*)dir)->fact)); 39 PetscCall(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 40 } 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(0); 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(0); 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(0); 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(PetscFree(pc->data)); 92 PetscFunctionReturn(0); 93 } 94 95 static PetscErrorCode PCApply_QR(PC pc,Vec x,Vec y) 96 { 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 { 108 PC_QR *dir = (PC_QR*)pc->data; 109 Mat fact; 110 111 PetscFunctionBegin; 112 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact; 113 PetscCall(MatMatSolve(fact,X,Y)); 114 PetscFunctionReturn(0); 115 } 116 117 static PetscErrorCode PCApplyTranspose_QR(PC pc,Vec x,Vec y) 118 { 119 PC_QR *dir = (PC_QR*)pc->data; 120 Mat fact; 121 122 PetscFunctionBegin; 123 fact = dir->hdr.inplace ? pc->pmat : ((PC_Factor*)dir)->fact; 124 PetscCall(MatSolveTranspose(fact,x,y)); 125 PetscFunctionReturn(0); 126 } 127 128 /* -----------------------------------------------------------------------------------*/ 129 130 /*MC 131 PCQR - Uses a direct solver, based on QR factorization, as a preconditioner 132 133 Level: beginner 134 135 Notes: 136 Usually this will compute an "exact" solution in one iteration and does 137 not need a Krylov method (i.e. you can use -ksp_type preonly, or 138 KSPSetType(ksp,KSPPREONLY) for the Krylov method 139 140 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 141 PCILU, PCLU, PCCHOLESKY, PCICC, PCFactorSetReuseOrdering(), PCFactorSetReuseFill(), PCFactorGetMatrix(), 142 PCFactorSetFill(), PCFactorSetUseInPlace(), PCFactorSetMatOrderingType(), PCFactorSetColumnPivot(), 143 PCFactorSetPivotingInBlocks(),PCFactorSetShiftType(),PCFactorSetShiftAmount() 144 PCFactorReorderForNonzeroDiagonal() 145 M*/ 146 147 PETSC_EXTERN PetscErrorCode PCCreate_QR(PC pc) 148 { 149 PC_QR *dir; 150 151 PetscFunctionBegin; 152 PetscCall(PetscNewLog(pc,&dir)); 153 pc->data = (void*)dir; 154 PetscCall(PCFactorInitialize(pc, MAT_FACTOR_QR)); 155 156 dir->col = NULL; 157 pc->ops->reset = PCReset_QR; 158 pc->ops->destroy = PCDestroy_QR; 159 pc->ops->apply = PCApply_QR; 160 pc->ops->matapply = PCMatApply_QR; 161 pc->ops->applytranspose = PCApplyTranspose_QR; 162 pc->ops->setup = PCSetUp_QR; 163 pc->ops->setfromoptions = PCSetFromOptions_Factor; 164 pc->ops->view = PCView_Factor; 165 pc->ops->applyrichardson = NULL; 166 PetscFunctionReturn(0); 167 } 168