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