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