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: `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;
35: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
37: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
39: MatDestroy(&ctx->B);
40: PetscObjectReference((PetscObject)B);
41: ctx->B = B;
42: return 0;
43: }
45: /*@
46: PCLMVMGetMatLMVM - Returns a pointer to the underlying `MATLMVM` matrix.
48: Input Parameter:
49: . pc - An `PCLMVM` preconditioner
51: Output Parameter:
52: . B - matrix inside the preconditioner, one of type `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`
54: Level: intermediate
56: .seealso: `PCLMVM`, `MATLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`, `PCLMVMSetMatLMVM()`
57: @*/
58: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
59: {
60: PC_LMVM *ctx = (PC_LMVM *)pc->data;
61: PetscBool same;
64: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
66: *B = ctx->B;
67: return 0;
68: }
70: /*@
71: PCLMVMSetIS - Sets the index sets that reduce the `PC` application.
73: Input Parameters:
74: + pc - An `PCLMVM` preconditioner
75: - inactive - Index set defining the variables removed from the problem
77: Level: intermediate
79: Developer Note:
80: Need to explain the purpose of this `IS`
82: .seealso: `PCLMVM`, `MatLMVMUpdate()`
83: @*/
84: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
85: {
86: PC_LMVM *ctx = (PC_LMVM *)pc->data;
87: PetscBool same;
91: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
93: PCLMVMClearIS(pc);
94: PetscObjectReference((PetscObject)inactive);
95: ctx->inactive = inactive;
96: return 0;
97: }
99: /*@
100: PCLMVMClearIS - Removes the inactive variable index set from a `PCLMVM`
102: Input Parameter:
103: . pc - An `PCLMVM` preconditioner
105: Level: intermediate
107: .seealso: `PCLMVM`, `MatLMVMUpdate()`
108: @*/
109: PetscErrorCode PCLMVMClearIS(PC pc)
110: {
111: PC_LMVM *ctx = (PC_LMVM *)pc->data;
112: PetscBool same;
115: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
117: if (ctx->inactive) ISDestroy(&ctx->inactive);
118: return 0;
119: }
121: static PetscErrorCode PCApply_LMVM(PC pc, Vec x, Vec y)
122: {
123: PC_LMVM *ctx = (PC_LMVM *)pc->data;
124: Vec xsub, ysub;
126: if (ctx->inactive) {
127: VecZeroEntries(ctx->xwork);
128: VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);
129: VecCopy(x, xsub);
130: VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);
131: } else {
132: VecCopy(x, ctx->xwork);
133: }
134: MatSolve(ctx->B, ctx->xwork, ctx->ywork);
135: if (ctx->inactive) {
136: VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);
137: VecCopy(ysub, y);
138: VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);
139: } else {
140: VecCopy(ctx->ywork, y);
141: }
142: return 0;
143: }
145: static PetscErrorCode PCReset_LMVM(PC pc)
146: {
147: PC_LMVM *ctx = (PC_LMVM *)pc->data;
149: if (ctx->xwork) VecDestroy(&ctx->xwork);
150: if (ctx->ywork) VecDestroy(&ctx->ywork);
151: return 0;
152: }
154: static PetscErrorCode PCSetUp_LMVM(PC pc)
155: {
156: PC_LMVM *ctx = (PC_LMVM *)pc->data;
157: PetscInt n, N;
158: PetscBool allocated;
159: const char *prefix;
161: if (pc->setupcalled) return 0;
162: MatLMVMIsAllocated(ctx->B, &allocated);
163: if (!allocated) {
164: MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
165: VecGetLocalSize(ctx->xwork, &n);
166: VecGetSize(ctx->xwork, &N);
167: MatSetSizes(ctx->B, n, n, N, N);
168: MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
169: } else {
170: MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
171: }
172: PCGetOptionsPrefix(pc, &prefix);
173: MatSetOptionsPrefix(ctx->B, prefix);
174: MatAppendOptionsPrefix(ctx->B, "pc_lmvm_");
175: return 0;
176: }
178: static PetscErrorCode PCView_LMVM(PC pc, PetscViewer viewer)
179: {
180: PC_LMVM *ctx = (PC_LMVM *)pc->data;
181: PetscBool iascii;
183: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
184: if (iascii && ctx->B->assembled) {
185: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
186: MatView(ctx->B, viewer);
187: PetscViewerPopFormat(viewer);
188: }
189: return 0;
190: }
192: static PetscErrorCode PCSetFromOptions_LMVM(PC pc, PetscOptionItems *PetscOptionsObject)
193: {
194: PC_LMVM *ctx = (PC_LMVM *)pc->data;
195: const char *prefix;
197: PCGetOptionsPrefix(pc, &prefix);
198: MatSetOptionsPrefix(ctx->B, prefix);
199: MatAppendOptionsPrefix(ctx->B, "pc_lmvm_");
200: MatSetFromOptions(ctx->B);
201: return 0;
202: }
204: static PetscErrorCode PCDestroy_LMVM(PC pc)
205: {
206: PC_LMVM *ctx = (PC_LMVM *)pc->data;
208: if (ctx->inactive) ISDestroy(&ctx->inactive);
209: if (pc->setupcalled) {
210: VecDestroy(&ctx->xwork);
211: VecDestroy(&ctx->ywork);
212: }
213: MatDestroy(&ctx->B);
214: PetscFree(pc->data);
215: return 0;
216: }
218: /*MC
219: PCLMVM - A a preconditioner constructed from a `MATLMVM` matrix. Options for the
220: underlying `MATLMVM` matrix can be access with the -pc_lmvm_ prefix.
222: Level: intermediate
224: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCLMVM`, `MATLDFP`, `MATLBFGS`, `MATLSR1`, `MATLBRDN`, `MATLMBRDN`, `MATLSBRDN`,
225: `PC`, `MATLMVM`, `PCLMVMUpdate()`, `PCLMVMSetMatLMVM()`, `PCLMVMGetMatLMVM()`
226: M*/
227: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
228: {
229: PC_LMVM *ctx;
231: PetscNew(&ctx);
232: pc->data = (void *)ctx;
234: pc->ops->reset = PCReset_LMVM;
235: pc->ops->setup = PCSetUp_LMVM;
236: pc->ops->destroy = PCDestroy_LMVM;
237: pc->ops->view = PCView_LMVM;
238: pc->ops->apply = PCApply_LMVM;
239: pc->ops->setfromoptions = PCSetFromOptions_LMVM;
240: pc->ops->applysymmetricleft = NULL;
241: pc->ops->applysymmetricright = NULL;
242: pc->ops->applytranspose = NULL;
243: pc->ops->applyrichardson = NULL;
244: pc->ops->presolve = NULL;
245: pc->ops->postsolve = NULL;
247: MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
248: MatSetType(ctx->B, MATLMVMBFGS);
249: PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
250: return 0;
251: }