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 189 // clang-format off 190 #define MATOP_TO_PCMATOP_CASE(var, OP) case MATOP_##OP: (var) = PCMATOP_##OP; break 191 switch (matop) { 192 MATOP_TO_PCMATOP_CASE(pcmatop, MULT); 193 MATOP_TO_PCMATOP_CASE(pcmatop, MULT_TRANSPOSE); 194 MATOP_TO_PCMATOP_CASE(pcmatop, MULT_HERMITIAN_TRANSPOSE); 195 MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE); 196 MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE_TRANSPOSE); 197 default: 198 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported MatOperation %d for PCMatSetApplyOperation()", (int)matop); 199 } 200 #undef MATOP_TO_PCMATOP_CASE 201 // clang-format on 202 203 pcmat->apply = pcmatop; 204 PetscFunctionReturn(PETSC_SUCCESS); 205 } 206 207 static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p) 208 { 209 PC_Mat *pcmat = (PC_Mat *)pc->data; 210 MatOperation matop = MATOP_MULT; 211 212 PetscFunctionBegin; 213 if (!pc->setupcalled) PetscCall(PCSetUp(pc)); 214 215 // clang-format off 216 #define PCMATOP_TO_MATOP_CASE(var, OP) case PCMATOP_##OP: (var) = MATOP_##OP; break 217 switch (pcmat->apply) { 218 PCMATOP_TO_MATOP_CASE(matop, MULT); 219 PCMATOP_TO_MATOP_CASE(matop, MULT_TRANSPOSE); 220 PCMATOP_TO_MATOP_CASE(matop, MULT_HERMITIAN_TRANSPOSE); 221 PCMATOP_TO_MATOP_CASE(matop, SOLVE); 222 PCMATOP_TO_MATOP_CASE(matop, SOLVE_TRANSPOSE); 223 default: 224 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]); 225 } 226 #undef PCMATOP_TO_MATOP_CASE 227 // clang-format on 228 229 *matop_p = matop; 230 PetscFunctionReturn(PETSC_SUCCESS); 231 } 232 233 static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer) 234 { 235 PC_Mat *pcmat = (PC_Mat *)pc->data; 236 PetscBool iascii; 237 238 PetscFunctionBegin; 239 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 240 if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply])); } 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /*MC 245 PCMAT - A preconditioner obtained by applying an operation of the `pmat` provided in 246 in `PCSetOperators(pc,amat,pmat)` or `KSPSetOperators(ksp,amat,pmat)`. By default, the operation is `MATOP_MULT`, 247 meaning that the `pmat` provides an approximate inverse of `amat`. 248 If some other operation of `pmat` implements the approximate inverse, 249 use `PCMatSetApplyOperation()` to select that operation. 250 251 Level: intermediate 252 253 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()` 254 M*/ 255 256 PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc) 257 { 258 PC_Mat *data; 259 260 PetscFunctionBegin; 261 PetscCall(PetscNew(&data)); 262 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", PCMatSetApplyOperation_Mat)); 263 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", PCMatGetApplyOperation_Mat)); 264 pc->data = data; 265 pc->ops->apply = PCApply_Mat; 266 pc->ops->matapply = PCMatApply_Mat; 267 pc->ops->applytranspose = PCApplyTranspose_Mat; 268 pc->ops->setup = PCSetUp_Mat; 269 pc->ops->destroy = PCDestroy_Mat; 270 pc->ops->setfromoptions = NULL; 271 pc->ops->view = PCView_Mat; 272 pc->ops->applyrichardson = NULL; 273 pc->ops->applysymmetricleft = NULL; 274 pc->ops->applysymmetricright = NULL; 275 PetscFunctionReturn(PETSC_SUCCESS); 276 } 277