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 if (pc->setupcalled) PetscFunctionReturn(0); 155 PetscCall(MatLMVMIsAllocated(ctx->B, &allocated)); 156 if (!allocated) { 157 PetscCall(MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork)); 158 PetscCall(VecGetLocalSize(ctx->xwork, &n)); 159 PetscCall(VecGetSize(ctx->xwork, &N)); 160 PetscCall(MatSetSizes(ctx->B, n, n, N, N)); 161 PetscCall(MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork)); 162 } else { 163 PetscCall(MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork)); 164 } 165 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 166 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 167 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 168 PetscFunctionReturn(0); 169 } 170 171 static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer) { 172 PC_LMVM *ctx = (PC_LMVM *)pc->data; 173 PetscBool iascii; 174 175 PetscFunctionBegin; 176 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 177 if (iascii && ctx->B->assembled) { 178 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 179 PetscCall(MatView(ctx->B, viewer)); 180 PetscCall(PetscViewerPopFormat(viewer)); 181 } 182 PetscFunctionReturn(0); 183 } 184 185 static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject) { 186 PC_LMVM *ctx = (PC_LMVM *)pc->data; 187 const char *prefix; 188 189 PetscFunctionBegin; 190 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 191 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 192 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 193 PetscCall(MatSetFromOptions(ctx->B)); 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode PCDestroy_LMVM(PC pc) { 198 PC_LMVM *ctx = (PC_LMVM *)pc->data; 199 200 PetscFunctionBegin; 201 if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive)); 202 if (pc->setupcalled) { 203 PetscCall(VecDestroy(&ctx->xwork)); 204 PetscCall(VecDestroy(&ctx->ywork)); 205 } 206 PetscCall(MatDestroy(&ctx->B)); 207 PetscCall(PetscFree(pc->data)); 208 PetscFunctionReturn(0); 209 } 210 211 /*MC 212 PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the 213 underlying LMVM matrix can be access with the "-pc_lmvm_" prefix. 214 215 Level: intermediate 216 217 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, 218 `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()` 219 M*/ 220 PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc) { 221 PC_LMVM *ctx; 222 223 PetscFunctionBegin; 224 PetscCall(PetscNewLog(pc, &ctx)); 225 pc->data = (void *)ctx; 226 227 pc->ops->reset = PCReset_LMVM; 228 pc->ops->setup = PCSetUp_LMVM; 229 pc->ops->destroy = PCDestroy_LMVM; 230 pc->ops->view = PCView_LMVM; 231 pc->ops->apply = PCApply_LMVM; 232 pc->ops->setfromoptions = PCSetFromOptions_LMVM; 233 pc->ops->applysymmetricleft = NULL; 234 pc->ops->applysymmetricright = NULL; 235 pc->ops->applytranspose = NULL; 236 pc->ops->applyrichardson = NULL; 237 pc->ops->presolve = NULL; 238 pc->ops->postsolve = NULL; 239 240 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B)); 241 PetscCall(MatSetType(ctx->B, MATLMVMBFGS)); 242 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1)); 243 PetscFunctionReturn(0); 244 } 245