ATOOLS/Math/Matrix.H

0001 #ifndef ATOOLS_Math_Matrix_H
0002 #define ATOOLS_Math_Matrix_H
0003
0004
0005 // (4x4) Matrix class...
0006
0007 #include "ATOOLS/Math/Vector.H"
0008 #include "ATOOLS/Math/MyComplex.H"
0009
0010 namespace ATOOLS {
0011
0012   template<int _rank>
0013   class Matrix {
0014   protected:
0015     double ** p_m;
0016     void NumRecipesNotation();
0017     void AmegicNotation();
0018     void Jacobi(double d[], Matrix<_rank>&, int *) const;
0019   public:
0020     Matrix();
0021     Matrix(const double ma[_rank][_rank]);
0022     Matrix(const Matrix<_rank>&);
0023     ~Matrix();
0024     Matrix<_rank>& operator=(const Matrix<_rank>&);
0025     Matrix<_rank>  operator*(const double);
0026     Matrix<_rank>  operator*(const Matrix<_rank>&);
0027     double       * operator[](int i);
0028     const double * operator[](int i) const;
0029     void MatrixOut() const;
0030     int Rank() const {return _rank;}
0031     void Diagonalize(double *,Matrix<_rank>&) const;
0032     void DiagonalizeSort(double *,Matrix<_rank>&) const;
0033     Matrix<_rank> Dagger();
0034     //Vec4D operator*(const Vec4D&);
0035   };
0036
0037   class CMatrix {
0038   protected:
0039     Complex ** p_m;
0040     int        m_rank;
0041   public:
0042     CMatrix() { m_rank = 0; p_m = NULL; }
0043     CMatrix(int);
0044     CMatrix(Complex**,int);
0045     CMatrix(const CMatrix&);
0046     ~CMatrix();
0047     CMatrix         Conjugate();
0048     int             Rank() const            { return m_rank; }
0049     Complex       * operator[](int i)       { return p_m[i]; }
0050     const Complex * operator[](int i) const { return p_m[i]; }
0051     Vec4C operator* (const Vec4C& cvec);
0052     Complex ** GetMatrix()      const { return p_m; }
0053   };
0054
0055
0056   //Global functions.
0057
0058   template<int _rank>
0059   inline Matrix<_rank> operator*(const double scal, const Matrix<_rank>& in) {
0060     Matrix<_rank> out;
0061     for(short int i=0; i<_rank; i++) {
0062       for(short int j=0; j<_rank; j++) {
0063     out[i][j]=scal*in[i][j];
0064       }
0065     }
0066     return out;
0067   }
0068
0069   template<int _rank>
0070   Matrix<_rank> operator*(const Matrix<_rank>& a,const Matrix<_rank>& b) {
0071     Matrix<_rank> out;
0072
0073     for(short int i=0; i<_rank; i++) {
0074       for(short int j=0; j<_rank; j++) {
0075     out[i][j] = 0.;
0076     for(short int k=0; k<_rank; k++) out[i][j] += a[i][k]*b[k][j];
0077       }
0078     }
0079     return out;
0080   }
0081
0082   inline Vec4D operator*(const Matrix<4>& a,const Vec4D& b) {
0083     Vec4D out;
0084
0085     for(short int i=0; i<4; i++) {
0086       out[i] = 0.;
0087       for(short int k=0; k<4; k++) out[i] += a[i][k]*b[k];
0088     }
0089     return out;
0090   }
0091
0092   inline Vec4D operator*(const Vec4D& a,const Matrix<4>& b) {
0093     Vec4D out;
0094
0095     for(short int i=0; i<4; i++) {
0096       out[i] = 0.;
0097       for(short int k=0; k<4; k++) out[i] += a[k]*b[k][i];
0098     }
0099     return out;
0100   }
0101
0102   CMatrix operator*(const Complex scal, const CMatrix& in);
0103   CMatrix operator*(const CMatrix& a,const CMatrix& b);
0104
0105   /*!
0106     \file
0107     \brief contains class Matrix
0108   */
0109
0110   /*!
0111     \class Matrix
0112     \brief is a template for a \f$n\times n\f$ matrix and some operations on it
0113
0114     This matrix template can be used to represent any (real) \f$n\times n\f$
0115     (i.e. quadratic) matrix.
0116
0117     The following operations are provided:
0118      - Multiplication with a scalar
0119      - Multiplication with another matrix of same rank
0120      - Transpose of a matrix
0121      - Diagonalization of a matrix
0122      - Diagonalization with elements sorted
0123      .
0124      For \f$4\times 4\f$ matrices there is the
0125      - Mulitplication with a 4-vector (Vec4D)
0126      .
0127      also implemented.
0128   */
0129
0130   /*!
0131     \fn void Matrix::NumRecipesNotation();
0132     \brief changes access mode to matrix
0133   */
0134
0135   /*!
0136     \fn void Matrix::AmegicNotation();
0137     \brief changes access mode to matrix
0138   */
0139
0140   /*!
0141     \fn void Matrix::Jacobi(double d[], Matrix<_rank>&, int *);
0142     \brief uses Jacobi method to obtain eigenvalues and eigenvectors for a matrix
0143   */
0144
0145   /*!
0146     \fn Matrix::Matrix();
0147     \brief Default constructor, initializes a matrix filed with zeros
0148   */
0149
0150   /*!
0151     \fn Matrix::Matrix(const double ma[_rank][_rank]);
0152     \brief Special constructor, initializes a matrix with a given C array
0153   */
0154
0155   /*!
0156     \fn Matrix::Matrix(const Matrix<_rank>&);
0157     \brief Copy constructor
0158   */
0159
0160   /*!
0161     \fn Matrix::~Matrix();
0162     \brief Destructor
0163   */
0164
0165   /*!
0166     \fn Matrix<_rank>& Matrix::operator=(const Matrix<_rank>&);
0167     \brief Assignment operator
0168   */
0169
0170   /*!
0171     \fn Matrix<_rank> Matrix::operator*(const double);
0172     \brief implements the multiplication with a scalar
0173   */
0174
0175   /*!
0176     \fn Matrix<_rank> Matrix::operator*(const Matrix<_rank>&);
0177     \brief implements the multiplication with a matrix of same rank
0178   */
0179
0180   /*!
0181     \fn double* Matrix::operator[](int i) {return p_m[i];}
0182     \brief provides write access matrix via p_m[i][j]
0183   */
0184
0185   /*!
0186     \fn const double* Matrix::operator[](int i) const {return p_m[i];}
0187     \brief provides read only acces to a constant matrix via p_m[i][j]
0188   */
0189
0190   /*!
0191     \fn void Matrix::MatrixOut() const;
0192     \brief prints a matrix on screen
0193   */
0194
0195   /*!
0196     \fn const int Matrix::Rank() const {return _rank;}
0197     \brief returns the rank of the matrix
0198   */
0199
0200   /*!
0201     \fn void Matrix::Diagonalize(double*,Matrix<_rank>&);
0202     \brief diagonizes the matrix
0203   */
0204
0205   /*!
0206     \fn void Matrix::DiagonalizeSort(double*,Matrix<_rank>&);
0207     \brief diagonizes the matrix (eigenvalues are sorted)
0208   */
0209
0210   /*!
0211     \fn Matrix<_rank> Matrix::Dagger();
0212     \brief transposes the matrix
0213   */
0214
0215   // other methods
0216
0217   /*!
0218     \fn inline Vec4D operator*(const Matrix<4>& a,const Vec4D& b)
0219     \brief multiplication of a \f$4\times 4\f$ matrix with a 4-vector
0220   */
0221
0222   /*!
0223     \fn inline Vec4D operator*(const Vec4D& a,const Matrix<4>& b)
0224     \brief multiplication of a 4-vector with a \f$4\times 4\f$ matrix
0225   */
0226 }
0227
0228 #endif
0229
0230