synaps/linalg/Eigen.h

Go to the documentation of this file.

/*********************************************************************
*      This file is part of the source code of SYNAPS kernel.        *
*      Author(s): B. Mourrain, GALAAD, INRIA                         *
**********************************************************************
History: 
$Id: Eigen.h,v 1.1 2005/07/11 08:33:33 mourrain Exp $
**********************************************************************/
#ifndef SYNAPS_LINALG_EIGEN_H
#define SYNAPS_LINALG_EIGEN_H
//----------------------------------------------------------------------
#include <synaps/init.h>
#include <synaps/linalg/VectDse.h>
#include <synaps/linalg/MatrDse.h>
#include <synaps/linalg/MATRIX.m>

//--------------------------------------------------------------------

__BEGIN_NAMESPACE_SYNAPS


namespace linalg
{
template<class T> struct rep1d;
}

namespace lapack
{
template<class C> struct rep2d;
}


//======================================================================
/*{ \latexonly \subsubsection*{Classical eigenproblem.} \endlatexonly }*/
//======================================================================
//--------------------------------------------------------------------
template<class VR_, class MR_> void eignval(VR_ & v, const MR_ & m);
template<class VR_, class MR_> void eignval(VR_ & v, const MR_ & m, const Real & mth);

template <class M> inline 
VectDse<typename algebraic_closure<typename M::value_type>::T> Eigenval(M & A)
{
  typedef VectDse<typename algebraic_closure<typename M::value_type>::T> V;
  V x(A.nbcol());
  eignval(x.rep(),A.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M> inline 
VectDse<typename M::value_type> Eigenval(M & A, Real)
{
  typedef VectDse<typename realof<typename M::value_type>::T> V;
  V x;
  eignval(x.rep(),A.rep(), Real());
  return x;
}
//----------------------------------------------------------------------
template <class M> inline 
VectDse<typename algebraic_closure<typename M::value_type>::T>
Eigenval(M & A, VectDse<typename M::value_type> & Er)
{
  typedef VectDse<typename algebraic_closure<typename M::value_type>::T> V;
  V x(A.nbcol());
  eignval(x.rep(),A.rep(), Er.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M>  inline 
MatrDse<typename algebraic_closure<typename M::value_type>::T,
  typename lapack::rep2d<typename algebraic_closure<typename M::value_type>::T> > 
  Eigenvct(const M & A)
{
  typedef MatrDse<typename algebraic_closure<typename M::value_type>::T, 
   typename lapack::rep2d<typename algebraic_closure<typename M::value_type>::T> > E;
  E x(A.nbcol(),A.nbrow());
  reigvct(x.rep(),A.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M> inline 
M Eigenvct(M & A, Real)
{
  M x(A.nbcol(),A.nbrow());
  reigvct(x.rep(), A.rep(), Real());
  return x;
}

//======================================================================
/*{ \latexonly \subsubsection*{Generalized eigenproblem.} \endlatexonly }*/
//======================================================================

template <class M>  inline 
VectDse<typename algebraic_closure<typename M::value_type>::T>
  Eigenval(const M & A, const M & B)
{
  typedef VectDse<typename algebraic_closure<typename M::value_type>::T> V;
  V x(A.nbcol());
  eignval(x.rep(),A.rep(),B.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M>  inline 
MatrDse<typename algebraic_closure<typename M::value_type>::T,
 typename lapack::rep2d<typename algebraic_closure<typename M::value_type>::T> >
  Eigenvct(const M & A, const M & B)
{
  typedef MatrDse<typename algebraic_closure<typename M::value_type>::T,
    lapack::rep2d<typename algebraic_closure<typename M::value_type>::T> > E;
  E x(A.nbcol(),A.nbrow());
  reigvct(x.rep(),A.rep(),B.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M>  inline 
VectDse<typename M::value_type> 
  Eigenval(const M & A, const M & B, Real reel)
{
  typedef VectDse<typename M::value_type> V;
  V x(A.nbcol());
  eignval(x.rep(),A.rep(),B.rep(), reel);
  return x;
}
//----------------------------------------------------------------------
template <class M> inline
M  Eigenvct(M & A, M & B, Real r)
{
  M x(A.nbcol(),A.nbrow());
  reigvct(x.rep(),A.rep(),B.rep(),r);
  return x;
}
//----------------------------------------------------------------------
template <typename C,class M,class V> inline
MatrDse<typename M::value_type, lapack::rep2d<typename M::value_type> > 
  Eigenvct(MatrDse<C,M> & A, 
           MatrDse<C,M> & B, Real r, 
           VectDse<V> & v)
{
  typedef MatrDse<C,lapack::rep2d<typename M::value_type> > E;
  E x(A.nbcol(),A.nbrow());
  reigvct(x.rep(),A.rep(),B.rep(),r,v.rep());
  return x;
}
//----------------------------------------------------------------------
template <class M> inline 
MatrDse<typename M::value_type, lapack::rep2d<typename M::value_type> >
  Eigenvct(MatrDse<M> & A, 
           MatrDse<M> & B, Real reel, 
           const typename M::size_type & m)
{
  typedef MatrDse<typename M::value_type, 
    lapack::rep2d<typename M::value_type> > E;
  E x(A.nbcol(),A.nbrow());
  reigvct(x.rep(),A.rep(),B.rep(),reel,m);
  return x;
}


__END_NAMESPACE_SYNAPS

#endif // SYNAPS_LINALG_EIGEN_H