#ifndef ROOT_TMatrixUtils #define ROOT_TMatrixUtils //+SEQ,CopyRight,T=NOINCLUDE. ////////////////////////////////////////////////////////////////////////// // // // Matrix utility classes. // // // // This file defines utility classes for the Linear Algebra Package. // // The following classes are defined here: // // TElementAction // // TElementPosAction // // TLazyMatrix // // THaarMatrix // // TMatrixRow // // TMatrixColumn // // TMatrixDiag // // TMatrixPivoting // // // ////////////////////////////////////////////////////////////////////////// #ifndef ROOT_TMatrix //*KEEP,TMatrix,T=C++. #include "TMatrix.h" //*KEND. #endif ////////////////////////////////////////////////////////////////////////// // // // TElementAction // // // // 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 TElementAction { friend class TMatrix; friend class TVector; private: virtual void Operation(Real_t &element) = 0; void operator=(const TElementAction &) { } }; ////////////////////////////////////////////////////////////////////////// // // // TElementPosAction // // // // 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 TElementPosAction { friend class TMatrix; friend class TVector; protected: Int_t fI; // i position of element being passed to Operation() Int_t fJ; // j position of element being passed to Operation() private: virtual void Operation(Real_t &element) = 0; void operator=(const TElementPosAction &) { } }; ////////////////////////////////////////////////////////////////////////// // // // TLazyMatrix // // // // Class used to make a lazy copy of a matrix, i.e. only copy matrix // // when really needed (when accessed). // // // ////////////////////////////////////////////////////////////////////////// class TLazyMatrix : public TObject { friend class TMatrix; protected: Int_t fRowUpb; Int_t fRowLwb; Int_t fColUpb; Int_t fColLwb; private: virtual void FillIn(TMatrix &m) const = 0; TLazyMatrix(const TLazyMatrix &) { } void operator=(const TLazyMatrix &) { } public: TLazyMatrix() { fRowUpb = fRowLwb = fColUpb = fColLwb = 0; } TLazyMatrix(Int_t nrows, Int_t ncols) : fRowUpb(nrows-1), fRowLwb(0), fColUpb(ncols-1), fColLwb(0) { } TLazyMatrix(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb) : fRowUpb(row_upb), fRowLwb(row_lwb), fColUpb(col_upb), fColLwb(col_lwb) { } ClassDef(TLazyMatrix,1) // Lazy matrix }; class THaarMatrix : public TLazyMatrix { private: void FillIn(TMatrix &m) const; public: THaarMatrix(Int_t n, Int_t no_cols = 0); }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixRow // // // // Class represents a row of a TMatrix. // // // ////////////////////////////////////////////////////////////////////////// class TMatrixRow : public TObject { friend class TMatrix; friend class TVector; private: const TMatrix *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] Real_t *fPtr; // pointer to the a[row,0] TMatrixRow() { fMatrix = 0; fInc = 0; fPtr = 0; } public: TMatrixRow(const TMatrix &matrix, Int_t row); void operator=(Real_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TVector &vec); const Real_t &operator()(Int_t i) const; Real_t &operator()(Int_t i); ClassDef(TMatrixRow,1) // One row of a matrix }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixColumn // // // // Class represents a column of a TMatrix. // // // ////////////////////////////////////////////////////////////////////////// class TMatrixColumn : public TObject { friend class TMatrix; friend class TVector; private: const TMatrix *fMatrix; // the matrix I am a column of Int_t fColInd; // effective column index Real_t *fPtr; // pointer to the a[0,i] column TMatrixColumn() { fMatrix = 0; fPtr = 0; } public: TMatrixColumn(const TMatrix &matrix, Int_t col); void operator=(Real_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TVector &vec); const Real_t &operator()(Int_t i) const; Real_t &operator()(Int_t i); ClassDef(TMatrixColumn,1) // One column of a matrix }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixDiag // // // // Class represents the diagonal of a matrix (for easy manipulation). // // // ////////////////////////////////////////////////////////////////////////// class TMatrixDiag : public TObject { friend class TMatrix; friend class TVector; private: const TMatrix *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) Real_t *fPtr; // pointer to the a[0,0] TMatrixDiag() { fMatrix = 0; fInc = 0; fNdiag = 0; fPtr = 0; } public: TMatrixDiag(const TMatrix &matrix); void operator=(Real_t val); void operator+=(Double_t val); void operator*=(Double_t val); void operator=(const TVector &vec); const Real_t &operator()(Int_t i) const; Real_t &operator()(Int_t i); Int_t GetNdiags() const { return fNdiag; } ClassDef(TMatrixDiag,1) // Diagonal of a matrix }; ////////////////////////////////////////////////////////////////////////// // // // TMatrixPivoting // // // // This class inherits from TMatrix and it keeps additional information // // about what is being/has been pivoted. // // // ////////////////////////////////////////////////////////////////////////// class TMatrixPivoting : public TMatrix { private: typedef Real_t* Index_t; // wanted to have typeof(index[0]) Index_t *const fRowIndex; // fRowIndex[i] = ptr to the i-th // matrix row, or 0 if the row // has been pivoted. Note, // pivoted columns are marked // by setting fIndex[j] to zero. // Information about the pivot that was // just picked up Double_t fPivotValue; // Value of the pivoting element Index_t fPivotRow; // pivot's location (ptrs) Index_t fPivotCol; Int_t fPivotOdd; // parity of the pivot // (0 for even, 1 for odd) void PickUpPivot(); // Pick up a pivot from // not-pivoted rows and cols public: TMatrixPivoting(const TMatrix &m); ~TMatrixPivoting(); Double_t PivotingAndElimination(); // Perform the pivoting, return // the pivot value times (-1)^(pi+pj) // (pi,pj - pivot el row & col) }; //----- inlines ---------------------------------------------------------------- #ifndef R__HPUX #ifndef __CINT__ inline TMatrixRow::TMatrixRow(const TMatrix &matrix, Int_t row) : fMatrix(&matrix), fInc(matrix.fNrows) { if (!matrix.IsValid()) { Error("TMatrixRow", "matrix is not initialized"); return; } fRowInd = row - matrix.fRowLwb; if (fRowInd >= matrix.fNrows || fRowInd < 0) { Error("TMatrixRow", "row #%d is not within the matrix", row); return; } fPtr = &(matrix.fIndex[0][fRowInd]); } inline const Real_t &TMatrixRow::operator()(Int_t i) const { // Get hold of the i-th row's element. static Real_t err; err = 0.0; if (!fMatrix->IsValid()) { Error("operator()", "matrix is not initialized"); return err; } Int_t acoln = i-fMatrix->fColLwb; // Effective index if (acoln >= fMatrix->fNcols || acoln < 0) { Error("operator()", "TMatrixRow index %d is out of row boundaries [%d,%d]", i, fMatrix->fColLwb, fMatrix->fNcols+fMatrix->fColLwb-1); return err; } return fMatrix->fIndex[acoln][fPtr-fMatrix->fElements]; } inline Real_t &TMatrixRow::operator()(Int_t i) { return (Real_t&)((*(const TMatrixRow *)this)(i)); } inline TMatrixColumn::TMatrixColumn(const TMatrix &matrix, Int_t col) : fMatrix(&matrix) { if (!matrix.IsValid()) { Error("TMatrixColumn", "matrix is not initialized"); return; } fColInd = col - matrix.fColLwb; if (fColInd >= matrix.fNcols || fColInd < 0) { Error("TMatrixColumn", "column #%d is not within the matrix", col); return; } fPtr = &(matrix.fIndex[fColInd][0]); } inline const Real_t &TMatrixColumn::operator()(Int_t i) const { // Access the i-th element of the column static Real_t err; err = 0.0; if (!fMatrix->IsValid()) { Error("operator()", "matrix is not initialized"); return err; } Int_t arown = i-fMatrix->fRowLwb; // Effective indices if (arown >= fMatrix->fNrows || arown < 0) { Error("operator()", "TMatrixColumn index %d is out of column boundaries [%d,%d]", i, fMatrix->fRowLwb, fMatrix->fNrows+fMatrix->fRowLwb-1); return err; } return fPtr[arown]; } inline Real_t &TMatrixColumn::operator()(Int_t i) { return (Real_t&)((*(const TMatrixColumn *)this)(i)); } inline TMatrixDiag::TMatrixDiag(const TMatrix &matrix) : fMatrix(&matrix), fInc(matrix.fNrows+1), fNdiag(TMath::Min(matrix.fNrows, matrix.fNcols)) { if (!matrix.IsValid()) { Error("TMatrixDiag", "matrix is not initialized"); return; } fPtr = &(matrix.fElements[0]); } inline const Real_t &TMatrixDiag::operator()(Int_t i) const { // Get hold of the i-th diag element (indexing always starts at 0, // regardless of matrix' col_lwb and row_lwb) static Real_t err; err = 0.0; if (!fMatrix->IsValid()) { Error("operator()", "matrix is not initialized"); return err; } if (i > fNdiag || i < 1) { Error("TMatrixDiag", "TMatrixDiag index %d is out of diag boundaries [1,%d]", i, fNdiag); return err; } return fMatrix->fIndex[i-1][i-1]; } inline Real_t &TMatrixDiag::operator()(Int_t i) { return (Real_t&)((*(const TMatrixDiag *)this)(i)); } #endif #endif #endif