#ifndef IntfLapack_H_SEEN #define IntfLapack_H_SEEN #include "machdefs.h" #include "tarray.h" #include "tvector.h" namespace SOPHYA { template class LapackServer { public: LapackServer(bool throw_on_error=false); virtual ~LapackServer(); virtual int LinSolve(TArray& a, TArray & b); virtual int LinSolveSym(TArray& a, TArray & b); virtual int LeastSquareSolve(TArray& a, TArray & b); virtual int LeastSquareSolveSVD_DC(TMatrix& a,TMatrix& b,TVector& s,int_4& rank,r_8 rcond=-1.); // Calcul de la matrice inverse en utilisant la resolution de syst. lineaire virtual int ComputeInverse(TMatrix& a, TMatrix& ainv); virtual int SVD(TArray& a, TArray & s); virtual int SVD(TArray& a, TArray & s, TArray & u, TArray & vt); virtual int SVD_DC(TMatrix& a, TVector& s, TMatrix& u, TMatrix& vt); virtual int LapackEigenSym(TArray& a, TVector& b, bool eigenvector=true); virtual int LapackEigen(TArray& a, TVector< complex >& eval, TMatrix& 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& a, TArray & s, TArray* up=NULL, TArray * 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; bool Throw_On_Error; }; /*! \ingroup LinAlg \fn LapackLinSolve(TArray&, TArray &) \brief Solves the linear system A*X = B using LapackServer. */ template inline int LapackLinSolve(TArray& a, TArray & b) { LapackServer lps; return( lps.LinSolve(a, b) ); } /*! \ingroup LinAlg \fn LapackLinSolveSym(TArray&, TArray &) \brief Solves the linear system A*X = B with A symetric using LapackServer. */ template inline int LapackLinSolveSym(TArray& a, TArray & b) { LapackServer lps; return( lps.LinSolveSym(a, b) ); } /*! \ingroup LinAlg \fn LapackLeastSquareSolve(TArray&, TArray &) \brief Solves the linear least squares problem A*X - B */ template inline int LapackLeastSquareSolve(TArray& a, TArray & b) { LapackServer lps; return( lps.LeastSquareSolve(a, b) ); } /*! \ingroup LinAlg \fn LapackInverse(TMatrix&) \brief Computes the inverse matrix using linear system solver LapackServer::LinSolve. */ template inline TMatrix LapackInverse(TMatrix& a) { LapackServer lps; TMatrix ainv; lps.ComputeInverse(a, ainv); return ainv; } /*! \ingroup LinAlg \fn LapackLeastSquareSolveSVD_DC \brief Solves the linear least squares problem A*X = B by SVD */ template inline int LapackLeastSquareSolveSVD_DC(TMatrix& a,TMatrix& b,TVector& s,int_4& rank,r_8 rcond=-1.) { LapackServer lps; return( lps.LeastSquareSolveSVD_DC(a,b,s,rank,rcond) ); } /*! \ingroup LinAlg \fn LapackSVD(TArray&, TArray &) \brief SVD decomposition using LapackServer. */ template inline int LapackSVD(TArray& a, TArray & s) { LapackServer lps; return( lps.SVD(a, s) ); } /*! \ingroup LinAlg \fn LapackSVD(TArray&, TArray &, TArray &, TArray &) \brief SVD decomposition using LapackServer. */ template inline int LapackSVD(TArray& a, TArray & s, TArray & u, TArray & vt) { LapackServer lps; return( lps.SVD(a, s, u, vt) ); } /*! \ingroup LinAlg \fn LapackSVD_DC(TMatrix&, TVector&, TMatrix&, TMatrix&) \brief SVD decomposition DC using LapackServer. */ template inline int LapackSVD_DC(TMatrix& a, TVector& s, TMatrix& u, TMatrix& vt) { LapackServer lps; return( lps.SVD_DC(a, s, u, vt) ); } /*! \ingroup LinAlg \fn LapackEigenSym(TArray&, TArray &) \brief Compute the eigenvalues and eigenvectors of A (symetric or hermitian). */ template inline int LapackEigenSym(TArray& a, TVector& b, bool eigenvector=true) { LapackServer lps; return( lps.LapackEigenSym(a,b,eigenvector) ); } /*! \ingroup LinAlg \fn LapackEigen(TArray&, TArray &) \brief Compute the eigenvalues and (right) eigenvectors of A (general square matrix). */ template inline int LapackEigen(TArray& a, TVector< complex >& eval, TMatrix& evec, bool eigenvector=true) { LapackServer lps; return( lps.LapackEigen(a,eval,evec,eigenvector) ); } } // Fin du namespace #endif