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 Vec X;
14 PetscObjectState Xstate;
15 PetscBool setfromoptionscalled;
16 } PC_LMVM;
17
18 /*@
19 PCLMVMSetUpdateVec - Set the vector to be used as solution update for the internal LMVM matrix.
20
21 Input Parameters:
22 + pc - The preconditioner
23 - X - Solution vector
24
25 Level: intermediate
26
27 Notes:
28 This is only needed if you want the preconditioner to automatically update the internal matrix.
29 It is called in some `SNES` implementations to update the preconditioner.
30 The right-hand side of the linear system is used as function vector.
31
32 .seealso: `MatLMVMUpdate()`, `PCLMVMSetMatLMVM()`
33 @*/
PCLMVMSetUpdateVec(PC pc,Vec X)34 PetscErrorCode PCLMVMSetUpdateVec(PC pc, Vec X)
35 {
36 PC_LMVM *ctx;
37 PetscBool same;
38
39 PetscFunctionBegin;
40 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
41 if (X) PetscValidHeaderSpecific(X, VEC_CLASSID, 2);
42 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
43 if (!same) PetscFunctionReturn(PETSC_SUCCESS);
44 ctx = (PC_LMVM *)pc->data;
45 PetscCall(PetscObjectReference((PetscObject)X));
46 PetscCall(VecDestroy(&ctx->X));
47 ctx->X = X;
48 ctx->Xstate = -1;
49 PetscFunctionReturn(PETSC_SUCCESS);
50 }
51
52 /*@
53 PCLMVMSetMatLMVM - Replaces the `MATLMVM` matrix inside the preconditioner with the one provided by the user.
54
55 Input Parameters:
56 + pc - An `PCLMVM` preconditioner
57 - B - An `MATLMVM` type matrix
58
59 Level: intermediate
60
61 .seealso: [](ch_ksp), `PCLMVMGetMatLMVM()`
62 @*/
PCLMVMSetMatLMVM(PC pc,Mat B)63 PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
64 {
65 PC_LMVM *ctx;
66 PetscBool same;
67
68 PetscFunctionBegin;
69 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
70 PetscValidHeaderSpecific(B, MAT_CLASSID, 2);
71 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
72 if (!same) PetscFunctionReturn(PETSC_SUCCESS);
73 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
74 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an MATLMVM.");
75 ctx = (PC_LMVM *)pc->data;
76 PetscCall(PetscObjectReference((PetscObject)B));
77 PetscCall(MatDestroy(&ctx->B));
78 ctx->B = B;
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81
82 /*@
83 PCLMVMGetMatLMVM - Returns a pointer to the underlying `MATLMVM` matrix.
84
85 Input Parameter:
86 . pc - An `PCLMVM` preconditioner
87
88 Output Parameter:
89 . B - `MATLMVM` matrix used by the preconditioner
90
91 Level: intermediate
92
93 .seealso: [](ch_ksp), `PCLMVMSetMatLMVM()`
94 @*/
PCLMVMGetMatLMVM(PC pc,Mat * B)95 PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
96 {
97 PC_LMVM *ctx;
98 PetscBool same;
99
100 PetscFunctionBegin;
101 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
102 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
103 PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
104 ctx = (PC_LMVM *)pc->data;
105 if (!ctx->B) {
106 Mat J;
107
108 if (pc->useAmat) J = pc->mat;
109 else J = pc->pmat;
110 PetscCall(PetscObjectBaseTypeCompare((PetscObject)J, MATLMVM, &same));
111 if (same) *B = J;
112 else {
113 const char *prefix;
114
115 PetscCall(PCGetOptionsPrefix(pc, &prefix));
116 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B));
117 PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
118 PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
119 PetscCall(MatSetType(ctx->B, MATLMVMBFGS));
120 PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1));
121 *B = ctx->B;
122 }
123 } else *B = ctx->B;
124 PetscFunctionReturn(PETSC_SUCCESS);
125 }
126
127 /*@
128 PCLMVMSetIS - Sets the index sets that reduce the `PC` application.
129
130 Input Parameters:
131 + pc - An `PCLMVM` preconditioner
132 - inactive - Index set defining the variables removed from the problem
133
134 Level: intermediate
135
136 Developer Notes:
137 Need to explain the purpose of this `IS`
138
139 .seealso: [](ch_ksp), `PCLMVMClearIS()`
140 @*/
PCLMVMSetIS(PC pc,IS inactive)141 PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
142 {
143 PC_LMVM *ctx;
144 PetscBool same;
145
146 PetscFunctionBegin;
147 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
148 PetscValidHeaderSpecific(inactive, IS_CLASSID, 2);
149 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
150 if (!same) PetscFunctionReturn(PETSC_SUCCESS);
151 ctx = (PC_LMVM *)pc->data;
152 PetscCall(PCLMVMClearIS(pc));
153 PetscCall(PetscObjectReference((PetscObject)inactive));
154 ctx->inactive = inactive;
155 PetscFunctionReturn(PETSC_SUCCESS);
156 }
157
158 /*@
159 PCLMVMClearIS - Removes the inactive variable index set from a `PCLMVM`
160
161 Input Parameter:
162 . pc - An `PCLMVM` preconditioner
163
164 Level: intermediate
165
166 .seealso: [](ch_ksp), `PCLMVMSetIS()`
167 @*/
PCLMVMClearIS(PC pc)168 PetscErrorCode PCLMVMClearIS(PC pc)
169 {
170 PC_LMVM *ctx;
171 PetscBool same;
172
173 PetscFunctionBegin;
174 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
175 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
176 if (!same) PetscFunctionReturn(PETSC_SUCCESS);
177 ctx = (PC_LMVM *)pc->data;
178 PetscCall(ISDestroy(&ctx->inactive));
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
PCApply_LMVM(PC pc,Vec x,Vec y)182 static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
183 {
184 PC_LMVM *ctx = (PC_LMVM *)pc->data;
185 Vec xsub, ysub, Bx = x, By = y;
186 Mat B = ctx->B ? ctx->B : (pc->useAmat ? pc->mat : pc->pmat);
187
188 PetscFunctionBegin;
189 if (ctx->inactive) {
190 if (!ctx->xwork) PetscCall(MatCreateVecs(B, &ctx->xwork, &ctx->ywork));
191 PetscCall(VecZeroEntries(ctx->xwork));
192 PetscCall(VecGetSubVector(ctx->xwork, ctx->inactive, &xsub));
193 PetscCall(VecCopy(x, xsub));
194 PetscCall(VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub));
195 Bx = ctx->xwork;
196 By = ctx->ywork;
197 }
198 PetscCall(MatSolve(B, Bx, By));
199 if (ctx->inactive) {
200 PetscCall(VecGetSubVector(ctx->ywork, ctx->inactive, &ysub));
201 PetscCall(VecCopy(ysub, y));
202 PetscCall(VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub));
203 }
204 PetscFunctionReturn(PETSC_SUCCESS);
205 }
206
PCReset_LMVM(PC pc)207 static PetscErrorCode PCReset_LMVM(PC pc)
208 {
209 PC_LMVM *ctx = (PC_LMVM *)pc->data;
210
211 PetscFunctionBegin;
212 PetscCall(VecDestroy(&ctx->xwork));
213 PetscCall(VecDestroy(&ctx->ywork));
214 PetscFunctionReturn(PETSC_SUCCESS);
215 }
216
PCSetUp_LMVM(PC pc)217 static PetscErrorCode PCSetUp_LMVM(PC pc)
218 {
219 PetscInt n, N;
220 PetscBool allocated;
221 Mat B;
222 PC_LMVM *ctx = (PC_LMVM *)pc->data;
223
224 PetscFunctionBegin;
225 PetscCall(PCLMVMGetMatLMVM(pc, &B));
226 PetscCall(MatLMVMIsAllocated(B, &allocated));
227 if (!allocated) {
228 Vec t1, t2;
229
230 PetscCall(MatCreateVecs(pc->mat, &t1, &t2));
231 PetscCall(VecGetLocalSize(t1, &n));
232 PetscCall(VecGetSize(t1, &N));
233 PetscCall(MatSetSizes(B, n, n, N, N));
234 PetscCall(MatLMVMAllocate(B, t1, t2));
235 PetscCall(VecDestroy(&t1));
236 PetscCall(VecDestroy(&t2));
237 }
238 /* Only call SetFromOptions if we internally handle the LMVM matrix */
239 if (B == ctx->B && ctx->setfromoptionscalled) PetscCall(MatSetFromOptions(ctx->B));
240 ctx->setfromoptionscalled = PETSC_FALSE;
241 PetscFunctionReturn(PETSC_SUCCESS);
242 }
243
PCView_LMVM(PC pc,PetscViewer viewer)244 static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
245 {
246 PC_LMVM *ctx = (PC_LMVM *)pc->data;
247 PetscBool isascii;
248
249 PetscFunctionBegin;
250 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
251 if (isascii && ctx->B && ctx->B->assembled) {
252 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
253 PetscCall(MatView(ctx->B, viewer));
254 PetscCall(PetscViewerPopFormat(viewer));
255 }
256 PetscFunctionReturn(PETSC_SUCCESS);
257 }
258
PCSetFromOptions_LMVM(PC pc,PetscOptionItems PetscOptionsObject)259 static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems PetscOptionsObject)
260 {
261 PC_LMVM *ctx = (PC_LMVM *)pc->data;
262
263 PetscFunctionBegin;
264 /* defer SetFromOptions calls to PCSetUp_LMVM */
265 ctx->setfromoptionscalled = PETSC_TRUE;
266 PetscFunctionReturn(PETSC_SUCCESS);
267 }
268
PCPreSolve_LMVM(PC pc,KSP ksp,Vec F,Vec X)269 static PetscErrorCode PCPreSolve_LMVM(PC pc, KSP ksp, Vec F, Vec X)
270 {
271 PC_LMVM *ctx = (PC_LMVM *)pc->data;
272
273 PetscFunctionBegin;
274 if (ctx->X && ctx->B) { /* Perform update only if requested. Otherwise we assume the user, e.g. TAO, has already taken care of it */
275 PetscObjectState Xstate;
276
277 PetscCall(PetscObjectStateGet((PetscObject)ctx->X, &Xstate));
278 if (ctx->Xstate != Xstate) PetscCall(MatLMVMUpdate(ctx->B, ctx->X, F));
279 ctx->Xstate = Xstate;
280 }
281 PetscFunctionReturn(PETSC_SUCCESS);
282 }
283
PCDestroy_LMVM(PC pc)284 static PetscErrorCode PCDestroy_LMVM(PC pc)
285 {
286 PC_LMVM *ctx = (PC_LMVM *)pc->data;
287
288 PetscFunctionBegin;
289 PetscCall(ISDestroy(&ctx->inactive));
290 PetscCall(VecDestroy(&ctx->xwork));
291 PetscCall(VecDestroy(&ctx->ywork));
292 PetscCall(VecDestroy(&ctx->X));
293 PetscCall(MatDestroy(&ctx->B));
294 PetscCall(PetscFree(pc->data));
295 PetscFunctionReturn(PETSC_SUCCESS);
296 }
297
298 /*MC
299 PCLMVM - A preconditioner constructed from a `MATLMVM` matrix.
300 If the matrix used to construct the preconditioner is not of type `MATLMVM`, an internal matrix is used.
301 Options for the internal `MATLMVM` matrix can be accessed with the `-pc_lmvm_` prefix.
302 Alternatively, the user can pass a suitable matrix with `PCLMVMSetMatLMVM()`.
303
304 Level: intermediate
305
306 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCLMVMSetUpdateVec()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
307 M*/
PCCreate_LMVM(PC pc)308 PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
309 {
310 PC_LMVM *ctx;
311
312 PetscFunctionBegin;
313 PetscCall(PetscNew(&ctx));
314 pc->data = (void *)ctx;
315
316 pc->ops->reset = PCReset_LMVM;
317 pc->ops->setup = PCSetUp_LMVM;
318 pc->ops->destroy = PCDestroy_LMVM;
319 pc->ops->view = PCView_LMVM;
320 pc->ops->apply = PCApply_LMVM;
321 pc->ops->setfromoptions = PCSetFromOptions_LMVM;
322 pc->ops->applysymmetricleft = NULL;
323 pc->ops->applysymmetricright = NULL;
324 pc->ops->applytranspose = NULL;
325 pc->ops->applyrichardson = NULL;
326 pc->ops->presolve = PCPreSolve_LMVM;
327 pc->ops->postsolve = NULL;
328 PetscFunctionReturn(PETSC_SUCCESS);
329 }
330