Actual source code: lmvmpc.c
petsc-3.10.5 2019-03-28
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>
8: typedef struct {
9: Vec xwork, ywork;
10: IS inactive;
11: Mat B;
12: PetscBool allocated;
13: } PC_LMVM;
15: /*@
16: PCLMVMSetMatLMVM - Replaces the LMVM matrix inside the preconditioner with
17: the one provided by the user.
18:
19: Input Parameters:
20: + pc - An LMVM preconditioner
21: - B - An LMVM-type matrix (MATLDFP, MATLBFGS, MATLSR1, MATLBRDN, MATLMBRDN, MATLSBRDN)
22:
23: Level: intermediate
24: @*/
25: PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B)
26: {
27: PC_LMVM *ctx = (PC_LMVM*)pc->data;
28: PetscErrorCode ierr;
29: PetscBool same;
30:
34: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
35: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
36: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
37: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
38: MatDestroy(&ctx->B);
39: PetscObjectReference((PetscObject)B);
40: ctx->B = B;
41: return(0);
42: }
44: /*@
45: PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix.
46:
47: Input Parameters:
48: . pc - An LMVM preconditioner
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: PetscErrorCode ierr;
59: PetscBool same;
60:
63: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
64: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
65: *B = ctx->B;
66: return(0);
67: }
69: /*@
70: PCLMVMSetIS - Sets the index sets that reduce the PC application.
71:
72: Input Parameters:
73: + pc - An LMVM preconditioner
74: - inactive - Index set defining the variables removed from the problem
75:
76: Level: intermediate
78: .seealso: MatLMVMUpdate()
79: @*/
80: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
81: {
82: PC_LMVM *ctx = (PC_LMVM*)pc->data;
83: PetscErrorCode ierr;
84: PetscBool same;
85:
89: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
90: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
91: PCLMVMClearIS(pc);
92: PetscObjectReference((PetscObject)inactive);
93: ctx->inactive = inactive;
94: return(0);
95: }
97: /*@
98: PCLMVMClearIS - Removes the inactive variable index set.
99:
100: Input Parameters:
101: . pc - An LMVM preconditioner
102:
103: Level: intermediate
105: .seealso: MatLMVMUpdate()
106: @*/
107: PetscErrorCode PCLMVMClearIS(PC pc)
108: {
109: PC_LMVM *ctx = (PC_LMVM*)pc->data;
110: PetscErrorCode ierr;
111: PetscBool same;
112:
115: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
116: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
117: if (ctx->inactive) {
118: ISDestroy(&ctx->inactive);
119: }
120: return(0);
121: }
123: static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y)
124: {
125: PC_LMVM *ctx = (PC_LMVM*)pc->data;
126: PetscErrorCode ierr;
127: Vec xsub, ysub;
130: if (ctx->inactive) {
131: VecZeroEntries(ctx->xwork);
132: VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);
133: VecCopy(x, xsub);
134: VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);
135: } else {
136: VecCopy(x, ctx->xwork);
137: }
138: MatSolve(ctx->B, ctx->xwork, ctx->ywork);
139: if (ctx->inactive) {
140: VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);
141: VecCopy(ysub, y);
142: VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);
143: } else {
144: VecCopy(ctx->ywork, y);
145: }
146: return(0);
147: }
149: static PetscErrorCode PCReset_LMVM(PC pc)
150: {
151: PC_LMVM *ctx = (PC_LMVM*)pc->data;
153:
155: if (ctx->xwork) {
156: VecDestroy(&ctx->xwork);
157: }
158: if (ctx->ywork) {
159: VecDestroy(&ctx->ywork);
160: }
161: return(0);
162: }
164: static PetscErrorCode PCSetUp_LMVM(PC pc)
165: {
166: PC_LMVM *ctx = (PC_LMVM*)pc->data;
168: PetscInt n, N;
169: PetscBool allocated;
172: MatLMVMIsAllocated(ctx->B, &allocated);
173: if (!allocated) {
174: MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
175: VecGetLocalSize(ctx->xwork, &n);
176: VecGetSize(ctx->xwork, &N);
177: MatSetSizes(ctx->B, n, n, N, N);
178: MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
179: } else {
180: MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
181: }
182: return(0);
183: }
185: static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer)
186: {
187: PC_LMVM *ctx = (PC_LMVM*)pc->data;
189: PetscBool iascii;
192: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
193: if (iascii) {
194: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
195: MatView(ctx->B, viewer);
196: PetscViewerPopFormat(viewer);
197: }
198: return(0);
199: }
201: static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc)
202: {
203: PC_LMVM *ctx = (PC_LMVM*)pc->data;
205:
207: MatSetFromOptions(ctx->B);
208: return(0);
209: }
211: static PetscErrorCode PCDestroy_LMVM(PC pc)
212: {
213: PC_LMVM *ctx = (PC_LMVM*)pc->data;
217: if (ctx->inactive) {
218: ISDestroy(&ctx->inactive);
219: }
220: if (pc->setupcalled) {
221: VecDestroy(&ctx->xwork);
222: VecDestroy(&ctx->ywork);
223: }
224: MatDestroy(&ctx->B);
225: PetscFree(pc->data);
226: return(0);
227: }
229: /*MC
230: PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the
231: underlying LMVM matrix can be access with the "-pc_lmvm_" prefix.
233: Level: intermediate
235: .seealso: PCCreate(), PCSetType(), PCType (for list of available types),
236: PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM()
237: M*/
238: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
239: {
241: PC_LMVM *ctx;
244: PetscNewLog(pc,&ctx);
245: pc->data = (void*)ctx;
246:
247: pc->ops->reset = PCReset_LMVM;
248: pc->ops->setup = PCSetUp_LMVM;
249: pc->ops->destroy = PCDestroy_LMVM;
250: pc->ops->view = PCView_LMVM;
251: pc->ops->apply = PCApply_LMVM;
252: pc->ops->setfromoptions = PCSetFromOptions_LMVM;
253: pc->ops->applysymmetricleft = 0;
254: pc->ops->applysymmetricright = 0;
255: pc->ops->applytranspose = 0;
256: pc->ops->applyrichardson = 0;
257: pc->ops->presolve = 0;
258: pc->ops->postsolve = 0;
259:
260: PCSetReusePreconditioner(pc, PETSC_TRUE);
261:
262: MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
263: MatSetType(ctx->B, MATLMVMBFGS);
264: PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
265: MatSetOptionsPrefix(ctx->B, "pc_lmvm_");
266: return(0);
267: }