xref: /petsc/src/tao/bound/impls/bqnk/bqnk.c (revision 2e4af2ae5ea14e06d4fbd1ab1c02cabcc89ffdd5)
1 #include <../src/tao/bound/impls/bqnk/bqnk.h> /*I "petsctao.h" I*/ /*I "petscmat.h" I*/
2 #include <petscksp.h>
3 
4 static const char *BQNK_INIT[64] = {"constant", "direction"};
5 static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"};
6 static const char *BNK_AS[64] = {"none", "bertsekas"};
7 
8 static PetscErrorCode TaoBQNKComputeHessian(Tao tao)
9 {
10   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
11   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
12   PetscErrorCode ierr;
13   PetscReal      gnorm2, delta;
14 
15   PetscFunctionBegin;
16   /* Alias the LMVM matrix into the TAO hessian */
17   if (tao->hessian) {
18     ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr);
19   }
20   if (tao->hessian_pre) {
21     ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr);
22   }
23   ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr);
24   tao->hessian = bqnk->B;
25   ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr);
26   tao->hessian_pre = bqnk->B;
27   /* Update the Hessian with the latest solution */
28   if (bqnk->is_spd) {
29     gnorm2 = bnk->gnorm*bnk->gnorm;
30     if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
31     if (bnk->f == 0.0) {
32       delta = 2.0 / gnorm2;
33     } else {
34       delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2;
35     }
36     ierr = MatLMVMSymBroydenSetDelta(bqnk->B, delta);CHKERRQ(ierr);
37   }
38   ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
39   ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr);
40   /* Prepare the reduced sub-matrices for the inactive set */
41   ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr);
42   if (bnk->active_idx) {
43     ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr);
44     ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr);
45   } else {
46     ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr);
47     bnk->H_inactive = tao->hessian;
48     ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr);
49   }
50   ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr);
51   ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr);
52   bnk->Hpre_inactive = bnk->H_inactive;
53   PetscFunctionReturn(0);
54 }
55 
56 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
57 {
58   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
59   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
60   PetscErrorCode ierr;
61 
62   PetscFunctionBegin;
63   ierr = TaoBNKComputeStep(tao, shift, ksp_reason, step_type);CHKERRQ(ierr);
64   if (*ksp_reason < 0) {
65     /* Krylov solver failed to converge so reset the LMVM matrix */
66     ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr);
67     ierr = MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr);
68   }
69   PetscFunctionReturn(0);
70 }
71 
72 PetscErrorCode TaoSolve_BQNK(Tao tao)
73 {
74   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
75   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
76   Mat_LMVM       *lmvm = (Mat_LMVM*)bqnk->B->data;
77   Mat_LMVM       *J0;
78   Mat_SymBrdn    *diag_ctx;
79   PetscBool      flg = PETSC_FALSE;
80   PetscErrorCode ierr;
81 
82   PetscFunctionBegin;
83   if (!tao->recycle) {
84     ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr);
85     lmvm->nresets = 0;
86     if (lmvm->J0) {
87       ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg);CHKERRQ(ierr);
88       if (flg) {
89         J0 = (Mat_LMVM*)lmvm->J0->data;
90         J0->nresets = 0;
91       }
92     }
93     flg = PETSC_FALSE;
94     ierr = PetscObjectTypeCompareAny((PetscObject)bqnk->B, &flg, MATLMVMSYMBROYDEN, MATLMVMSYMBADBROYDEN, MATLMVMBFGS, MATLMVMDFP, "");
95     if (flg) {
96       diag_ctx = (Mat_SymBrdn*)lmvm->ctx;
97       J0 = (Mat_LMVM*)diag_ctx->D->data;
98       J0->nresets = 0;
99     }
100   }
101   ierr = (*bqnk->solve)(tao);CHKERRQ(ierr);
102   PetscFunctionReturn(0);
103 }
104 
105 PetscErrorCode TaoSetUp_BQNK(Tao tao)
106 {
107   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
108   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
109   PetscErrorCode ierr;
110   PetscInt       n, N;
111   PetscBool      is_lmvm, is_sym, is_spd;
112 
113   PetscFunctionBegin;
114   ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr);
115   ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
116   ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
117   ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr);
118   ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr);
119   ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr);
120   if (!is_lmvm) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type");
121   ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr);
122   if (!is_sym) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric");
123   ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr);
124   ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr);
125   ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr);
126   ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr);
127   PetscFunctionReturn(0);
128 }
129 
130 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao)
131 {
132   TAO_BNK        *bnk = (TAO_BNK *)tao->data;
133   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
134   PetscErrorCode ierr;
135   KSPType        ksp_type;
136 
137   PetscFunctionBegin;
138   ierr = PetscOptionsHead(PetscOptionsObject,"Quasi-Newton-Krylov method for bound constrained optimization");CHKERRQ(ierr);
139   ierr = PetscOptionsEList("-tao_bqnk_init_type", "radius initialization type", "", BQNK_INIT, BQNK_INIT_TYPES, BQNK_INIT[bnk->init_type], &bnk->init_type, NULL);CHKERRQ(ierr);
140   ierr = PetscOptionsEList("-tao_bqnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL);CHKERRQ(ierr);
141   ierr = PetscOptionsEList("-tao_bqnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL);CHKERRQ(ierr);
142   ierr = PetscOptionsReal("-tao_bqnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);CHKERRQ(ierr);
143   ierr = PetscOptionsReal("-tao_bqnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);CHKERRQ(ierr);
144   ierr = PetscOptionsReal("-tao_bqnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);CHKERRQ(ierr);
145   ierr = PetscOptionsReal("-tao_bqnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);CHKERRQ(ierr);
146   ierr = PetscOptionsReal("-tao_bqnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);CHKERRQ(ierr);
147   ierr = PetscOptionsReal("-tao_bqnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);CHKERRQ(ierr);
148   ierr = PetscOptionsReal("-tao_bqnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);CHKERRQ(ierr);
149   ierr = PetscOptionsReal("-tao_bqnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);CHKERRQ(ierr);
150   ierr = PetscOptionsReal("-tao_bqnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);CHKERRQ(ierr);
151   ierr = PetscOptionsReal("-tao_bqnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);CHKERRQ(ierr);
152   ierr = PetscOptionsReal("-tao_bqnk_eta1", "(developer) threshold for rejecting step (-tao_bqnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);CHKERRQ(ierr);
153   ierr = PetscOptionsReal("-tao_bqnk_eta2", "(developer) threshold for accepting marginal step (-tao_bqnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);CHKERRQ(ierr);
154   ierr = PetscOptionsReal("-tao_bqnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bqnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);CHKERRQ(ierr);
155   ierr = PetscOptionsReal("-tao_bqnk_eta4", "(developer) threshold for accepting good step (-tao_bqnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);CHKERRQ(ierr);
156   ierr = PetscOptionsReal("-tao_bqnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bqnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);CHKERRQ(ierr);
157   ierr = PetscOptionsReal("-tao_bqnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bqnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);CHKERRQ(ierr);
158   ierr = PetscOptionsReal("-tao_bqnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);CHKERRQ(ierr);
159   ierr = PetscOptionsReal("-tao_bqnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);CHKERRQ(ierr);
160   ierr = PetscOptionsReal("-tao_bqnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bqnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);CHKERRQ(ierr);
161   ierr = PetscOptionsReal("-tao_bqnk_nu1", "(developer) threshold for small line-search step length (-tao_bqnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);CHKERRQ(ierr);
162   ierr = PetscOptionsReal("-tao_bqnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bqnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);CHKERRQ(ierr);
163   ierr = PetscOptionsReal("-tao_bqnk_nu3", "(developer) threshold for large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);CHKERRQ(ierr);
164   ierr = PetscOptionsReal("-tao_bqnk_nu4", "(developer) threshold for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);CHKERRQ(ierr);
165   ierr = PetscOptionsReal("-tao_bqnk_omega1", "(developer) radius reduction factor for very small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega1, &bnk->omega1,NULL);CHKERRQ(ierr);
166   ierr = PetscOptionsReal("-tao_bqnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bqnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);CHKERRQ(ierr);
167   ierr = PetscOptionsReal("-tao_bqnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bqnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);CHKERRQ(ierr);
168   ierr = PetscOptionsReal("-tao_bqnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);CHKERRQ(ierr);
169   ierr = PetscOptionsReal("-tao_bqnk_omega5", "(developer) radius increase factor for very large line-search step length (-tao_bqnk_update_type step)", "", bnk->omega5, &bnk->omega5,NULL);CHKERRQ(ierr);
170   ierr = PetscOptionsReal("-tao_bqnk_mu1", "(developer) threshold for accepting very good step (-tao_bqnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);CHKERRQ(ierr);
171   ierr = PetscOptionsReal("-tao_bqnk_mu2", "(developer) threshold for accepting good step (-tao_bqnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);CHKERRQ(ierr);
172   ierr = PetscOptionsReal("-tao_bqnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);CHKERRQ(ierr);
173   ierr = PetscOptionsReal("-tao_bqnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bqnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);CHKERRQ(ierr);
174   ierr = PetscOptionsReal("-tao_bqnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);CHKERRQ(ierr);
175   ierr = PetscOptionsReal("-tao_bqnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bqnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);CHKERRQ(ierr);
176   ierr = PetscOptionsReal("-tao_bqnk_theta", "(developer) trust region interpolation factor (-tao_bqnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);CHKERRQ(ierr);
177   ierr = PetscOptionsReal("-tao_bqnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);CHKERRQ(ierr);
178   ierr = PetscOptionsReal("-tao_bqnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);CHKERRQ(ierr);
179   ierr = PetscOptionsReal("-tao_bqnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);CHKERRQ(ierr);
180   ierr = PetscOptionsReal("-tao_bqnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);CHKERRQ(ierr);
181   ierr = PetscOptionsReal("-tao_bqnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);CHKERRQ(ierr);
182   ierr = PetscOptionsInt("-tao_bqnk_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);CHKERRQ(ierr);
183   ierr = PetscOptionsTail();CHKERRQ(ierr);
184   ierr = TaoSetFromOptions(bnk->bncg);CHKERRQ(ierr);
185   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
186   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
187   ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr);
188   ierr = PetscStrcmp(ksp_type,KSPNASH,&bnk->is_nash);CHKERRQ(ierr);
189   ierr = PetscStrcmp(ksp_type,KSPSTCG,&bnk->is_stcg);CHKERRQ(ierr);
190   ierr = PetscStrcmp(ksp_type,KSPGLTR,&bnk->is_gltr);CHKERRQ(ierr);
191   if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION;
192   ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr);
193   ierr = MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);CHKERRQ(ierr);
194   PetscFunctionReturn(0);
195 }
196 
197 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer)
198 {
199   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
200   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
201   PetscErrorCode ierr;
202   PetscBool      isascii;
203 
204   PetscFunctionBegin;
205   ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr);
206   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
207   if (isascii) {
208     ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
209     ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr);
210     ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
211   }
212   PetscFunctionReturn(0);
213 }
214 
215 static PetscErrorCode TaoDestroy_BQNK(Tao tao)
216 {
217   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
218   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
219   PetscErrorCode ierr;
220 
221   PetscFunctionBegin;
222   ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr);
223   ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr);
224   ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr);
225   ierr = PetscFree(bnk->ctx);CHKERRQ(ierr);
226   ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr);
227   PetscFunctionReturn(0);
228 }
229 
230 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao)
231 {
232   TAO_BNK        *bnk;
233   TAO_BQNK       *bqnk;
234   PetscErrorCode ierr;
235 
236   PetscFunctionBegin;
237   ierr = TaoCreate_BNK(tao);CHKERRQ(ierr);
238   ierr = KSPSetOptionsPrefix(tao->ksp, "tao_bqnk_");CHKERRQ(ierr);
239   tao->ops->solve = TaoSolve_BQNK;
240   tao->ops->setfromoptions = TaoSetFromOptions_BQNK;
241   tao->ops->destroy = TaoDestroy_BQNK;
242   tao->ops->view = TaoView_BQNK;
243   tao->ops->setup = TaoSetUp_BQNK;
244 
245   bnk = (TAO_BNK *)tao->data;
246   bnk->computehessian = TaoBQNKComputeHessian;
247   bnk->computestep = TaoBQNKComputeStep;
248   bnk->init_type = BNK_INIT_DIRECTION;
249 
250   ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr);
251   bnk->ctx = (void*)bqnk;
252   bqnk->is_spd = PETSC_TRUE;
253 
254   ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr);
255   ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr);
256   ierr = MatSetOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr);
257   ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr);
258   PetscFunctionReturn(0);
259 }
260 
261 /*@
262    TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid
263    only for quasi-Newton family of methods.
264 
265    Input Parameters:
266 .  tao - Tao solver context
267 
268    Output Parameters:
269 .  B - LMVM matrix
270 
271    Level: advanced
272 
273 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoSetLMVMMatrix()
274 @*/
275 PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B)
276 {
277   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
278   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
279   PetscErrorCode ierr;
280   PetscBool      flg = PETSC_FALSE;
281 
282   PetscFunctionBegin;
283   ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr);
284   if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
285   *B = bqnk->B;
286   PetscFunctionReturn(0);
287 }
288 
289 /*@
290    TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid
291    only for quasi-Newton family of methods.
292 
293    QN family of methods create their own LMVM matrices and users who wish to
294    manipulate this matrix should use TaoGetLMVMMatrix() instead.
295 
296    Input Parameters:
297 +  tao - Tao solver context
298 -  B - LMVM matrix
299 
300    Level: advanced
301 
302 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoGetLMVMMatrix()
303 @*/
304 PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B)
305 {
306   TAO_BNK        *bnk = (TAO_BNK*)tao->data;
307   TAO_BQNK       *bqnk = (TAO_BQNK*)bnk->ctx;
308   PetscErrorCode ierr;
309   PetscBool      flg = PETSC_FALSE;
310 
311   PetscFunctionBegin;
312   ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr);
313   if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms");
314   ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg);CHKERRQ(ierr);
315   if (!flg) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix");
316   if (bqnk->B) {
317     ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr);
318   }
319   ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr);
320   bqnk->B = B;
321   PetscFunctionReturn(0);
322 }
323