STK_ISymEigen.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2016  Serge Iovleff, Université Lille 1, Inria
 
    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
 * Purpose:  Define The Interface SymEigen Class.
 * Author:   Serge Iovleff, S..._Dot_I..._At_stkpp_Dot_org (see copyright for ...)
 **/
 
#ifndef STK_ISYMEIGEN_H
#define STK_ISYMEIGEN_H
 
#include <Sdk.h>
#include "STK_Algebra_Util.h"
 
#include <Arrays/include/STK_CArraySquare.h>
#include <Arrays/include/STK_CArrayVector.h>
 
namespace STK
{
template<class Derived>
class ISymEigen: public IRunnerBase, public IRecursiveTemplate<Derived>
{
  protected:
    typedef IRunnerBase Base;
    typedef typename hidden::AlgebraTraits<Derived>::SquareArray SquareArray;
 
    typedef typename hidden::Traits< SquareArray >::Type Type;
    enum
    {
      structure_ = hidden::Traits< SquareArray >::structure_,
      orient_    = hidden::Traits< SquareArray >::orient_,
      sizeRows_  = hidden::Traits< SquareArray >::sizeRows_,
      sizeCols_  = hidden::Traits< SquareArray >::sizeCols_,
      size_ = (sizeRows_<sizeCols_) ? sizeRows_ : sizeCols_
 
    };
    ISymEigen();
    ISymEigen( SquareArray const& data, bool ref =false);
    template<class OtherDerived>
    ISymEigen( ExprBase<OtherDerived> const& data);
    ISymEigen( ISymEigen const& eigen);
 
  public:
    inline ~ISymEigen() {}
    ISymEigen& operator=( ISymEigen const& eigen)
    {
      range_ = eigen.range_;
      trace_ = eigen.trace_;    // trace of the matrix
      norm_  = eigen.norm_;     // norm of the matrix
      rank_  = eigen.rank_;     // rank of the matrix
      det_   = eigen.det_;      // determinant of the matrix
      eigenVectors_        = eigen.eigenVectors_;
      eigenValues_         = eigen.eigenValues_;
      SupportEigenVectors_ = eigen.SupportEigenVectors_;
      return *this;
    }
    // getters
    inline Range const& range()  const { return range_;}
    inline Type const& norm()  const { return norm_;}
    inline int const& rank()  const { return rank_;}
    inline Type const& det()  const { return det_;}
    inline Type const& trace()  const { return trace_;}
    inline CArraySquare<Type, size_> const& rotation() const{ return eigenVectors_;}
    inline CArraySquare<Type, size_> const& eigenVectors() const{ return eigenVectors_;}
    inline CArrayVector<Type, size_> const& eigenValues() const { return eigenValues_;}
 
    template<class OtherDerived>
    void setData( ExprBase<OtherDerived> const& data)
    {
      STK_STATIC_ASSERT(OtherDerived::structure_==(int)Arrays::square_,YOU_HAVE_TO_USE_A_SQUARE_MATRIX_IN_THIS_METHOD)
      range_ = data.range(); trace_ = 0; norm_ = 0.; rank_ = 0; det_ = 0.;
      eigenVectors_ = data;
      eigenValues_.resize(range_) = 0;
      SupportEigenVectors_.resize(2*data.size());
      this->hasRun_ = false;
    }
    // computation
    template<class ArraySquare>
    ArraySquare& ginv(ArraySquare& res)  const
    {
      Type tol = Arithmetic<Type>::epsilon() * norm_;
      if(tol==0) { tol = Arithmetic<Type>::min();}
      // compute and return PD^{-1}P'
      return (res = eigenVectors_ * eigenValues_.diagonalize().safeInverse(tol) * eigenVectors_.transpose());
    }
    template<class ArraySquare>
    ArraySquare& ginvsqrt(ArraySquare& res) const
    {
      Type tol = Arithmetic<Type>::epsilon() * norm_;
      if(tol==0) { tol = Arithmetic<Type>::min();}
      // compute and return PD^{-1/2}P'
      return(res = eigenVectors_ * eigenValues_.diagonalize().sqrt().safeInverse(tol) * eigenVectors_.transpose());
    }
    template<class ArraySquare>
    ArraySquare& gsqrt(ArraySquare& res) const
    {
      Type tol = Arithmetic<Type>::epsilon() * norm_;
      if(tol==0) { tol = Arithmetic<Type>::min();}
      // compute and return PD^{1/2}P'
      return(res = eigenVectors_ * eigenValues_.diagonalize().sqrt() * eigenVectors_.transpose());
    }
    virtual bool run()
    {
      if (eigenVectors_.empty()) { return true;}
      return this->asDerived().runImpl();
    }
 
  protected:
    Range range_;
    Type trace_;
    Type norm_;
    int rank_;
    Type det_;
    CArraySquare<Type, size_> eigenVectors_;
    CArrayVector<Type, size_> eigenValues_;
    CVectorXi SupportEigenVectors_;
    void finalizeStep();
};
 
/* @brief Default constructor */
template<class Derived>
ISymEigen<Derived>::ISymEigen()
                             : Base()
                             , range_(), trace_(Type(0)), norm_(Type(0)), rank_(Type(0)), det_(Type(0))
                             , eigenVectors_()
                             , eigenValues_()
                             , SupportEigenVectors_()
{
  STK_STATIC_ASSERT(   SquareArray::structure_==(int)Arrays::square_
                    || SquareArray::structure_==(int)Arrays::symmetric_
                    || SquareArray::structure_==(int)Arrays::upper_symmetric_
                    || SquareArray::structure_==(int)Arrays::lower_symmetric_
                   ,YOU_HAVE_TO_USE_A_SQUARE_MATRIX_IN_THIS_METHOD)
}
/* @brief Constructor
 *  The original data set can be overwritten by the eigenvectors if it is
 *  stored in a CSquareXd. Observe that in this case the base index have
 *  to be 0.
 *  @param data reference on a symmetric square matrix
 *  @param ref @c true if we overwrite the data set, @c false otherwise
 */
template<class Derived>
ISymEigen<Derived>::ISymEigen( SquareArray const& data, bool ref)
                             : Base()
                             , range_(data.range()), trace_(Type(0)), norm_(Type(0)), rank_(Type(0)), det_(Type(0))
                             , eigenVectors_(data, ref)
                             , eigenValues_(data.size(), 0.)
                             , SupportEigenVectors_(2*data.size(), 0)
{
    STK_STATIC_ASSERT(   SquareArray::structure_==(int)Arrays::square_
                      || SquareArray::structure_==(int)Arrays::symmetric_
                      || SquareArray::structure_==(int)Arrays::upper_symmetric_
                      || SquareArray::structure_==(int)Arrays::lower_symmetric_
                     ,YOU_HAVE_TO_USE_A_SQUARE_MATRIX_IN_THIS_METHOD)
}
/* @brief template constructor
 *  @param data reference on a symmetric square expression
 */
template<class Derived>
template<class OtherDerived>
ISymEigen<Derived>::ISymEigen( ExprBase<OtherDerived> const& data)
                             : Base()
                             , range_(data.range()), trace_(Type(0)), norm_(0.), rank_(0), det_(0.)
                             , eigenVectors_(data.asDerived())
                             , eigenValues_(data.size(), 0.)
                             , SupportEigenVectors_(2*data.size(), 0)
{
    STK_STATIC_ASSERT(   SquareArray::structure_==(int)Arrays::square_
                      || SquareArray::structure_==(int)Arrays::symmetric_
                      || SquareArray::structure_==(int)Arrays::upper_symmetric_
                      || SquareArray::structure_==(int)Arrays::lower_symmetric_
                     ,YOU_HAVE_TO_USE_A_SQUARE_MATRIX_IN_THIS_METHOD)
}
/* Copy constructor.
 *  @param eigen the EigenValue to copy
 **/
template<class Derived>
ISymEigen<Derived>::ISymEigen( ISymEigen const& eigen)
                             : Base(eigen)
                             , range_(eigen.range_), trace_(eigen.trace_), norm_(eigen.norm_), rank_(eigen.rank_), det_(eigen.det_)
                             , eigenVectors_(eigen.eigenVectors_)
                             , eigenValues_(eigen.eigenValues_)
                             , SupportEigenVectors_(eigen.SupportEigenVectors_)
{
    STK_STATIC_ASSERT(   SquareArray::structure_==(int)Arrays::square_
                      || SquareArray::structure_==(int)Arrays::symmetric_
                      || SquareArray::structure_==(int)Arrays::upper_symmetric_
                      || SquareArray::structure_==(int)Arrays::lower_symmetric_
                     ,YOU_HAVE_TO_USE_A_SQUARE_MATRIX_IN_THIS_METHOD)
}
 
/* finalize the computation by computing the rank, the trace norm and the
 * determinant of the matrix.
 **/
template<class Derived>
void ISymEigen<Derived>::finalizeStep()
{
  // compute sign and trace
  trace_ = 0.;
  int s = 1;
  for (int i=eigenValues_.begin(); i< eigenValues_.end(); ++i )
  {
    Type value = eigenValues_[i];
    s *= sign(value);
    trace_ += value;
  }
 
  // compute norm_ (2-norm) and determinant
  norm_ = eigenValues_.maxElt();
  det_  = 0;
  if (eigenValues_.abs().minElt() > 0)
  { det_ = s * std::exp(eigenValues_.abs().log().sum());}
 
  // compute norm_ rank_
  rank_ = eigenValues_.size(); // full rank
  Type tol = norm_ * Arithmetic<Type>::epsilon();
  for (int i=eigenValues_.begin(); i< eigenValues_.end(); ++i )
  { if (std::abs(eigenValues_[i]) < tol) { rank_--;}}
}
 
} // namespace STK
 
#endif //STK_ISYMEIGEN_H