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 CHKERRQ(MatSetErrorIfFailure(pc->pmat,pc->erroriffailure)); 20 if (dir->hdr.inplace) { 21 MatFactorType ftype; 22 23 CHKERRQ(MatGetFactorType(pc->pmat, &ftype)); 24 if (ftype == MAT_FACTOR_NONE) { 25 CHKERRQ(MatQRFactor(pc->pmat,dir->col,&((PC_Factor*)dir)->info)); 26 CHKERRQ(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 CHKERRQ(MatGetFactor(pc->pmat,((PC_Factor*)dir)->solvertype,MAT_FACTOR_QR,&((PC_Factor*)dir)->fact)); 39 CHKERRQ(PetscLogObjectParent((PetscObject)pc,(PetscObject)((PC_Factor*)dir)->fact)); 40 } 41 CHKERRQ(MatQRFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->col,&((PC_Factor*)dir)->info)); 42 CHKERRQ(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 CHKERRQ(MatQRFactorSymbolic(((PC_Factor*)dir)->fact,pc->pmat,dir->col,&((PC_Factor*)dir)->info)); 46 CHKERRQ(MatGetInfo(((PC_Factor*)dir)->fact,MAT_LOCAL,&info)); 47 dir->hdr.actualfill = info.fill_ratio_needed; 48 } else { 49 CHKERRQ(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 50 } 51 CHKERRQ(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 52 if (err) { /* FactorSymbolic() fails */ 53 pc->failedreason = (PCFailedReason)err; 54 PetscFunctionReturn(0); 55 } 56 57 CHKERRQ(MatQRFactorNumeric(((PC_Factor*)dir)->fact,pc->pmat,&((PC_Factor*)dir)->info)); 58 CHKERRQ(MatFactorGetError(((PC_Factor*)dir)->fact,&err)); 59 if (err) { /* FactorNumeric() fails */ 60 pc->failedreason = (PCFailedReason)err; 61 } 62 } 63 64 CHKERRQ(PCFactorGetMatSolverType(pc,&stype)); 65 if (!stype) { 66 MatSolverType solverpackage; 67 CHKERRQ(MatFactorGetSolverType(((PC_Factor*)dir)->fact,&solverpackage)); 68 CHKERRQ(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) CHKERRQ(MatDestroy(&((PC_Factor*)dir)->fact)); 79 CHKERRQ(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 CHKERRQ(PCReset_QR(pc)); 89 CHKERRQ(PetscFree(((PC_Factor*)dir)->ordering)); 90 CHKERRQ(PetscFree(((PC_Factor*)dir)->solvertype)); 91 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscNewLog(pc,&dir)); 153 pc->data = (void*)dir; 154 CHKERRQ(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