STK_lapack_SymEigen.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff
 
    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA
 
    Contact : S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
*/
 
/*
 * Project:  stkpp::Algebra
 * created on: 20 nov. 2013
 * Author:   iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef STK_LAPACK_SYMEIGEN_H
#define STK_LAPACK_SYMEIGEN_H
 
#include "STK_ISymEigen.h"
#include "STK_lapack_Util.h"
 
namespace STK
{
 
// forward declaration
namespace lapack
{
template<class SquareArray> class SymEigen;
}
 
namespace hidden
{
template<class SquareArray_>
struct AlgebraTraits< lapack::SymEigen<SquareArray_> >
{
  typedef SquareArray_ SquareArray;
};
 
} // namespace hidden
 
 
namespace lapack
{
template<class SquareArray>
class SymEigen: public ISymEigen< SymEigen<SquareArray> >
{
  public:
    typedef ISymEigen< SymEigen<SquareArray> > Base;
    using Base::eigenValues_;
    using Base::eigenVectors_;
    using Base::range_;
    using Base::norm_;
    using Base::rank_;
    using Base::det_;
 
    SymEigen();
    SymEigen( SquareArray const& data, bool ref =false);
    template<class Derived>
    SymEigen( ArrayBase<Derived> const& data);
    SymEigen( SymEigen const& eigen);
    virtual ~SymEigen() {}
 
    // setters
    inline void setJobz(char jobz) { jobz_ = jobz;}
    inline void setRange(char range) { RANGE_ = range;}
    inline void setUplo(char uplo) { UPLO_ = uplo;}
    inline void setVlAndVu(Real const& vl, Real const& vu) { VL_ = vl; VU_ = vu;}
    inline void setIlAndIu(int il, int iu) { IL_ = il; IU_ = iu;}
 
    inline virtual SymEigen* clone() const { return new SymEigen(*this);}
    bool runImpl();
 
  private:
    char jobz_, RANGE_, UPLO_;
    Real VL_, VU_;
    int IL_, IU_;
};
/* @brief Default Constructor */
template<class SquareArray>
SymEigen<SquareArray>::SymEigen(): Base()
                                 , jobz_('V'), RANGE_('A'), UPLO_('U')
                                 , VL_(0.0), VU_(0.0), IL_(0), IU_(0)
{}
/* @brief Constructor
 *  @param data reference on a symmetric square matrix
 *  @param ref @c true if we overwrite the data set, @c false otherwise
 */
template<class SquareArray>
SymEigen<SquareArray>::SymEigen( SquareArray const& data, bool ref)
                              : Base(data)
                               , jobz_('V'), RANGE_('A'), UPLO_('U')
                               , VL_(0.0), VU_(0.0), IL_(0), IU_(0)
{}
/* @brief Constructor
 *  @param data reference on a symmetric square matrix
 *  @param ref @c true if we overwrite the data set, @c false otherwise
 */
template<>
inline SymEigen<CSquareX>::SymEigen( CSquareX const& data, bool ref)
                                  : Base(data, ref)
                                   , jobz_('V'), RANGE_('A'), UPLO_('U')
                                   , VL_(0.0), VU_(0.0), IL_(0), IU_(0)
{}
 
template<class SquareArray>
template<class Derived>
SymEigen<SquareArray>::SymEigen( ArrayBase<Derived> const& data)
                              : Base(data)
                               , jobz_('V'), RANGE_('A'), UPLO_('U')
                               , VL_(0.0), VU_(0.0), IL_(0), IU_(0)
{}
template<class SquareArray>
SymEigen<SquareArray>::SymEigen( SymEigen const& eigen)
                              : Base(eigen)
                               , jobz_(eigen.jobz_), RANGE_(eigen.RANGE_), UPLO_(eigen.UPLO_)
                               , VL_(eigen.VL_), VU_(eigen.VU_), IL_(eigen.IL_), IU_(eigen.IU_)
 {}
 
/* @brief Run eigen decomposition
 *  Launch SYEVR LAPACK routine to perform the eigenvalues decomposition.
 *  @return @c true if no error occur, @c false otherwise
 */
template<class SquareArray>
bool SymEigen<SquareArray>::runImpl()
{
#ifdef STK_ALGEBRA_VERY_VERBOSE
  stk_cout << _T("Enter in SymEigen::run\n");
#endif
  /* copy square matrix with the original data set. */
  CSquareX data_ = eigenVectors_;
  // shift data sets
  data_.shift(0);
  eigenVectors_.shift(0);
  eigenValues_.shift(0);
  this->SupportEigenVectors_.shift(0);
  /* set default behavior */
  Real absTol = 0.0; // let Lapack chose the correct tolerance
  // get optimal size necessary for work
  Real work; // work is just one place, get the optimal size for work
  // iwork is just on place, get the optimal size for iWork
  int iwork, lWork =-1, liwork =-1; // workspace variable
 
  int info = 1;
#ifdef STK_ALGEBRA_DEBUG
  stk_cout << _T("Data dimensions: ") << data_.rows() << " " << data_.cols() << "\n";
  stk_cout << _T("eigenValues_ dimensions: ") << eigenValues_.rows() << " " << eigenValues_.cols() << "\n";
  stk_cout << _T("eigenVectors_ dimensions: ") << eigenVectors_.rows() << " " << eigenVectors_.cols() << "\n";
  stk_cout << _T("Options: ") << jobz_ << " " << RANGE_ << " " << UPLO_ << "\n";
#endif
  info = syevr( jobz_, RANGE_, UPLO_
              , range_.size(), data_.p_data(), range_.size()
              , VL_, VU_, IL_, IU_
              , absTol, &rank_,  eigenValues_.p_data()
              , eigenVectors_.p_data(), range_.size(), this->SupportEigenVectors_.p_data()
              , &work, lWork, &iwork, liwork);
  // check any error
  if (info!=0)
  {
    if (info>0)
    { this->msg_error_ = STKERROR_NO_ARG(SymEigen::run,internal error);
      return false;
    }
    this->msg_error_= STKERROR_1ARG(SymEigen::run,-info,error parameter);
    return false;
  }
#ifdef STK_ALGEBRA_DEBUG
  stk_cout << _T("Size needed:") << (int)work << " " << iwork << " " << "\n";
#endif
  // get results and allocate space
  lWork = (int)work;
  liwork = iwork;
  Real* p_work = new Real[lWork];
  int* p_iwork = new int[liwork];
 
  // Call SYEVR with the optimal block size
  info = syevr( jobz_, RANGE_, UPLO_
              , range_.size(), data_.p_data(), range_.size()
              , VL_, VU_, IL_, IU_
              , absTol, &rank_, eigenValues_.p_data()
              , eigenVectors_.p_data(), range_.size(), this->SupportEigenVectors_.p_data()
              , p_work, lWork, p_iwork, liwork);
  // recover memory
  delete[] p_work;
  delete[] p_iwork;
 
  // finalize
  data_.shift(range_.begin());
  eigenVectors_.shift(range_.begin());
  eigenValues_.shift(range_.begin());
  this->SupportEigenVectors_.shift(range_.begin());
  this->finalizeStep();
  this->hasRun_ = true;
  // return the result of the computation
  if (!info) return true;
  if (info>0)
  { this->msg_error_ = STKERROR_NO_ARG(SymEigen::run,internal error);
    return false;
  }
  this->msg_error_= STKERROR_1ARG(SymEigen::run,-info,error parameter);
  return false;
}
 
} // namespace lapack
 
} // namespace STK
 
#endif /* STK_LAPACK_SYMEIGEN_H */