Actual source code: lmvmpc.c
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: */
6: #include <petsc/private/pcimpl.h>
7: #include <petsc/private/matimpl.h>
9: typedef struct {
10: Vec xwork, ywork;
11: IS inactive;
12: Mat B;
13: PetscBool allocated;
14: } PC_LMVM;
16: /*@
17: PCLMVMSetMatLMVM - Replaces the `MATLMVM` matrix inside the preconditioner with
18: the one provided by the user.
20: Input Parameters:
21: + pc - An `PCLMVM` preconditioner
22: - B - An LMVM-type matrix (`MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`)
24: Level: intermediate
26: .seealso: [](ch_ksp), `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMGetMatLMVM()`
27: @*/
28: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
29: {
30: PC_LMVM *ctx = (PC_LMVM *)pc->data;
31: PetscBool same;
33: PetscFunctionBegin;
36: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
37: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
38: PetscCall(PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same));
39: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
40: PetscCall(MatDestroy(&ctx->B));
41: PetscCall(PetscObjectReference((PetscObject)B));
42: ctx->B = B;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: PCLMVMGetMatLMVM - Returns a pointer to the underlying `MATLMVM` matrix.
49: Input Parameter:
50: . pc - An `PCLMVM` preconditioner
52: Output Parameter:
53: . B - matrix inside the preconditioner, one of type `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`
55: Level: intermediate
57: .seealso: [](ch_ksp), `PCLMVM`, `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMSetMatLMVM()`
58: @*/
59: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
60: {
61: PC_LMVM *ctx = (PC_LMVM *)pc->data;
62: PetscBool same;
64: PetscFunctionBegin;
66: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
67: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
68: *B = ctx->B;
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: /*@
73: PCLMVMSetIS - Sets the index sets that reduce the `PC` application.
75: Input Parameters:
76: + pc - An `PCLMVM` preconditioner
77: - inactive - Index set defining the variables removed from the problem
79: Level: intermediate
81: Developer Notes:
82: Need to explain the purpose of this `IS`
84: .seealso: [](ch_ksp), `PCLMVM`, `MatLMVMUpdate()`
85: @*/
86: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
87: {
88: PC_LMVM *ctx = (PC_LMVM *)pc->data;
89: PetscBool same;
91: PetscFunctionBegin;
94: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
95: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
96: PetscCall(PCLMVMClearIS(pc));
97: PetscCall(PetscObjectReference((PetscObject)inactive));
98: ctx->inactive = inactive;
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*@
103: PCLMVMClearIS - Removes the inactive variable index set from a `PCLMVM`
105: Input Parameter:
106: . pc - An `PCLMVM` preconditioner
108: Level: intermediate
110: .seealso: [](ch_ksp), `PCLMVM`, `MatLMVMUpdate()`
111: @*/
112: PetscErrorCode PCLMVMClearIS(PC pc)
113: {
114: PC_LMVM *ctx = (PC_LMVM *)pc->data;
115: PetscBool same;
117: PetscFunctionBegin;
119: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same));
120: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
121: if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
126: {
127: PC_LMVM *ctx = (PC_LMVM *)pc->data;
128: Vec xsub, ysub;
130: PetscFunctionBegin;
131: if (ctx->inactive) {
132: PetscCall(VecZeroEntries(ctx->xwork));
133: PetscCall(VecGetSubVector(ctx->xwork, ctx->inactive, &xsub));
134: PetscCall(VecCopy(x, xsub));
135: PetscCall(VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub));
136: } else {
137: PetscCall(VecCopy(x, ctx->xwork));
138: }
139: PetscCall(MatSolve(ctx->B, ctx->xwork, ctx->ywork));
140: if (ctx->inactive) {
141: PetscCall(VecGetSubVector(ctx->ywork, ctx->inactive, &ysub));
142: PetscCall(VecCopy(ysub, y));
143: PetscCall(VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub));
144: } else {
145: PetscCall(VecCopy(ctx->ywork, y));
146: }
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode PCReset_LMVM(PC pc)
151: {
152: PC_LMVM *ctx = (PC_LMVM *)pc->data;
154: PetscFunctionBegin;
155: if (ctx->xwork) PetscCall(VecDestroy(&ctx->xwork));
156: if (ctx->ywork) PetscCall(VecDestroy(&ctx->ywork));
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: static PetscErrorCode PCSetUp_LMVM(PC pc)
161: {
162: PC_LMVM *ctx = (PC_LMVM *)pc->data;
163: PetscInt n, N;
164: PetscBool allocated;
165: const char *prefix;
167: PetscFunctionBegin;
168: if (pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
169: PetscCall(MatLMVMIsAllocated(ctx->B, &allocated));
170: if (!allocated) {
171: PetscCall(MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork));
172: PetscCall(VecGetLocalSize(ctx->xwork, &n));
173: PetscCall(VecGetSize(ctx->xwork, &N));
174: PetscCall(MatSetSizes(ctx->B, n, n, N, N));
175: PetscCall(MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork));
176: } else {
177: PetscCall(MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork));
178: }
179: PetscCall(PCGetOptionsPrefix(pc, &prefix));
180: PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
181: PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
186: {
187: PC_LMVM *ctx = (PC_LMVM *)pc->data;
188: PetscBool iascii;
190: PetscFunctionBegin;
191: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
192: if (iascii && ctx->B->assembled) {
193: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
194: PetscCall(MatView(ctx->B, viewer));
195: PetscCall(PetscViewerPopFormat(viewer));
196: }
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject)
201: {
202: PC_LMVM *ctx = (PC_LMVM *)pc->data;
203: const char *prefix;
205: PetscFunctionBegin;
206: PetscCall(PCGetOptionsPrefix(pc, &prefix));
207: PetscCall(MatSetOptionsPrefix(ctx->B, prefix));
208: PetscCall(MatAppendOptionsPrefix(ctx->B, "pc_lmvm_"));
209: PetscCall(MatSetFromOptions(ctx->B));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: static PetscErrorCode PCDestroy_LMVM(PC pc)
214: {
215: PC_LMVM *ctx = (PC_LMVM *)pc->data;
217: PetscFunctionBegin;
218: if (ctx->inactive) PetscCall(ISDestroy(&ctx->inactive));
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(PETSC_SUCCESS);
226: }
228: /*MC
229: PCLMVM - A a preconditioner constructed from a `MATLMVM` matrix. Options for the
230: underlying `MATLMVM` matrix can be access with the -pc_lmvm_ prefix.
232: Level: intermediate
234: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`,
235: `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
236: M*/
237: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
238: {
239: PC_LMVM *ctx;
241: PetscFunctionBegin;
242: PetscCall(PetscNew(&ctx));
243: pc->data = (void *)ctx;
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;
258: PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B));
259: PetscCall(MatSetType(ctx->B, MATLMVMBFGS));
260: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }