source: Sophya/trunk/SophyaExt/LinAlg/intflapack.h@ 2572

#ifndef  IntfLapack_H_SEEN
#define  IntfLapack_H_SEEN

#include "machdefs.h"
#include "tarray.h"
#include "tvector.h"

namespace SOPHYA {

template <class T>
class LapackServer { 
public:
  LapackServer();
  virtual ~LapackServer();

  virtual int LinSolve(TArray<T>& a, TArray<T> & b); 
  virtual int LinSolveSym(TArray<T>& a, TArray<T> & b); 
  virtual int LeastSquareSolve(TArray<T>& a, TArray<T> & b); 
  virtual int LeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond=-1.); 

  virtual int SVD(TArray<T>& a, TArray<T> & s); 
  virtual int SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt);
  virtual int SVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt);
 
  virtual int LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector=true);
  virtual int LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector);

  //! Set the workspace size factor for LAPACK routines
  inline void SetWorkSpaceSizeFactor(int f = 2)
  { wspace_size_factor = (f > 1) ? f : 1; }

  //! Returns the workspace size factor
  inline int  GetWorkSpaceSizeFactor() 
  { return wspace_size_factor; }

private:
  int SVDDriver(TArray<T>& a, TArray<T> & s, 
                TArray<T>* up=NULL, TArray<T> * vtp=NULL);
  int_4 ilaenv_en_C(int_4 ispec,char *name,char *opts,int_4 n1,int_4 n2,int_4 n3,int_4 n4);
  int_4 type2i4(void *val,int nbytes);

  int wspace_size_factor;
};

/*! \ingroup LinAlg
    \fn LapackLinSolve(TArray<T>&, TArray<T> &)
    \brief Solves the linear system A*X = B using LapackServer. 
*/
template <class T>
inline int LapackLinSolve(TArray<T>& a, TArray<T> & b)
{ LapackServer<T> lps; return( lps.LinSolve(a, b) );  }

/*! \ingroup LinAlg
    \fn LapackLinSolveSym(TArray<T>&, TArray<T> &)
    \brief Solves the linear system A*X = B with A symetric using LapackServer. 
*/
template <class T>
inline int LapackLinSolveSym(TArray<T>& a, TArray<T> & b)
{ LapackServer<T> lps; return( lps.LinSolveSym(a, b) );  }

/*! \ingroup LinAlg
    \fn LapackLeastSquareSolve(TArray<T>&, TArray<T> &)
    \brief Solves the linear least squares problem A*X - B
*/
template <class T>
inline int LapackLeastSquareSolve(TArray<T>& a, TArray<T> & b)
{ LapackServer<T> lps; return( lps.LeastSquareSolve(a, b) );  }

/*! \ingroup LinAlg
    \fn LapackLeastSquareSolveSVD_DC
    \brief Solves the linear least squares problem A*X = B by SVD
*/
template <class T>
inline int LapackLeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond=-1.)
{ LapackServer<T> lps; return( lps.LeastSquareSolveSVD_DC(a,b,s,rank,rcond) );  }

/*! \ingroup LinAlg
    \fn LapackSVD(TArray<T>&, TArray<T> &)
    \brief SVD decomposition using LapackServer. 
*/
template <class T>
inline int LapackSVD(TArray<T>& a, TArray<T> & s)
{ LapackServer<T> lps; return( lps.SVD(a, s) ); }


/*! \ingroup LinAlg
    \fn LapackSVD(TArray<T>&, TArray<T> &, TArray<T> &, TArray<T> &)
    \brief SVD decomposition using LapackServer. 
*/
template <class T>
inline int LapackSVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt)
{ LapackServer<T> lps; return( lps.SVD(a, s, u, vt) ); }


/*! \ingroup LinAlg
    \fn LapackSVD_DC(TMatrix<T>&, TVector<r_8>&, TMatrix<T>&, TMatrix<T>&)
    \brief SVD decomposition DC using LapackServer. 
*/
template <class T>
inline int LapackSVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt)
{ LapackServer<T> lps; return( lps.SVD_DC(a, s, u, vt) ); }


/*! \ingroup LinAlg
    \fn LapackEigenSym(TArray<T>&, TArray<T> &)
    \brief Compute the eigenvalues and eigenvectors of A (symetric or hermitian). 
*/
template <class T>
inline int LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector=true)
{ LapackServer<T> lps; return( lps.LapackEigenSym(a,b,eigenvector) );  }

/*! \ingroup LinAlg
    \fn LapackEigen(TArray<T>&, TArray<T> &)
    \brief Compute the eigenvalues and (right) eigenvectors of A (general square matrix). 
*/
template <class T>
inline int LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector=true)
{ LapackServer<T> lps; return( lps.LapackEigen(a,eval,evec,eigenvector) );  }

} // Fin du namespace

#endif