ROL
ROL_HelperFunctions.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #ifndef ROL_HELPERFUNCTIONS_HPP
50 #define ROL_HELPERFUNCTIONS_HPP
51 
52 #include "ROL_Vector.hpp"
53 #include "ROL_Objective.hpp"
54 #include "ROL_BoundConstraint.hpp"
55 #include "ROL_Secant.hpp"
56 #include "Teuchos_SerialDenseMatrix.hpp"
57 #include "Teuchos_SerialDenseVector.hpp"
58 #include "Teuchos_LAPACK.hpp"
59 
60 namespace ROL {
61 
62  template<class Real>
63  Teuchos::SerialDenseMatrix<int, Real> computeDenseHessian(Objective<Real> &obj, const Vector<Real> &x) {
64 
65  Real tol = std::sqrt(ROL_EPSILON<Real>());
66 
67  int dim = x.dimension();
68  Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
69 
70  ROL::Ptr<Vector<Real> > e = x.clone();
71  ROL::Ptr<Vector<Real> > h = x.dual().clone();
72 
73  for (int i=0; i<dim; i++) {
74  e = x.basis(i);
75  obj.hessVec(*h, *e, x, tol);
76  for (int j=0; j<dim; j++) {
77  e = x.basis(j);
78  //H(j,i) = e->dot(h->dual());
79  H(j,i) = e->apply(*h);
80  }
81  }
82 
83  return H;
84 
85  }
86 
87 
88  template<class Real>
89  Teuchos::SerialDenseMatrix<int, Real> computeScaledDenseHessian(Objective<Real> &obj, const Vector<Real> &x) {
90 
91  Real tol = std::sqrt(ROL_EPSILON<Real>());
92 
93  int dim = x.dimension();
94  Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
95 
96  ROL::Ptr<Vector<Real> > ei = x.clone();
97  ROL::Ptr<Vector<Real> > ej = x.dual().clone();
98  ROL::Ptr<Vector<Real> > h = x.dual().clone();
99 
100  for (int i=0; i<dim; i++) {
101  ei = ei->basis(i);
102  obj.hessVec(*h, *ei, x, tol);
103  for (int j=0; j<dim; j++) {
104  ej = ej->basis(j);
105  H(j,i) = ej->dot(*h);
106  }
107  }
108 
109  return H;
110 
111  }
112 
113 
114  template<class Real>
115  Teuchos::SerialDenseMatrix<int, Real> computeDotMatrix(const Vector<Real> &x) {
116 
117  int dim = x.dimension();
118  Teuchos::SerialDenseMatrix<int, Real> M(dim, dim);
119 
120  ROL::Ptr<Vector<Real> > ei = x.clone();
121  ROL::Ptr<Vector<Real> > ej = x.clone();
122 
123  for (int i=0; i<dim; i++) {
124  ei = x.basis(i);
125  for (int j=0; j<dim; j++) {
126  ej = x.basis(j);
127  M(j,i) = ej->dot(*ei);
128  }
129  }
130 
131  return M;
132 
133  }
134 
135  template<class Real>
136  std::vector<std::vector<Real> > computeEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
137 
138  Teuchos::LAPACK<int, Real> lapack;
139 
140  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
141 
142  char jobvl = 'N';
143  char jobvr = 'N';
144 
145  int n = mat.numRows();
146 
147  std::vector<Real> real(n, 0);
148  std::vector<Real> imag(n, 0);
149  std::vector<std::vector<Real> > eigenvals;
150 
151  Real* vl = 0;
152  Real* vr = 0;
153 
154  int ldvl = 1;
155  int ldvr = 1;
156 
157  int lwork = 4*n;
158 
159  std::vector<Real> work(lwork, 0);
160 
161  int info = 0;
162 
163  lapack.GEEV(jobvl, jobvr, n, &mymat(0,0), n, &real[0], &imag[0], vl, ldvl, vr, ldvr, &work[0], lwork, &info);
164 
165  eigenvals.push_back(real);
166  eigenvals.push_back(imag);
167 
168  return eigenvals;
169 
170  }
171 
172 
173  template<class Real>
174  std::vector<std::vector<Real> > computeGenEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & A,
175  const Teuchos::SerialDenseMatrix<int, Real> & B) {
176 
177  Teuchos::LAPACK<int, Real> lapack;
178 
179  Teuchos::SerialDenseMatrix<int, Real> myA(Teuchos::Copy, A);
180  Teuchos::SerialDenseMatrix<int, Real> myB(Teuchos::Copy, B);
181 
182  char jobvl = 'N';
183  char jobvr = 'N';
184 
185  int n = A.numRows();
186 
187  std::vector<Real> real(n, 0);
188  std::vector<Real> imag(n, 0);
189  std::vector<Real> beta(n, 0);
190  std::vector<std::vector<Real> > eigenvals;
191 
192  Real* vl = 0;
193  Real* vr = 0;
194 
195  int ldvl = 1;
196  int ldvr = 1;
197 
198  int lwork = 10*n;
199 
200  std::vector<Real> work(lwork, 0);
201 
202  int info = 0;
203 
204  lapack.GGEV(jobvl, jobvr, n, &myA(0,0), n, &myB(0,0), n, &real[0], &imag[0], &beta[0],
205  vl, ldvl, vr, ldvr, &work[0], lwork, &info);
206 
207  for (int i=0; i<n; i++) {
208  real[i] /= beta[i];
209  imag[i] /= beta[i];
210  }
211 
212  eigenvals.push_back(real);
213  eigenvals.push_back(imag);
214 
215  return eigenvals;
216 
217  }
218 
219 
220  template<class Real>
221  Teuchos::SerialDenseMatrix<int, Real> computeInverse(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
222 
223  Teuchos::LAPACK<int, Real> lapack;
224 
225  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
226 
227  int n = mat.numRows();
228 
229  std::vector<int> ipiv(n, 0);
230 
231  int lwork = 5*n;
232 
233  std::vector<Real> work(lwork, 0);
234 
235  int info = 0;
236 
237  lapack.GETRF(n, n, &mymat(0,0), n, &ipiv[0], &info);
238  lapack.GETRI(n, &mymat(0,0), n, &ipiv[0], &work[0], lwork, &info);
239 
240  return mymat;
241 
242  }
243 
244 
245 
246  template<class Real>
247  class ProjectedObjective : public Objective<Real> {
248  private:
249  ROL::Ptr<Objective<Real> > obj_;
250  ROL::Ptr<BoundConstraint<Real> > con_;
251  ROL::Ptr<Secant<Real> > secant_;
252 
253  ROL::Ptr<ROL::Vector<Real> > primalV_;
254  ROL::Ptr<ROL::Vector<Real> > dualV_;
256 
259  Real eps_;
260 
261  public:
263  ROL::Ptr<Secant<Real> > &secant,
264  bool useSecantPrecond = false,
265  bool useSecantHessVec = false,
266  Real eps = 0.0 )
267  : isInitialized_(false), useSecantPrecond_(useSecantPrecond),
268  useSecantHessVec_(useSecantHessVec), eps_(eps) {
269  obj_ = ROL::makePtrFromRef(obj);
270  con_ = ROL::makePtrFromRef(con);
271  secant_ = secant;
272  }
273 
274  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
275  obj_->update(x,flag,iter);
276  con_->update(x,flag,iter);
277  }
278 
279  Real value( const Vector<Real> &x, Real &tol ) {
280  return obj_->value(x,tol);
281  }
282 
283  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
284  obj_->gradient(g,x,tol);
285  }
286 
287  Real dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol ) {
288  return obj_->dirDeriv(x,d,tol);
289  }
290 
291  void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
292  if ( useSecantHessVec_ ) {
293  secant_->applyB( Hv, v );
294  }
295  else {
296  obj_->hessVec( Hv, v, x, tol );
297  }
298  }
299 
300  void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
301  if ( useSecantHessVec_ ) {
302  secant_->applyH(Hv,v);
303  }
304  else {
305  obj_->invHessVec(Hv,v,x,tol);
306  }
307  }
308 
309  void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
310  if ( useSecantPrecond_ ) {
311  secant_->applyH( Mv, v );
312  }
313  else {
314  obj_->precond( Mv, v, x, tol );
315  }
316  }
317 
329  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
330  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
331  if ( con_->isActivated() ) {
332  if (!isInitialized_) {
333  primalV_ = x.clone();
334  dualV_ = x.dual().clone();
335  isInitialized_ = true;
336  }
337  // Set vnew to v
338  primalV_->set(v);
339  // Remove elements of vnew corresponding to binding set
340  con_->pruneActive(*primalV_,d,p,eps_);
341  // Apply full Hessian to reduced vector
342  hessVec(Hv,*primalV_,x,tol);
343  // Remove elements of Hv corresponding to binding set
344  con_->pruneActive(Hv,d,p,eps_);
345  // Set vnew to v
346  primalV_->set(v);
347  // Remove Elements of vnew corresponding to complement of binding set
348  con_->pruneInactive(*primalV_,d,p,eps_);
349  dualV_->set(primalV_->dual());
350  con_->pruneInactive(*dualV_,d,p,eps_);
351  // Fill complement of binding set elements in Hp with v
352  Hv.plus(*dualV_);
353  }
354  else {
355  hessVec(Hv,v,x,tol);
356  }
357  }
358 
369  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
370  const Vector<Real> &x, Real &tol ) {
371  if ( con_->isActivated() ) {
372  if (!isInitialized_) {
373  primalV_ = x.clone();
374  dualV_ = x.dual().clone();
375  isInitialized_ = true;
376  }
377  // Set vnew to v
378  primalV_->set(v);
379  // Remove elements of vnew corresponding to binding set
380  con_->pruneActive(*primalV_,p,eps_);
381  // Apply full Hessian to reduced vector
382  hessVec(Hv,*primalV_,x,tol);
383  // Remove elements of Hv corresponding to binding set
384  con_->pruneActive(Hv,p,eps_);
385  // Set vnew to v
386  primalV_->set(v);
387  // Remove Elements of vnew corresponding to complement of binding set
388  con_->pruneInactive(*primalV_,p,eps_);
389  dualV_->set(primalV_->dual());
390  con_->pruneInactive(*dualV_,p,eps_);
391  // Fill complement of binding set elements in Hp with v
392  Hv.plus(*dualV_);
393  }
394  else {
395  hessVec(Hv,v,x,tol);
396  }
397  }
398 
410  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
411  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
412  if ( con_->isActivated() ) {
413  if (!isInitialized_) {
414  primalV_ = x.clone();
415  dualV_ = x.dual().clone();
416  isInitialized_ = true;
417  }
418  // Set vnew to v
419  dualV_->set(v);
420  // Remove elements of vnew corresponding to binding set
421  con_->pruneActive(*dualV_,d,p,eps_);
422  // Apply full Hessian to reduced vector
423  invHessVec(Hv,*dualV_,x,tol);
424  // Remove elements of Hv corresponding to binding set
425  con_->pruneActive(Hv,d,p,eps_);
426  // Set vnew to v
427  dualV_->set(v);
428  // Remove Elements of vnew corresponding to complement of binding set
429  con_->pruneInactive(*dualV_,d,p,eps_);
430  primalV_->set(dualV_->dual());
431  con_->pruneInactive(*primalV_,d,p,eps_);
432  // Fill complement of binding set elements in Hv with v
433  Hv.plus(*primalV_);
434  }
435  else {
436  invHessVec(Hv,v,x,tol);
437  }
438  }
439 
450  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
451  const Vector<Real> &x, Real &tol ) {
452  if ( con_->isActivated() ) {
453  if (!isInitialized_) {
454  primalV_ = x.clone();
455  dualV_ = x.dual().clone();
456  isInitialized_ = true;
457  }
458  // Set vnew to v
459  dualV_->set(v);
460  // Remove elements of vnew corresponding to binding set
461  con_->pruneActive(*dualV_,p,eps_);
462  // Apply full Hessian to reduced vector
463  invHessVec(Hv,*dualV_,x,tol);
464  // Remove elements of Hv corresponding to binding set
465  con_->pruneActive(Hv,p,eps_);
466  // Set vnew to v
467  dualV_->set(v);
468  // Remove Elements of vnew corresponding to complement of binding set
469  con_->pruneInactive(*dualV_,p,eps_);
470  primalV_->set(dualV_->dual());
471  con_->pruneInactive(*primalV_,p,eps_);
472  // Fill complement of binding set elements in Hv with v
473  Hv.plus(*primalV_);
474  }
475  else {
476  invHessVec(Hv,v,x,tol);
477  }
478  }
479 
491  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
492  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
493  if ( con_->isActivated() ) {
494  if (!isInitialized_) {
495  primalV_ = x.clone();
496  dualV_ = x.dual().clone();
497  isInitialized_ = true;
498  }
499  // Set vnew to v
500  dualV_->set(v);
501  // Remove elements of vnew corresponding to binding set
502  con_->pruneActive(*dualV_,d,p,eps_);
503  // Apply full Hessian to reduced vector
504  precond(Mv,*dualV_,x,tol);
505  // Remove elements of Mv corresponding to binding set
506  con_->pruneActive(Mv,d,p,eps_);
507  // Set vnew to v
508  dualV_->set(v);
509  // Remove Elements of vnew corresponding to complement of binding set
510  con_->pruneInactive(*dualV_,d,p,eps_);
511  primalV_->set(dualV_->dual());
512  con_->pruneInactive(*primalV_,d,p,eps_);
513  // Fill complement of binding set elements in Mv with v
514  Mv.plus(*primalV_);
515  }
516  else {
517  precond(Mv,v,x,tol);
518  }
519  }
520 
531  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
532  const Vector<Real> &x, Real &tol ) {
533  if ( con_->isActivated() ) {
534  if (!isInitialized_) {
535  primalV_ = x.clone();
536  dualV_ = x.dual().clone();
537  isInitialized_ = true;
538  }
539  // Set vnew to v
540  dualV_->set(v);
541  // Remove elements of vnew corresponding to binding set
542  con_->pruneActive(*dualV_,p,eps_);
543  // Apply full Hessian to reduced vector
544  precond(Mv,*dualV_,x,tol);
545  // Remove elements of Mv corresponding to binding set
546  con_->pruneActive(Mv,p,eps_);
547  // Set vnew to v
548  dualV_->set(v);
549  // Remove Elements of vnew corresponding to complement of binding set
550  con_->pruneInactive(*dualV_,p,eps_);
551  primalV_->set(dualV_->dual());
552  con_->pruneInactive(*primalV_,p,eps_);
553  // Fill complement of binding set elements in Mv with v
554  Mv.plus(*primalV_);
555  }
556  else {
557  precond(Mv,v,x,tol);
558  }
559  }
560 
561  void project( Vector<Real> &x ) {
562  con_->project(x);
563  }
564 
565  void pruneActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
566  con_->pruneActive(v,g,x,eps_);
567  }
568 
569  void pruneActive( Vector<Real> &v, const Vector<Real> &x ) {
570  con_->pruneActive(v,x,eps_);
571  }
572 
573  void pruneInactive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
574  con_->pruneInactive(v,g,x,eps_);
575  }
576 
577  void pruneInactive( Vector<Real> &v, const Vector<Real> &x ) {
578  con_->pruneInactive(v,x,eps_);
579  }
580 
581  bool isFeasible( const Vector<Real> &v ) {
582  return con_->isFeasible(v);
583  }
584 
585  bool isConActivated(void) {
586  return con_->isActivated();
587  }
588 
590  con_->computeProjectedStep(v,x);
591  }
592  };
593 
594 } // namespace ROL
595 
596 #endif
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:226
Real value(const Vector< Real > &x, Real &tol)
Compute value.
ROL::Ptr< Secant< Real > > secant_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
ROL::Ptr< ROL::Vector< Real > > primalV_
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:182
virtual void plus(const Vector &x)=0
Compute , where .
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Real dirDeriv(const Vector< Real > &x, const Vector< Real > &d, Real &tol)
Compute directional derivative.
void pruneInactive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
Teuchos::SerialDenseMatrix< int, Real > computeScaledDenseHessian(Objective< Real > &obj, const Vector< Real > &x)
ProjectedObjective(Objective< Real > &obj, BoundConstraint< Real > &con, ROL::Ptr< Secant< Real > > &secant, bool useSecantPrecond=false, bool useSecantHessVec=false, Real eps=0.0)
ROL::Ptr< Objective< Real > > obj_
bool isFeasible(const Vector< Real > &v)
void pruneActive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
void reducedPrecond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced preconditioner to a vector, v. The reduced preconditioner first removes elements of...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
void hessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Teuchos::SerialDenseMatrix< int, Real > computeDenseHessian(Objective< Real > &obj, const Vector< Real > &x)
std::vector< std::vector< Real > > computeGenEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &A, const Teuchos::SerialDenseMatrix< int, Real > &B)
void reducedInvHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced inverse Hessian to a vector, v. The reduced inverse Hessian first removes elements ...
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
ROL::Ptr< ROL::Vector< Real > > dualV_
std::vector< std::vector< Real > > computeEigenvalues(const Teuchos::SerialDenseMatrix< int, Real > &mat)
void reducedInvHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced inverse Hessian to a vector, v. The reduced inverse Hessian first removes elements ...
void pruneInactive(Vector< Real > &v, const Vector< Real > &x)
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
Provides the interface to apply upper and lower bound constraints.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Teuchos::SerialDenseMatrix< int, Real > computeDotMatrix(const Vector< Real > &x)
ROL::Ptr< BoundConstraint< Real > > con_
void computeProjectedStep(Vector< Real > &v, const Vector< Real > &x)
void reducedPrecond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced preconditioner to a vector, v. The reduced preconditioner first removes elements of...
void project(Vector< Real > &x)
void invHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply inverse Hessian approximation to vector.
Teuchos::SerialDenseMatrix< int, Real > computeInverse(const Teuchos::SerialDenseMatrix< int, Real > &mat)
void precond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply preconditioner to vector.
constexpr auto dim
void pruneActive(Vector< Real > &v, const Vector< Real > &x)