From a0363eb8c5440334c964b179266ac4e7f1f918f3 Mon Sep 17 00:00:00 2001 From: Giulio Eulisse Date: Fri, 25 Apr 2014 21:17:20 +0200 Subject: [PATCH] Drop spurious Math directory. --- Math/MatrixRepresentationsStatic.h | 492 ---------------- Math/SMatrix.h | 748 ------------------------ Math/SMatrix.icc | 907 ----------------------------- 3 files changed, 2147 deletions(-) delete mode 100644 Math/MatrixRepresentationsStatic.h delete mode 100644 Math/SMatrix.h delete mode 100644 Math/SMatrix.icc diff --git a/Math/MatrixRepresentationsStatic.h b/Math/MatrixRepresentationsStatic.h deleted file mode 100644 index b71872f1abcbf..0000000000000 --- a/Math/MatrixRepresentationsStatic.h +++ /dev/null @@ -1,492 +0,0 @@ -// @(#)root/smatrix:$Id: MatrixRepresentationsStatic.h 34964 2010-08-24 13:58:51Z moneta $ -// Author: L. Moneta, J. Palacios 2006 - -#ifndef ROOT_Math_MatrixRepresentationsStatic -#define ROOT_Math_MatrixRepresentationsStatic 1 - -#warning "using MatrixRepresentationsStatic by vin" - - -//#ifndef SMATRIX_USE_COMPUTATION -//#define SMATRIX_USE_CONSTEXPR -//#endif - - -#ifndef SMATRIX_USE_CONSTEXPR -#define SMATRIX_USE_COMPUTATION -#endif - - -// Include files - -/** - @defgroup MatRep SMatrix Storage Representation - @ingroup SMatrixGroup - - @author Juan Palacios - @date 2006-01-15 - - Classes MatRepStd and MatRepSym for generic and symmetric matrix - data storage and manipulation. Define data storage and access, plus - operators =, +=, -=, ==. - - */ - -#ifndef ROOT_Math_StaticCheck -#include "Math/StaticCheck.h" -#endif - -#ifdef SMATRIX_USE_CONSTEXPR -#warning "using MatrixRepresentationsStatic constexpr" - -#include -#include -#include -#include -#endif - - -namespace ROOT { - -namespace Math { - - //________________________________________________________________________________ - /** - MatRepStd - Standard Matrix representation for a general D1 x D2 matrix. - This class is itself a template on the contained type T, the number of rows and the number of columns. - Its data member is an array T[nrows*ncols] containing the matrix data. - The data are stored in the row-major C convention. - For example, for a matrix, M, of size 3x3, the data \f$ \left[a_0,a_1,a_2,.......,a_7,a_8 \right] \f$d are stored in the following order: - \f[ - M = \left( \begin{array}{ccc} - a_0 & a_1 & a_2 \\ - a_3 & a_4 & a_5 \\ - a_6 & a_7 & a_8 \end{array} \right) - \f] - - @ingroup MatRep - */ - - - template - class MatRepStd { - - public: - - typedef T value_type; - - inline const T& operator()(unsigned int i, unsigned int j) const { - return fArray[i*D2+j]; - } - inline T& operator()(unsigned int i, unsigned int j) { - return fArray[i*D2+j]; - } - inline T& operator[](unsigned int i) { return fArray[i]; } - - inline const T& operator[](unsigned int i) const { return fArray[i]; } - - inline T apply(unsigned int i) const { return fArray[i]; } - - inline T* Array() { return fArray; } - - inline const T* Array() const { return fArray; } - - template - inline MatRepStd& operator+=(const R& rhs) { - for(unsigned int i=0; i - inline MatRepStd& operator-=(const R& rhs) { - for(unsigned int i=0; i - inline MatRepStd& operator=(const R& rhs) { - for(unsigned int i=0; i - inline bool operator==(const R& rhs) const { - bool rc = true; - for(unsigned int i=0; i -// struct Creator { -// static const RowOffsets & Offsets() { -// static RowOffsets off; -// return off; -// } - - /** - Static structure to keep the conversion from (i,j) to offsets in the storage data for a - symmetric matrix - */ - - - - template - struct RowOffsets { - inline RowOffsets() { - int v[D]; - v[0]=0; - for (unsigned int i=1; i struct indices{}; - - template - struct make_indices_impl; - - template - struct make_indices_impl, N> - { - typedef typename make_indices_impl, - N>::type type; - }; - - template - struct make_indices_impl, N> { - typedef indices type; - }; - - template - struct make_indices : make_indices_impl<0, indices<>, N> {}; - // end of stuff - - - - template - constexpr std::array()(std::declval())), sizeof...(I)> - do_make(F f, indices) - { - return std::array()(std::declval())), - sizeof...(I)>{{ f(I0 + I)... }}; - } - - template - constexpr std::array()(std::declval())), N> - make(F f) { - return do_make(f, typename make_indices::type()); - } - - } // namespace rowOffsetsUtils - -#elif defined(SMATRIX_USE_COMPUTATION) - -#else - - - -// Make the lookup tables available at compile time: -// Add them to a namespace? -static const int fOff1x1[] = {0}; -static const int fOff2x2[] = {0, 1, 1, 2}; -static const int fOff3x3[] = {0, 1, 3, 1, 2, 4, 3, 4, 5}; -static const int fOff4x4[] = {0, 1, 3, 6, 1, 2, 4, 7, 3, 4, 5, 8, 6, 7, 8, 9}; -static const int fOff5x5[] = {0, 1, 3, 6, 10, 1, 2, 4, 7, 11, 3, 4, 5, 8, 12, 6, 7, 8, 9, 13, 10, 11, 12, 13, 14}; -static const int fOff6x6[] = {0, 1, 3, 6, 10, 15, 1, 2, 4, 7, 11, 16, 3, 4, 5, 8, 12, 17, 6, 7, 8, 9, 13, 18, 10, 11, 12, 13, 14, 19, 15, 16, 17, 18, 19, 20}; - -static const int fOff7x7[] = {0, 1, 3, 6, 10, 15, 21, 1, 2, 4, 7, 11, 16, 22, 3, 4, 5, 8, 12, 17, 23, 6, 7, 8, 9, 13, 18, 24, 10, 11, 12, 13, 14, 19, 25, 15, 16, 17, 18, 19, 20, 26, 21, 22, 23, 24, 25, 26, 27}; - -static const int fOff8x8[] = {0, 1, 3, 6, 10, 15, 21, 28, 1, 2, 4, 7, 11, 16, 22, 29, 3, 4, 5, 8, 12, 17, 23, 30, 6, 7, 8, 9, 13, 18, 24, 31, 10, 11, 12, 13, 14, 19, 25, 32, 15, 16, 17, 18, 19, 20, 26, 33, 21, 22, 23, 24, 25, 26, 27, 34, 28, 29, 30, 31, 32, 33, 34, 35}; - -static const int fOff9x9[] = {0, 1, 3, 6, 10, 15, 21, 28, 36, 1, 2, 4, 7, 11, 16, 22, 29, 37, 3, 4, 5, 8, 12, 17, 23, 30, 38, 6, 7, 8, 9, 13, 18, 24, 31, 39, 10, 11, 12, 13, 14, 19, 25, 32, 40, 15, 16, 17, 18, 19, 20, 26, 33, 41, 21, 22, 23, 24, 25, 26, 27, 34, 42, 28, 29, 30, 31, 32, 33, 34, 35, 43, 36, 37, 38, 39, 40, 41, 42, 43, 44}; - -static const int fOff10x10[] = {0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 1, 2, 4, 7, 11, 16, 22, 29, 37, 46, 3, 4, 5, 8, 12, 17, 23, 30, 38, 47, 6, 7, 8, 9, 13, 18, 24, 31, 39, 48, 10, 11, 12, 13, 14, 19, 25, 32, 40, 49, 15, 16, 17, 18, 19, 20, 26, 33, 41, 50, 21, 22, 23, 24, 25, 26, 27, 34, 42, 51, 28, 29, 30, 31, 32, 33, 34, 35, 43, 52, 36, 37, 38, 39, 40, 41, 42, 43, 44, 53, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54}; - -template<> - struct RowOffsets<1> { - RowOffsets() {} - int operator()(unsigned int , unsigned int ) const { return 0; } // Just one element - int apply(unsigned int ) const { return 0; } - }; - -template<> - struct RowOffsets<2> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return i+j; /*fOff2x2[i*2+j];*/ } - int apply(unsigned int i) const { return fOff2x2[i]; } - }; - -template<> - struct RowOffsets<3> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff3x3[i*3+j]; } - int apply(unsigned int i) const { return fOff3x3[i]; } - }; - -template<> - struct RowOffsets<4> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff4x4[i*4+j]; } - int apply(unsigned int i) const { return fOff4x4[i]; } - }; - - template<> - struct RowOffsets<5> { - inline RowOffsets() {} - inline int operator()(unsigned int i, unsigned int j) const { return fOff5x5[i*5+j]; } -// int operator()(unsigned int i, unsigned int j) const { -// if(j <= i) return (i * (i + 1)) / 2 + j; -// else return (j * (j + 1)) / 2 + i; -// } - inline int apply(unsigned int i) const { return fOff5x5[i]; } - }; - -template<> - struct RowOffsets<6> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff6x6[i*6+j]; } - int apply(unsigned int i) const { return fOff6x6[i]; } - }; - -template<> - struct RowOffsets<7> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff7x7[i*7+j]; } - int apply(unsigned int i) const { return fOff7x7[i]; } - }; - -template<> - struct RowOffsets<8> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff8x8[i*8+j]; } - int apply(unsigned int i) const { return fOff8x8[i]; } - }; - -template<> - struct RowOffsets<9> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff9x9[i*9+j]; } - int apply(unsigned int i) const { return fOff9x9[i]; } - }; - -template<> - struct RowOffsets<10> { - RowOffsets() {} - int operator()(unsigned int i, unsigned int j) const { return fOff10x10[i*10+j]; } - int apply(unsigned int i) const { return fOff10x10[i]; } - }; - -#endif - -//_________________________________________________________________________________ - /** - MatRepSym - Matrix storage representation for a symmetric matrix of dimension NxN - This class is a template on the contained type and on the symmetric matrix size, N. - It has as data member an array of type T of size N*(N+1)/2, - containing the lower diagonal block of the matrix. - The order follows the lower diagonal block, still in a row-major convention. - For example for a symmetric 3x3 matrix the order of the 6 elements - \f$ \left[a_0,a_1.....a_5 \right]\f$ is: - \f[ - M = \left( \begin{array}{ccc} - a_0 & a_1 & a_3 \\ - a_1 & a_2 & a_4 \\ - a_3 & a_4 & a_5 \end{array} \right) - \f] - - @ingroup MatRep - */ - template - class MatRepSym { - - public: - -#if defined(SMATRIX_USE_CONSTEXPR) || defined(SMATRIX_USE_COMPUTATION) - /* constexpr */ MatRepSym(){} -#else - MatRepSym() :fOff(0) { CreateOffsets(); } -#endif - typedef T value_type; - - - inline T & operator()(unsigned int i, unsigned int j) - { return fArray[offset(i, j)]; } - - inline /* constexpr */ T const & operator()(unsigned int i, unsigned int j) const - { return fArray[offset(i, j)]; } - - inline T& operator[](unsigned int i) { - return fArray[off(i)]; - } - - inline /* constexpr */ T const & operator[](unsigned int i) const { - return fArray[off(i)]; - } - - inline /* constexpr */ T apply(unsigned int i) const { - return fArray[off(i)]; - } - - inline T* Array() { return fArray; } - - inline const T* Array() const { return fArray; } - - /** - assignment : only symmetric to symmetric allowed - */ - template - inline MatRepSym& operator=(const R&) { - STATIC_CHECK(0==1, - Cannot_assign_general_to_symmetric_matrix_representation); - return *this; - } - inline MatRepSym& operator=(const MatRepSym& rhs) { - for(unsigned int i=0; i - inline MatRepSym& operator+=(const R&) { - STATIC_CHECK(0==1, - Cannot_add_general_to_symmetric_matrix_representation); - return *this; - } - inline MatRepSym& operator+=(const MatRepSym& rhs) { - for(unsigned int i=0; i - inline MatRepSym& operator-=(const R&) { - STATIC_CHECK(0==1, - Cannot_substract_general_to_symmetric_matrix_representation); - return *this; - } - inline MatRepSym& operator-=(const MatRepSym& rhs) { - for(unsigned int i=0; i - inline bool operator==(const R& rhs) const { - bool rc = true; - for(unsigned int i=0; i(off1); - return v[i]; - } - - static inline constexpr unsigned int - offset(unsigned int i, unsigned int j) - { - //if (j > i) std::swap(i, j); - return off(i*D+j); - // return (i>j) ? (i * (i+1) / 2) + j : (j * (j+1) / 2) + i; - } -#elif defined(SMATRIX_USE_COMPUTATION) - static inline unsigned int - offset(unsigned int i, unsigned int j) - { - return (i>j) ? (i * (i+1) / 2) + j : (j * (j+1) / 2) + i; - } - static inline unsigned int - off(unsigned int i) { - return offset(i/D, i%D); - } -#else - unsigned int - offset(unsigned int i, unsigned int j) const - { - return Offsets()(i,j); - } - unsigned int - off(unsigned int i) const { - return Offsets().apply(i); - } - - - void CreateOffsets() { - const static RowOffsets off; - fOff = &off; - } - - inline const RowOffsets & Offsets() const { - return *fOff; - } -#endif - - - private: - //T __attribute__ ((aligned (16))) fArray[kSize]; - T fArray[kSize]; - -#if !defined(SMATRIX_USE_CONSTEXPR) && !defined(SMATRIX_USE_COMPUTATION) - const RowOffsets * fOff; //! transient -#endif - }; - - - -} // namespace Math -} // namespace ROOT - - -#endif // MATH_MATRIXREPRESENTATIONSSTATIC_H diff --git a/Math/SMatrix.h b/Math/SMatrix.h deleted file mode 100644 index c23f9db6dc6c7..0000000000000 --- a/Math/SMatrix.h +++ /dev/null @@ -1,748 +0,0 @@ -// @(#)root/smatrix:$Id$ -// Author: T. Glebe, L. Moneta, J. Palacios 2005 - -#ifndef ROOT_Math_SMatrix -#define ROOT_Math_SMatrix - -/********************************************************************************* -// -// source: -// -// type: source code -// -// created: 20. Mar 2001 -// -// author: Thorsten Glebe -// HERA-B Collaboration -// Max-Planck-Institut fuer Kernphysik -// Saupfercheckweg 1 -// 69117 Heidelberg -// Germany -// E-mail: T.Glebe@mpi-hd.mpg.de -// -// Description: A fixed size two dimensional Matrix class -// -// changes: -// 20 Mar 2001 (TG) creation -// 21 Mar 2001 (TG) added operators +=, -=, *=, /= -// 26 Mar 2001 (TG) place_in_row(), place_in_col() added -// 02 Apr 2001 (TG) non-const Array() added -// 03 Apr 2001 (TG) invert() added -// 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added -// 09 Apr 2001 (TG) CTOR from array added -// 11 Apr 2001 (TG) rows(), cols(), size() replaced by rows, cols, size -// 25 Mai 2001 (TG) row(), col() added -// 04 Sep 2001 (TG) moved inlined functions to .icc file -// 11 Jan 2002 (TG) added operator==(), operator!=() -// 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<() -// -***************************************************************************/ -// for platform specific configurations - -#ifndef ROOT_Math_MnConfig -#include "Math/MConfig.h" -#endif - -#include - - - - -//doxygen tag - -/** - @defgroup SMatrixGroup SMatrix - - \ref SMatrix for high performance vector and matrix computations. - Classes representing Matrices and Vectors of arbitrary type and dimension and related functions. - For a detailed description and usage examples see: -
    -
  • \ref SMatrix home page -
  • \ref SVectorDoc -
  • \ref SMatrixDoc -
  • \ref MatVecFunctions -
- -*/ - -/** - @defgroup SMatrixSVector Matrix and Vector classes - - @ingroup SMatrixGroup - - Classes representing Matrices and Vectors of arbitrary type and dimension. - For a detailed description and usage examples see: -
    -
  • \ref SVectorDoc -
  • \ref SMatrixDoc -
  • \ref MatVecFunctions -
- -*/ - - -#ifndef ROOT_Math_Expression -#include "Math/Expression.h" -#endif -#ifndef ROOT_Math_MatrixRepresentationsStatic -#include "Math/MatrixRepresentationsStatic.h" -#endif - - -namespace ROOT { - -namespace Math { - - -template class SVector; - -struct SMatrixIdentity { }; -struct SMatrixNoInit { }; - -//__________________________________________________________________________ -/** - SMatrix: a generic fixed size D1 x D2 Matrix class. - The class is template on the scalar type, on the matrix sizes: - D1 = number of rows and D2 = number of columns - amd on the representation storage type. - By default the representation is MatRepStd (standard D1xD2 of type T), - but it can be of type MatRepSym for symmetric matrices DxD, where the storage is only - D*(D+1)/2. - - See \ref SMatrixDoc. - - Original author is Thorsten Glebe - HERA-B Collaboration, MPI Heidelberg (Germany) - - @ingroup SMatrixSVector - - @authors T. Glebe, L. Moneta and J. Palacios -*/ -//============================================================================== -// SMatrix: column-wise storage -//============================================================================== -template > -class SMatrix { -public: - /** @name --- Typedefs --- */ - - /** contained scalar type */ - typedef T value_type; - - /** storage representation type */ - typedef R rep_type; - - /** STL iterator interface. */ - typedef T* iterator; - - /** STL const_iterator interface. */ - typedef const T* const_iterator; - - - - /** @name --- Constructors and Assignment --- */ - - /** - Default constructor: - */ - SMatrix(); - /// - /** - construct from an without initialization - */ - SMatrix( SMatrixNoInit ){} - - /** - construct from an identity matrix - */ - SMatrix( SMatrixIdentity ); - /** - copy constructor (from a matrix of the same representation - */ - SMatrix(const SMatrix& rhs); - /** - construct from a matrix with different representation. - Works only from symmetric to general and not viceversa. - */ - template - SMatrix(const SMatrix& rhs); - - /** - Construct from an expression. - In case of symmetric matrices does not work if expression is of type general - matrices. In case one needs to force the assignment from general to symmetric, one can use the - ROOT::Math::AssignSym::Evaluate function. - */ - template - SMatrix(const Expr& rhs); - - - /** - Constructor with STL iterator interface. The data will be copied into the matrix - \param begin start iterator position - \param end end iterator position - \param triang if true only the triangular lower/upper part of the matrix is filled from the iterators - \param lower if true the lower triangular part is filled - - Size of the matrix must match size of the iterators, if triang is false, otherwise the size of the - triangular block. In the case of symmetric matrices triang is considered always to be true - (what-ever the user specifies) and the size of the iterators must be equal to the size of the - triangular block, which is the number of independent elements of a symmetric matrix: N*(N+1)/2 - - */ - template - SMatrix(InputIterator begin, InputIterator end, bool triang = false, bool lower = true); - - /** - Constructor with STL iterator interface. The data will be copied into the matrix - \param begin start iterator position - \param size iterator size - \param triang if true only the triangular lower/upper part of the matrix is filled from the iterators - \param lower if true the lower triangular part is filled - - Size of the iterators must not be larger than the size of the matrix representation. - In the case of symmetric matrices the size is N*(N+1)/2. - - */ - template - SMatrix(InputIterator begin, unsigned int size, bool triang = false, bool lower = true); - - /** - constructor of a symmetrix a matrix from a SVector containing the lower (upper) - triangular part. - */ -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION - SMatrix(const SVector & v, bool lower = true ); -#else - template - SMatrix(const SVector & v, bool lower = true ); -#endif - - - /** - Construct from a scalar value (only for size 1 matrices) - */ - explicit SMatrix(const T& rhs); - - /** - Assign from another compatible matrix. - Possible Symmetirc to general but NOT vice-versa - */ - template - SMatrix& operator=(const M& rhs); - - /** - Assign from a matrix expression - */ - template - SMatrix& operator=(const Expr& rhs); - - /** - Assign from an identity matrix - */ - SMatrix & operator=(SMatrixIdentity ); - - /** - Assign from a scalar value (only for size 1 matrices) - */ - SMatrix& operator=(const T& rhs); - - /** @name --- Matrix dimension --- */ - - /** - Enumeration defining the matrix dimension, - number of rows, columns and size = rows*columns) - */ - enum { - /// return no. of matrix rows - kRows = D1, - /// return no. of matrix columns - kCols = D2, - /// return no of elements: rows*columns - kSize = D1*D2 - }; - - /** @name --- Access functions --- */ - - /** access the parse tree with the index starting from zero and - following the C convention for the order in accessing - the matrix elements. - Same convention for general and symmetric matrices. - */ - T apply(unsigned int i) const; - - /// return read-only pointer to internal array - const T* Array() const; - /// return pointer to internal array - T* Array(); - - /** @name --- STL-like interface --- - The iterators access the matrix element in the order how they are - stored in memory. The C (row-major) convention is used, and in the - case of symmetric matrices the iterator spans only the lower diagonal - block. For example for a symmetric 3x3 matrices the order of the 6 - elements \f${a_0,...a_5}\f$ is: - \f[ - M = \left( \begin{array}{ccc} - a_0 & a_1 & a_3 \\ - a_1 & a_2 & a_4 \\ - a_3 & a_4 & a_5 \end{array} \right) - \f] - */ - - /** STL iterator interface. */ - iterator begin(); - - /** STL iterator interface. */ - iterator end(); - - /** STL const_iterator interface. */ - const_iterator begin() const; - - /** STL const_iterator interface. */ - const_iterator end() const; - - /** - Set matrix elements with STL iterator interface. The data will be copied into the matrix - \param begin start iterator position - \param end end iterator position - \param triang if true only the triangular lower/upper part of the matrix is filled from the iterators - \param lower if true the lower triangular part is filled - - Size of the matrix must match size of the iterators, if triang is false, otherwise the size of the - triangular block. In the case of symmetric matrices triang is considered always to be true - (what-ever the user specifies) and the size of the iterators must be equal to the size of the - triangular block, which is the number of independent elements of a symmetric matrix: N*(N+1)/2 - - */ - template - void SetElements(InputIterator begin, InputIterator end, bool triang = false, bool lower = true); - - /** - Constructor with STL iterator interface. The data will be copied into the matrix - \param begin start iterator position - \param size iterator size - \param triang if true only the triangular lower/upper part of the matrix is filled from the iterators - \param lower if true the lower triangular part is filled - - Size of the iterators must not be larger than the size of the matrix representation. - In the case of symmetric matrices the size is N*(N+1)/2. - - */ - template - void SetElements(InputIterator begin, unsigned int size, bool triang = false, bool lower = true); - - - /** @name --- Operators --- */ - /// element wise comparison - bool operator==(const T& rhs) const; - /// element wise comparison - bool operator!=(const T& rhs) const; - /// element wise comparison - template - bool operator==(const SMatrix& rhs) const; - /// element wise comparison - bool operator!=(const SMatrix& rhs) const; - /// element wise comparison - template - bool operator==(const Expr& rhs) const; - /// element wise comparison - template - bool operator!=(const Expr& rhs) const; - - /// element wise comparison - bool operator>(const T& rhs) const; - /// element wise comparison - bool operator<(const T& rhs) const; - /// element wise comparison - template - bool operator>(const SMatrix& rhs) const; - /// element wise comparison - template - bool operator<(const SMatrix& rhs) const; - /// element wise comparison - template - bool operator>(const Expr& rhs) const; - /// element wise comparison - template - bool operator<(const Expr& rhs) const; - - /** - read only access to matrix element, with indices starting from 0 - */ - const T& operator()(unsigned int i, unsigned int j) const; - /** - read/write access to matrix element with indices starting from 0 - */ - T& operator()(unsigned int i, unsigned int j); - - /** - read only access to matrix element, with indices starting from 0. - Function will check index values and it will assert if they are wrong - */ - const T& At(unsigned int i, unsigned int j) const; - /** - read/write access to matrix element with indices starting from 0. - Function will check index values and it will assert if they are wrong - */ - T& At(unsigned int i, unsigned int j); - - - // helper class for implementing the m[i][j] operator - - class SMatrixRow { - public: - SMatrixRow ( SMatrix & rhs, unsigned int i ) : - fMat(&rhs), fRow(i) - {} - T & operator[](int j) { return (*fMat)(fRow,j); } - private: - SMatrix * fMat; - unsigned int fRow; - }; - - class SMatrixRow_const { - public: - SMatrixRow_const ( const SMatrix & rhs, unsigned int i ) : - fMat(&rhs), fRow(i) - {} - - const T & operator[](int j) const { return (*fMat)(fRow, j); } - - private: - const SMatrix * fMat; - unsigned int fRow; - }; - - /** - read only access to matrix element, with indices starting from 0 : m[i][j] - */ - SMatrixRow_const operator[](unsigned int i) const { return SMatrixRow_const(*this, i); } - /** - read/write access to matrix element with indices starting from 0 : m[i][j] - */ - SMatrixRow operator[](unsigned int i) { return SMatrixRow(*this, i); } - - - /** - addition with a scalar - */ - SMatrix&operator+=(const T& rhs); - - /** - addition with another matrix of any compatible representation - */ - template - SMatrix&operator+=(const SMatrix& rhs); - - /** - addition with a compatible matrix expression - */ - template - SMatrix& operator+=(const Expr& rhs); - - /** - subtraction with a scalar - */ - SMatrix& operator-=(const T& rhs); - - /** - subtraction with another matrix of any compatible representation - */ - template - SMatrix&operator-=(const SMatrix& rhs); - - /** - subtraction with a compatible matrix expression - */ - template - SMatrix& operator-=(const Expr& rhs); - - /** - multiplication with a scalar - */ - SMatrix& operator*=(const T& rhs); - -#ifndef __CINT__ - - - /** - multiplication with another compatible matrix (it is a real matrix multiplication) - Note that this operation does not avid to create a temporary to store intermidiate result - */ - template - SMatrix& operator*=(const SMatrix& rhs); - - /** - multiplication with a compatible matrix expression (it is a real matrix multiplication) - */ - template - SMatrix& operator*=(const Expr& rhs); - -#endif - - /** - division with a scalar - */ - SMatrix& operator/=(const T& rhs); - - - - /** @name --- Linear Algebra Functions --- */ - - /** - Invert a square Matrix ( this method changes the current matrix). - Return true if inversion is successfull. - The method used for general square matrices is the LU factorization taken from Dinv routine - from the CERNLIB (written in C++ from CLHEP authors) - In case of symmetric matrices Bunch-Kaufman diagonal pivoting method is used - (The implementation is the one written by the CLHEP authors) - */ - bool Invert(); - - /** - Invert a square Matrix and returns a new matrix. In case the inversion fails - the current matrix is returned. - \param ifail . ifail will be set to 0 when inversion is successfull. - See ROOT::Math::SMatrix::Invert for the inversion algorithm - */ - SMatrix Inverse(int & ifail ) const; - - /** - Fast Invertion of a square Matrix ( this method changes the current matrix). - Return true if inversion is successfull. - The method used is based on direct inversion using the Cramer rule for - matrices upto 5x5. Afterwards the same defult algorithm of Invert() is used. - Note that this method is faster but can suffer from much larger numerical accuracy - when the condition of the matrix is large - */ - bool InvertFast(); - - /** - Invert a square Matrix and returns a new matrix. In case the inversion fails - the current matrix is returned. - \param ifail . ifail will be set to 0 when inversion is successfull. - See ROOT::Math::SMatrix::InvertFast for the inversion algorithm - */ - SMatrix InverseFast(int & ifail ) const; - - /** - Invertion of a symmetric positive defined Matrix using Choleski decomposition. - ( this method changes the current matrix). - Return true if inversion is successfull. - The method used is based on Choleski decomposition - A compile error is given if the matrix is not of type symmetric and a run-time failure if the - matrix is not positive defined. - For solving a linear system, it is possible to use also the function - ROOT::Math::SolveChol(matrix, vector) which will be faster than performing the inversion - */ - bool InvertChol(); - - /** - Invert of a symmetric positive defined Matrix using Choleski decomposition. - A compile error is given if the matrix is not of type symmetric and a run-time failure if the - matrix is not positive defined. - In case the inversion fails the current matrix is returned. - \param ifail . ifail will be set to 0 when inversion is successfull. - See ROOT::Math::SMatrix::InvertChol for the inversion algorithm - */ - SMatrix InverseChol(int & ifail ) const; - - /** - determinant of square Matrix via Dfact. - Return true when the calculation is successfull. - \param det will contain the calculated determinant value - \b Note: this will destroy the contents of the Matrix! - */ - bool Det(T& det); - - /** - determinant of square Matrix via Dfact. - Return true when the calculation is successfull. - \param det will contain the calculated determinant value - \b Note: this will preserve the content of the Matrix! - */ - bool Det2(T& det) const; - - - /** @name --- Matrix Slice Functions --- */ - - /// place a vector in a Matrix row - template - SMatrix& Place_in_row(const SVector& rhs, - unsigned int row, - unsigned int col); - /// place a vector expression in a Matrix row - template - SMatrix& Place_in_row(const VecExpr& rhs, - unsigned int row, - unsigned int col); - /// place a vector in a Matrix column - template - SMatrix& Place_in_col(const SVector& rhs, - unsigned int row, - unsigned int col); - /// place a vector expression in a Matrix column - template - SMatrix& Place_in_col(const VecExpr& rhs, - unsigned int row, - unsigned int col); - /// place a matrix in this matrix - template - SMatrix& Place_at(const SMatrix& rhs, - unsigned int row, - unsigned int col); - /// place a matrix expression in this matrix - template - SMatrix& Place_at(const Expr& rhs, - unsigned int row, - unsigned int col); - - /** - return a full Matrix row as a vector (copy the content in a new vector) - */ - SVector Row(unsigned int therow) const; - - /** - return a full Matrix column as a vector (copy the content in a new vector) - */ - SVector Col(unsigned int thecol) const; - - /** - return a slice of therow as a vector starting at the colum value col0 until col0+N, - where N is the size of the vector (SubVector::kSize ) - Condition col0+N <= D2 - */ - template - SubVector SubRow(unsigned int therow, unsigned int col0 = 0 ) const; - - /** - return a slice of the column as a vector starting at the row value row0 until row0+Dsub. - where N is the size of the vector (SubVector::kSize ) - Condition row0+N <= D1 - */ - template - SubVector SubCol(unsigned int thecol, unsigned int row0 = 0) const; - - /** - return a submatrix with the upper left corner at the values (row0, col0) and with sizes N1, N2 - where N1 and N2 are the dimension of the sub-matrix (SubMatrix::kRows and SubMatrix::kCols ) - Condition row0+N1 <= D1 && col0+N2 <=D2 - */ - template - SubMatrix Sub(unsigned int row0, unsigned int col0) const; - - /** - return diagonal elements of a matrix as a Vector. - It works only for squared matrices D1 == D2, otherwise it will produce a compile error - */ - SVector Diagonal() const; - - /** - Set the diagonal elements from a Vector - Require that vector implements ::kSize since a check (statically) is done on - diagonal size == vector size - */ - template - void SetDiagonal(const Vector & v); - - /** - return the trace of a matrix - Sum of the diagonal elements - */ - T Trace() const; - - - /** - return the upper Triangular block of the matrices (including the diagonal) as - a vector of sizes N = D1 * (D1 + 1)/2. - It works only for square matrices with D1==D2, otherwise it will produce a compile error - */ -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION - SVector UpperBlock() const; -#else - template - SubVector UpperBlock() const; -#endif - - /** - return the lower Triangular block of the matrices (including the diagonal) as - a vector of sizes N = D1 * (D1 + 1)/2. - It works only for square matrices with D1==D2, otherwise it will produce a compile error - */ -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION - SVector LowerBlock() const; -#else - template - SubVector LowerBlock() const; -#endif - - - /** @name --- Other Functions --- */ - - /** - Function to check if a matrix is sharing same memory location of the passed pointer - This function is used by the expression templates to avoid the alias problem during - expression evaluation. When the matrix is in use, for example in operations - like A = B * A, a temporary object storing the intermediate result is automatically - created when evaluating the expression. - - */ - bool IsInUse(const T* p) const; - - // submatrices - - /// Print: used by operator<<() - std::ostream& Print(std::ostream& os) const; - - - - -public: - - /** @name --- Data Member --- */ - - /** - Matrix Storage Object containing matrix data - */ - R fRep; - -}; // end of class SMatrix - - - - -//============================================================================== -// operator<< -//============================================================================== -template -inline std::ostream& operator<<(std::ostream& os, const ROOT::Math::SMatrix& rhs) { - return rhs.Print(os); -} - - - } // namespace Math - -} // namespace ROOT - - - - - - -#ifndef __CINT__ - -#ifndef ROOT_Math_SMatrix_icc -#include "Math/SMatrix.icc" -#endif - -#ifndef ROOT_Math_MatrixFunctions -#include "Math/MatrixFunctions.h" -#endif - -#endif //__CINT__ - -#endif /* ROOT_Math_SMatrix */ diff --git a/Math/SMatrix.icc b/Math/SMatrix.icc deleted file mode 100644 index 6f19dcb6d42b1..0000000000000 --- a/Math/SMatrix.icc +++ /dev/null @@ -1,907 +0,0 @@ -// @(#)root/smatrix:$Id$ -// Authors: T. Glebe, L. Moneta 2005 - -#ifndef ROOT_Math_SMatrix_icc -#define ROOT_Math_SMatrix_icc -// ******************************************************************** -// -// source: -// -// type: source code -// -// created: 21. Mar 2001 -// -// author: Thorsten Glebe -// HERA-B Collaboration -// Max-Planck-Institut fuer Kernphysik -// Saupfercheckweg 1 -// 69117 Heidelberg -// Germany -// E-mail: T.Glebe@mpi-hd.mpg.de -// -// Description: SMatrix implementation file -// -// changes: -// 21 Mar 2001 (TG) creation -// 26 Mar 2001 (TG) place_in_row(), place_in_col() added -// 03 Apr 2001 (TG) invert() added -// 07 Apr 2001 (TG) CTOR from SVertex (dyadic product) added -// 09 Apr 2001 (TG) CTOR from array added -// 25 Mai 2001 (TG) row(), col() added -// 11 Jul 2001 (TG) added #include Functions.hh -// 11 Jan 2002 (TG) added operator==(), operator!=() -// 14 Jan 2002 (TG) added more operator==(), operator!=(), operator>(), operator<() -// -// ******************************************************************** -#include -#include -#include -//#ifndef ROOT_Math_Dsinv -//#include "Math/Dsinv.h" -//#endif -//#include "Math/Dsinv_array.h" -//#include "Math/Dsfact.h" - -#ifndef ROOT_Math_Dfact -#include "Math/Dfact.h" -#endif -#ifndef ROOT_Math_Dinv -#include "Math/Dinv.h" -#endif -#ifndef ROOT_Math_Functions -#include "Math/Functions.h" -#endif -#ifndef ROOT_Math_HelperOps -#include "Math/HelperOps.h" -#endif -#ifndef ROOT_Math_StaticCheck -#include "Math/StaticCheck.h" -#endif - - - - - - - -namespace ROOT { - -namespace Math { - - - -//============================================================================== -// Constructors -//============================================================================== -template -SMatrix::SMatrix() { - // operator=(0); // if operator=(T ) is defined - for(unsigned int i=0; i -SMatrix::SMatrix( SMatrixIdentity ) { - for(unsigned int i=0; i -SMatrix::SMatrix(const SMatrix& rhs) { - fRep = rhs.fRep; -} - - -template -template -SMatrix::SMatrix(const SMatrix& rhs) { - operator=(rhs); -} - - -template -template -SMatrix::SMatrix(const Expr& rhs) { - operator=(rhs); -} - - -//============================================================================= -// New Constructors from STL interfaces -//============================================================================= -template -template -SMatrix::SMatrix(InputIterator ibegin, InputIterator iend, bool triang, bool lower) { - // assume iterator size == matrix size - for(unsigned int i=0; i::Evaluate(*this,ibegin,iend,triang,lower); -} - -template -template -SMatrix::SMatrix(InputIterator ibegin, unsigned int size, bool triang, bool lower) { - // assume iterator size <= matrix size (no check needed in AssignItr) - assert( size <= R::kSize); - for(unsigned int i=0; i::Evaluate(*this,ibegin,ibegin+size,triang,lower,false); -} - - -//============================================================================== -// Assignment and operator= for scalar types for matrices of size 1 -// compiles only for matrices of size 1 -//============================================================================== -template -SMatrix::SMatrix(const T& rhs) { - STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 ); - fRep[0] = rhs; -} - -template -SMatrix& SMatrix::operator=(const T& rhs) { - STATIC_CHECK( kSize == 1,SMatrix_NOT_of_size_1 ); - fRep[0] = rhs; - return *this; -} - -//============================================================================= -//============================================================================= - -template -template -SMatrix& SMatrix::operator=(const M& rhs) { - fRep = rhs.fRep; - return *this; -} - -template -template -SMatrix& SMatrix::operator=(const Expr& rhs) { - - Assign::Evaluate(*this, rhs); - return *this; -} - -//============================================================================= -// assign from an identity -template -SMatrix& SMatrix::operator= ( SMatrixIdentity ) { - for(unsigned int i=0; i -SMatrix& SMatrix::operator+=(const T& rhs) { - // self-addition with a scalar value - for(unsigned int i=0; i -template -SMatrix& SMatrix::operator+=(const SMatrix& rhs) { - // self-addition with another matrix (any representation) - // use operator+= of the representation object - fRep += rhs.fRep; - return *this; -} - - -template -template -SMatrix& SMatrix::operator+=(const Expr& rhs) { - // self-addition with an expression - PlusEquals::Evaluate(*this, rhs); - return *this; -} - - -//============================================================================== -// operator-= -//============================================================================== -template -SMatrix& SMatrix::operator-=(const T& rhs) { - // self-subtraction with a scalar value - for(unsigned int i=0; i -template -SMatrix& SMatrix::operator-=(const SMatrix& rhs) { - // self-subtraction with another matrix (any representation) - // use operator-= of the representation object - fRep -= rhs.fRep; - return *this; -} - - -template -template -SMatrix& SMatrix::operator-=(const Expr& rhs) { - // self-subtraction with an expression - MinusEquals::Evaluate(*this, rhs); - return *this; -} - -//============================================================================== -// operator*= -//============================================================================== -template -SMatrix& SMatrix::operator*=(const T& rhs) { - // case of multiplication with a scalar - for(unsigned int i=0; i -template -SMatrix& SMatrix::operator*=(const SMatrix& rhs) { - // self-multiplication with another matrix (will work only for square matrices) - // a temporary is needed and will be created automatically to store intermediate result - return operator=(*this * rhs); -} - -template -template -SMatrix& SMatrix::operator*=(const Expr& rhs) { - // self-multiplication with an expression (will work only for square matrices) - // a temporary is needed and will be created automatically to store intermediate result - return operator=(*this * rhs); -} - - -//============================================================================== -// operator/= (only for scalar values) -//============================================================================== -template -SMatrix& SMatrix::operator/=(const T& rhs) { - // division with a scalar - for(unsigned int i=0; i -bool SMatrix::operator==(const T& rhs) const { - bool rc = true; - for(unsigned int i=0; i -template -bool SMatrix::operator==(const SMatrix& rhs) const { - return fRep == rhs.fRep; -} - -template -template -bool SMatrix::operator==(const Expr& rhs) const { - bool rc = true; - for(unsigned int i=0; i -inline bool SMatrix::operator!=(const T& rhs) const { - return !operator==(rhs); -} - -template -inline bool SMatrix::operator!=(const SMatrix& rhs) const { - return !operator==(rhs); -} - -template -template -inline bool SMatrix::operator!=(const Expr& rhs) const { - return !operator==(rhs); -} - - -//============================================================================== -// operator> (element wise comparison) -//============================================================================== -template -bool SMatrix::operator>(const T& rhs) const { - bool rc = true; - for(unsigned int i=0; i rhs); - } - return rc; -} - -template -template -bool SMatrix::operator>(const SMatrix& rhs) const { - bool rc = true; - for(unsigned int i=0; i rhs.fRep[i]); - } - return rc; -} - -template -template -bool SMatrix::operator>(const Expr& rhs) const { - bool rc = true; - for(unsigned int i=0; i rhs.apply(i)); - } - return rc; -} - -//============================================================================== -// operator< (element wise comparison) -//============================================================================== -template -bool SMatrix::operator<(const T& rhs) const { - bool rc = true; - for(unsigned int i=0; i -template -bool SMatrix::operator<(const SMatrix& rhs) const { - bool rc = true; - for(unsigned int i=0; i -template -bool SMatrix::operator<(const Expr& rhs) const { - bool rc = true; - for(unsigned int i=0; i -inline bool SMatrix::Invert() { - STATIC_CHECK( D1 == D2,SMatrix_not_square); - return Inverter::Dinv((*this).fRep); -} - -// invert returning a new matrix -template -inline SMatrix SMatrix::Inverse(int & ifail) const { - SMatrix tmp(*this); - bool ok = tmp.Invert(); - ifail = 0; - if (!ok) ifail = 1; - return tmp; -} - -// fast inversion -template -inline bool SMatrix::InvertFast() { - STATIC_CHECK( D1 == D2,SMatrix_not_square); - return FastInverter::Dinv((*this).fRep); -} - -// fast inversion returning a new matrix -template -inline SMatrix SMatrix::InverseFast(int & ifail) const { - SMatrix tmp(*this); - bool ok = tmp.InvertFast(); - ifail = 0; - if (!ok) ifail = 1; - return tmp; -} - -// Choleski inversion for symmetric and positive defined matrices -template -inline bool SMatrix::InvertChol() { - STATIC_CHECK( D1 == D2,SMatrix_not_square); - return CholInverter::Dinv((*this).fRep); -} - -template -inline SMatrix SMatrix::InverseChol(int &ifail) const { - SMatrix tmp(*this); - bool ok = tmp.InvertChol(); - ifail = 0; - if (!ok) ifail = 1; - return tmp; -} - - - -//============================================================================== -// det -//============================================================================== -template -inline bool SMatrix::Det(T& det) { - STATIC_CHECK( D1 == D2,SMatrix_not_square); - // return Dfact, D1, D1>(*this,det); - //return Dfact(fRep, det); - return Determinant::Dfact(fRep, det); -} -template -inline bool SMatrix::Det2(T& det) const { - SMatrix tmp(*this); - return tmp.Det(det); -} - - -//============================================================================== -// place_in_row -//============================================================================== -template -template -SMatrix& SMatrix::Place_in_row(const SVector& rhs, - unsigned int row, - unsigned int col) { - - assert(col+D <= D2); - - for(unsigned int i=row*D2+col, j=0; j -template -SMatrix& SMatrix::Place_in_row(const VecExpr& rhs, - unsigned int row, - unsigned int col) { - - assert(col+D <= D2); - - for(unsigned int i=row*D2+col, j=0; j -template -SMatrix& SMatrix::Place_in_col(const SVector& rhs, - unsigned int row, - unsigned int col) { - - assert(row+D <= D1); - - for(unsigned int i=row*D2+col, j=0; j -template -SMatrix& SMatrix::Place_in_col(const VecExpr& rhs, - unsigned int row, - unsigned int col) { - - assert(row+D <= D1); - - for(unsigned int i=row*D2+col, j=0; j -template -SMatrix& SMatrix::Place_at(const SMatrix& rhs, - unsigned int row, - unsigned int col) { - PlaceMatrix::Evaluate(*this,rhs,row,col); - return *this; -} - -//============================================================================== -// place_at -//============================================================================== -template -template -SMatrix& SMatrix::Place_at(const Expr& rhs, - unsigned int row, - unsigned int col) { - PlaceExpr::Evaluate(*this,rhs,row,col); - return *this; -} - -//============================================================================== -// row -//============================================================================== -template -SVector SMatrix::Row(unsigned int therow) const { - - const unsigned int offset = therow*D2; - - /*static*/ SVector tmp; - for(unsigned int i=0; i -SVector SMatrix::Col(unsigned int thecol) const { - - /*static */ SVector tmp; - for(unsigned int i=0; i -std::ostream& SMatrix::Print(std::ostream& os) const { - const std::ios_base::fmtflags prevFmt = os.setf(std::ios::right,std::ios::adjustfield); - // os.setf(ios::fixed); - - os << "[ "; - for (unsigned int i=0; i < D1; ++i) { - for (unsigned int j=0; j < D2; ++j) { - os << std::setw(12) << fRep[i*D2+j]; - if ((!((j+1)%12)) && (j < D2-1)) - os << std::endl << " ..."; - } - if (i != D1 - 1) - os << std::endl << " "; - } - os << " ]"; - - if (prevFmt != os.flags() ) os.setf(prevFmt, std::ios::adjustfield); - return os; -} - -//============================================================================== -// Access functions -//============================================================================== -template -inline T SMatrix::apply(unsigned int i) const { return fRep[i]; } - -template -inline const T* SMatrix::Array() const { return fRep.Array(); } - -template -inline T* SMatrix::Array() { return fRep.Array(); } - -//============================================================================== -// Operators -//============================================================================== -template -inline const T& SMatrix::operator()(unsigned int i, unsigned int j) const { - return fRep(i,j); -} - -template -inline T& SMatrix::operator()(unsigned int i, unsigned int j) { - return fRep(i,j); -} - - -//============================================================================== -// Element access with At() -//============================================================================== -template -inline const T& SMatrix::At(unsigned int i, unsigned int j) const { - assert(i < D1); - assert(j < D2); - return fRep(i,j); -} - -template -inline T& SMatrix::At(unsigned int i, unsigned int j) { - assert(i < D1); - assert(j < D2); - return fRep(i,j); -} - -//============================================================================== -// STL interface -//============================================================================== -template -inline T * SMatrix::begin() { - return fRep.Array(); -} - -template -inline T * SMatrix::end() { - return fRep.Array() + R::kSize; -} - -template -inline const T * SMatrix::begin() const { - return fRep.Array(); -} - -template -inline const T * SMatrix::end() const { - return fRep.Array() + R::kSize; -} - - -template -template -void SMatrix::SetElements(InputIterator ibegin, InputIterator iend, bool triang, bool lower) { - // assume iterator size == matrix size when filling full matrix - AssignItr::Evaluate(*this,ibegin,iend,triang,lower); -} - -template -template -void SMatrix::SetElements(InputIterator ibegin, unsigned int size, bool triang, bool lower) { - // assume iterator size <= matrix size (no check to be done in AssignItr) - assert( size <= R::kSize); - AssignItr::Evaluate(*this,ibegin,ibegin+size,triang,lower,false); -} - - - -//============================================================================== -// SubMatrices and slices of columns and rows -//============================================================================== -template -template -SubVector SMatrix::SubRow(unsigned int therow, unsigned int col0 ) const { - - const unsigned int offset = therow*D2 + col0; - - STATIC_CHECK( SubVector::kSize <= D2,SVector_dimension_too_small); - assert(col0 + SubVector::kSize <= D2); - - SubVector tmp; - for(unsigned int i=0; i -template -SubVector SMatrix::SubCol(unsigned int thecol, unsigned int row0 ) const { - - const unsigned int offset = thecol + row0*D1; - - STATIC_CHECK( SubVector::kSize <= D1,SVector_dimension_too_small); - assert(row0 + SubVector::kSize <= D1); - - SubVector tmp; - for(unsigned int i=0; i -template -SubMatrix SMatrix::Sub(unsigned int row0, unsigned int col0) const { - - SubMatrix tmp; - RetrieveMatrix::Evaluate - (tmp,*this,row0,col0); - return tmp; -} - -//diagonal -template -SVector SMatrix::Diagonal( ) const { - - // only for squared matrices - STATIC_CHECK( D1 == D2,SMatrix_NOT_square ); - - SVector tmp; - for(unsigned int i=0; i -template -void SMatrix::SetDiagonal( const Vector & v) { - - // check size that size of vector is correct - STATIC_CHECK( ( ( D1 <= D2) && Vector::kSize == D1 ) || - ( ( D2 < D1 ) && Vector::kSize == D2 ), SVector_size_NOT_correct ); - - - for(unsigned int i=0; i -T SMatrix::Trace( ) const { - unsigned int diagSize = D1; - if (D2 < D1) diagSize = D2; - T trace = 0; - for(unsigned int i=0; i< diagSize; ++i) { - trace += fRep[ i*D2 + i] ; - } - return trace; -} - -//upper block -template -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION -SVector SMatrix::UpperBlock( ) const { -#else -template -SubVector SMatrix::UpperBlock( ) const { -#endif - // only for squared matrices - STATIC_CHECK( D1 == D2,SMatrix_NOT_square ); - -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION - SVector tmp; -#else - // N must be equal D1 *(D1 +1)/2 - STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size ); - SubVector tmp; -#endif - - int k = 0; - for(unsigned int i=0; i -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION -SVector SMatrix::LowerBlock( ) const { -#else -template -SubVector SMatrix::LowerBlock( ) const { -#endif - - // only for squared matrices - STATIC_CHECK( D1 == D2,SMatrix_NOT_square ); - -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION - SVector tmp; -#else - // N must be equal D1 *(D1 +1)/2 - STATIC_CHECK( SubVector::kSize == D1*(D1+1)/2,SVector_Wrong_Size ); - SubVector tmp; -#endif - - int k = 0; - for(unsigned int i=0; i -#ifndef UNSUPPORTED_TEMPLATE_EXPRESSION -SMatrix::SMatrix(const SVector & v, bool lower ) { -#else -template -SMatrix::SMatrix(const SVector & v, bool lower ) { -#endif - - // only for squared matrices - STATIC_CHECK( D1 == D2,SMatrix_NOT_square ); - -#ifdef UNSUPPORTED_TEMPLATE_EXPRESSION - STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size ); -#endif - - int k = 0; - if (lower) { - // case of lower block - for(unsigned int i=0; i -bool SMatrix::IsInUse( const T * p) const { - return p == fRep.Array(); -} - - - - - } // namespace Math - -} // namespace ROOT - - - -#endif