1 2 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 3 4 typedef enum { 5 PCMATOP_MULT, 6 PCMATOP_MULT_TRANSPOSE, 7 PCMATOP_MULT_HERMITIAN_TRANSPOSE, 8 PCMATOP_SOLVE, 9 PCMATOP_SOLVE_TRANSPOSE, 10 } PCMatOperation; 11 12 const char *const PCMatOpTypes[] = {"Mult", "MultTranspose", "MultHermitianTranspose", "Solve", "SolveTranspose", NULL}; 13 14 typedef struct _PCMAT { 15 PCMatOperation apply; 16 } PC_Mat; 17 18 static PetscErrorCode PCApply_Mat(PC pc, Vec x, Vec y) 19 { 20 PC_Mat *pcmat = (PC_Mat *)pc->data; 21 22 PetscFunctionBegin; 23 switch (pcmat->apply) { 24 case PCMATOP_MULT: 25 PetscCall(MatMult(pc->pmat, x, y)); 26 break; 27 case PCMATOP_MULT_TRANSPOSE: 28 PetscCall(MatMultTranspose(pc->pmat, x, y)); 29 break; 30 case PCMATOP_SOLVE: 31 PetscCall(MatSolve(pc->pmat, x, y)); 32 break; 33 case PCMATOP_SOLVE_TRANSPOSE: 34 PetscCall(MatSolveTranspose(pc->pmat, x, y)); 35 break; 36 case PCMATOP_MULT_HERMITIAN_TRANSPOSE: 37 PetscCall(MatMultHermitianTranspose(pc->pmat, x, y)); 38 break; 39 } 40 PetscFunctionReturn(PETSC_SUCCESS); 41 } 42 43 static PetscErrorCode PCMatApply_Mat(PC pc, Mat X, Mat Y) 44 { 45 PC_Mat *pcmat = (PC_Mat *)pc->data; 46 47 PetscFunctionBegin; 48 switch (pcmat->apply) { 49 case PCMATOP_MULT: 50 PetscCall(MatMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y)); 51 break; 52 case PCMATOP_MULT_TRANSPOSE: 53 PetscCall(MatTransposeMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y)); 54 break; 55 case PCMATOP_SOLVE: 56 PetscCall(MatMatSolve(pc->pmat, X, Y)); 57 break; 58 case PCMATOP_SOLVE_TRANSPOSE: 59 PetscCall(MatMatSolveTranspose(pc->pmat, X, Y)); 60 break; 61 case PCMATOP_MULT_HERMITIAN_TRANSPOSE: { 62 Mat W; 63 64 PetscCall(MatDuplicate(X, MAT_COPY_VALUES, &W)); 65 PetscCall(MatConjugate(W)); 66 PetscCall(MatTransposeMatMult(pc->pmat, W, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y)); 67 PetscCall(MatConjugate(Y)); 68 PetscCall(MatDestroy(&W)); 69 break; 70 } 71 } 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 75 static PetscErrorCode PCApplyTranspose_Mat(PC pc, Vec x, Vec y) 76 { 77 PC_Mat *pcmat = (PC_Mat *)pc->data; 78 79 PetscFunctionBegin; 80 switch (pcmat->apply) { 81 case PCMATOP_MULT: 82 PetscCall(MatMultTranspose(pc->pmat, x, y)); 83 break; 84 case PCMATOP_MULT_TRANSPOSE: 85 PetscCall(MatMult(pc->pmat, x, y)); 86 break; 87 case PCMATOP_SOLVE: 88 PetscCall(MatSolveTranspose(pc->pmat, x, y)); 89 break; 90 case PCMATOP_SOLVE_TRANSPOSE: 91 PetscCall(MatSolve(pc->pmat, x, y)); 92 break; 93 case PCMATOP_MULT_HERMITIAN_TRANSPOSE: { 94 Vec w; 95 96 PetscCall(VecDuplicate(x, &w)); 97 PetscCall(VecCopy(x, w)); 98 PetscCall(VecConjugate(w)); 99 PetscCall(MatMult(pc->pmat, w, y)); 100 PetscCall(VecConjugate(y)); 101 PetscCall(VecDestroy(&w)); 102 break; 103 } 104 } 105 PetscFunctionReturn(PETSC_SUCCESS); 106 } 107 108 static PetscErrorCode PCDestroy_Mat(PC pc) 109 { 110 PetscFunctionBegin; 111 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", NULL)); 112 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", NULL)); 113 PetscCall(PetscFree(pc->data)); 114 PetscFunctionReturn(PETSC_SUCCESS); 115 } 116 117 /*@ 118 PCMatSetApplyOperation - Set which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`. 119 120 Logically collective 121 122 Input Parameters: 123 + pc - An instance of `PCMAT` 124 - matop - The selected `MatOperation` 125 126 Level: intermediate 127 128 Note: 129 If you have a matrix type that implements an exact inverse that isn't a factorization, 130 you can use `PCMatSetApplyOperation(pc, MATOP_SOLVE)`. 131 132 .seealso: `PCMAT`, `PCMatGetApplyOperation()`, `PCApply()`, `MatOperation` 133 @*/ 134 PetscErrorCode PCMatSetApplyOperation(PC pc, MatOperation matop) 135 { 136 PetscFunctionBegin; 137 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 138 PetscTryMethod((PetscObject)pc, "PCMatSetApplyOperation_C", (PC, MatOperation), (pc, matop)); 139 PetscFunctionReturn(PETSC_SUCCESS); 140 } 141 142 /*@ 143 PCMatGetApplyOperation - Get which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`. 144 145 Logically collective 146 147 Input Parameter: 148 . pc - An instance of `PCMAT` 149 150 Output Parameter: 151 . matop - The `MatOperation` 152 153 Level: intermediate 154 155 .seealso: `PCMAT`, `PCMatSetApplyOperation()`, `PCApply()`, `MatOperation` 156 @*/ 157 PetscErrorCode PCMatGetApplyOperation(PC pc, MatOperation *matop) 158 { 159 PetscFunctionBegin; 160 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 161 PetscValidPointer(matop, 2); 162 PetscUseMethod((PetscObject)pc, "PCMatGetApplyOperation_C", (PC, MatOperation *), (pc, matop)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 static PetscErrorCode PCMatSetApplyOperation_Mat(PC pc, MatOperation matop) 167 { 168 PC_Mat *pcmat = (PC_Mat *)pc->data; 169 PCMatOperation pcmatop; 170 171 PetscFunctionBegin; 172 173 // clang-format off 174 #define MATOP_TO_PCMATOP_CASE(var, OP) case MATOP_##OP: (var) = PCMATOP_##OP; break 175 switch (matop) { 176 MATOP_TO_PCMATOP_CASE(pcmatop, MULT); 177 MATOP_TO_PCMATOP_CASE(pcmatop, MULT_TRANSPOSE); 178 MATOP_TO_PCMATOP_CASE(pcmatop, MULT_HERMITIAN_TRANSPOSE); 179 MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE); 180 MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE_TRANSPOSE); 181 default: 182 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported MatOperation %d for PCMatSetApplyOperation()", (int)matop); 183 } 184 #undef MATOP_TO_PCMATOP_CASE 185 // clang-format on 186 187 pcmat->apply = pcmatop; 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p) 192 { 193 PC_Mat *pcmat = (PC_Mat *)pc->data; 194 MatOperation matop = MATOP_MULT; 195 196 PetscFunctionBegin; 197 198 // clang-format off 199 #define PCMATOP_TO_MATOP_CASE(var, OP) case PCMATOP_##OP: (var) = MATOP_##OP; break 200 switch (pcmat->apply) { 201 PCMATOP_TO_MATOP_CASE(matop, MULT); 202 PCMATOP_TO_MATOP_CASE(matop, MULT_TRANSPOSE); 203 PCMATOP_TO_MATOP_CASE(matop, MULT_HERMITIAN_TRANSPOSE); 204 PCMATOP_TO_MATOP_CASE(matop, SOLVE); 205 PCMATOP_TO_MATOP_CASE(matop, SOLVE_TRANSPOSE); 206 } 207 #undef PCMATOP_TO_MATOP_CASE 208 // clang-format on 209 210 *matop_p = matop; 211 PetscFunctionReturn(PETSC_SUCCESS); 212 } 213 214 static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer) 215 { 216 PC_Mat *pcmat = (PC_Mat *)pc->data; 217 PetscBool iascii; 218 219 PetscFunctionBegin; 220 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 221 if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply])); } 222 PetscFunctionReturn(PETSC_SUCCESS); 223 } 224 225 /*MC 226 PCMAT - A preconditioner obtained by applying an operation of the `pmat` provided in 227 in `PCSetOperators(pc,amat,pmat)` or `KSPSetOperators(ksp,amat,pmat)`. By default, the operation is `MATOP_MULT`, 228 meaning that the `pmat` provides an approximate inverse of `amat`. 229 If some other operation of `pmat` implements the approximate inverse, 230 use `PCMatSetApplyOperation()` to select that operation. 231 232 Level: intermediate 233 234 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()` 235 M*/ 236 237 PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc) 238 { 239 PC_Mat *data; 240 241 PetscFunctionBegin; 242 PetscCall(PetscNew(&data)); 243 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", PCMatSetApplyOperation_Mat)); 244 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", PCMatGetApplyOperation_Mat)); 245 data->apply = PCMATOP_MULT; 246 pc->data = data; 247 pc->ops->apply = PCApply_Mat; 248 pc->ops->matapply = PCMatApply_Mat; 249 pc->ops->applytranspose = PCApplyTranspose_Mat; 250 pc->ops->setup = NULL; 251 pc->ops->destroy = PCDestroy_Mat; 252 pc->ops->setfromoptions = NULL; 253 pc->ops->view = PCView_Mat; 254 pc->ops->applyrichardson = NULL; 255 pc->ops->applysymmetricleft = NULL; 256 pc->ops->applysymmetricright = NULL; 257 PetscFunctionReturn(PETSC_SUCCESS); 258 } 259