orthogonal transformation, where
each transformation eat away the off-diagonal elements, except the
inner most.
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()
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
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).
Allocate new matrix. Arguments are number of rows, columns, row lowerbound (0 default) and column lowerbound (0 default).
TMatrixD destructor.
Draw this matrix using an intermediate histogram The histogram is named "TMatrixD" by default and no title
Erase the old matrix and create a new one according to new boundaries with indexation starting at 0.
Erase the old matrix and create a new one according to new boudaries.
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.
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).
symmetrize matrix (matrix needs to be a square one).
make a unit matrix (matrix need not be a square one). The matrix is traversed in the natural (that is, column by column) order.
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.
Take an absolute value of a matrix, i.e. apply Abs() to each element.
Square each element of the matrix.
Take square root of all elements.
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.
Row matrix norm, MAX{ SUM{ |M(i,j)|, over j}, over i}.
The norm is induced by the infinity vector norm.
Column matrix norm, MAX{ SUM{ |M(i,j)|, over i}, over j}.
The norm is induced by the 1 vector norm.
Square of the Euclidian norm, SUM{ m(i,j)^2 }.
b(i,j) = a(i,j)/sqrt(abs*(v(i)*v(j)))
Multiply/divide a matrix columns with a vector: matrix(i,j) *= v(i)
Multiply/divide a matrix row with a vector: matrix(i,j) *= v(j)
Print the matrix as a table of elements (zeros are printed as dots).
Transpose a matrix.
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.
Allocate new matrix and set it to inv(m).
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
Allocate new matrix and set it to inv(m).
Return a matrix containing the eigen-vectors; also fill the supplied vector with the eigen values.
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:
The basic idea is to perform orthogonal transformation, where
each transformation eat away the off-diagonal elements, except the
inner most.
*/
/*
PRIVATE METHOD:
The basic idea is to find matrices 
and 
so that
, where is orthogonal and
is lower triangular. The QL algorithm
consist of a
sequence of orthogonal transformations
have eigenvalues with different absolute value
, then
[lower triangular form] as
. The eigenvalues appear on the diagonal in
increasing order of absolute magnitude. (2) If If has an
eigenvalue of multiplicty of order 
,
[lower triangular form] as
, except for a diagona block matrix of order ,
whose eigenvalues
.
*/
/*
PRIVATE METHOD:
*/
General matrix multiplication. Create a matrix C such that C = A * B. Note, matrix C needs to be allocated.
Compute C = A*B. The same as AMultB(), only matrix C is already allocated, and it is *this.
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.
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.
Stream an object of class TMatrixD.
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]
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)