ROL_MINRES.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#pragma once
#ifndef ROL_MINRES_HPP
#define ROL_MINRES_HPP

#include <array>
#include "ROL_Krylov.hpp"
#include "ROL_VectorClone.hpp"

namespace ROL {

namespace details {

template<typename Real>
class MINRES : public Krylov<Real> {

  using V  = Vector<Real>;
  using OP = LinearOperator<Real>;

private:

  // Givens rotation matrix elements
  Real resnorm_;
  int maxiter_;  
  bool useInexact_;
  std::array<Real,4> H_;
  std::array<Real,2> rhs_;

  VectorCloneMap<Real> clones_;      

  void givens( Real& c, Real& s, Real& r, Real a, Real b ) const {

    Real zero(0), one(1);

    if( b == zero ) {
      c = ( a >= zero ? one : -one );
      s = zero;
      r = std::abs(a);
    }
    else if( a == zero ) {
      c = zero;
      s = ( b >= zero ? one : -one );
      r = std::abs(b);
    }
    else if( std::abs(a) > std::abs(b) ) {
      auto t = b/a;
      auto u = std::copysign(std::sqrt(one+t*t),a);
      c = one/u;
      s = c*t;
      r = a*u;
    }
    else {
      auto t = a/b;
      auto u = std::copysign(std::sqrt(one+t*t),b);
      s = 1/u;
      c = s*t;
      r = b*u;
    }
  } // givens()

public:

  MINRES( Real absTol = 1.e-4, Real relTol = 1.e-2, unsigned maxit = 100, bool useInexact = false) :
    Krylov<Real>(absTol,relTol,maxit), useInexact_(useInexact),
    clones_("v_prev","v_curr","v_next","w_prev","w_curr","w_next") { }

  // Note: Preconditioner is not implemented
  virtual Real run( V &x, OP &A, const V &b, OP &M, int &iter, int &flag ) override {

    auto v_prev = clones_( x, "v_prev" );  v_prev->zero(); 
    auto v_curr = clones_( x, "v_curr" );  v_curr->set(b);
    auto v_next = clones_( x, "v_next" );  v_next->zero();
    auto w_prev = clones_( x, "w_prev" );  w_prev->zero();
    auto w_curr = clones_( x, "w_curr" );  w_curr->zero();
    auto w_next = clones_( x, "w_next" );  w_next->zero();

    Real c_prev{0}, s_prev{0}, c_curr{0}, s_curr{0}, c_next{0}, s_next{0};

    resnorm_ = v_curr->norm();    
    Real rtol = std::min(Krylov<Real>::getAbsoluteTolerance(),Krylov<Real>::getRelativeTolerance()*resnorm_);
    Real itol = std::sqrt(ROL_EPSILON<Real>());

    for( auto &e: H_ ) e = 0;

    rhs_[0] = resnorm_; rhs_[1] = 0;

    v_curr->scale(1.0/resnorm_);

    for( iter=0;  iter < (int)Krylov<Real>::getMaximumIteration(); iter++) {
      if ( useInexact_ ) {
        itol = rtol/((Real)Krylov<Real>::getMaximumIteration() * resnorm_);
      }

      if( resnorm_ < rtol ) break;

      A.apply( *v_next, *v_curr, itol );

      if( iter>0 ) v_next->axpy(-H_[1],*v_prev);    

      H_[2] = v_next->dot(*v_curr);

      v_next->axpy(-H_[2],*v_curr);

      H_[3] = v_next->norm();

      v_next->scale(1.0/H_[3]);

      // Rotation on H_[0] and H_[1]. 
      if( iter>1 ) {
        H_[0]  = s_prev*H_[1];
        H_[1] *= c_prev;
      }     

      // Rotation on H_[1] and H_[2]
      if( iter>0 ) {
        auto tmp = c_curr*H_[2]-s_curr*H_[1];
        H_[1] = c_curr*H_[1] + s_curr*H_[2];
        H_[2] = tmp;
      }

      // Form rotation coefficients
      givens( c_next, s_next, H_[2], H_[2], H_[3] );
      rhs_[1]  = -s_next*rhs_[0];
      rhs_[0] *=  c_next;

      w_next->set( *v_curr );

      if( iter>0 )  w_next->axpy( -H_[1], *w_curr );
      if( iter>1 )  w_next->axpy( -H_[0], *w_prev );

      w_next->scale(1.0/H_[2]);

      x.axpy( rhs_[0], *w_next );

      v_prev->set( *v_curr ); 
      v_curr->set( *v_next ); 
      w_prev->set( *w_curr ); 
      w_curr->set( *w_next ); 

      c_prev = c_curr;
      c_curr = c_next;
      s_prev = s_curr;
      s_curr = s_next;

      rhs_[0] = rhs_[1];

      H_[1] = H_[3];

      resnorm_ = std::abs( rhs_[1] );

    } // for (iter)

    if ( iter == (int)Krylov<Real>::getMaximumIteration() ) flag = 1;
    else iter++;

    return resnorm_;
  } // run()

}; // class MINRES

} // namespace details


using details::MINRES;


} // namespace ROL


#endif // ROL_MINRES_HPP