#ifndef ATOOLS_Math_Matrix_H #define ATOOLS_Math_Matrix_H // (4x4) Matrix class... #include "ATOOLS/Math/Vector.H" #include "ATOOLS/Math/MyComplex.H" namespace ATOOLS { template class Matrix { protected: double ** p_m; void NumRecipesNotation(); void AmegicNotation(); void Jacobi(double d[], Matrix<_rank>&, int *) const; public: Matrix(); Matrix(const double ma[_rank][_rank]); Matrix(const Matrix<_rank>&); ~Matrix(); Matrix<_rank>& operator=(const Matrix<_rank>&); Matrix<_rank> operator*(const double); Matrix<_rank> operator*(const Matrix<_rank>&); double * operator[](int i); const double * operator[](int i) const; void MatrixOut() const; int Rank() const {return _rank;} void Diagonalize(double *,Matrix<_rank>&) const; void DiagonalizeSort(double *,Matrix<_rank>&) const; Matrix<_rank> Dagger(); //Vec4D operator*(const Vec4D&); }; class CMatrix { protected: Complex ** p_m; int m_rank; public: CMatrix() { m_rank = 0; p_m = NULL; } CMatrix(int); CMatrix(Complex**,int); CMatrix(const CMatrix&); ~CMatrix(); CMatrix Conjugate(); int Rank() const { return m_rank; } Complex * operator[](int i) { return p_m[i]; } const Complex * operator[](int i) const { return p_m[i]; } Vec4C operator* (const Vec4C& cvec); Complex ** GetMatrix() const { return p_m; } }; //Global functions. template inline Matrix<_rank> operator*(const double scal, const Matrix<_rank>& in) { Matrix<_rank> out; for(short int i=0; i<_rank; i++) { for(short int j=0; j<_rank; j++) { out[i][j]=scal*in[i][j]; } } return out; } template Matrix<_rank> operator*(const Matrix<_rank>& a,const Matrix<_rank>& b) { Matrix<_rank> out; for(short int i=0; i<_rank; i++) { for(short int j=0; j<_rank; j++) { out[i][j] = 0.; for(short int k=0; k<_rank; k++) out[i][j] += a[i][k]*b[k][j]; } } return out; } inline Vec4D operator*(const Matrix<4>& a,const Vec4D& b) { Vec4D out; for(short int i=0; i<4; i++) { out[i] = 0.; for(short int k=0; k<4; k++) out[i] += a[i][k]*b[k]; } return out; } inline Vec4D operator*(const Vec4D& a,const Matrix<4>& b) { Vec4D out; for(short int i=0; i<4; i++) { out[i] = 0.; for(short int k=0; k<4; k++) out[i] += a[k]*b[k][i]; } return out; } CMatrix operator*(const Complex scal, const CMatrix& in); CMatrix operator*(const CMatrix& a,const CMatrix& b); /*! \file \brief contains class Matrix */ /*! \class Matrix \brief is a template for a \f$n\times n\f$ matrix and some operations on it This matrix template can be used to represent any (real) \f$n\times n\f$ (i.e. quadratic) matrix. The following operations are provided: - Multiplication with a scalar - Multiplication with another matrix of same rank - Transpose of a matrix - Diagonalization of a matrix - Diagonalization with elements sorted . For \f$4\times 4\f$ matrices there is the - Mulitplication with a 4-vector (Vec4D) . also implemented. */ /*! \fn void Matrix::NumRecipesNotation(); \brief changes access mode to matrix */ /*! \fn void Matrix::AmegicNotation(); \brief changes access mode to matrix */ /*! \fn void Matrix::Jacobi(double d[], Matrix<_rank>&, int *); \brief uses Jacobi method to obtain eigenvalues and eigenvectors for a matrix */ /*! \fn Matrix::Matrix(); \brief Default constructor, initializes a matrix filed with zeros */ /*! \fn Matrix::Matrix(const double ma[_rank][_rank]); \brief Special constructor, initializes a matrix with a given C array */ /*! \fn Matrix::Matrix(const Matrix<_rank>&); \brief Copy constructor */ /*! \fn Matrix::~Matrix(); \brief Destructor */ /*! \fn Matrix<_rank>& Matrix::operator=(const Matrix<_rank>&); \brief Assignment operator */ /*! \fn Matrix<_rank> Matrix::operator*(const double); \brief implements the multiplication with a scalar */ /*! \fn Matrix<_rank> Matrix::operator*(const Matrix<_rank>&); \brief implements the multiplication with a matrix of same rank */ /*! \fn double* Matrix::operator[](int i) {return p_m[i];} \brief provides write access matrix via p_m[i][j] */ /*! \fn const double* Matrix::operator[](int i) const {return p_m[i];} \brief provides read only acces to a constant matrix via p_m[i][j] */ /*! \fn void Matrix::MatrixOut() const; \brief prints a matrix on screen */ /*! \fn const int Matrix::Rank() const {return _rank;} \brief returns the rank of the matrix */ /*! \fn void Matrix::Diagonalize(double*,Matrix<_rank>&); \brief diagonizes the matrix */ /*! \fn void Matrix::DiagonalizeSort(double*,Matrix<_rank>&); \brief diagonizes the matrix (eigenvalues are sorted) */ /*! \fn Matrix<_rank> Matrix::Dagger(); \brief transposes the matrix */ // other methods /*! \fn inline Vec4D operator*(const Matrix<4>& a,const Vec4D& b) \brief multiplication of a \f$4\times 4\f$ matrix with a 4-vector */ /*! \fn inline Vec4D operator*(const Vec4D& a,const Matrix<4>& b) \brief multiplication of a 4-vector with a \f$4\times 4\f$ matrix */ } #endif