// @(#)root/matrix:$Name: $:$Id: TMatrix.h,v 1.25 2003/09/05 09:21:54 brun Exp $ // Authors: Oleg E. Kiselyov, Fons Rademakers 03/11/97 /************************************************************************* * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. * * All rights reserved. * * * * For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see $ROOTSYS/README/CREDITS. * *************************************************************************/ #ifndef ROOT_TMatrix #define ROOT_TMatrix ////////////////////////////////////////////////////////////////////////// // // // 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. // // // // The present package provides all facilities to completely AVOID // // returning matrices. Use "TMatrix A(TMatrix::kTransposed,B);" and // // other fancy constructors as much as possible. If one really needs // // to return a matrix, return a TLazyMatrix object instead. The // // conversion is completely transparent to the end user, e.g. // // "TMatrix m = THaarMatrix(5);" and _is_ efficient. // // // // For usage examples see $ROOTSYS/test/vmatrix.cxx and vvector.cxx // // and also: // // http://root.cern.ch/root/html/TMatrix.html#TMatrix:description // // // // The implementation is based on original code by // // Oleg E. Kiselyov (oleg@pobox.com). // // Several additions/optimisations by Eddy Offermann // // // ////////////////////////////////////////////////////////////////////////// #ifndef ROOT_TVector #include "TVector.h" #endif class TMatrix; class TLazyMatrix; class TMatrixRow; class TMatrixColumn; class TMatrixDiag; class TMatrixFlat; class TMatrixPivoting; TMatrix &operator+=(TMatrix &target, const TMatrix &source); TMatrix &operator-=(TMatrix &target, const TMatrix &source); TMatrix operator+(const TMatrix &source1, const TMatrix &source2); TMatrix operator-(const TMatrix &source1, const TMatrix &source2); TMatrix operator*(const TMatrix &source1, const TMatrix &source2); TMatrix &Add(TMatrix &target, Double_t scalar, const TMatrix &source); TMatrix &ElementMult(TMatrix &target, const TMatrix &source); TMatrix &ElementDiv(TMatrix &target, const TMatrix &source); Bool_t operator==(const TMatrix &im1, const TMatrix &im2); void Compare(const TMatrix &im1, const TMatrix &im2); Bool_t AreCompatible(const TMatrix &im1, const TMatrix &im2); Double_t E2Norm(const TMatrix &m1, const TMatrix &m2); class TMatrix : public TObject { friend class TVector; friend class TMatrixRow; friend class TMatrixColumn; friend class TMatrixDiag; friend class TMatrixFlat; friend class TMatrixPivoting; 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 Real_t *fElements; //[fNelems] elements themselves Real_t **fIndex; //! index[i] = &matrix(0,i) (col index) void Allocate(Int_t nrows, Int_t ncols, Int_t row_lwb = 0, Int_t col_lwb = 0); void Invalidate() { fNrows = fNcols = fNelems = -1; fElements = 0; fIndex = 0; } static Int_t Pdcholesky(const Real_t *a, Real_t *u, const Int_t n); static void MakeTridiagonal(TMatrix &a,TVector &d,TVector &e); static void MakeEigenVectors(TVector &d,TVector &e,TMatrix &z); static void EigenSort(TMatrix &eigenVectors,TVector &eigenValues); // Elementary constructors void Transpose(const TMatrix &m); void Invert(const TMatrix &m); void InvertPosDef(const TMatrix &m); void AMultB(const TMatrix &a, const TMatrix &b); void AtMultB(const TMatrix &a, const TMatrix &b); friend void MakeHaarMatrix(TMatrix &m); friend void MakeHilbertMatrix(TMatrix &m); public: enum EMatrixCreatorsOp1 { kZero, kUnit, kTransposed, kInverted, kInvertedPosDef }; enum EMatrixCreatorsOp2 { kMult, kTransposeMult, kInvMult, kInvPosDefMult, kAtBA }; TMatrix() { Invalidate(); } TMatrix(Int_t nrows, Int_t ncols); TMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb); TMatrix(Int_t nrows, Int_t ncols, const Real_t *elements, Option_t *option=""); TMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, const Real_t *elements, Option_t *option=""); TMatrix(const TMatrix &another); TMatrix(EMatrixCreatorsOp1 op, const TMatrix &prototype); TMatrix(const TMatrix &a, EMatrixCreatorsOp2 op, const TMatrix &b); TMatrix(const TLazyMatrix &lazy_constructor); virtual ~TMatrix(); void Clear(Option_t *option=""); void Draw(Option_t *option=""); // *MENU* 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 TMatrix &m); Bool_t IsValid() const; Bool_t IsSymmetric() const; Int_t GetRowLwb() const { return fRowLwb; } Int_t GetRowUpb() const { return fNrows+fRowLwb-1; } Int_t GetNrows() const { return fNrows; } Int_t GetColLwb() const { return fColLwb; } Int_t GetColUpb() const { return fNcols+fColLwb-1; } Int_t GetNcols() const { return fNcols; } Int_t GetNoElements() const { return fNelems; } const Real_t *GetElements() const { return fElements; } Real_t *GetElements() { return fElements; } void GetElements(Real_t *elements, Option_t *option="") const; void SetElements(const Real_t *elements, Option_t *option=""); const Real_t &operator()(Int_t rown, Int_t coln) const; Real_t &operator()(Int_t rown, Int_t coln); const TMatrixRow operator[](Int_t rown) const; TMatrixRow operator[](Int_t rown); TMatrix &operator=(const TMatrix &source); TMatrix &operator=(const TLazyMatrix &source); TMatrix &operator=(Real_t val); TMatrix &operator-=(Double_t val); TMatrix &operator+=(Double_t val); TMatrix &operator*=(Double_t val); Bool_t operator==(Real_t val) const; Bool_t operator!=(Real_t val) const; Bool_t operator<(Real_t val) const; Bool_t operator<=(Real_t val) const; Bool_t operator>(Real_t val) const; Bool_t operator>=(Real_t val) const; TMatrix &Zero(); TMatrix &Abs(); TMatrix &Sqr(); TMatrix &Sqrt(); TMatrix &Apply(const TElementAction &action); TMatrix &Apply(const TElementPosAction &action); TMatrix &Invert(Double_t *determ_ptr = 0); TMatrix &InvertPosDef(); const TMatrix EigenVectors(TVector &eigenValues) const; TMatrix &MakeSymmetric(); TMatrix &UnitMatrix(); TMatrix &operator*=(const TMatrix &source); TMatrix &operator*=(const TMatrixDiag &diag); TMatrix &operator/=(const TMatrixDiag &diag); TMatrix &operator*=(const TMatrixRow &row); TMatrix &operator/=(const TMatrixRow &row); TMatrix &operator*=(const TMatrixColumn &col); TMatrix &operator/=(const TMatrixColumn &col); void Mult(const TMatrix &a, const TMatrix &b); Double_t RowNorm() const; Double_t NormInf() const { return RowNorm(); } Double_t ColNorm() const; Double_t Norm1() const { return ColNorm(); } Double_t E2Norm() const; TMatrix &NormByDiag(const TVector &v, Option_t *option="D"); TMatrix &NormByColumn(const TVector &v, Option_t *option="D"); TMatrix &NormByRow(const TVector &v, Option_t *option="D"); Double_t Determinant() const; void Print(Option_t *option="") const; // *MENU* friend TMatrix &operator+=(TMatrix &target, const TMatrix &source); friend TMatrix &operator-=(TMatrix &target, const TMatrix &source); friend TMatrix &Add(TMatrix &target, Double_t scalar, const TMatrix &source); friend TMatrix &ElementMult(TMatrix &target, const TMatrix &source); friend TMatrix &ElementDiv(TMatrix &target, const TMatrix &source); friend TMatrix operator+(const TMatrix &source1, const TMatrix &source2); friend TMatrix operator-(const TMatrix &source1, const TMatrix &source2); friend TMatrix operator*(const TMatrix &source1, const TMatrix &source2); friend Bool_t operator==(const TMatrix &im1, const TMatrix &im2); friend void Compare(const TMatrix &im1, const TMatrix &im2); friend Bool_t AreCompatible(const TMatrix &im1, const TMatrix &im2); friend Double_t E2Norm(const TMatrix &m1, const TMatrix &m2); ClassDef(TMatrix,2) // Matrix class }; // Service functions (useful in the verification code). // They print some detail info if the validation condition fails void VerifyElementValue(const TMatrix &m, Real_t val); void VerifyMatrixIdentity(const TMatrix &m1, const TMatrix &m2); inline Bool_t TMatrix::IsValid() const { if (fNrows == -1) return kFALSE; return kTRUE; } #ifndef ROOT_TMatrixUtils #include "TMatrixUtils.h" #endif #endif