1 /* 2 This provides a thin wrapper around LMVM matrices in order to use their MatLMVMSolve 3 methods as preconditioner applications in KSP solves. 4 */ 5 6 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 7 #include <petsc/private/matimpl.h> 8 9 typedef struct { 10 Vec xwork, ywork; 11 IS inactive; 12 Mat B; 13 PetscBool allocated; 14 } PC_LMVM; 15 16 /*@ 17 PCLMVMSetMatLMVM - Replaces the `MATLMVM` matrix inside the preconditioner with 18 the one provided by the user. 19 20 Input Parameters: 21 + pc - An `PCLMVM` preconditioner 22 - B - An LMVM-type matrix (`MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`) 23 24 Level: intermediate 25 26 .seealso: `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMGetMatLMVM()` 27 @*/ 28 PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B) { 29 PC_LMVM *ctx = (PC_LMVM *)pc->data; 30 PetscBool same; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 34 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 35 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 36 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 37 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same)); 38 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type."); 39 PetscCall(MatDestroy(&ctx->B)); 40 PetscCall(PetscObjectReference((PetscObject)B)); 41 ctx->B = B; 42 PetscFunctionReturn(0); 43 } 44 45 /*@ 46 PCLMVMGetMatLMVM - Returns a pointer to the underlying `MATLMVM` matrix. 47 48 Input Parameter: 49 . pc - An `PCLMVM` preconditioner 50 51 Output Parameter: 52 . B - matrix inside the preconditioner, one of type `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN` 53 54 Level: intermediate 55 56 .seealso: `PCLMVM`, `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMSetMatLMVM()` 57 @*/ 58 PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B) { 59 PC_LMVM *ctx = (PC_LMVM *)pc->data; 60 PetscBool same; 61 62 PetscFunctionBegin; 63 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 64 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 65 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 66 *B = ctx->B; 67 PetscFunctionReturn(0); 68 } 69 70 /*@ 71 PCLMVMSetIS - Sets the index sets that reduce the `PC` application. 72 73 Input Parameters: 74 + pc - An `PCLMVM` preconditioner 75 - inactive - Index set defining the variables removed from the problem 76 77 Level: intermediate 78 79 Developer Note: 80 Need to explain the purpose of this `IS` 81 82 .seealso: `PCLMVM`, `MatLMVMUpdate()` 83 @*/ 84 PetscErrorCode PCLMVMSetIS(PC pc, IS inactive) { 85 PC_LMVM *ctx = (PC_LMVM *)pc->data; 86 PetscBool same; 87 88 PetscFunctionBegin; 89 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 90 PetscValidHeaderSpecific(inactive, IS_CLASSID, 2); 91 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 92 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 93 PetscCall(PCLMVMClearIS(pc)); 94 PetscCall(PetscObjectReference((PetscObject)inactive)); 95 ctx->inactive = inactive; 96 PetscFunctionReturn(0); 97 } 98 99 /*@ 100 PCLMVMClearIS - Removes the inactive variable index set from a `PCLMVM` 101 102 Input Parameter: 103 . pc - An `PCLMVM` preconditioner 104 105 Level: intermediate 106 107 .seealso: `PCLMVM`, `MatLMVMUpdate()` 108 @*/ 109 PetscErrorCode PCLMVMClearIS(PC pc) { 110 PC_LMVM *ctx = (PC_LMVM *)pc->data; 111 PetscBool same; 112 113 PetscFunctionBegin; 114 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 115 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 116 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 117 if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive)); 118 PetscFunctionReturn(0); 119 } 120 121 static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y) { 122 PC_LMVM *ctx = (PC_LMVM *)pc->data; 123 Vec xsub, ysub; 124 125 PetscFunctionBegin; 126 if (ctx->inactive) { 127 PetscCall(VecZeroEntries(ctx->xwork)); 128 PetscCall(VecGetSubVector(ctx->xwork, ctx->inactive, &xsub)); 129 PetscCall(VecCopy(x, xsub)); 130 PetscCall(VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub)); 131 } else { 132 PetscCall(VecCopy(x, ctx->xwork)); 133 } 134 PetscCall(MatSolve(ctx->B, ctx->xwork, ctx->ywork)); 135 if (ctx->inactive) { 136 PetscCall(VecGetSubVector(ctx->ywork, ctx->inactive, &ysub)); 137 PetscCall(VecCopy(ysub, y)); 138 PetscCall(VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub)); 139 } else { 140 PetscCall(VecCopy(ctx->ywork, y)); 141 } 142 PetscFunctionReturn(0); 143 } 144 145 static PetscErrorCode PCReset_LMVM(PC pc) { 146 PC_LMVM *ctx = (PC_LMVM *)pc->data; 147 148 PetscFunctionBegin; 149 if (ctx->xwork) PetscCall(VecDestroy(&ctx->xwork)); 150 if (ctx->ywork) PetscCall(VecDestroy(&ctx->ywork)); 151 PetscFunctionReturn(0); 152 } 153 154 static PetscErrorCode PCSetUp_LMVM(PC pc) { 155 PC_LMVM *ctx = (PC_LMVM *)pc->data; 156 PetscInt n, N; 157 PetscBool allocated; 158 const char *prefix; 159 160 PetscFunctionBegin; 161 if (pc->setupcalled) PetscFunctionReturn(0); 162 PetscCall(MatLMVMIsAllocated(ctx->B, &allocated)); 163 if (!allocated) { 164 PetscCall(MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork)); 165 PetscCall(VecGetLocalSize(ctx->xwork, &n)); 166 PetscCall(VecGetSize(ctx->xwork, &N)); 167 PetscCall(MatSetSizes(ctx->B, n, n, N, N)); 168 PetscCall(MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork)); 169 } else { 170 PetscCall(MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork)); 171 } 172 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 173 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 174 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 175 PetscFunctionReturn(0); 176 } 177 178 static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer) { 179 PC_LMVM *ctx = (PC_LMVM *)pc->data; 180 PetscBool iascii; 181 182 PetscFunctionBegin; 183 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 184 if (iascii && ctx->B->assembled) { 185 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 186 PetscCall(MatView(ctx->B, viewer)); 187 PetscCall(PetscViewerPopFormat(viewer)); 188 } 189 PetscFunctionReturn(0); 190 } 191 192 static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject) { 193 PC_LMVM *ctx = (PC_LMVM *)pc->data; 194 const char *prefix; 195 196 PetscFunctionBegin; 197 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 198 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 199 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 200 PetscCall(MatSetFromOptions(ctx->B)); 201 PetscFunctionReturn(0); 202 } 203 204 static PetscErrorCode PCDestroy_LMVM(PC pc) { 205 PC_LMVM *ctx = (PC_LMVM *)pc->data; 206 207 PetscFunctionBegin; 208 if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive)); 209 if (pc->setupcalled) { 210 PetscCall(VecDestroy(&ctx->xwork)); 211 PetscCall(VecDestroy(&ctx->ywork)); 212 } 213 PetscCall(MatDestroy(&ctx->B)); 214 PetscCall(PetscFree(pc->data)); 215 PetscFunctionReturn(0); 216 } 217 218 /*MC 219 PCLMVM - A a preconditioner constructed from a `MATLMVM` matrix. Options for the 220 underlying `MATLMVM` matrix can be access with the -pc_lmvm_ prefix. 221 222 Level: intermediate 223 224 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, 225 `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()` 226 M*/ 227 PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc) { 228 PC_LMVM *ctx; 229 230 PetscFunctionBegin; 231 PetscCall(PetscNew(&ctx)); 232 pc->data = (void *)ctx; 233 234 pc->ops->reset = PCReset_LMVM; 235 pc->ops->setup = PCSetUp_LMVM; 236 pc->ops->destroy = PCDestroy_LMVM; 237 pc->ops->view = PCView_LMVM; 238 pc->ops->apply = PCApply_LMVM; 239 pc->ops->setfromoptions = PCSetFromOptions_LMVM; 240 pc->ops->applysymmetricleft = NULL; 241 pc->ops->applysymmetricright = NULL; 242 pc->ops->applytranspose = NULL; 243 pc->ops->applyrichardson = NULL; 244 pc->ops->presolve = NULL; 245 pc->ops->postsolve = NULL; 246 247 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B)); 248 PetscCall(MatSetType(ctx->B, MATLMVMBFGS)); 249 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1)); 250 PetscFunctionReturn(0); 251 } 252