// smatrix.h
 
#ifndef smatrix_H
#define smatrix_H

#undef  trf_info
#define trf_error

#include "nvector.h"

template <class T> class matrix;
template <class T> class smatrix;
template <class T>
smatrix<T> operator% (const  matrix<T>&  mat, const smatrix<T>&  sma);

template <class T>
class smatrix : public array<T> {
private:
  int _nrow;
  int ij( int i, int j ) const { return i>j ? i*(i+1)/2+j : j*(j+1)/2+i ; };
  void ostr(std::ostream& stream) const;
public:
  explicit smatrix(int N) : array<T>( N*(N+1)/2 ) { _nrow = N; };
  smatrix(const smatrix<T>& sma);
  virtual ~smatrix(void) { };
  smatrix<T> &operator=(const smatrix<T>& sma);
  int nrow() const { return _nrow; };
  T& operator()( int i, int j ) const { return data[ij(i,j)]; };
  smatrix<T> &operator+=(const smatrix<T>& sma);
  smatrix<T> &operator-=(const smatrix<T>& sma);
  smatrix<T> &operator*=(const T& tval);
#if defined(DEFECT_GUIDING_DECLS)
  friend smatrix<T> operator+(const smatrix<T>& sma1, const smatrix<T>& sma2);
  friend smatrix<T> operator-(const smatrix<T>& sma1, const smatrix<T>& sma2);
  friend smatrix<T> operator%(const  matrix<T>&  mat, const smatrix<T>&  sma);
#elif defined(DEFECT_NO_EXPLICIT_QUALIFICATION)
  template<class S>
  friend smatrix<S> operator+(const smatrix<S>& sma1, const smatrix<S>& sma2);
  template<class S>
  friend smatrix<S> operator-(const smatrix<S>& sma1, const smatrix<S>& sma2);
  template<class S>
  friend smatrix<S> operator%(const  matrix<S>&  mat, const smatrix<S>&  sma);
#else
  friend smatrix<T> operator+ <> (const smatrix<T>& sma1, const smatrix<T>& sma2);
  friend smatrix<T> operator- <> (const smatrix<T>& sma1, const smatrix<T>& sma2);
  friend smatrix<T> operator% <> (const  matrix<T>&  mat, const smatrix<T>&  sma);
#endif
};

template <class T>
void normalize(smatrix<T>& sma);

// copy
template <class T>
inline smatrix<T>::smatrix(const smatrix<T>& sma)
: array<T>(sma)
{
#ifdef trf_info
 std::cout << "copying a symmetric matrix\n";
#endif
 _nrow = sma._nrow;
}
 
// assignment
template <class T>
inline smatrix<T> &smatrix<T>::operator=(const smatrix<T>& sma)
{
#ifdef trf_info
  std::cout << "assigning a symmetric matrix\n";
#endif
 (array<T>&)*this = sma;
 _nrow = sma._nrow;
 return *this;
}
 
// inversion
extern int invert( smatrix<double>& sma );

// need to add code for dim>1
/*
template <class T>
smatrix<T> smatrix<T>::invert( ) {
 smatrix<T> sma(_nrow);
 if ( _nrow == 1 ) sma(0,0) = 1.0/data[0];
 else fill( T(0) );
 return sma;
}
*/

// +=
template <class T>
inline smatrix<T> &smatrix<T>::operator+=(const smatrix<T>& sma)
{
 (array<T>&)*this += sma;
 return *this;
}

// -=
template <class T>
inline smatrix<T> &smatrix<T>::operator-=(const smatrix<T>& sma)
{
 (array<T>&)*this -= sma;
 return *this;
}

// *=
template <class T>
inline smatrix<T> &smatrix<T>::operator*=(const T& t)
{
 (array<T>&)*this *= t;
 return *this;
}

template <class T>
inline smatrix<T> operator*(const T& val,const smatrix<T>& sma)
{
 smatrix<T> lhs = sma;
 lhs *= val;
 return lhs;
}

template <class T>
inline smatrix<T> operator*(const smatrix<T>& sma,const T& val)
{
 smatrix<T> lhs = sma;
 lhs *= val;
 return lhs;
}

// +
template <class T>
inline smatrix<T> operator+(const smatrix<T>& sma1,const smatrix<T>& sma2)
{
  smatrix<T> sum = sma1;
  sum += sma2;
  return sum;
}

// -
template <class T>
inline smatrix<T> operator-(const smatrix<T>& sma1,const smatrix<T>& sma2)
{
  smatrix<T> sum = sma1;
  sum -= sma2;
  return sum;
}

//

//
// multiplication: t = S % v == vT * S * v
//
template <class T>
T operator%(const smatrix<T>& sma, const nvector<T>& vec);

#ifndef __xlC__
#include "smatrix.c"
#endif

#endif