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