SMatrix.icc

// @(#)root/smatrix:$Name:  $:$Id: SMatrix.iccv 1.0 2005/11/24 12:00:00 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"


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);
}
 
#ifdef OLD_IMPL
   
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() {
  return Dsinv<T,D1,D1>(fArray);
}

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

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

//==============================================================================
// det
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline bool SMatrix<T,D1,D2>::Det(T& det) {
  return Dfact<SMatrix<T,D1,D1>, D1, D1>(*this,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];
}


  }  // namespace Math

}  // namespace ROOT
          


#endif