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: PetscErrorCode ierr;
30: PetscBool same;
35: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
36: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
37: PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);
38: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type.");
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 LMVM matrix.
48: Input Parameters:
49: . pc - An LMVM preconditioner
51: Output Parameters:
52: . B - LMVM matrix inside the preconditioner
54: Level: intermediate
55: @*/
56: PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B)
57: {
58: PC_LMVM *ctx = (PC_LMVM*)pc->data;
59: PetscErrorCode ierr;
60: PetscBool same;
64: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
65: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
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 LMVM preconditioner
75: - inactive - Index set defining the variables removed from the problem
77: Level: intermediate
79: .seealso: MatLMVMUpdate()
80: @*/
81: PetscErrorCode PCLMVMSetIS(PC pc, IS inactive)
82: {
83: PC_LMVM *ctx = (PC_LMVM*)pc->data;
84: PetscErrorCode ierr;
85: PetscBool same;
90: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
91: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
92: PCLMVMClearIS(pc);
93: PetscObjectReference((PetscObject)inactive);
94: ctx->inactive = inactive;
95: return(0);
96: }
98: /*@
99: PCLMVMClearIS - Removes the inactive variable index set.
101: Input Parameters:
102: . pc - An LMVM preconditioner
104: Level: intermediate
106: .seealso: MatLMVMUpdate()
107: @*/
108: PetscErrorCode PCLMVMClearIS(PC pc)
109: {
110: PC_LMVM *ctx = (PC_LMVM*)pc->data;
111: PetscErrorCode ierr;
112: PetscBool same;
116: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);
117: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type.");
118: if (ctx->inactive) {
119: ISDestroy(&ctx->inactive);
120: }
121: return(0);
122: }
124: static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y)
125: {
126: PC_LMVM *ctx = (PC_LMVM*)pc->data;
127: PetscErrorCode ierr;
128: Vec xsub, ysub;
131: if (ctx->inactive) {
132: VecZeroEntries(ctx->xwork);
133: VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);
134: VecCopy(x, xsub);
135: VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);
136: } else {
137: VecCopy(x, ctx->xwork);
138: }
139: MatSolve(ctx->B, ctx->xwork, ctx->ywork);
140: if (ctx->inactive) {
141: VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);
142: VecCopy(ysub, y);
143: VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);
144: } else {
145: VecCopy(ctx->ywork, y);
146: }
147: return(0);
148: }
150: static PetscErrorCode PCReset_LMVM(PC pc)
151: {
152: PC_LMVM *ctx = (PC_LMVM*)pc->data;
156: if (ctx->xwork) {
157: VecDestroy(&ctx->xwork);
158: }
159: if (ctx->ywork) {
160: VecDestroy(&ctx->ywork);
161: }
162: return(0);
163: }
165: static PetscErrorCode PCSetUp_LMVM(PC pc)
166: {
167: PC_LMVM *ctx = (PC_LMVM*)pc->data;
169: PetscInt n, N;
170: PetscBool allocated;
173: MatLMVMIsAllocated(ctx->B, &allocated);
174: if (!allocated) {
175: MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);
176: VecGetLocalSize(ctx->xwork, &n);
177: VecGetSize(ctx->xwork, &N);
178: MatSetSizes(ctx->B, n, n, N, N);
179: MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);
180: } else {
181: MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);
182: }
183: return(0);
184: }
186: static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer)
187: {
188: PC_LMVM *ctx = (PC_LMVM*)pc->data;
190: PetscBool iascii;
193: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
194: if (iascii && ctx->B->assembled) {
195: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
196: MatView(ctx->B, viewer);
197: PetscViewerPopFormat(viewer);
198: }
199: return(0);
200: }
202: static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc)
203: {
204: PC_LMVM *ctx = (PC_LMVM*)pc->data;
208: MatSetFromOptions(ctx->B);
209: return(0);
210: }
212: static PetscErrorCode PCDestroy_LMVM(PC pc)
213: {
214: PC_LMVM *ctx = (PC_LMVM*)pc->data;
218: if (ctx->inactive) {
219: ISDestroy(&ctx->inactive);
220: }
221: if (pc->setupcalled) {
222: VecDestroy(&ctx->xwork);
223: VecDestroy(&ctx->ywork);
224: }
225: MatDestroy(&ctx->B);
226: PetscFree(pc->data);
227: return(0);
228: }
230: /*MC
231: PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the
232: underlying LMVM matrix can be access with the "-pc_lmvm_" prefix.
234: Level: intermediate
236: .seealso: PCCreate(), PCSetType(), PCType (for list of available types),
237: PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM()
238: M*/
239: PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc)
240: {
242: PC_LMVM *ctx;
245: PetscNewLog(pc,&ctx);
246: pc->data = (void*)ctx;
248: pc->ops->reset = PCReset_LMVM;
249: pc->ops->setup = PCSetUp_LMVM;
250: pc->ops->destroy = PCDestroy_LMVM;
251: pc->ops->view = PCView_LMVM;
252: pc->ops->apply = PCApply_LMVM;
253: pc->ops->setfromoptions = PCSetFromOptions_LMVM;
254: pc->ops->applysymmetricleft = NULL;
255: pc->ops->applysymmetricright = NULL;
256: pc->ops->applytranspose = NULL;
257: pc->ops->applyrichardson = NULL;
258: pc->ops->presolve = NULL;
259: pc->ops->postsolve = NULL;
261: PCSetReusePreconditioner(pc, PETSC_TRUE);
263: MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);
264: MatSetType(ctx->B, MATLMVMBFGS);
265: PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);
266: MatSetOptionsPrefix(ctx->B, "pc_lmvm_");
267: return(0);
268: }