// @(#)root/matrix:$Name: $:$Id: TMatrixDUtils.h,v 1.30 2005/04/07 14:43:35 rdm Exp $ // Authors: Fons Rademakers, Eddy Offermann Nov 2003 /************************************************************************* * 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_TMatrixDUtils #define ROOT_TMatrixDUtils ////////////////////////////////////////////////////////////////////////// // // // Matrix utility classes. // // // // This file defines utility classes for the Linear Algebra Package. // // The following classes are defined here: // // // // Different matrix views without copying data elements : // // TMatrixDRow_const TMatrixDRow // // TMatrixDColumn_const TMatrixDColumn // // TMatrixDDiag_const TMatrixDDiag // // TMatrixDFlat_const TMatrixDFlat // // TMatrixDSub_const TMatrixDSub // // TMatrixDSparseRow_const TMatrixDSparseRow // // TMatrixDSparseDiag_const TMatrixDSparseDiag // // // // TElementActionD // // TElementPosActionD // // // ////////////////////////////////////////////////////////////////////////// #ifndef ROOT_TMatrixDBase #include "TMatrixDBase.h" #endif class TVectorD; class TMatrixDBase; class TMatrixD; class TMatrixDSym; ////////////////////////////////////////////////////////////////////////// // // // TElementActionD // // // // A class to do a specific operation on every vector or matrix element // // (regardless of it position) as the object is being traversed. // // This is an abstract class. Derived classes need to implement the // // action function Operation(). // // // ////////////////////////////////////////////////////////////////////////// class TElementActionD { friend class TMatrixDBase; friend class TMatrixD; friend class TMatrixDSym; friend class TMatrixDSparse; friend class TVectorD; protected: virtual ~TElementActionD() { } virtual void Operation(Double_t &element) const = 0; private: void operator=(const TElementActionD &) { } }; ////////////////////////////////////////////////////////////////////////// // // // TElementPosActionD // // // // A class to do a specific operation on every vector or matrix element // // as the object is being traversed. This is an abstract class. // // Derived classes need to implement the action function Operation(). // // In the action function the location of the current element is // // known (fI=row, fJ=columns). // // // ////////////////////////////////////////////////////////////////////////// class TElementPosActionD { friend class TMatrixDBase; friend class TMatrixD; friend class TMatrixDSym; friend class TMatrixDSparse; friend class TVectorD; protected: mutable Int_t fI; // i position of element being passed to Operation() mutable Int_t fJ; // j position of element being passed to Operation() virtual ~TElementPosActionD() { } virtual void Operation(Double_t &element) const = 0; private: void operator=(const TElementPosActionD &) { } }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDRow_const // // // // Class represents a row of a TMatrixDBase // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDRow_const { protected: const TMatrixDBase *fMatrix; // the matrix I am a row of Int_t fRowInd; // effective row index Int_t fInc; // if ptr = @a[row,i], then ptr+inc = @a[row,i+1] const Double_t *fPtr; // pointer to the a[row,0] public: TMatrixDRow_const() { fMatrix = 0; fInc = 0; fPtr = 0; } TMatrixDRow_const(const TMatrixD &matrix,Int_t row); TMatrixDRow_const(const TMatrixDSym &matrix,Int_t row); virtual ~TMatrixDRow_const() { } inline const TMatrixDBase *GetMatrix () const { return fMatrix; } inline Int_t GetRowIndex() const { return fRowInd; } inline Int_t GetInc () const { return fInc; } inline const Double_t *GetPtr () const { return fPtr; } inline const Double_t &operator ()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t acoln = i-fMatrix->GetColLwb(); Assert(acoln < fMatrix->GetNcols() && acoln >= 0); return fPtr[acoln]; } inline const Double_t &operator [](Int_t i) const { return (*(const TMatrixDRow_const *)this)(i); } ClassDef(TMatrixDRow_const,0) // One row of a dense matrix (double precision) }; class TMatrixDRow : public TMatrixDRow_const { public: TMatrixDRow() {} TMatrixDRow(TMatrixD &matrix,Int_t row); TMatrixDRow(TMatrixDSym &matrix,Int_t row); TMatrixDRow(const TMatrixDRow &mr); inline Double_t *GetPtr() const { return const_cast(fPtr); } inline const Double_t &operator()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t acoln = i-fMatrix->GetColLwb(); Assert(acoln < fMatrix->GetNcols() && acoln >= 0); return fPtr[acoln]; } inline Double_t &operator()(Int_t i) { Assert(fMatrix->IsValid()); const Int_t acoln = i-fMatrix->GetColLwb(); Assert(acoln < fMatrix->GetNcols() && acoln >= 0); return (const_cast(fPtr))[acoln]; } inline const Double_t &operator[](Int_t i) const { return (*(const TMatrixDRow *)this)(i); } inline Double_t &operator[](Int_t i) { return (*( TMatrixDRow *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDRow_const &r); void operator=(const TMatrixDRow &r) { operator=((TMatrixDRow_const &)r); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDRow_const &r); void operator*=(const TMatrixDRow_const &r); ClassDef(TMatrixDRow,0) // One row of a dense matrix (double precision) }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDColumn_const // // // // Class represents a column of a TMatrixDBase // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDColumn_const { protected: const TMatrixDBase *fMatrix; // the matrix I am a column of Int_t fColInd; // effective column index Int_t fInc; // if ptr = @a[i,col], then ptr+inc = @a[i+1,col] const Double_t *fPtr; // pointer to the a[0,col] column public: TMatrixDColumn_const() { fMatrix = 0; fInc = 0; fPtr = 0; } TMatrixDColumn_const(const TMatrixD &matrix,Int_t col); TMatrixDColumn_const(const TMatrixDSym &matrix,Int_t col); virtual ~TMatrixDColumn_const() { } inline const TMatrixDBase *GetMatrix () const { return fMatrix; } inline Int_t GetColIndex() const { return fColInd; } inline Int_t GetInc () const { return fInc; } inline const Double_t *GetPtr () const { return fPtr; } inline const Double_t &operator ()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t arown = i-fMatrix->GetRowLwb(); Assert(arown < fMatrix->GetNrows() && arown >= 0); return fPtr[arown*fInc]; } inline const Double_t &operator [](Int_t i) const { return (*(const TMatrixDColumn_const *)this)(i); } ClassDef(TMatrixDColumn_const,0) // One column of a dense matrix (double precision) }; class TMatrixDColumn : public TMatrixDColumn_const { public: TMatrixDColumn() {} TMatrixDColumn(TMatrixD &matrix,Int_t col); TMatrixDColumn(TMatrixDSym &matrix,Int_t col); TMatrixDColumn(const TMatrixDColumn &mc); inline Double_t *GetPtr() const { return const_cast(fPtr); } inline const Double_t &operator()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t arown = i-fMatrix->GetRowLwb(); Assert(arown < fMatrix->GetNrows() && arown >= 0); return fPtr[arown]; } inline Double_t &operator()(Int_t i) { Assert(fMatrix->IsValid()); const Int_t arown = i-fMatrix->GetRowLwb(); Assert(arown < fMatrix->GetNrows() && arown >= 0); return (const_cast(fPtr))[arown*fInc]; } inline const Double_t &operator[](Int_t i) const { return (*(const TMatrixDColumn *)this)(i); } inline Double_t &operator[](Int_t i) { return (*( TMatrixDColumn *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDColumn_const &c); void operator=(const TMatrixDColumn &c) { operator=((TMatrixDColumn_const &)c); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDColumn_const &c); void operator*=(const TMatrixDColumn_const &c); ClassDef(TMatrixDColumn,0) // One column of a dense matrix (double precision) }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDDiag_const // // // // Class represents the diagonal of a matrix (for easy manipulation). // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDDiag_const { protected: const TMatrixDBase *fMatrix; // the matrix I am the diagonal of Int_t fInc; // if ptr=@a[i,i], then ptr+inc = @a[i+1,i+1] Int_t fNdiag; // number of diag elems, min(nrows,ncols) const Double_t *fPtr; // pointer to the a[0,0] public: TMatrixDDiag_const() { fMatrix = 0; fInc = 0; fNdiag = 0; fPtr = 0; } TMatrixDDiag_const(const TMatrixD &matrix); TMatrixDDiag_const(const TMatrixDSym &matrix); virtual ~TMatrixDDiag_const() { } inline const TMatrixDBase *GetMatrix() const { return fMatrix; } inline const Double_t *GetPtr () const { return fPtr; } inline Int_t GetInc () const { return fInc; } inline const Double_t &operator ()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i < fNdiag && i >= 0); return fPtr[i*fInc]; } inline const Double_t &operator [](Int_t i) const { return (*(const TMatrixDDiag_const *)this)(i); } Int_t GetNdiags() const { return fNdiag; } ClassDef(TMatrixDDiag_const,0) // Diagonal of a dense matrix (double precision) }; class TMatrixDDiag : public TMatrixDDiag_const { public: TMatrixDDiag() {} TMatrixDDiag(TMatrixD &matrix); TMatrixDDiag(TMatrixDSym &matrix); TMatrixDDiag(const TMatrixDDiag &md); inline Double_t *GetPtr() const { return const_cast(fPtr); } inline const Double_t &operator()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i < fNdiag && i >= 0); return fPtr[i*fInc]; } inline Double_t &operator()(Int_t i) { Assert(fMatrix->IsValid()); Assert(i < fNdiag && i >= 0); return (const_cast(fPtr))[i*fInc]; } inline const Double_t &operator[](Int_t i) const { return (*(const TMatrixDDiag *)this)(i); } inline Double_t &operator[](Int_t i) { return (*( TMatrixDDiag *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDDiag_const &d); void operator=(const TMatrixDDiag &d) { operator=((TMatrixDDiag_const &)d); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDDiag_const &d); void operator*=(const TMatrixDDiag_const &d); ClassDef(TMatrixDDiag,0) // Diagonal of a dense matrix (double precision) }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDFlat_const // // // // Class represents a flat matrix (for easy manipulation). // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDFlat_const { protected: const TMatrixDBase *fMatrix; // the matrix I am the diagonal of Int_t fNelems; // const Double_t *fPtr; // pointer to the a[0,0] public: TMatrixDFlat_const() { fMatrix = 0; fNelems = 0; fPtr = 0; } TMatrixDFlat_const(const TMatrixD &matrix); TMatrixDFlat_const(const TMatrixDSym &matrix); virtual ~TMatrixDFlat_const() { } inline const TMatrixDBase *GetMatrix() const { return fMatrix; } inline const Double_t *GetPtr () const { return fPtr; } inline const Double_t &operator ()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i >=0 && i < fNelems); return fPtr[i]; } inline const Double_t &operator [](Int_t i) const { return (*(const TMatrixDFlat_const *)this)(i); } ClassDef(TMatrixDFlat_const,0) // Flat representation of a dense matrix }; class TMatrixDFlat : public TMatrixDFlat_const { public: TMatrixDFlat() {} TMatrixDFlat(TMatrixD &matrix); TMatrixDFlat(TMatrixDSym &matrix); TMatrixDFlat(const TMatrixDFlat &mf); inline Double_t *GetPtr() const { return const_cast(fPtr); } inline const Double_t &operator()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i >=0 && i < fNelems); return fPtr[i]; } inline Double_t &operator()(Int_t i) { Assert(fMatrix->IsValid()); Assert(i >=0 && i < fNelems); return (const_cast(fPtr))[i]; } inline const Double_t &operator[](Int_t i) const { return (*(const TMatrixDFlat *)this)(i); } inline Double_t &operator[](Int_t i) { return (*( TMatrixDFlat *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDFlat_const &f); void operator=(const TMatrixDFlat &f) { operator=((TMatrixDFlat_const &)f); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDFlat_const &f); void operator*=(const TMatrixDFlat_const &f); ClassDef(TMatrixDFlat,0) // Flat representation of a dense matrix }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDSub_const // // // // Class represents a sub matrix of a TMatrixDBase // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDSub_const { protected: const TMatrixDBase *fMatrix; // the matrix I am a submatrix of Int_t fRowOff; // Int_t fColOff; // Int_t fNrowsSub; // Int_t fNcolsSub; // enum {kWorkMax = 100}; public: TMatrixDSub_const() { fRowOff = fColOff = fNrowsSub = fNcolsSub = 0; fMatrix = 0; } TMatrixDSub_const(const TMatrixD &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb); TMatrixDSub_const(const TMatrixDSym &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb); virtual ~TMatrixDSub_const() { } inline const TMatrixDBase *GetMatrix() const { return fMatrix; } inline Int_t GetRowOff() const { return fRowOff; } inline Int_t GetColOff() const { return fColOff; } inline Int_t GetNrows () const { return fNrowsSub; } inline Int_t GetNcols () const { return fNcolsSub; } inline const Double_t &operator ()(Int_t rown,Int_t coln) const { Assert(fMatrix->IsValid()); Assert(rown < fNrowsSub && rown >= 0); Assert(coln < fNcolsSub && coln >= 0); const Int_t index = (rown+fRowOff)*fMatrix->GetNcols()+coln+fColOff; const Double_t *ptr = fMatrix->GetMatrixArray(); return ptr[index]; } ClassDef(TMatrixDSub_const,0) // sub matrix (double precision) }; class TMatrixDSub : public TMatrixDSub_const { public: TMatrixDSub() {} TMatrixDSub(TMatrixD &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb); TMatrixDSub(TMatrixDSym &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb); TMatrixDSub(const TMatrixDSub &ms); inline Double_t &operator()(Int_t rown,Int_t coln) { Assert(fMatrix->IsValid()); Assert(rown < fNrowsSub && rown >= 0); Assert(coln < fNcolsSub && coln >= 0); const Int_t index = (rown+fRowOff)*fMatrix->GetNcols()+coln+fColOff; const Double_t *ptr = fMatrix->GetMatrixArray(); return (const_cast(ptr))[index]; } void Rank1Update(const TVectorD &vec,Double_t alpha=1.0); void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDSub_const &s); void operator=(const TMatrixDSub &s) { operator=((TMatrixDSub_const &)s); } void operator=(const TMatrixDBase &m); void operator+=(const TMatrixDSub_const &s); void operator*=(const TMatrixDSub_const &s); void operator+=(const TMatrixDBase &m); void operator*=(const TMatrixD &m); void operator*=(const TMatrixDSym &m); ClassDef(TMatrixDSub,0) // sub matrix (double precision) }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDSparseRow_const // // // // Class represents a row of a TMatrixDSparse // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDSparse; class TMatrixDSparseRow_const { protected: const TMatrixDBase *fMatrix; // the matrix I am a row of Int_t fRowInd; // effective row index Int_t fNindex; // index range const Int_t *fColPtr; // column index pointer const Double_t *fDataPtr; // data pointer public: TMatrixDSparseRow_const() { fMatrix = 0; fRowInd = 0; fNindex = 0; fColPtr = 0; fDataPtr = 0; } TMatrixDSparseRow_const(const TMatrixDSparse &matrix,Int_t row); virtual ~TMatrixDSparseRow_const() { } inline const TMatrixDBase *GetMatrix () const { return fMatrix; } inline const Double_t *GetDataPtr () const { return fDataPtr; } inline const Int_t *GetColPtr () const { return fColPtr; } inline Int_t GetRowIndex() const { return fRowInd; } inline Int_t GetNindex () const { return fNindex; } inline Double_t operator()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t acoln = i-fMatrix->GetColLwb(); Assert(acoln < fMatrix->GetNcols() && acoln >= 0); const Int_t index = TMath::BinarySearch(fNindex,fColPtr,acoln); if (index >= 0 && fColPtr[index] == acoln) return fDataPtr[index]; else return 0.0; } inline Double_t operator[](Int_t i) const { return (*(const TMatrixDSparseRow_const *)this)(i); } ClassDef(TMatrixDSparseRow_const,0) // One row of a sparse matrix (double precision) }; class TMatrixDSparseRow : public TMatrixDSparseRow_const { public: TMatrixDSparseRow() {} TMatrixDSparseRow(TMatrixDSparse &matrix,Int_t row); TMatrixDSparseRow(const TMatrixDSparseRow &mr); inline Double_t *GetDataPtr() const { return const_cast(fDataPtr); } inline Double_t operator()(Int_t i) const { Assert(fMatrix->IsValid()); const Int_t acoln = i-fMatrix->GetColLwb(); Assert(acoln < fMatrix->GetNcols() && acoln >= 0); const Int_t index = TMath::BinarySearch(fNindex,fColPtr,acoln); if (index >= 0 && fColPtr[index] == acoln) return fDataPtr[index]; else return 0.0; } Double_t &operator()(Int_t i); inline Double_t operator[](Int_t i) const { return (*(const TMatrixDSparseRow *)this)(i); } inline Double_t &operator[](Int_t i) { return (*(TMatrixDSparseRow *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDSparseRow_const &r); void operator=(const TMatrixDSparseRow &r) { operator=((TMatrixDSparseRow_const &)r); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDSparseRow_const &r); void operator*=(const TMatrixDSparseRow_const &r); ClassDef(TMatrixDSparseRow,0) // One row of a sparse matrix (double precision) }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDSparseDiag_const // // // // Class represents the diagonal of a matrix (for easy manipulation). // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDSparseDiag_const { protected: const TMatrixDBase *fMatrix; // the matrix I am the diagonal of Int_t fNdiag; // number of diag elems, min(nrows,ncols) const Double_t *fDataPtr; // data pointer public: TMatrixDSparseDiag_const() { fMatrix = 0; fNdiag = 0; fDataPtr = 0; } TMatrixDSparseDiag_const(const TMatrixDSparse &matrix); virtual ~TMatrixDSparseDiag_const() { } inline const TMatrixDBase *GetMatrix () const { return fMatrix; } inline const Double_t *GetDataPtr() const { return fDataPtr; } inline Int_t GetNdiags () const { return fNdiag; } inline Double_t operator ()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i < fNdiag && i >= 0); const Int_t * const pR = fMatrix->GetRowIndexArray(); const Int_t * const pC = fMatrix->GetColIndexArray(); const Double_t * const pD = fMatrix->GetMatrixArray(); const Int_t sIndex = pR[i]; const Int_t eIndex = pR[i+1]; const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex; if (index >= sIndex && pC[index] == i) return pD[index]; else return 0.0; } inline Double_t operator [](Int_t i) const { return (*(const TMatrixDSparseRow_const *)this)(i); } ClassDef(TMatrixDSparseDiag_const,0) // Diagonal of a sparse matrix (double precision) }; class TMatrixDSparseDiag : public TMatrixDSparseDiag_const { public: TMatrixDSparseDiag() {} TMatrixDSparseDiag(TMatrixDSparse &matrix); TMatrixDSparseDiag(const TMatrixDSparseDiag &md); inline Double_t *GetDataPtr() const { return const_cast(fDataPtr); } inline Double_t operator()(Int_t i) const { Assert(fMatrix->IsValid()); Assert(i < fNdiag && i >= 0); const Int_t * const pR = fMatrix->GetRowIndexArray(); const Int_t * const pC = fMatrix->GetColIndexArray(); const Double_t * const pD = fMatrix->GetMatrixArray(); const Int_t sIndex = pR[i]; const Int_t eIndex = pR[i+1]; const Int_t index = TMath::BinarySearch(eIndex-sIndex,pC+sIndex,i)+sIndex; if (index >= sIndex && pC[index] == i) return pD[index]; else return 0.0; } Double_t &operator()(Int_t i); inline Double_t operator[](Int_t i) const { return (*(const TMatrixDSparseDiag *)this)(i); } inline Double_t &operator[](Int_t i) { return (*(TMatrixDSparseDiag *)this)(i); } void operator= (Double_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TMatrixDSparseDiag_const &d); void operator=(const TMatrixDSparseDiag &d) { operator=((TMatrixDSparseDiag_const &)d); } void operator=(const TVectorD &vec); void operator+=(const TMatrixDSparseDiag_const &d); void operator*=(const TMatrixDSparseDiag_const &d); ClassDef(TMatrixDSparseDiag,0) // Diagonal of a dense matrix (double precision) }; Double_t Drand(Double_t &ix); #endif