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 { 28 PC_LMVM *ctx = (PC_LMVM*)pc->data; 29 PetscBool same; 30 31 PetscFunctionBegin; 32 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 33 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 34 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 35 PetscCheck(same,PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 36 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same)); 37 PetscCheck(same,PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type."); 38 PetscCall(MatDestroy(&ctx->B)); 39 PetscCall(PetscObjectReference((PetscObject)B)); 40 ctx->B = B; 41 PetscFunctionReturn(0); 42 } 43 44 /*@ 45 PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix. 46 47 Input Parameters: 48 . pc - An LMVM preconditioner 49 50 Output Parameters: 51 . B - LMVM matrix inside the preconditioner 52 53 Level: intermediate 54 @*/ 55 PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B) 56 { 57 PC_LMVM *ctx = (PC_LMVM*)pc->data; 58 PetscBool same; 59 60 PetscFunctionBegin; 61 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 62 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 63 PetscCheck(same,PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 64 *B = ctx->B; 65 PetscFunctionReturn(0); 66 } 67 68 /*@ 69 PCLMVMSetIS - Sets the index sets that reduce the PC application. 70 71 Input Parameters: 72 + pc - An LMVM preconditioner 73 - inactive - Index set defining the variables removed from the problem 74 75 Level: intermediate 76 77 .seealso: `MatLMVMUpdate()` 78 @*/ 79 PetscErrorCode PCLMVMSetIS(PC pc, IS inactive) 80 { 81 PC_LMVM *ctx = (PC_LMVM*)pc->data; 82 PetscBool same; 83 84 PetscFunctionBegin; 85 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 86 PetscValidHeaderSpecific(inactive, IS_CLASSID, 2); 87 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 88 PetscCheck(same,PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 89 PetscCall(PCLMVMClearIS(pc)); 90 PetscCall(PetscObjectReference((PetscObject)inactive)); 91 ctx->inactive = inactive; 92 PetscFunctionReturn(0); 93 } 94 95 /*@ 96 PCLMVMClearIS - Removes the inactive variable index set. 97 98 Input Parameters: 99 . pc - An LMVM preconditioner 100 101 Level: intermediate 102 103 .seealso: `MatLMVMUpdate()` 104 @*/ 105 PetscErrorCode PCLMVMClearIS(PC pc) 106 { 107 PC_LMVM *ctx = (PC_LMVM*)pc->data; 108 PetscBool same; 109 110 PetscFunctionBegin; 111 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 112 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same)); 113 PetscCheck(same,PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 114 if (ctx->inactive) { 115 PetscCall(ISDestroy(&ctx->inactive)); 116 } 117 PetscFunctionReturn(0); 118 } 119 120 static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y) 121 { 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 { 147 PC_LMVM *ctx = (PC_LMVM*)pc->data; 148 149 PetscFunctionBegin; 150 if (ctx->xwork) { 151 PetscCall(VecDestroy(&ctx->xwork)); 152 } 153 if (ctx->ywork) { 154 PetscCall(VecDestroy(&ctx->ywork)); 155 } 156 PetscFunctionReturn(0); 157 } 158 159 static PetscErrorCode PCSetUp_LMVM(PC pc) 160 { 161 PC_LMVM *ctx = (PC_LMVM*)pc->data; 162 PetscInt n, N; 163 PetscBool allocated; 164 const char *prefix; 165 166 PetscFunctionBegin; 167 PetscCall(MatLMVMIsAllocated(ctx->B, &allocated)); 168 if (!allocated) { 169 PetscCall(MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork)); 170 PetscCall(VecGetLocalSize(ctx->xwork, &n)); 171 PetscCall(VecGetSize(ctx->xwork, &N)); 172 PetscCall(MatSetSizes(ctx->B, n, n, N, N)); 173 PetscCall(MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork)); 174 } else { 175 PetscCall(MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork)); 176 } 177 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 178 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 179 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 180 PetscFunctionReturn(0); 181 } 182 183 static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer) 184 { 185 PC_LMVM *ctx = (PC_LMVM*)pc->data; 186 PetscBool iascii; 187 188 PetscFunctionBegin; 189 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 190 if (iascii && ctx->B->assembled) { 191 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 192 PetscCall(MatView(ctx->B, viewer)); 193 PetscCall(PetscViewerPopFormat(viewer)); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc) 199 { 200 PC_LMVM *ctx = (PC_LMVM*)pc->data; 201 const char *prefix; 202 203 PetscFunctionBegin; 204 PetscCall(PCGetOptionsPrefix(pc, &prefix)); 205 PetscCall(MatSetOptionsPrefix(ctx->B, prefix)); 206 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_")); 207 PetscCall(MatSetFromOptions(ctx->B)); 208 PetscFunctionReturn(0); 209 } 210 211 static PetscErrorCode PCDestroy_LMVM(PC pc) 212 { 213 PC_LMVM *ctx = (PC_LMVM*)pc->data; 214 215 PetscFunctionBegin; 216 if (ctx->inactive) { 217 PetscCall(ISDestroy(&ctx->inactive)); 218 } 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(0); 226 } 227 228 /*MC 229 PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the 230 underlying LMVM matrix can be access with the "-pc_lmvm_" prefix. 231 232 Level: intermediate 233 234 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, 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(PetscNewLog(pc,&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(PCSetReusePreconditioner(pc, PETSC_TRUE)); 259 260 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B)); 261 PetscCall(MatSetType(ctx->B, MATLMVMBFGS)); 262 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1)); 263 PetscFunctionReturn(0); 264 } 265