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