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 LMVM matrix inside the preconditioner with
18: the one provided by the user.
20: Input Parameters:
21: + pc - An LMVM preconditioner
22: - B - An LMVM-type matrix (MATLDFP, MATLBFGS, MATLSR1, MATLBRDN, MATLMBRDN, MATLSBRDN)
24: Level: intermediate
25: @*/
26: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
27: {
28: PC_LMVM *ctx = (PC_LMVM*)pc->data;
29: PetscBool same;
33: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
35: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
37: MatDestroy(&ctx->B);
38: PetscObjectReference((PetscObject)B);
39: ctx->B = B;
40: return 0;
41: }
43: /*@
44: PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix.
46: Input Parameters:
47: . pc - An LMVM preconditioner
49: Output Parameters:
50: . B - LMVM matrix inside the preconditioner
52: Level: intermediate
53: @*/
54: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
55: {
56: PC_LMVM *ctx = (PC_LMVM*)pc->data;
57: PetscBool same;
60: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
62: *B = ctx->B;
63: return 0;
64: }
66: /*@
67: PCLMVMSetIS - Sets the index sets that reduce the PC application.
69: Input Parameters:
70: + pc - An LMVM preconditioner
71: - inactive - Index set defining the variables removed from the problem
73: Level: intermediate
75: .seealso: MatLMVMUpdate()
76: @*/
77: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
78: {
79: PC_LMVM *ctx = (PC_LMVM*)pc->data;
80: PetscBool same;
84: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
86: PCLMVMClearIS(pc);
87: PetscObjectReference((PetscObject)inactive);
88: ctx->inactive = inactive;
89: return 0;
90: }
92: /*@
93: PCLMVMClearIS - Removes the inactive variable index set.
95: Input Parameters:
96: . pc - An LMVM preconditioner
98: Level: intermediate
100: .seealso: MatLMVMUpdate()
101: @*/
102: PetscErrorCode PCLMVMClearIS(PC pc)
103: {
104: PC_LMVM *ctx = (PC_LMVM*)pc->data;
105: PetscBool same;
108: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
110: if (ctx->inactive) {
111: ISDestroy(&ctx->inactive);
112: }
113: return 0;
114: }
116: static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y)
117: {
118: PC_LMVM *ctx = (PC_LMVM*)pc->data;
119: Vec xsub, ysub;
121: if (ctx->inactive) {
122: VecZeroEntries(ctx->xwork);
123: VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);
124: VecCopy(x, xsub);
125: VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);
126: } else {
127: VecCopy(x, ctx->xwork);
128: }
129: MatSolve(ctx->B, ctx->xwork, ctx->ywork);
130: if (ctx->inactive) {
131: VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);
132: VecCopy(ysub, y);
133: VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);
134: } else {
135: VecCopy(ctx->ywork, y);
136: }
137: return 0;
138: }
140: static PetscErrorCode PCReset_LMVM(PC pc)
141: {
142: PC_LMVM *ctx = (PC_LMVM*)pc->data;
144: if (ctx->xwork) {
145: VecDestroy(&ctx->xwork);
146: }
147: if (ctx->ywork) {
148: VecDestroy(&ctx->ywork);
149: }
150: return 0;
151: }
153: static PetscErrorCode PCSetUp_LMVM(PC pc)
154: {
155: PC_LMVM *ctx = (PC_LMVM*)pc->data;
156: PetscInt n, N;
157: PetscBool allocated;
159: MatLMVMIsAllocated(ctx->B, &allocated);
160: if (!allocated) {
161: MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
162: VecGetLocalSize(ctx->xwork, &n);
163: VecGetSize(ctx->xwork, &N);
164: MatSetSizes(ctx->B, n, n, N, N);
165: MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
166: } else {
167: MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
168: }
169: return 0;
170: }
172: static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer)
173: {
174: PC_LMVM *ctx = (PC_LMVM*)pc->data;
175: PetscBool iascii;
177: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
178: if (iascii && ctx->B->assembled) {
179: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
180: MatView(ctx->B, viewer);
181: PetscViewerPopFormat(viewer);
182: }
183: return 0;
184: }
186: static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc)
187: {
188: PC_LMVM *ctx = (PC_LMVM*)pc->data;
190: MatSetFromOptions(ctx->B);
191: return 0;
192: }
194: static PetscErrorCode PCDestroy_LMVM(PC pc)
195: {
196: PC_LMVM *ctx = (PC_LMVM*)pc->data;
198: if (ctx->inactive) {
199: ISDestroy(&ctx->inactive);
200: }
201: if (pc->setupcalled) {
202: VecDestroy(&ctx->xwork);
203: VecDestroy(&ctx->ywork);
204: }
205: MatDestroy(&ctx->B);
206: PetscFree(pc->data);
207: return 0;
208: }
210: /*MC
211: PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the
212: underlying LMVM matrix can be access with the "-pc_lmvm_" prefix.
214: Level: intermediate
216: .seealso: PCCreate(), PCSetType(), PCType (for list of available types),
217: PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM()
218: M*/
219: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
220: {
221: PC_LMVM *ctx;
223: PetscNewLog(pc,&ctx);
224: pc->data = (void*)ctx;
226: pc->ops->reset = PCReset_LMVM;
227: pc->ops->setup = PCSetUp_LMVM;
228: pc->ops->destroy = PCDestroy_LMVM;
229: pc->ops->view = PCView_LMVM;
230: pc->ops->apply = PCApply_LMVM;
231: pc->ops->setfromoptions = PCSetFromOptions_LMVM;
232: pc->ops->applysymmetricleft = NULL;
233: pc->ops->applysymmetricright = NULL;
234: pc->ops->applytranspose = NULL;
235: pc->ops->applyrichardson = NULL;
236: pc->ops->presolve = NULL;
237: pc->ops->postsolve = NULL;
239: PCSetReusePreconditioner(pc, PETSC_TRUE);
241: MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
242: MatSetType(ctx->B, MATLMVMBFGS);
243: PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
244: MatSetOptionsPrefix(ctx->B, "pc_lmvm_");
245: return 0;
246: }