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