SMatrix.icc

// @(#)root/smatrix:$Name:  $:$Id: SMatrix.icc,v 1.2 2005/12/05 16:33:47 moneta Exp $
// 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 <iostream>
#include <iomanip>
#include <assert.h>
#include "Math/Dsinv.h"
#include "Math/Dsfact.h"
#include "Math/Dfact.h"
#include "Math/Dinv.h"
#include "Math/Functions.h"

#include "Math/StaticCheck.h"


namespace ROOT { 

  namespace Math { 



//==============================================================================
// Constructors
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>::SMatrix() {
  // operator=(0);   // if operator=(T ) is defined
    for(unsigned int i=0; i<D1*D2; ++i) {
      fArray[i] = 0;
    }
}

template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>::SMatrix(const SMatrix<T,D1,D2>& rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] = rhs.fArray[i];
  }
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>::SMatrix(const Expr<A,T,D1,D2>& rhs) {
  operator=(rhs);
}


//==============================================================================
// New Constructors from STL interfaces
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <class InputIterator>
SMatrix<T,D1,D2>::SMatrix(InputIterator begin, InputIterator end) {
  assert( std::distance(begin,end) <= kSize);
  std::copy(begin,end, fArray);
}

template <class T, unsigned int D1, unsigned int D2>
template <class InputIterator>
SMatrix<T,D1,D2>::SMatrix(InputIterator begin, unsigned int size) {
  assert( size <=  kSize);
  std::copy(begin,begin+size,fArray);
}


#ifdef OLD_IMPL
//=============================================================================
// old constructors (disabled)
//============================================================================
   
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>::SMatrix(const T& rhs, const bool diagonal) {
  if(diagonal) {
    // symmetric matrix!
    assert(D1==D2);
    // 50% faster than with nested loops
    for(unsigned int i=0; i<D1*D2; ++i)
      fArray[i] = 0;
    for(unsigned int i=0; i<D1; ++i)
      fArray[i*D1+i] = rhs;
  } else {
    for(unsigned int i=0; i<D1*D2; ++i) {
      fArray[i] = rhs;
    }
  }
}

template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>::SMatrix(const SVector<T,D1>& rhs) {
  for(unsigned int i=0; i<D1; ++i) {
    const unsigned int row = i*D1;
    for(unsigned int j=0; j<D1; ++j) {
      fArray[row+j] = rhs[i]*rhs[j];
    }
  }
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>::SMatrix(const Expr<A,T,D1>& rhs) {
  for(unsigned int i=0; i<D1; ++i) {
    const unsigned int row = i*D1;
    for(unsigned int j=0; j<D1; ++j) {
      fArray[row+j] = rhs.apply(i)*rhs.apply(j);
    }
  }
}

template <class T, unsigned int D1, unsigned int D2>
template <class T1>
SMatrix<T,D1,D2>::SMatrix(const T1* a, const bool triang, const unsigned int len) {
  // fill from array with lower triangular matrix?
  if(triang==false) {
    for(unsigned int i=0; i<len; ++i) {
      fArray[i] = a[i];
    }
    if(len<D1*D2)
      for(unsigned int i=len; i<D1*D2; ++i)
	fArray[i] = 0.;
  } else {
    // symmetric matrix required!
    assert(D1==D2);

    unsigned int k=0;
    if(len == (Square(D1)+D1)/2) {
      for(unsigned int i=0; i<D1; ++i) {
	const unsigned int row = i*D1;
	for(unsigned int j=0; j<=i && k<len; ++j) {
	  fArray[row+j] = fArray[j*D1+i] = a[k++];
	}
      }
    } else {
      unsigned int k=0;
      for(unsigned int i=0; i<D1; ++i) {
	const unsigned int row = i*D1;
	for(unsigned int j=0; j<=i; ++j) {
	  if(k<len)
	    fArray[row+j] = fArray[j*D1+i] = a[k++];
	  else
	    fArray[row+j] = fArray[j*D1+i] = 0.;
	} // for j
      } // for i
    } // if len == 
  } // if triang
}


//==============================================================================
// operator=
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator=(const T&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] = rhs;
  }
  return *this;
}

#endif


template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator=(const Expr<A,T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] = rhs.apply(i);
  }
  return *this;
}


//==============================================================================
// operator+=
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator+=(const SMatrix<T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] += rhs.apply(i);
  }
  return *this;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator+=(const Expr<A,T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] += rhs.apply(i);
  }
  return *this;
}


//==============================================================================
// operator-=
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator-=(const SMatrix<T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] -= rhs.apply(i);
  }
  return *this;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator-=(const Expr<A,T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] -= rhs.apply(i);
  }
  return *this;
}


//==============================================================================
// operator*=
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator*=(const SMatrix<T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] *= rhs.apply(i);
  }
  return *this;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator*=(const Expr<A,T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] *= rhs.apply(i);
  }
  return *this;
}


//==============================================================================
// operator/=
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator/=(const SMatrix<T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] /= rhs.apply(i);
  }
  return *this;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::operator/=(const Expr<A,T,D1,D2>&  rhs) {
  for(unsigned int i=0; i<D1*D2; ++i) {
    fArray[i] /= rhs.apply(i);
  }
  return *this;
}


//==============================================================================
// operator== (element wise comparison)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator==(const T& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] == rhs);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator==(const SMatrix<T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] == rhs.fArray[i]);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
bool SMatrix<T,D1,D2>::operator==(const Expr<A,T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] == rhs.apply(i));
  }
  return rc;
}

//==============================================================================
// operator!= (element wise comparison)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::operator!=(const T& rhs) const {
  return !operator==(rhs);
}

template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::operator!=(const SMatrix<T,D1,D2>& rhs) const {
  return !operator==(rhs);
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
inline bool SMatrix<T,D1,D2>::operator!=(const Expr<A,T,D1,D2>& rhs) const {
  return !operator==(rhs);
}


//==============================================================================
// operator> (element wise comparison)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator>(const T& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] > rhs);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator>(const SMatrix<T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] > rhs.fArray[i]);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
bool SMatrix<T,D1,D2>::operator>(const Expr<A,T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] > rhs.apply(i));
  }
  return rc;
}

//==============================================================================
// operator< (element wise comparison)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator<(const T& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] < rhs);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
bool SMatrix<T,D1,D2>::operator<(const SMatrix<T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] < rhs.fArray[i]);
  }
  return rc;
}

template <class T, unsigned int D1, unsigned int D2>
template <class A>
bool SMatrix<T,D1,D2>::operator<(const Expr<A,T,D1,D2>& rhs) const {
  bool rc = true;
  for(unsigned int i=0; i<D1*D2; ++i) {
    rc = rc && (fArray[i] < rhs.apply(i));
  }
  return rc;
}

//==============================================================================
// sinvert
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Sinvert() {
  STATIC_CHECK( D1 == D2,SMatrix_not_square); 
  return Dsinv<T,D1,D1>(fArray);
}

template <class T, unsigned int D1, unsigned int D2>
inline SMatrix<T,D1,D2> SMatrix<T,D1,D2>::Sinverse() const {
  SMatrix<T,D1,D2> tmp(*this);
  tmp.Sinvert();
  return tmp;
}



//==============================================================================
// sdet
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Sdet(T& det) {
  STATIC_CHECK( D1 == D2,SMatrix_not_square); 
  return Dsfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
}
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Sdet2(T& det) const {
  SMatrix<T,D1,D2> tmp(*this);
  return tmp.Sdet(det);
}


//==============================================================================
// invert
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Invert() {
  STATIC_CHECK( D1 == D2,SMatrix_not_square); 
  return Inverter<D2>::Dinv(*this);
}

template <class T, unsigned int D1, unsigned int D2>
inline SMatrix<T,D1,D2> SMatrix<T,D1,D2>::Inverse() const {
  SMatrix<T,D1,D2> tmp(*this);
  tmp.Invert();
  return tmp;
}

//==============================================================================
// det
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Det(T& det) {
  STATIC_CHECK( D1 == D2,SMatrix_not_square); 
  return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,det);
}
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Det2(T& det) const {
  SMatrix<T,D1,D2> tmp(*this);
  return tmp.Det(det);
}


//==============================================================================
// place_in_row
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <unsigned int D>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_row(const SVector<T,D>& rhs,
						 const unsigned int row,
						 const unsigned int col) {

  assert(col+D <= D2);

  for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
    fArray[i] = rhs.apply(j);
  }
  return *this;
}

//==============================================================================
// place_in_row
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <class A, unsigned int D>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_row(const Expr<A,T,D>& rhs,
						 const unsigned int row,
						 const unsigned int col) {

  assert(col+D <= D2);

  for(unsigned int i=row*D2+col, j=0; j<D; ++i, ++j) {
    fArray[i] = rhs.apply(j);
  }
  return *this;
}

//==============================================================================
// place_in_col
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <unsigned int D>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_col(const SVector<T,D>& rhs,
						 const unsigned int row,
						 const unsigned int col) {

  assert(row+D <= D1);

  for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
    fArray[i] = rhs.apply(j);
  }
  return *this;
}

//==============================================================================
// place_in_col
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <class A, unsigned int D>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_in_col(const Expr<A,T,D>& rhs,
						 const unsigned int row,
						 const unsigned int col) {

  assert(row+D <= D1);

  for(unsigned int i=row*D2+col, j=0; j<D; i+=D2, ++j) {
    fArray[i] = rhs.apply(j);
  }
  return *this;
}

//==============================================================================
// place_at
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <unsigned int D3, unsigned int D4>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_at(const SMatrix<T,D3,D4>& rhs,
					     const unsigned int row,
					     const unsigned int col) {
  assert(row+D3 <= D1 && col+D4 <= D2);

  const unsigned int offset = row*D2+col;

  for(unsigned int i=0; i<D3*D4; ++i) {
    fArray[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
  }

  return *this;
}

//==============================================================================
// place_at
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <class A, unsigned int D3, unsigned int D4>
SMatrix<T,D1,D2>& SMatrix<T,D1,D2>::Place_at(const Expr<A,T,D3,D4>& rhs,
					     const unsigned int row,
					     const unsigned int col) {
  assert(row+D3 <= D1 && col+D4 <= D2);

  const unsigned int offset = row*D2+col;

  for(unsigned int i=0; i<D3*D4; ++i) {
    fArray[offset+(i/D4)*D2+i%D4] = rhs.apply(i);
  }

  return *this;
}

//==============================================================================
// row
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SVector<T,D2> SMatrix<T,D1,D2>::Row(const unsigned int therow) const {

  const unsigned int offset = therow*D2;

  /*static*/  SVector<T,D2> tmp;
  for(unsigned int i=0; i<D2; ++i) {
    tmp[i] = fArray[offset+i];
  }
  return tmp;
}

//==============================================================================
// col
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
SVector<T,D1> SMatrix<T,D1,D2>::Col(const unsigned int thecol) const {

  /*static */ SVector<T,D1> tmp;
  for(unsigned int i=0; i<D1; ++i) {
    tmp[i] = fArray[thecol+i*D2];
  }
  return tmp;
}

//==============================================================================
// print
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
std::ostream& SMatrix<T,D1,D2>::Print(std::ostream& os) const {
  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) << fArray[i*D2+j];
      if ((!((j+1)%12)) && (j < D2-1))
	os << std::endl << "         ...";
    }
    if (i != D1 - 1)
      os << std::endl  << "  ";
  }
  os << " ]";
  
  return os;
}

//==============================================================================
// Access functions
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline T SMatrix<T,D1,D2>::apply(unsigned int i) const { return fArray[i]; }

template <class T, unsigned int D1, unsigned int D2>
inline const T* SMatrix<T,D1,D2>::Array() const { return fArray; }

template <class T, unsigned int D1, unsigned int D2>
inline T* SMatrix<T,D1,D2>::Array() { return fArray; }

//==============================================================================
// Operators
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline const T& SMatrix<T,D1,D2>::operator()(unsigned int i, unsigned int j) const {
  return fArray[i*D2+j];
}

template <class T, unsigned int D1, unsigned int D2>
inline T& SMatrix<T,D1,D2>::operator()(unsigned int i, unsigned int j) {
  return fArray[i*D2+j];
}

//==============================================================================
// STL interface
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline T * SMatrix<T,D1,D2>::begin() {
  return fArray;
}

template <class T, unsigned int D1, unsigned int D2>
inline T * SMatrix<T,D1,D2>::end() {
  return fArray + kSize;
}

template <class T, unsigned int D1, unsigned int D2>
inline const T * SMatrix<T,D1,D2>::begin() const {
  return fArray;
}

template <class T, unsigned int D1, unsigned int D2>
inline const T * SMatrix<T,D1,D2>::end() const {
  return fArray + kSize;
}



//==============================================================================
// SubMatrices and slices of columns and rows
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
template <unsigned int N>  
SVector<T,N> SMatrix<T,D1,D2>::SubRow(const unsigned int therow, const unsigned int col0 ) const { 

  const unsigned int offset = therow*D2 + col0;

  STATIC_CHECK( N <= D2,SVector_dimension_too_small); 
  assert(col0 + N <= D2);

  SVector<T,N> tmp;
  for(unsigned int i=0; i<N; ++i) {
    tmp[i] = fArray[offset+i];
  }
  return tmp;
}

template <class T, unsigned int D1, unsigned int D2>
template <unsigned int N>  
SVector<T,N> SMatrix<T,D1,D2>::SubCol(const unsigned int thecol, const unsigned int row0 ) const { 

  const unsigned int offset = thecol + row0*D1;

  STATIC_CHECK( N <= D1,SVector_dimension_too_small); 
  assert(row0 + N <= D1);

  SVector<T,N> tmp;
  for(unsigned int i=0; i<N; ++i) {
    tmp[i] = fArray[offset+i*D1];
  }
  return tmp;
}

// submatrix
template <class T, unsigned int D1, unsigned int D2>
template <unsigned int N1, unsigned int N2>  
SMatrix<T,N1,N2> SMatrix<T,D1,D2>::SubMatrix(const unsigned int row0, const unsigned int col0) const { 

  STATIC_CHECK( N1 <= D1,Smatrix_dimension1_too_small); 
  STATIC_CHECK( N2 <= D2,Smatrix_dimension2_too_small); 
  assert(row0 + N1 <= D1);
  assert(col0 + N2 <= D2);

  SMatrix<T,N1,N2> tmp;
  for(unsigned int i=0; i<N1; ++i) {
    for(unsigned int j=0; j<N2; ++j) {
      tmp(i,j) = fArray[(i+row0)*D2+j+col0];
    }
  }
  return tmp;
}

//diagonal
template <class T, unsigned int D1, unsigned int D2>
SVector<T,D1> SMatrix<T,D1,D2>::Diagonal( ) const { 

  // only for squared matrices
  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );

  SVector<T,D1> tmp;
  for(unsigned int i=0; i<D1; ++i) {
    tmp[i] = fArray[ i*D2 + i];
  }
  return tmp;
}

//upper block
template <class T, unsigned int D1, unsigned int D2>
#ifndef UNSUPPORTED
SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2>::UpperBlock( ) const { 
#else
template <unsigned int N>  
SVector<T,N> SMatrix<T,D1,D2>::UpperBlock( ) const { 
#endif
  // only for squared matrices
  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );

#ifndef UNSUPPORTED
  SVector<T, D1 * (D2 +1)/2 >  tmp;
#else
  // N must be equal D1 *(D1 +1)/2
  STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
  SVector<T,N> tmp;
#endif

  int k = 0;
  for(unsigned int i=0; i<D1; ++i) {
    for(unsigned int j=i; j<D2; ++j) {
      tmp[k] = fArray[ i*D2 + j];
      k++;
    }
  }
  return tmp;
}

//lower block
template <class T, unsigned int D1, unsigned int D2>
#ifndef UNSUPPORTED
SVector<T, D1 * (D2 +1)/2 >  SMatrix<T,D1,D2>::LowerBlock( ) const { 
#else
template <unsigned int N>  
SVector<T,N> SMatrix<T,D1,D2>::LowerBlock( ) const { 
#endif

  // only for squared matrices
  STATIC_CHECK( D1 == D2,SMatrix_NOT_square );

#ifndef UNSUPPORTED
  SVector<T, D1 * (D2 +1)/2 >  tmp;
#else
  // N must be equal D1 *(D1 +1)/2
  STATIC_CHECK( N == D1*(D1+1)/2,SVector_Wrong_Size );
  SVector<T,N> tmp;
#endif

  int k = 0;
  for(unsigned int i=0; i<D1; ++i) {
    for(unsigned int j=0; j<=i; ++j) {
      tmp[k] = fArray[ i*D2 + j];
      k++;
    }
  }
  return tmp;
}




  }  // namespace Math

}  // namespace ROOT
          


#endif