CLHEP/Matrix/DiagMatrix.h

0001 // -*- C++ -*-
0002 // CLASSDOC OFF
0003 // ---------------------------------------------------------------------------
0004 // CLASSDOC ON
0005 //
0006 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
0007 //
0008 // This software written by Nobu Katayama and Mike Smyth, Cornell University.
0009 //
0010 // DiagMatrix is a class for diagonal matrix. This is useful for a covariance
0011 // matrix of measured quantities since they are uncorrelated to each other
0012 // and therefore diagonal. It is obviously smaller and faster to manipulate.
0013 //
0014
0015 #ifndef _DIAGMatrix_H_
0016 #define _DIAGMatrix_H_
0017
0018 #include <vector>
0019
0020 #include "CLHEP/Matrix/defs.h"
0021 #include "CLHEP/Matrix/GenMatrix.h"
0022
0023 namespace CLHEP {
0024
0025 class HepRandom;
0026
0027 class HepMatrix;
0028 class HepSymMatrix;
0029 class HepVector;
0030
0031 /**
0032  * @author
0033  * @ingroup matrix
0034  */
0035 class HepDiagMatrix: public HepGenMatrix {
0036 public:
0037    inline HepDiagMatrix();
0038    // Default constructor. Gives 0x0 matrix. Another Matrix can be assigned
0039    // to it.
0040
0041    explicit HepDiagMatrix(int p);
0042    HepDiagMatrix(int p, int);
0043    // Constructor. Gives p x p diagonal matrix.
0044    // With a second argument, either 0 or 1, the matrix is initialized.
0045    // 0 means a zero matrix, 1 means the identity matrix.
0046
0047    HepDiagMatrix(int p, HepRandom &r);
0048
0049    HepDiagMatrix(const HepDiagMatrix &hm1);
0050    // Copy constructor.
0051
0052    virtual ~HepDiagMatrix();
0053    // Destructor.
0054
0055    inline int num_row() const;
0056    inline int num_col() const;
0057    // Returns the number of rows/columns. (Should be equal.)
0058
0059    double &operator()(int row, int col);
0060    const double &operator()(int row, int col) const;
0061    // Read or write a matrix element. row must be equal to col.
0062    // ** Note that indexing starts from (1,1). **
0063
0064    double &fast(int row, int col);
0065    const double &fast(int row, int col) const;
0066    // fast element access.
0067    // Must be row>=col;
0068    // ** Note that indexing starts from (1,1). **
0069
0070    void assign(const HepMatrix &hm2);
0071    // Assigns hm2 to d, assuming hm2 is a diagnal matrix.
0072
0073    void assign(const HepSymMatrix &hm2);
0074    // Assigns hm2 to d, assuming hm2 is a diagnal matrix.
0075
0076    void assign(const HepDiagMatrix &hm2);
0077    // Another form of assignment. For consistency.
0078
0079    HepDiagMatrix & operator*=(double t);
0080    // Multiply a DiagMatrix by a floating number
0081
0082    HepDiagMatrix & operator/=(double t);
0083    // Divide a DiagMatrix by a floating number
0084
0085    HepDiagMatrix & operator+=( const HepDiagMatrix &hm2);
0086    HepDiagMatrix & operator-=( const HepDiagMatrix &hm2);
0087    // Add or subtract a DiagMatrix.
0088
0089    HepDiagMatrix & operator=( const HepDiagMatrix &hm2);
0090    // Assignment operator. To assign SymMatrix to DiagMatrix, use d<<s.
0091
0092    HepDiagMatrix operator- () const;
0093    // unary minus, ie. flip the sign of each element.
0094
0095    HepDiagMatrix T() const;
0096    // Returns the transpose of a DiagMatrix (which is itself).
0097
0098    HepDiagMatrix apply(double (*f)(double,
0099                            int, int)) const;
0100    // Apply a function to all elements of the matrix.
0101
0102    HepSymMatrix similarity(const HepMatrix &hm1) const;
0103    // Returns hm1*s*hm1.T().
0104    HepSymMatrix similarityT(const HepMatrix &hm1) const;
0105    // Returns hm1.T()*s*hm1.
0106
0107    double similarity(const HepVector &) const;
0108    // Returns v.T()*s*v (This is a scaler).
0109
0110    HepDiagMatrix sub(int min_row, int max_row) const;
0111    // Returns a sub matrix of a SymMatrix.
0112    HepDiagMatrix sub(int min_row, int max_row);
0113    // SGI CC bug. I have to have both with/without const. I should not need
0114    // one without const.
0115
0116    void sub(int row, const HepDiagMatrix &hm1);
0117    // Sub matrix of this SymMatrix is replaced with hm1.
0118
0119    HepDiagMatrix inverse(int&ierr) const;
0120    // Invert a Matrix. The matrix is not changed
0121    // Returns 0 when successful, otherwise non-zero.
0122
0123    void invert(int&ierr);
0124    // Invert a Matrix.
0125    // N.B. the contents of the matrix are replaced by the inverse.
0126    // Returns ierr = 0 when successful, otherwise non-zero.
0127    // This method has less overhead then inverse().
0128
0129    inline void invert();
0130    // Invert a matrix. Throw std::runtime_error on failure.
0131
0132    inline HepDiagMatrix inverse() const;
0133    // Invert a matrix. Throw std::runtime_error on failure.
0134
0135    double determinant() const;
0136    // calculate the determinant of the matrix.
0137
0138    double trace() const;
0139    // calculate the trace of the matrix (sum of diagonal elements).
0140
0141    class HepDiagMatrix_row {
0142    public:
0143       inline HepDiagMatrix_row(HepDiagMatrix&,int);
0144       inline double & operator[](int);
0145    private:
0146       HepDiagMatrix& _a;
0147       int _r;
0148    };
0149    class HepDiagMatrix_row_const {
0150    public:
0151       inline HepDiagMatrix_row_const(const HepDiagMatrix&,int);
0152       inline const double & operator[](int) const;
0153    private:
0154       const HepDiagMatrix& _a;
0155       int _r;
0156    };
0157    // helper classes to implement m[i][j]
0158
0159    inline HepDiagMatrix_row operator[] (int);
0160    inline HepDiagMatrix_row_const operator[] (int) const;
0161    // Read or write a matrix element.
0162    // While it may not look like it, you simply do m[i][j] to get an
0163    // element.
0164    // ** Note that the indexing starts from [0][0]. **
0165
0166 protected:
0167    inline int num_size() const;
0168
0169 private:
0170    friend class HepDiagMatrix_row;
0171    friend class HepDiagMatrix_row_const;
0172    friend class HepMatrix;
0173    friend class HepSymMatrix;
0174
0175    friend HepDiagMatrix operator*(const HepDiagMatrix &hm1,
0176                    const HepDiagMatrix &hm2);
0177    friend HepDiagMatrix operator+(const HepDiagMatrix &hm1,
0178                    const HepDiagMatrix &hm2);
0179    friend HepDiagMatrix operator-(const HepDiagMatrix &hm1,
0180                    const HepDiagMatrix &hm2);
0181    friend HepMatrix operator*(const HepDiagMatrix &hm1, const HepMatrix &hm2);
0182    friend HepMatrix operator*(const HepMatrix &hm1, const HepDiagMatrix &hm2);
0183    friend HepVector operator*(const HepDiagMatrix &hm1, const HepVector &hm2);
0184
0185 #ifdef DISABLE_ALLOC
0186    std::vector<double > m;
0187 #else
0188    std::vector<double,Alloc<double,25> > m;
0189 #endif
0190    int nrow;
0191 #if defined(__sun) || !defined(__GNUG__)
0192 //
0193 // Sun CC 4.0.1 has this bug.
0194 //
0195    static double zero;
0196 #else
0197    static const double zero;
0198 #endif
0199 };
0200
0201 std::ostream& operator<<(std::ostream &s, const HepDiagMatrix &q);
0202 // Write out Matrix, SymMatrix, DiagMatrix and Vector into ostream.
0203
0204 HepMatrix operator*(const HepMatrix &hm1, const HepDiagMatrix &hm2);
0205 HepMatrix operator*(const HepDiagMatrix &hm1, const HepMatrix &hm2);
0206 HepDiagMatrix operator*(double t, const HepDiagMatrix &d1);
0207 HepDiagMatrix operator*(const HepDiagMatrix &d1, double t);
0208 // Multiplication operators
0209 // Note that m *= hm1 is always faster than m = m * hm1
0210
0211 HepDiagMatrix operator/(const HepDiagMatrix &hm1, double t);
0212 // d = d1 / t. (d /= t is faster if you can use it.)
0213
0214 HepMatrix operator+(const HepMatrix &hm1, const HepDiagMatrix &d2);
0215 HepMatrix operator+(const HepDiagMatrix &d1, const HepMatrix &hm2);
0216 HepDiagMatrix operator+(const HepDiagMatrix &hm1, const HepDiagMatrix &d2);
0217 HepSymMatrix operator+(const HepSymMatrix &s1, const HepDiagMatrix &d2);
0218 HepSymMatrix operator+(const HepDiagMatrix &d1, const HepSymMatrix &s2);
0219 // Addition operators
0220
0221 HepMatrix operator-(const HepMatrix &hm1, const HepDiagMatrix &d2);
0222 HepMatrix operator-(const HepDiagMatrix &d1, const HepMatrix &hm2);
0223 HepDiagMatrix operator-(const HepDiagMatrix &d1, const HepDiagMatrix &d2);
0224 HepSymMatrix operator-(const HepSymMatrix &s1, const HepDiagMatrix &d2);
0225 HepSymMatrix operator-(const HepDiagMatrix &d1, const HepSymMatrix &s2);
0226 // Subtraction operators
0227
0228 HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2);
0229 // Direct sum of two diagonal matricies;
0230
0231 }  // namespace CLHEP
0232
0233 #ifdef ENABLE_BACKWARDS_COMPATIBILITY
0234 //  backwards compatibility will be enabled ONLY in CLHEP 1.9
0235 using namespace CLHEP;
0236 #endif
0237
0238 #ifndef HEP_DEBUG_INLINE
0239 #include "CLHEP/Matrix/DiagMatrix.icc"
0240 #endif
0241
0242 #endif /*!_DIAGMatrix_H*/