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  }
80  }
81 
82  return H;
83 
84  }
85 
86 
87  template<class Real>
88  Teuchos::SerialDenseMatrix<int, Real> computeScaledDenseHessian(Objective<Real> &obj, const Vector<Real> &x) {
89 
90  Real tol = std::sqrt(ROL_EPSILON<Real>());
91 
92  int dim = x.dimension();
93  Teuchos::SerialDenseMatrix<int, Real> H(dim, dim);
94 
95  ROL::Ptr<Vector<Real> > ei = x.clone();
96  ROL::Ptr<Vector<Real> > ej = x.dual().clone();
97  ROL::Ptr<Vector<Real> > h = x.dual().clone();
98 
99  for (int i=0; i<dim; i++) {
100  ei = ei->basis(i);
101  obj.hessVec(*h, *ei, x, tol);
102  for (int j=0; j<dim; j++) {
103  ej = ej->basis(j);
104  H(j,i) = ej->dot(*h);
105  }
106  }
107 
108  return H;
109 
110  }
111 
112 
113  template<class Real>
114  Teuchos::SerialDenseMatrix<int, Real> computeDotMatrix(const Vector<Real> &x) {
115 
116  int dim = x.dimension();
117  Teuchos::SerialDenseMatrix<int, Real> M(dim, dim);
118 
119  ROL::Ptr<Vector<Real> > ei = x.clone();
120  ROL::Ptr<Vector<Real> > ej = x.clone();
121 
122  for (int i=0; i<dim; i++) {
123  ei = x.basis(i);
124  for (int j=0; j<dim; j++) {
125  ej = x.basis(j);
126  M(j,i) = ej->dot(*ei);
127  }
128  }
129 
130  return M;
131 
132  }
133 
134  template<class Real>
135  std::vector<std::vector<Real> > computeEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
136 
137  Teuchos::LAPACK<int, Real> lapack;
138 
139  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
140 
141  char jobvl = 'N';
142  char jobvr = 'N';
143 
144  int n = mat.numRows();
145 
146  std::vector<Real> real(n, 0);
147  std::vector<Real> imag(n, 0);
148  std::vector<std::vector<Real> > eigenvals;
149 
150  Real* vl = 0;
151  Real* vr = 0;
152 
153  int ldvl = 1;
154  int ldvr = 1;
155 
156  int lwork = 4*n;
157 
158  std::vector<Real> work(lwork, 0);
159 
160  int info = 0;
161 
162  lapack.GEEV(jobvl, jobvr, n, &mymat(0,0), n, &real[0], &imag[0], vl, ldvl, vr, ldvr, &work[0], lwork, &info);
163 
164  eigenvals.push_back(real);
165  eigenvals.push_back(imag);
166 
167  return eigenvals;
168 
169  }
170 
171 
172  template<class Real>
173  std::vector<std::vector<Real> > computeGenEigenvalues(const Teuchos::SerialDenseMatrix<int, Real> & A,
174  const Teuchos::SerialDenseMatrix<int, Real> & B) {
175 
176  Teuchos::LAPACK<int, Real> lapack;
177 
178  Teuchos::SerialDenseMatrix<int, Real> myA(Teuchos::Copy, A);
179  Teuchos::SerialDenseMatrix<int, Real> myB(Teuchos::Copy, B);
180 
181  char jobvl = 'N';
182  char jobvr = 'N';
183 
184  int n = A.numRows();
185 
186  std::vector<Real> real(n, 0);
187  std::vector<Real> imag(n, 0);
188  std::vector<Real> beta(n, 0);
189  std::vector<std::vector<Real> > eigenvals;
190 
191  Real* vl = 0;
192  Real* vr = 0;
193 
194  int ldvl = 1;
195  int ldvr = 1;
196 
197  int lwork = 10*n;
198 
199  std::vector<Real> work(lwork, 0);
200 
201  int info = 0;
202 
203  lapack.GGEV(jobvl, jobvr, n, &myA(0,0), n, &myB(0,0), n, &real[0], &imag[0], &beta[0],
204  vl, ldvl, vr, ldvr, &work[0], lwork, &info);
205 
206  for (int i=0; i<n; i++) {
207  real[i] /= beta[i];
208  imag[i] /= beta[i];
209  }
210 
211  eigenvals.push_back(real);
212  eigenvals.push_back(imag);
213 
214  return eigenvals;
215 
216  }
217 
218 
219  template<class Real>
220  Teuchos::SerialDenseMatrix<int, Real> computeInverse(const Teuchos::SerialDenseMatrix<int, Real> & mat) {
221 
222  Teuchos::LAPACK<int, Real> lapack;
223 
224  Teuchos::SerialDenseMatrix<int, Real> mymat(Teuchos::Copy, mat);
225 
226  int n = mat.numRows();
227 
228  std::vector<int> ipiv(n, 0);
229 
230  int lwork = 5*n;
231 
232  std::vector<Real> work(lwork, 0);
233 
234  int info = 0;
235 
236  lapack.GETRF(n, n, &mymat(0,0), n, &ipiv[0], &info);
237  lapack.GETRI(n, &mymat(0,0), n, &ipiv[0], &work[0], lwork, &info);
238 
239  return mymat;
240 
241  }
242 
243 
244 
245  template<class Real>
246  class ProjectedObjective : public Objective<Real> {
247  private:
248  ROL::Ptr<Objective<Real> > obj_;
249  ROL::Ptr<BoundConstraint<Real> > con_;
250  ROL::Ptr<Secant<Real> > secant_;
251 
252  ROL::Ptr<ROL::Vector<Real> > primalV_;
253  ROL::Ptr<ROL::Vector<Real> > dualV_;
255 
258  Real eps_;
259 
260  public:
262  ROL::Ptr<Secant<Real> > &secant,
263  bool useSecantPrecond = false,
264  bool useSecantHessVec = false,
265  Real eps = 0.0 )
266  : isInitialized_(false), useSecantPrecond_(useSecantPrecond),
267  useSecantHessVec_(useSecantHessVec), eps_(eps) {
268  obj_ = ROL::makePtrFromRef(obj);
269  con_ = ROL::makePtrFromRef(con);
270  secant_ = secant;
271  }
272 
273  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
274  obj_->update(x,flag,iter);
275  con_->update(x,flag,iter);
276  }
277 
278  Real value( const Vector<Real> &x, Real &tol ) {
279  return obj_->value(x,tol);
280  }
281 
282  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
283  obj_->gradient(g,x,tol);
284  }
285 
286  Real dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol ) {
287  return obj_->dirDeriv(x,d,tol);
288  }
289 
290  void hessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
291  if ( useSecantHessVec_ ) {
292  secant_->applyB( Hv, v );
293  }
294  else {
295  obj_->hessVec( Hv, v, x, tol );
296  }
297  }
298 
299  void invHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
300  if ( useSecantHessVec_ ) {
301  secant_->applyH(Hv,v);
302  }
303  else {
304  obj_->invHessVec(Hv,v,x,tol);
305  }
306  }
307 
308  void precond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
309  if ( useSecantPrecond_ ) {
310  secant_->applyH( Mv, v );
311  }
312  else {
313  obj_->precond( Mv, v, x, tol );
314  }
315  }
316 
328  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
329  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
330  if ( con_->isActivated() ) {
331  if (!isInitialized_) {
332  primalV_ = x.clone();
333  dualV_ = x.dual().clone();
334  isInitialized_ = true;
335  }
336  // Set vnew to v
337  primalV_->set(v);
338  // Remove elements of vnew corresponding to binding set
339  con_->pruneActive(*primalV_,d,p,eps_);
340  // Apply full Hessian to reduced vector
341  hessVec(Hv,*primalV_,x,tol);
342  // Remove elements of Hv corresponding to binding set
343  con_->pruneActive(Hv,d,p,eps_);
344  // Set vnew to v
345  primalV_->set(v);
346  // Remove Elements of vnew corresponding to complement of binding set
347  con_->pruneInactive(*primalV_,d,p,eps_);
348  dualV_->set(primalV_->dual());
349  con_->pruneInactive(*dualV_,d,p,eps_);
350  // Fill complement of binding set elements in Hp with v
351  Hv.plus(*dualV_);
352  }
353  else {
354  hessVec(Hv,v,x,tol);
355  }
356  }
357 
368  void reducedHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
369  const Vector<Real> &x, Real &tol ) {
370  if ( con_->isActivated() ) {
371  if (!isInitialized_) {
372  primalV_ = x.clone();
373  dualV_ = x.dual().clone();
374  isInitialized_ = true;
375  }
376  // Set vnew to v
377  primalV_->set(v);
378  // Remove elements of vnew corresponding to binding set
379  con_->pruneActive(*primalV_,p,eps_);
380  // Apply full Hessian to reduced vector
381  hessVec(Hv,*primalV_,x,tol);
382  // Remove elements of Hv corresponding to binding set
383  con_->pruneActive(Hv,p,eps_);
384  // Set vnew to v
385  primalV_->set(v);
386  // Remove Elements of vnew corresponding to complement of binding set
387  con_->pruneInactive(*primalV_,p,eps_);
388  dualV_->set(primalV_->dual());
389  con_->pruneInactive(*dualV_,p,eps_);
390  // Fill complement of binding set elements in Hp with v
391  Hv.plus(*dualV_);
392  }
393  else {
394  hessVec(Hv,v,x,tol);
395  }
396  }
397 
409  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
410  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
411  if ( con_->isActivated() ) {
412  if (!isInitialized_) {
413  primalV_ = x.clone();
414  dualV_ = x.dual().clone();
415  isInitialized_ = true;
416  }
417  // Set vnew to v
418  dualV_->set(v);
419  // Remove elements of vnew corresponding to binding set
420  con_->pruneActive(*dualV_,d,p,eps_);
421  // Apply full Hessian to reduced vector
422  invHessVec(Hv,*dualV_,x,tol);
423  // Remove elements of Hv corresponding to binding set
424  con_->pruneActive(Hv,d,p,eps_);
425  // Set vnew to v
426  dualV_->set(v);
427  // Remove Elements of vnew corresponding to complement of binding set
428  con_->pruneInactive(*dualV_,d,p,eps_);
429  primalV_->set(dualV_->dual());
430  con_->pruneInactive(*primalV_,d,p,eps_);
431  // Fill complement of binding set elements in Hv with v
432  Hv.plus(*primalV_);
433  }
434  else {
435  invHessVec(Hv,v,x,tol);
436  }
437  }
438 
449  void reducedInvHessVec( Vector<Real> &Hv, const Vector<Real> &v, const Vector<Real> &p,
450  const Vector<Real> &x, Real &tol ) {
451  if ( con_->isActivated() ) {
452  if (!isInitialized_) {
453  primalV_ = x.clone();
454  dualV_ = x.dual().clone();
455  isInitialized_ = true;
456  }
457  // Set vnew to v
458  dualV_->set(v);
459  // Remove elements of vnew corresponding to binding set
460  con_->pruneActive(*dualV_,p,eps_);
461  // Apply full Hessian to reduced vector
462  invHessVec(Hv,*dualV_,x,tol);
463  // Remove elements of Hv corresponding to binding set
464  con_->pruneActive(Hv,p,eps_);
465  // Set vnew to v
466  dualV_->set(v);
467  // Remove Elements of vnew corresponding to complement of binding set
468  con_->pruneInactive(*dualV_,p,eps_);
469  primalV_->set(dualV_->dual());
470  con_->pruneInactive(*primalV_,p,eps_);
471  // Fill complement of binding set elements in Hv with v
472  Hv.plus(*primalV_);
473  }
474  else {
475  invHessVec(Hv,v,x,tol);
476  }
477  }
478 
490  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
491  const Vector<Real> &d, const Vector<Real> &x, Real &tol ) {
492  if ( con_->isActivated() ) {
493  if (!isInitialized_) {
494  primalV_ = x.clone();
495  dualV_ = x.dual().clone();
496  isInitialized_ = true;
497  }
498  // Set vnew to v
499  dualV_->set(v);
500  // Remove elements of vnew corresponding to binding set
501  con_->pruneActive(*dualV_,d,p,eps_);
502  // Apply full Hessian to reduced vector
503  precond(Mv,*dualV_,x,tol);
504  // Remove elements of Mv corresponding to binding set
505  con_->pruneActive(Mv,d,p,eps_);
506  // Set vnew to v
507  dualV_->set(v);
508  // Remove Elements of vnew corresponding to complement of binding set
509  con_->pruneInactive(*dualV_,d,p,eps_);
510  primalV_->set(dualV_->dual());
511  con_->pruneInactive(*primalV_,d,p,eps_);
512  // Fill complement of binding set elements in Mv with v
513  Mv.plus(*primalV_);
514  }
515  else {
516  precond(Mv,v,x,tol);
517  }
518  }
519 
530  void reducedPrecond( Vector<Real> &Mv, const Vector<Real> &v, const Vector<Real> &p,
531  const Vector<Real> &x, Real &tol ) {
532  if ( con_->isActivated() ) {
533  if (!isInitialized_) {
534  primalV_ = x.clone();
535  dualV_ = x.dual().clone();
536  isInitialized_ = true;
537  }
538  // Set vnew to v
539  dualV_->set(v);
540  // Remove elements of vnew corresponding to binding set
541  con_->pruneActive(*dualV_,p,eps_);
542  // Apply full Hessian to reduced vector
543  precond(Mv,*dualV_,x,tol);
544  // Remove elements of Mv corresponding to binding set
545  con_->pruneActive(Mv,p,eps_);
546  // Set vnew to v
547  dualV_->set(v);
548  // Remove Elements of vnew corresponding to complement of binding set
549  con_->pruneInactive(*dualV_,p,eps_);
550  primalV_->set(dualV_->dual());
551  con_->pruneInactive(*primalV_,p,eps_);
552  // Fill complement of binding set elements in Mv with v
553  Mv.plus(*primalV_);
554  }
555  else {
556  precond(Mv,v,x,tol);
557  }
558  }
559 
560  void project( Vector<Real> &x ) {
561  con_->project(x);
562  }
563 
564  void pruneActive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
565  con_->pruneActive(v,g,x,eps_);
566  }
567 
568  void pruneActive( Vector<Real> &v, const Vector<Real> &x ) {
569  con_->pruneActive(v,x,eps_);
570  }
571 
572  void pruneInactive( Vector<Real> &v, const Vector<Real> &g, const Vector<Real> &x ) {
573  con_->pruneInactive(v,g,x,eps_);
574  }
575 
576  void pruneInactive( Vector<Real> &v, const Vector<Real> &x ) {
577  con_->pruneInactive(v,x,eps_);
578  }
579 
580  bool isFeasible( const Vector<Real> &v ) {
581  return con_->isFeasible(v);
582  }
583 
584  bool isConActivated(void) {
585  return con_->isActivated();
586  }
587 
589  con_->computeProjectedStep(v,x);
590  }
591  };
592 
593 } // namespace ROL
594 
595 #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:70
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)