// This may look like C code, but it is really -*- C++ -*- #ifndef SOpeMatrix_SEEN #define SOpeMatrix_SEEN #include "machdefs.h" #include "tmatrix.h" #include "tvector.h" //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// //------------------------------------------------------------// // La classe de calcul simple sur les TMatrix // //------------------------------------------------------------// //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// namespace SOPHYA { /*! \class SimpleMatrixOperation \ingroup TArray \brief Class for simple operation on TMatrix This class provides presently a basic implementation of the Gauss pivoting for linear system solving and matrix inversion. \warning For large or ill conditioned matrices, use the methods from LinAlg module. \sa SOPHYA::LapackServer SOPHYA::TMatrix */ //! Class for simple operation on TMatrix template class SimpleMatrixOperation { public: //! To define the type of data scaling for Gauss pivoting method (type T) enum GausPivScal { PIV_NO_SCALE = 0, //!< do not scale datas for gauss pivoting PIV_GLOB_SCALE = 1, //!< do a global scaling of datas for gauss pivoting (default) PIV_ROW_SCALE = 2 //!< do a scaling by row of datas for gauss pivoting }; static inline int SetGausPivScal(int gps = PIV_GLOB_SCALE) { GausPivScaling = gps; return GausPivScaling; } static TMatrix Inverse(TMatrix const & A); static T GausPiv(TMatrix& A, TMatrix& B, bool computedet=false); protected: static int GausPivScaling; }; } // Fin du namespace //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// //------------------------------------------------------------// // Resolution de systemes lineaires // //------------------------------------------------------------// //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// namespace SOPHYA { //------------------------------------------------------------ // Resolution du systeme A*C = B //------------------------------------------------------------ /*! \ingroup TArray \fn LinSolveInPlace(TMatrix&,TVector&) \brief Solve A*C = B for C in place and return 1 if OK, 0 if not. */ template inline T LinSolveInPlace(TMatrix& a, TVector& b) { if(a.NCols() != b.NRows() || a.NCols() != a.NRows()) throw(SzMismatchError("LinSolveInPlace(TMatrix,TVector) size mismatch")); return SimpleMatrixOperation::GausPiv(a,b); } //------------------------------------------------------------ // Resolution du systeme A*C = B, avec C retourne dans B //------------------------------------------------------------ /*! \ingroup TArray \fn LinSolve(const TMatrix&, const TVector&, TVector&) \brief Solve A*C = B and return C. Return 1 if OK, 0 if not. */ template inline T LinSolve(const TMatrix& a, const TVector& b, TVector& c) { if(a.NCols()!=b.NRows() || a.NCols()!=a.NRows()) throw(SzMismatchError("LinSolve(TMatrix,TVector) size mismatch")); c = b; TMatrix a1(a); return SimpleMatrixOperation::GausPiv(a1,c); } } // Fin du namespace //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// //------------------------------------------------------------// // Inverse d'une matrice // //------------------------------------------------------------// //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// namespace SOPHYA { /*! \ingroup TArray \fn Inverse(TMatrix const &) \brief To inverse a TMatrix. */ template inline TMatrix Inverse(TMatrix const & A) {return SimpleMatrixOperation::Inverse(A);} /*! \ingroup TArray \fn Determinant(TMatrix const &) \brief Return the TMatrix determinant */ template inline T Determinant(TMatrix const & A) { TMatrix a(A,false); TMatrix b(a.NCols(),a.NRows()); b = IdentityMatrix(1.); return SimpleMatrixOperation::GausPiv(a,b,true); } /*! \ingroup TArray \fn GausPiv(TMatrix const &,TMatrix &, bool) \brief Gauss pivoting on matrix \b A, doing the same operations on matrix \b B. \param computedet = true : return the determinant of \b A (beware of over/underfloat), default is false. \return determinant if requested, or 1 if OK. \sa SOPHYA::SimpleMatrixOperation::GausPiv(TMatrix&,TMatrix&,bool) */ template inline T GausPiv(TMatrix const & A,TMatrix & B, bool computedet=false) { TMatrix a(A,false); return SimpleMatrixOperation::GausPiv(a,B,computedet); } } // Fin du namespace //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// //------------------------------------------------------------// // Linear fitting // //------------------------------------------------------------// //////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// namespace SOPHYA { /*! \class LinFitter \ingroup TArray Class for linear fitting \sa TMatrix TArray */ //! Class for linear fitting template class LinFitter { public : LinFitter(); virtual ~LinFitter(); //! Linear fitting r_8 LinFit(const TVector& x, const TVector& y, sa_size_t nf, T (*f)(sa_size_t,T), TVector& c); //! Linear fitting r_8 LinFit(const TMatrix& fx, const TVector& y, TVector& c); //! Linear fitting with errors r_8 LinFit(const TVector& x, const TVector& y, const TVector& errY2, sa_size_t nf,T (*f)(sa_size_t,T), TVector& c, TVector& errC); //! Linear fitting with errors r_8 LinFit(const TMatrix& fx, const TVector& y, const TVector& errY2, TVector& c, TVector& errC); }; } // Fin du namespace #endif