text: #ifndef VIENNACL_LINALG_LANCZOS_HPP_
  #define VIENNACL_LINALG_LANCZOS_HPP_
  
  /* =========================================================================
     Copyright (c) 2010-2014, Institute for Microelectronics,
                              Institute for Analysis and Scientific Computing,
                              TU Wien.
     Portions of this software are copyright by UChicago Argonne, LLC.
  
                              -----------------
                    ViennaCL - The Vienna Computing Library
                              -----------------
  
     Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
  
     (A list of authors and contributors can be found in the PDF manual)
  
     License:         MIT (X11), see file LICENSE in the base directory
  ============================================================================= */
  
  /** @file viennacl/linalg/lanczos.hpp
  *   @brief Generic interface for the Lanczos algorithm.
  *
  *   Contributed by Guenther Mader and Astrid Rupp.
  */
  
  #include <cmath>
  #include <vector>
  #include "viennacl/vector.hpp"
  #include "viennacl/compressed_matrix.hpp"
  #include "viennacl/linalg/prod.hpp"
  #include "viennacl/linalg/inner_prod.hpp"
  #include "viennacl/linalg/norm_2.hpp"
  #include "viennacl/io/matrix_market.hpp"
  #include "viennacl/linalg/bisect.hpp"
  #include <boost/random.hpp>
  #include <boost/random/mersenne_twister.hpp>
  #include <boost/numeric/ublas/matrix.hpp>
  #include <boost/numeric/ublas/matrix_proxy.hpp>
  #include <boost/numeric/ublas/matrix_expression.hpp>
  #include <boost/numeric/ublas/matrix_sparse.hpp>
  #include <boost/numeric/ublas/vector.hpp>
  #include <boost/numeric/ublas/operation.hpp>
  #include <boost/numeric/ublas/vector_expression.hpp>
  #include <boost/numeric/ublas/io.hpp>
  
  namespace viennacl
  {
  namespace linalg
  {
  
  /** @brief A tag for the lanczos algorithm.
  */
  class lanczos_tag
  {
  public:
  
    enum
    {
      partial_reorthogonalization = 0,
      full_reorthogonalization,
      no_reorthogonalization
    };
  
    /** @brief The constructor
    *
    * @param factor                 Exponent of epsilon - tolerance for batches of Reorthogonalization
    * @param numeig                 Number of eigenvalues to be returned
    * @param met                    Method for Lanczos-Algorithm: 0 for partial Reorthogonalization, 1 for full Reorthogonalization and 2 for Lanczos without Reorthogonalization
    * @param krylov                 Maximum krylov-space size
    */
  
    lanczos_tag(double factor = 0.75,
                vcl_size_t numeig = 10,
                int met = 0,
                vcl_size_t krylov = 100) : factor_(factor), num_eigenvalues_(numeig), method_(met), krylov_size_(krylov) {}
  
    /** @brief Sets the number of eigenvalues */
    void num_eigenvalues(vcl_size_t numeig){ num_eigenvalues_ = numeig; }
  
      /** @brief Returns the number of eigenvalues */
    vcl_size_t num_eigenvalues() const { return num_eigenvalues_; }
  
      /** @brief Sets the exponent of epsilon */
    void factor(double fct) { factor_ = fct; }
  
    /** @brief Returns the exponent */
    double factor() const { return factor_; }
  
    /** @brief Sets the size of the kylov space */
    void krylov_size(vcl_size_t max) { krylov_size_ = max; }
  
    /** @brief Returns the size of the kylov space */
    vcl_size_t  krylov_size() const { return krylov_size_; }
  
    /** @brief Sets the reorthogonalization method */
    void method(int met){ method_ = met; }
  
    /** @brief Returns the reorthogonalization method */
    int method() const { return method_; }
  
  
  private:
    double factor_;
    vcl_size_t num_eigenvalues_;
    int method_; // see enum defined above for possible values
    vcl_size_t krylov_size_;
  };
  
  
  namespace detail
  {
    /**
    *   @brief Implementation of the Lanczos PRO algorithm
    *
    *   @param A            The system matrix
    *   @param r            Random start vector
    *   @param size         Size of krylov-space
    *   @param tag          Lanczos_tag with several options for the algorithm
    *   @return             Returns the eigenvalues (number of eigenvalues equals size of krylov-space)
    */
  
    template< typename MatrixT, typename VectorT >
    std::vector<
            typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
            >
    lanczosPRO (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag const & tag)
    {
      typedef typename viennacl::result_of::value_type<MatrixT>::type        ScalarType;
      typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
  
  
      // generation of some random numbers, used for lanczos PRO algorithm
      boost::mt11213b mt;
      boost::normal_distribution<CPU_ScalarType> N(0, 1);
      boost::bernoulli_distribution<CPU_ScalarType> B(0.5);
      boost::triangle_distribution<CPU_ScalarType> T(-1, 0, 1);
  
      boost::variate_generator<boost::mt11213b&, boost::normal_distribution<CPU_ScalarType> >     get_N(mt, N);
      boost::variate_generator<boost::mt11213b&, boost::bernoulli_distribution<CPU_ScalarType> >  get_B(mt, B);
      boost::variate_generator<boost::mt11213b&, boost::triangle_distribution<CPU_ScalarType> >   get_T(mt, T);
  
  
      long i, k, retry, reorths;
      std::vector<long> l_bound(size/2), u_bound(size/2);
      bool second_step;
      CPU_ScalarType squ_eps, eta, temp, eps, retry_th;
      vcl_size_t n = r.size();
      std::vector< std::vector<CPU_ScalarType> > w(2, std::vector<CPU_ScalarType>(size));
      CPU_ScalarType cpu_beta;
  
      boost::numeric::ublas::vector<CPU_ScalarType> s(n);
  
      VectorT t(n);
      CPU_ScalarType inner_rt;
      ScalarType vcl_beta;
      ScalarType vcl_alpha;
      std::vector<CPU_ScalarType> alphas, betas;
      boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
  
      second_step = false;
      eps = std::numeric_limits<CPU_ScalarType>::epsilon();
      squ_eps = std::sqrt(eps);
      retry_th = 1e-2;
      eta = std::exp(std::log(eps) * tag.factor());
      reorths = 0;
      retry = 0;
  
      vcl_beta = viennacl::linalg::norm_2(r);
  
      r /= vcl_beta;
  
      detail::copy_vec_to_vec(r,s);
      boost::numeric::ublas::column(Q, 0) = s;
  
      VectorT u = viennacl::linalg::prod(A, r);
      vcl_alpha = viennacl::linalg::inner_prod(u, r);
      alphas.push_back(vcl_alpha);
      w[0][0] = 1;
      betas.push_back(vcl_beta);
  
      long batches = 0;
      for (i = 1;i < static_cast<long>(size); i++)
      {
        r = u - vcl_alpha * r;
        vcl_beta = viennacl::linalg::norm_2(r);
  
        betas.push_back(vcl_beta);
        r = r / vcl_beta;
  
        vcl_size_t index = vcl_size_t(i % 2);
        w[index][vcl_size_t(i)] = 1;
        k = (i + 1) % 2;
        w[index][0] = (betas[1] * w[vcl_size_t(k)][1] + (alphas[0] - vcl_alpha) * w[vcl_size_t(k)][0] - betas[vcl_size_t(i) - 1] * w[index][0]) / vcl_beta + eps * 0.3 * get_N() * (betas[1] + vcl_beta);
  
        for (vcl_size_t j = 1; j < vcl_size_t(i - 1); j++)
        {
                w[index][j] = (betas[j + 1] * w[vcl_size_t(k)][j + 1] + (alphas[j] - vcl_alpha) * w[vcl_size_t(k)][j] + betas[j] * w[vcl_size_t(k)][j - 1] - betas[vcl_size_t(i) - 1] * w[index][j]) / vcl_beta + eps * 0.3 * get_N() * (betas[j + 1] + vcl_beta);
        }
        w[index][vcl_size_t(i) - 1] = 0.6 * eps * CPU_ScalarType(n) * get_N() * betas[1] / vcl_beta;
  
        if (second_step)
        {
          for (vcl_size_t j = 0; j < vcl_size_t(batches); j++)
          {
            l_bound[vcl_size_t(j)]++;
            u_bound[vcl_size_t(j)]--;
  
            for (k = l_bound[j];k < u_bound[j];k++)
            {
              detail::copy_vec_to_vec(boost::numeric::ublas::column(Q, vcl_size_t(k)), t);
              inner_rt = viennacl::linalg::inner_prod(r,t);
              r = r - inner_rt * t;
              w[index][vcl_size_t(k)] = 1.5 * eps * get_N();
              reorths++;
            }
          }
          temp = viennacl::linalg::norm_2(r);
          r = r / temp;
          vcl_beta = vcl_beta * temp;
          second_step = false;
        }
        batches = 0;
  
        for (vcl_size_t j = 0; j < vcl_size_t(i); j++)
        {
          if (std::fabs(w[index][j]) >= squ_eps)
          {
            detail::copy_vec_to_vec(boost::numeric::ublas::column(Q, j), t);
            inner_rt = viennacl::linalg::inner_prod(r,t);
            r = r - inner_rt * t;
            w[index][j] = 1.5 * eps * get_N();
            k = long(j) - 1;
            reorths++;
            while (k >= 0 && std::fabs(w[index][vcl_size_t(k)]) > eta)
            {
              detail::copy_vec_to_vec(boost::numeric::ublas::column(Q, vcl_size_t(k)), t);
              inner_rt = viennacl::linalg::inner_prod(r,t);
              r = r - inner_rt * t;
              w[index][vcl_size_t(k)] = 1.5 * eps * get_N();
              k--;
              reorths++;
            }
            l_bound[vcl_size_t(batches)] = k + 1;
            k = long(j) + 1;
  
            while (k < i && std::fabs(w[index][vcl_size_t(k)]) > eta)
            {
              detail::copy_vec_to_vec(boost::numeric::ublas::column(Q, vcl_size_t(k)), t);
              inner_rt = viennacl::linalg::inner_prod(r,t);
              r = r - inner_rt * t;
              w[index][vcl_size_t(k)] = 1.5 * eps * get_N();
              k++;
              reorths++;
            }
            u_bound[vcl_size_t(batches)] = k - 1;
            batches++;
            j = vcl_size_t(k);
          }
        }
  
        if (batches > 0)
        {
          temp = viennacl::linalg::norm_2(r);
          r = r / temp;
          vcl_beta = vcl_beta * temp;
          second_step = true;
  
          while (temp < retry_th)
          {
            for (vcl_size_t j = 0; j < vcl_size_t(i); j++)
            {
              detail::copy_vec_to_vec(boost::numeric::ublas::column(Q, vcl_size_t(k)), t);
              inner_rt = viennacl::linalg::inner_prod(r,t);
              r = r - inner_rt * t;
              reorths++;
            }
            retry++;
            temp = viennacl::linalg::norm_2(r);
            r = r / temp;
            vcl_beta = vcl_beta * temp;
          }
        }
  
        detail::copy_vec_to_vec(r,s);
        boost::numeric::ublas::column(Q, vcl_size_t(i)) = s;
  
        cpu_beta = vcl_beta;
        s = - cpu_beta * boost::numeric::ublas::column(Q, vcl_size_t(i - 1));
        detail::copy_vec_to_vec(s, u);
        u += viennacl::linalg::prod(A, r);
        vcl_alpha = viennacl::linalg::inner_prod(u, r);
        alphas.push_back(vcl_alpha);
      }
  
      return bisect(alphas, betas);
    }
  
  
    /**
    *   @brief Implementation of the lanczos algorithm without reorthogonalization
    *
    *   @param A            The system matrix
    *   @param r            Random start vector
    *   @param size         Size of krylov-space
    *   @return             Returns the eigenvalues (number of eigenvalues equals size of krylov-space)
    */
    template<typename MatrixT, typename VectorT>
    std::vector<
            typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
            >
    lanczos (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag)
    {
      typedef typename viennacl::result_of::value_type<MatrixT>::type        ScalarType;
      typedef typename viennacl::result_of::cpu_value_type<ScalarType>::type    CPU_ScalarType;
  
      ScalarType vcl_beta;
      ScalarType vcl_alpha;
      std::vector<CPU_ScalarType> alphas, betas;
      CPU_ScalarType norm;
      vcl_size_t n = r.size();
      VectorT u(n), t(n);
      boost::numeric::ublas::vector<CPU_ScalarType> s(r.size()), u_zero(n), q(n);
      boost::numeric::ublas::matrix<CPU_ScalarType> Q(n, size);
  
      u_zero = boost::numeric::ublas::zero_vector<CPU_ScalarType>(n);
      detail::copy_vec_to_vec(u_zero, u);
      norm = norm_2(r);
  
      for (vcl_size_t i = 0;i < size; i++)
      {
        r /= norm;
        vcl_beta = norm;
  
        detail::copy_vec_to_vec(r,s);
        boost::numeric::ublas::column(Q, i) = s;
  
        u += prod(A, r);
        vcl_alpha = inner_prod(u, r);
        r = u - vcl_alpha * r;
        norm = norm_2(r);
  
        q = boost::numeric::ublas::column(Q, i);
        detail::copy_vec_to_vec(q, t);
  
        u = - norm * t;
        alphas.push_back(vcl_alpha);
        betas.push_back(vcl_beta);
        s.clear();
      }
  
      return bisect(alphas, betas);
    }
  
    /**
    *   @brief Implementation of the Lanczos FRO algorithm
    *
    *   @param A            The system matrix
    *   @param r            Random start vector
    *   @param size         Size of krylov-space
    *   @return             Returns the eigenvalues (number of eigenvalues equals size of krylov-space)
    */
    template< typename MatrixT, typename VectorT >
    std::vector<
            typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type
            >
    lanczosFRO (MatrixT const& A, VectorT & r, vcl_size_t size, lanczos_tag)
    {
      typedef typename viennacl::result_of::value_type<MatrixT>::type            NumericType;
      typedef typename viennacl::result_of::cpu_value_type<NumericType>::type    CPU_NumericType;
  
      CPU_NumericType temp;
      CPU_NumericType norm;
      NumericType vcl_beta;
      NumericType vcl_alpha;
      std::vector<CPU_NumericType> alphas, betas;
      vcl_size_t n = r.size();
      VectorT u(n), t(n);
      NumericType inner_rt;
      boost::numeric::ublas::vector<CPU_NumericType> u_zero(n), s(r.size()), q(n);
      boost::numeric::ublas::matrix<CPU_NumericType> Q(n, size);
  
      long reorths = 0;
      norm = norm_2(r);
  
  
      for (vcl_size_t i = 0; i < size; i++)
      {
        r /= norm;
  
        for (vcl_size_t j = 0; j < i; j++)
        {
          q = boost::numeric::ublas::column(Q, j);
          detail::copy_vec_to_vec(q, t);
          inner_rt = viennacl::linalg::inner_prod(r,t);
          r = r - inner_rt * t;
          reorths++;
        }
        temp = viennacl::linalg::norm_2(r);
        r = r / temp;
        vcl_beta = temp * norm;
        detail::copy_vec_to_vec(r,s);
        boost::numeric::ublas::column(Q, i) = s;
  
        u += viennacl::linalg::prod(A, r);
        vcl_alpha = viennacl::linalg::inner_prod(u, r);
        r = u - vcl_alpha * r;
        norm = viennacl::linalg::norm_2(r);
        q = boost::numeric::ublas::column(Q, i);
        detail::copy_vec_to_vec(q, t);
        u = - norm * t;
        alphas.push_back(vcl_alpha);
        betas.push_back(vcl_beta);
      }
  
      return bisect(alphas, betas);
    }
  
  } // end namespace detail
  
  /**
  *   @brief Implementation of the calculation of eigenvalues using lanczos
  *
  *   @param matrix        The system matrix
  *   @param tag           Tag with several options for the lanczos algorithm
  *   @return              Returns the n largest eigenvalues (n defined in the lanczos_tag)
  */
  template<typename MatrixT>
  std::vector< typename viennacl::result_of::cpu_value_type<typename MatrixT::value_type>::type >
  eig(MatrixT const & matrix, lanczos_tag const & tag)
  {
    typedef typename viennacl::result_of::value_type<MatrixT>::type           NumericType;
    typedef typename viennacl::result_of::cpu_value_type<NumericType>::type   CPU_NumericType;
    typedef typename viennacl::result_of::vector_for_matrix<MatrixT>::type    VectorT;
  
    boost::mt11213b mt;
    boost::normal_distribution<CPU_NumericType>    N(0, 1);
    boost::bernoulli_distribution<CPU_NumericType> B(0.5);
    boost::triangle_distribution<CPU_NumericType>  T(-1, 0, 1);
  
    boost::variate_generator<boost::mt11213b&, boost::normal_distribution<CPU_NumericType> >     get_N(mt, N);
    boost::variate_generator<boost::mt11213b&, boost::bernoulli_distribution<CPU_NumericType> >  get_B(mt, B);
    boost::variate_generator<boost::mt11213b&, boost::triangle_distribution<CPU_NumericType> >   get_T(mt, T);
  
    std::vector<CPU_NumericType> eigenvalues;
    vcl_size_t matrix_size = matrix.size1();
    VectorT r(matrix_size);
    std::vector<CPU_NumericType> s(matrix_size);
  
    for (vcl_size_t i=0; i<s.size(); ++i)
      s[i] = 3.0 * get_B() + get_T() - 1.5;
  
    detail::copy_vec_to_vec(s,r);
  
    vcl_size_t size_krylov = (matrix_size < tag.krylov_size()) ? matrix_size
                                                                : tag.krylov_size();
  
    switch (tag.method())
    {
      case lanczos_tag::partial_reorthogonalization:
        eigenvalues = detail::lanczosPRO(matrix, r, size_krylov, tag);
        break;
      case lanczos_tag::full_reorthogonalization:
        eigenvalues = detail::lanczosFRO(matrix, r, size_krylov, tag);
        break;
      case lanczos_tag::no_reorthogonalization:
        eigenvalues = detail::lanczos(matrix, r, size_krylov, tag);
        break;
    }
  
    std::vector<CPU_NumericType> largest_eigenvalues;
  
    for (vcl_size_t i = 1; i<=tag.num_eigenvalues(); i++)
      largest_eigenvalues.push_back(eigenvalues[size_krylov-i]);
  
  
    return largest_eigenvalues;
  }
  
  
  
  
  } // end namespace linalg
  } // end namespace viennacl
  #endif
  