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