TMatrixD


class description - source file - inheritance tree

class TMatrixD : public TObject


    protected:
void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb = 0, Int_t col_lwb = 0) void AMultB(const TMatrixD& a, const TMatrixD& b) void AtMultB(const TMatrixD& a, const TMatrixD& b) void EigenSort(TMatrixD& eigenVectors, TVectorD& eigenValues) void Invalidate() void Invert(const TMatrixD& m) void InvertPosDef(const TMatrixD& m) void MakeEigenVectors(TVectorD& d, TVectorD& e, TMatrixD& z) void MakeTridiagonal(TMatrixD& a, TVectorD& d, TVectorD& e) Int_t Pdcholesky(const Double_t* a, Double_t* u, const Int_t n) void Transpose(const TMatrixD& m) public:
TMatrixD TMatrixD() TMatrixD TMatrixD(Int_t nrows, Int_t ncols) TMatrixD TMatrixD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb) TMatrixD TMatrixD(Int_t nrows, Int_t ncols, const Double_t* elements, Option_t* option) TMatrixD TMatrixD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, const Double_t* elements, Option_t* option) TMatrixD TMatrixD(const TMatrixD& another) TMatrixD TMatrixD(TMatrixD::EMatrixCreatorsOp1 op, const TMatrixD& prototype) TMatrixD TMatrixD(const TMatrixD& a, TMatrixD::EMatrixCreatorsOp2 op, const TMatrixD& b) TMatrixD TMatrixD(const TLazyMatrixD& lazy_constructor) TMatrixD EigenVectors(TVectorD& eigenValues) virtual void ~TMatrixD() TMatrixD& Abs() TMatrixD& Apply(TElementActionD& action) TMatrixD& Apply(TElementPosActionD& action) static TClass* Class() Double_t ColNorm() const Double_t Determinant() const virtual void Draw(Option_t* option) Double_t E2Norm() const Int_t GetColLwb() const Int_t GetColUpb() const Double_t* GetElements() void GetElements(Double_t* elements, Option_t* option) const Int_t GetNcols() const Int_t GetNoElements() const Int_t GetNrows() const Int_t GetRowLwb() const Int_t GetRowUpb() const TMatrixD& HilbertMatrix() TMatrixD& Invert(Double_t* determ_ptr = 0) TMatrixD& InvertPosDef() virtual TClass* IsA() const Bool_t IsSymmetric() const Bool_t IsValid() const TMatrixD& MakeSymmetric() void Mult(const TMatrixD& a, const TMatrixD& b) Double_t Norm1() const TMatrixD& NormByColumn(const TVectorD& v, Option_t* option = "D") TMatrixD& NormByDiag(const TVectorD& v, Option_t* option = "D") TMatrixD& NormByRow(const TVectorD& v, Option_t* option = "D") Double_t NormInf() const Bool_t operator!=(Double_t val) const const Double_t& operator()(Int_t rown, Int_t coln) const Double_t& operator()(Int_t rown, Int_t coln) TMatrixD& operator*=(Double_t val) TMatrixD& operator*=(const TMatrixD& source) TMatrixD& operator*=(const TMatrixDDiag& diag) TMatrixD& operator*=(const TMatrixDRow& diag) TMatrixD& operator*=(const TMatrixDColumn& diag) TMatrixD& operator+=(Double_t val) TMatrixD& operator-=(Double_t val) TMatrixD& operator/=(const TMatrixDDiag& diag) TMatrixD& operator/=(const TMatrixDRow& diag) TMatrixD& operator/=(const TMatrixDColumn& diag) Bool_t operator<(Double_t val) const Bool_t operator<=(Double_t val) const TMatrixD& operator=(const TMatrixD& source) TMatrixD& operator=(const TLazyMatrixD& source) TMatrixD& operator=(Double_t val) Bool_t operator==(Double_t val) const Bool_t operator>(Double_t val) const Bool_t operator>=(Double_t val) const const TMatrixDRow operator[](Int_t rown) const TMatrixDRow operator[](Int_t rown) virtual void Print(Option_t* option) const void ResizeTo(Int_t nrows, Int_t ncols) void ResizeTo(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb) void ResizeTo(const TMatrixD& m) Double_t RowNorm() const void SetElements(const Double_t* elements, Option_t* option) virtual void ShowMembers(TMemberInspector& insp, char* parent) TMatrixD& Sqr() TMatrixD& Sqrt() virtual void Streamer(TBuffer& b) void StreamerNVirtual(TBuffer& b) TMatrixD& UnitMatrix() TMatrixD& Zero()

Data Members


    protected:
Int_t fNrows number of rows Int_t fNcols number of columns Int_t fNelems number of elements in matrix Int_t fRowLwb lower bound of the row index Int_t fColLwb lower bound of the col index Double_t* fElements [fNelems] elements themselves Double_t** fIndex ! index[i] = &matrix(0,i) (col index) public:
static const TMatrixD::EMatrixCreatorsOp1 kZero static const TMatrixD::EMatrixCreatorsOp1 kUnit static const TMatrixD::EMatrixCreatorsOp1 kTransposed static const TMatrixD::EMatrixCreatorsOp1 kInverted static const TMatrixD::EMatrixCreatorsOp1 kInvertedPosDef static const TMatrixD::EMatrixCreatorsOp2 kMult static const TMatrixD::EMatrixCreatorsOp2 kTransposeMult static const TMatrixD::EMatrixCreatorsOp2 kInvMult static const TMatrixD::EMatrixCreatorsOp2 kInvPosDefMult static const TMatrixD::EMatrixCreatorsOp2 kAtBA

Class Description

                                                                      
 Linear Algebra Package                                               
                                                                      
 The present package implements all the basic algorithms dealing      
 with vectors, matrices, matrix columns, rows, diagonals, etc.        
                                                                      
 Matrix elements are arranged in memory in a COLUMN-wise              
 fashion (in FORTRAN's spirit). In fact, it makes it very easy to     
 feed the matrices to FORTRAN procedures, which implement more        
 elaborate algorithms.                                                
                                                                      
 Unless otherwise specified, matrix and vector indices always start   
 with 0, spanning up to the specified limit-1. However, there are     
 constructors to which one can specify aribtrary lower and upper      
 bounds, e.g. TMatrixD m(1,10,1,5) defines a matrix that ranges, and  
 that can be addresses, from 1..10, 1..5 (a(1,1)..a(10,5)).           
                                                                      
 The present package provides all facilities to completely AVOID      
 returning matrices. Use "TMatrixD A(TMatrixD::kTransposed,B);" and   
 other fancy constructors as much as possible. If one really needs    
 to return a matrix, return a TLazyMatrixD object instead. The        
 conversion is completely transparent to the end user, e.g.           
 "TMatrixD m = THaarMatrix(5);" and _is_ efficient.                   
                                                                      
 Since TMatrixD et al. are fully integrated in ROOT they of course    
 can be stored in a ROOT database.                                    
                                                                      
                                                                      
 How to efficiently use this package                                  
 -----------------------------------                                  
                                                                      
 1. Never return complex objects (matrices or vectors)                
    Danger: For example, when the following snippet:                  
        TMatrixD foo(int n)                                           
        {                                                             
           TMatrixD foom(n,n); fill_in(foom); return foom;            
        }                                                             
        TMatrixD m = foo(5);                                          
    runs, it constructs matrix foo:foom, copies it onto stack as a    
    return value and destroys foo:foom. Return value (a matrix)       
    from foo() is then copied over to m (via a copy constructor),     
    and the return value is destroyed. So, the matrix constructor is  
    called 3 times and the destructor 2 times. For big matrices,      
    the cost of multiple constructing/copying/destroying of objects   
    may be very large. *Some* optimized compilers can cut down on 1   
    copying/destroying, but still it leaves at least two calls to     
    the constructor. Note, TLazyMatrices (see below) can construct    
    TMatrixD m "inplace", with only a _single_ call to the            
    constructor.                                                      
                                                                      
 2. Use "two-address instructions"                                    
        "void TMatrixD::operator += (const TMatrixD &B);"             
    as much as possible.                                              
    That is, to add two matrices, it's much more efficient to write   
        A += B;                                                       
    than                                                              
        TMatrixD C = A + B;                                           
    (if both operand should be preserved,                             
        TMatrixD C = A; C += B;                                       
    is still better).                                                 
                                                                      
 3. Use glorified constructors when returning of an object seems      
    inevitable:                                                       
        "TMatrixD A(TMatrixD::kTransposed,B);"                        
        "TMatrixD C(A,TMatrixD::kTransposeMult,B);"                   
                                                                      
    like in the following snippet (from $ROOTSYS/test/vmatrix.cxx)    
    that verifies that for an orthogonal matrix T, T'T = TT' = E.     
                                                                      
    TMatrixD haar = THaarMatrix(5);                                   
    TMatrixD unit(TMatrixD::kUnit,haar);                              
    TMatrixD haar_t(TMatrixD::kTransposed,haar);                      
    TMatrixD hth(haar,TMatrixD::kTransposeMult,haar);                 
    TMatrixD hht(haar,TMatrixD::kMult,haar_t);                        
    TMatrixD hht1 = haar; hht1 *= haar_t;                             
    VerifyMatrixIdentity(unit,hth);                                   
    VerifyMatrixIdentity(unit,hht);                                   
    VerifyMatrixIdentity(unit,hht1);                                  
                                                                      
 4. Accessing row/col/diagonal of a matrix without much fuss          
    (and without moving a lot of stuff around):                       
                                                                      
        TMatrixD m(n,n); TVectorD v(n); TMatrixDDiag(m) += 4;         
        v = TMatrixDRow(m,0);                                         
        TMatrixDColumn m1(m,1); m1(2) = 3; // the same as m(2,1)=3;   
    Note, constructing of, say, TMatrixDDiag does *not* involve any   
    copying of any elements of the source matrix.                     
                                                                      
 5. It's possible (and encouraged) to use "nested" functions          
    For example, creating of a Hilbert matrix can be done as follows: 
                                                                      
    void foo(const TMatrixD &m)                                       
    {                                                                 
      TMatrixD m1(TMatrixD::kZero,m);                                 
      struct MakeHilbert : public TElementPosAction {                 
        void Operation(Double_t &element) { element = 1./(fI+fJ-1); } 
      };                                                              
      m1.Apply(MakeHilbert());                                        
    }                                                                 
                                                                      
    of course, using a special method TMatrixD::HilbertMatrix() is    
    still more optimal, but not by a whole lot. And that's right,     
    class MakeHilbert is declared *within* a function and local to    
    that function. It means one can define another MakeHilbert class  
    (within another function or outside of any function, that is, in  
    the global scope), and it still will be OK. Note, this currently  
    is not yet supported by the interpreter CINT.                     
                                                                      
    Another example is applying of a simple function to each matrix   
    element:                                                          
                                                                      
    void foo(TMatrixD &m, TMatrixD &m1)                               
    {                                                                 
      typedef  double (*dfunc_t)(double);                             
      class ApplyFunction : public TElementActionD {                  
        dfunc_t fFunc;                                                
        void Operation(Double_t &element) { element=fFunc(element); } 
      public:                                                         
        ApplyFunction(dfunc_t func):fFunc(func) {}                    
      };                                                              
      ApplyFunction x(TMath::Sin);                                    
      m.Apply(x);                                                     
    }                                                                 
                                                                      
    Validation code $ROOTSYS/test/vmatrix.cxx and vvector.cxx contain 
    a few more examples of that kind.                                 
                                                                      
 6. Lazy matrices: instead of returning an object return a "recipe"   
    how to make it. The full matrix would be rolled out only when     
    and where it's needed:                                            
       TMatrixD haar = THaarMatrix(5);                                
    THaarMatrix() is a *class*, not a simple function. However        
    similar this looks to a returning of an object (see note #1       
    above), it's dramatically different. THaarMatrix() constructs a   
    TLazyMatrixD, an object of just a few bytes long. A
    "TMatrixD(const TLazyMatrixD &recipe)" constructor follows the    
    recipe and makes the matrix haar() right in place. No matrix      
    element is moved whatsoever!                                      
                                                                      
 The implementation is based on original code by                      
 Oleg E. Kiselyov (oleg@pobox.com).                                   
                                                                      


void Allocate(Int_t no_rows, Int_t no_cols, Int_t row_lwb, Int_t col_lwb)
 Allocate new matrix. Arguments are number of rows, columns, row
 lowerbound (0 default) and column lowerbound (0 default).

~TMatrixD()
 TMatrixD destructor.

void Draw(Option_t *option)
 Draw this matrix using an intermediate histogram
 The histogram is named "TMatrixD" by default and no title

void ResizeTo(Int_t nrows, Int_t ncols)
 Erase the old matrix and create a new one according to new boundaries
 with indexation starting at 0.

void ResizeTo(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)
 Erase the old matrix and create a new one according to new boudaries.

TMatrixD(EMatrixCreatorsOp1 op, const TMatrixD &prototype)
 Create a matrix applying a specific operation to the prototype.
 Example: TMatrixD a(10,12); ...; TMatrixD b(TMatrixD::kTransposed, a);
 Supported operations are: kZero, kUnit, kTransposed, kInverted and kInvertedPosDef.

TMatrixD(const TMatrixD &a, EMatrixCreatorsOp2 op, const TMatrixD &b)
 Create a matrix applying a specific operation to two prototypes.
 Example: TMatrixD a(10,12), b(12,5); ...; TMatrixD c(a, TMatrixD::kMult, b);
 Supported operations are: kMult (a*b), kTransposeMult (a'*b),
 kInvMult,kInvPosDefMult (a^(-1)*b) and kAtBA (a'*b*a).

TMatrixD MakeSymmetric()
 symmetrize matrix (matrix needs to be a square one).

TMatrixD UnitMatrix()
 make a unit matrix (matrix need not be a square one). The matrix
 is traversed in the natural (that is, column by column) order.

TMatrixD HilbertMatrix()
 Make a Hilbert matrix. Hilb[i,j] = 1/(i+j-1), i,j=1...max, OR
 Hilb[i,j] = 1/(i+j+1), i,j=0...max-1 (matrix need not be a square one).
 The matrix is traversed in the natural (that is, column by column) order.

TMatrixD Abs()
 Take an absolute value of a matrix, i.e. apply Abs() to each element.

TMatrixD Sqr()
 Square each element of the matrix.

TMatrixD Sqrt()
 Take square root of all elements.

TMatrixD Apply(TElementPosActionD &action)
 Apply action to each element of the matrix. In action the location
 of the current element is known. The matrix is traversed in the
 natural (that is, column by column) order.

Double_t RowNorm() const
 Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
 The norm is induced by the infinity vector norm.

Double_t ColNorm() const
 Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
 The norm is induced by the 1 vector norm.

Double_t E2Norm() const
 Square of the Euclidian norm, SUM{ m(i,j)^2 }.

TMatrixD NormByDiag(const TVectorD &v, Option_t *option)
 b(i,j) = a(i,j)/sqrt(abs*(v(i)*v(j)))

TMatrixD NormByColumn(const TVectorD &v, Option_t *option)
 Multiply/divide a matrix columns with a vector:
 matrix(i,j) *= v(i)

TMatrixD NormByRow(const TVectorD &v, Option_t *option)
 Multiply/divide a matrix row with a vector:
 matrix(i,j) *= v(j)

void Print(Option_t *) const
 Print the matrix as a table of elements (zeros are printed as dots).

void Transpose(const TMatrixD &prototype)
 Transpose a matrix.

TMatrixD Invert(Double_t *determ_ptr)
 The most general (Gauss-Jordan) matrix inverse

 This method works for any matrix (which of course must be square and
 non-singular). Use this method only if none of specialized algorithms
 (for symmetric, tridiagonal, etc) matrices isn't applicable/available.
 Also, the matrix to invert has to be _well_ conditioned:
 Gauss-Jordan eliminations (even with pivoting) perform poorly for
 near-singular matrices (e.g., Hilbert matrices).

 The method inverts matrix inplace and returns the determinant if
 determ_ptr was specified as not 0. Determinant will be exactly zero
 if the matrix turns out to be (numerically) singular. If determ_ptr is
 0 and matrix happens to be singular, throw up.

 The algorithm perform inplace Gauss-Jordan eliminations with
 full pivoting. It was adapted from my Algol-68 "translation" (ca 1986)
 of the FORTRAN code described in
 Johnson, K. Jeffrey, "Numerical methods in chemistry", New York,
 N.Y.: Dekker, c1980, 503 pp, p.221

 Note, since it's much more efficient to perform operations on matrix
 columns rather than matrix rows (due to the layout of elements in the
 matrix), the present method implements a "transposed" algorithm.

Bool_t IsSymmetric() const

void Invert(const TMatrixD &m)
 Allocate new matrix and set it to inv(m).

TMatrixD InvertPosDef()

Int_t Pdcholesky( const Double_t *a, Double_t *u, const Int_t n)
  Program Pdcholesky inverts a positiv definite (n x n) - matrix A,
  using the Cholesky decomposition

  Input:	a	- (n x n)- Matrix A
  		n	- dimensions n of matrices

  Output:	u	- (n x n)- Matrix U so that U^T . U = A
		return	- 0 decomposition succesful
			- 1 decomposition failed

void InvertPosDef(const TMatrixD &m)
 Allocate new matrix and set it to inv(m).

TMatrixD EigenVectors( TVectorD &eigenValues)
 Return a matrix containing the eigen-vectors; also fill the
 supplied vector with the eigen values.

void MakeTridiagonal( TMatrixD &a, TVectorD &d, TVectorD &e)
 The comments in this algorithm are modified version of those in
 "Numerical ...". Please refer to that book (web-page) for more on
 the algorithm.
 
  /*
    
PRIVATE METHOD:
Tridiagonalise the covariance matrix according to the Householder method as described in Numerical Recipes in C section 11.2.

The basic idea is to perform $P-2$ orthogonal transformation, where each transformation eat away the off-diagonal elements, except the inner most.

   */
 

void MakeEigenVectors( TVectorD &d, TVectorD &e, TMatrixD &z)
 
  /*
    
PRIVATE METHOD:
Find eigenvalues and vectors of tridiagonalised covariance matrix according to the QL with implicit shift algorithm from Numerical Recipes in C section 11.3.

The basic idea is to find matrices $\mathsf{Q}$ and $\mathsf{L}$ so that $\mathsf{C} = \mathsf{Q} \cdot \mathsf{L}$, where $\mathsf{Q}$ is orthogonal and $\mathsf{L}$ is lower triangular. The QL algorithm consist of a sequence of orthogonal transformations

\begin{displaymath}
    \mathsf{C}_s = \mathsf{Q}_s \cdot \mathsf{L}_s
    \end{displaymath}


\begin{displaymath}
    \mathsf{C}_{s+1} = \mathsf{L}_s \cdot \mathsf{Q}_s
    = \mathsf{Q}_s^T \cdot \mathsf{C}_s \cdot \mathsf{Q}_s
    \end{displaymath}

(1) If $\mathsf{C}$ have eigenvalues with different absolute value $\vert l_i\vert$, then $\mathsf{C}_s \rightarrow$ [lower triangular form] as $s\rightarrow\infty$. The eigenvalues appear on the diagonal in increasing order of absolute magnitude. (2) If If $\mathsf{C}$ has an eigenvalue $\vert l_i\vert$ of multiplicty of order $p$, $\mathsf{C}_s \rightarrow$ [lower triangular form] as $s\rightarrow\infty$, except for a diagona block matrix of order $p$, whose eigenvalues $\rightarrow l_i$.
  */
 

void EigenSort( TMatrixD &eigenVectors, TVectorD &eigenValues)
 
  /*
    
PRIVATE METHOD:
Order the eigenvalues and vectors by ascending eigenvalue. The algorithm is a straight insertion. It's taken from Numerical Recipes in C section 11.1.
  */
 

void AMultB(const TMatrixD &a, const TMatrixD &b)
 General matrix multiplication. Create a matrix C such that C = A * B.
 Note, matrix C needs to be allocated.

void Mult(const TMatrixD &a, const TMatrixD &b)
 Compute C = A*B. The same as AMultB(), only matrix C is already
 allocated, and it is *this.

void AtMultB(const TMatrixD &a, const TMatrixD &b)
 Create a matrix C such that C = A' * B. In other words,
 c[i,j] = SUM{ a[k,i] * b[k,j] }. Note, matrix C needs to be allocated.

Double_t Determinant() const
 Compute the determinant of a general square matrix.
 Example: Matrix A; Double_t A.Determinant();

 Gauss-Jordan transformations of the matrix with a slight
 modification to take advantage of the *column*-wise arrangement
 of Matrix elements. Thus we eliminate matrix's columns rather than
 rows in the Gauss-Jordan transformations. Note that determinant
 is invariant to matrix transpositions.
 The matrix is copied to a special object of type TMatrixDPivoting,
 where all Gauss-Jordan eliminations with full pivoting are to
 take place.

void Streamer(TBuffer &R__b)
 Stream an object of class TMatrixD.

TMatrixD(Int_t no_rows, Int_t no_cols)

TMatrixD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb)

Bool_t IsValid() const

void SetElements(const Double_t *elements, Option_t *option)

TMatrixD(Int_t no_rows, Int_t no_cols, const Double_t *elements, Option_t *option)
 option="F": array elements contains the matrix stored column-wise
             like in Fortran, so a[i,j] = elements[i+no_rows*j],
 else        it is supposed that array elements are stored row-wise
             a[i,j] = elements[i*no_cols+j]

TMatrixD(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, const Double_t *elements, Option_t *option)

void GetElements(Double_t *elements, Option_t *option) const

TMatrixD(const TLazyMatrixD &lazy_constructor)

TMatrixD(const TMatrixD &another) : TObject(another)

void ResizeTo(const TMatrixD &m)

TMatrixD Zero()

TMatrixD Apply(TElementActionD &action)



Inline Functions


                     void Invalidate()
                 TMatrixD TMatrixD(const TLazyMatrixD& lazy_constructor)
                    Int_t GetRowLwb() const
                    Int_t GetRowUpb() const
                    Int_t GetNrows() const
                    Int_t GetColLwb() const
                    Int_t GetColUpb() const
                    Int_t GetNcols() const
                    Int_t GetNoElements() const
                     void GetElements(Double_t* elements, Option_t* option) const
          const Double_t& operator()(Int_t rown, Int_t coln) const
                Double_t& operator()(Int_t rown, Int_t coln)
        const TMatrixDRow operator[](Int_t rown) const
              TMatrixDRow operator[](Int_t rown)
                TMatrixD& operator=(const TMatrixD& source)
                TMatrixD& operator=(const TLazyMatrixD& source)
                TMatrixD& operator=(Double_t val)
                TMatrixD& operator-=(Double_t val)
                TMatrixD& operator+=(Double_t val)
                TMatrixD& operator*=(Double_t val)
                   Bool_t operator==(Double_t val) const
                   Bool_t operator!=(Double_t val) const
                   Bool_t operator<(Double_t val) const
                   Bool_t operator<=(Double_t val) const
                   Bool_t operator>(Double_t val) const
                   Bool_t operator>=(Double_t val) const
                TMatrixD& operator*=(const TMatrixD& source)
                TMatrixD& operator*=(const TMatrixDDiag& diag)
                TMatrixD& operator/=(const TMatrixDDiag& diag)
                TMatrixD& operator*=(const TMatrixDRow& diag)
                TMatrixD& operator/=(const TMatrixDRow& diag)
                TMatrixD& operator*=(const TMatrixDColumn& diag)
                TMatrixD& operator/=(const TMatrixDColumn& diag)
                 Double_t NormInf() const
                 Double_t Norm1() const
                  TClass* Class()
                  TClass* IsA() const
                     void ShowMembers(TMemberInspector& insp, char* parent)
                     void StreamerNVirtual(TBuffer& b)


Author: Fons Rademakers 03/11/97
Last update: root/matrix:$Name: $:$Id: TMatrixD.cxx,v 1.26 2002/09/15 10:16:44 brun Exp $
Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *


ROOT page - Class index - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.