#ifndef FILE_MATRIX_EXPR
#define FILE_MATRIX_EXPR

/**************************************************************************/
/* File:   matrix.hpp                                                     */
/* Author: Joachim Schoeberl                                              */
/* Date:   01. Jan. 02                                                    */
/**************************************************************************/


template <int H, int W, typename T> class Mat;

///
extern void CheckMatRange(int h, int w, int i);
///
extern void CheckMatRange(int h, int w, int i, int j);



/**
   A simple matrix.
   Has height, width and data-pointer. 
   No memory allocation/deallocation. User must provide memory.
 */
template <typename T = double>
class FlatMatrix : public MatExpr<FlatMatrix<T> >
{
protected:
  /// the height
  int h;
  /// the width
  int w;
  /// the data
  T * data;
public:

  /// element type
  typedef T TELEM;
  /// scalar type of elements (double or Complex)
  typedef typename mat_traits<T>::TSCAL TSCAL;


  /// nothing done in default constructor
  FlatMatrix () throw () { ; }
  
  /// set height, width, and mem
  FlatMatrix (int ah, int aw, T * adata) throw ()
    : h(ah), w(aw), data(adata) { ; }
  
  /// set height = width, and mem
  FlatMatrix (int ah, T * adata) throw()
    : h(ah), w(ah), data(adata) { ; }

  /// allocates at local heap
  FlatMatrix (int ah, int aw, LocalHeap & lh) throw (LocalHeapOverflow)
    : h(ah), w(aw), data((T*)lh.Alloc(ah*aw*sizeof(T))) { ; }
  
  /// allocates at local heap
  FlatMatrix (int ah, LocalHeap & lh) throw (LocalHeapOverflow)
    : h(ah), w(ah), data((T*)lh.Alloc(ah*ah*sizeof(T))) { ; }
  
  /// copy constructor. copies pointers, not contents
  FlatMatrix (const FlatMatrix & m) throw () 
    : MatExpr<FlatMatrix> (), h(m.h), w(m.w) , data(m.data)
  { ; }
  

  /// useful to put FlatMatrix over other matrix
  template <typename T2>
  explicit FlatMatrix (const MatExpr<T2> & m)
    : h(m.Height()), w(m.Width()),
      data(const_cast<T*>(&m.Spec()(0,0)))  
  { ; }
  
  /// useful to put FlatMatrix over other Mat
  template <int H, int W>
  FlatMatrix (const Mat<H,W,TSCAL> & m) throw()
    : h(H), w(W), data(const_cast<T*>(&m(0,0)))
  { ; }

  /// do nothing 
  ~FlatMatrix () throw() { ; }

  /// set size, and assign mem
  void AssignMemory (int ah, int aw, LocalHeap & lh)  throw (LocalHeapOverflow)
  {
    h = ah;
    w = aw;
    data = (T*)lh.Alloc(h*w*sizeof(T));
  }
  
  /// set size, and assign mem
  void AssignMemory (int ah, int aw, T * mem) throw()
  {
    h = ah;
    w = aw;
    data = mem;
  }
  

  /// assign contents
  template<typename TB>
  FlatMatrix & operator= (const Expr<TB> & m)
  {
    return MatExpr<FlatMatrix>::operator= (m);
  }

  /// copy contents
  FlatMatrix & operator= (const FlatMatrix & m) throw()
  {
    for (int i = 0; i < h*w; i++)
      data[i] = m(i);
    return *this;
  }

  /// assign constant
  FlatMatrix & operator= (TSCAL s) throw()
  {
    for (int i = 0; i < h*w; i++)
      data[i] = s; 
    return *this;
  }

  /// copy size and pointers
  FlatMatrix & Assign (const FlatMatrix & m) throw()
  {
    h = m.h;
    w = m.w;
    data = m.data;
    return *this;
  }


  /// access operator, linear access
  TELEM & operator() (int i)
  { 
#ifdef CHECK_RANGE
    CheckMatRange(h,w,i);
#endif
    return data[i]; 
  }

  /// access operator
  TELEM & operator() (int i, int j) 
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,w,i,j);
#endif
    return data[i*w+j]; 
  }


  /// access operator, linear access
  const TELEM & operator() (int i) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,w,i);
#endif
    return data[i]; 
  }

  /// access operator
  const TELEM & operator() (int i, int j) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,w,i,j);
#endif
    return data[i*w+j]; 
  }

  /// the height
  int Height () const throw() { return h; }

  /// the width
  int Width () const throw() { return w; }
};



/**
   A Matrix class with memory allocation/deallocation
 */
template <typename T = double>
class Matrix : public FlatMatrix<T>
{
public:

  /// element type
  typedef T TELEM;
  /// scalar type of elements (double or Complex)
  typedef typename mat_traits<T>::TSCAL TSCAL;

  /// default constructor
  Matrix () throw () : FlatMatrix<T> (0, 0) { ; }

  /// allocate matrix of size ah * ah
  Matrix (int ah) : FlatMatrix<T> (ah, new T[ah*ah]) { ; }

  /// allocate matrix of size ah * aw
  Matrix (int ah, int aw) : FlatMatrix<T> (ah, aw, new T[ah*aw]) { ; }

  /// delete memory
  ~Matrix() { delete [] this->data; }

  /// sets new size of matrix
  void SetSize(int ah, int aw)
  {
    if (this->h == ah && this->w == aw) return;
    delete [] this->data;
    this->h = ah;
    this->w = aw;
    this->data = new T[this->h * this->w];
  }

  /// sets new size of matrix
  void SetSize(int ah)
  {
    if (this->h == ah && this->w == ah) return;
    delete [] this->data;
    this->h = ah;
    this->w = ah;
    this->data = new T[this->h * this->w];
  }

  /// assign matrix, sizes must match
  template<typename TB>
  Matrix & operator= (const Expr<TB> & m)
  {
    MatExpr<FlatMatrix<T> >::operator= (m);
    return *this;
  }

  /// fill matrix with scalar
  Matrix & operator= (TSCAL s)
  {
    FlatMatrix<T>::operator= (s);
    return *this;
  }
};




/**
   A matrix of fixed size.
   Useful as entry type in system matrices,...
 */
template <int H, int W = H, typename T = double>
class Mat : public MatExpr<Mat<H,W,T> >
{
  T data[H*W];
public:
  typedef T TELEM;
  typedef typename mat_traits<T>::TSCAL TSCAL;
  typedef Vec<H, typename mat_traits<T>::TV_COL> TV_COL;
  typedef Vec<W, typename mat_traits<T>::TV_ROW> TV_ROW;
  enum { HEIGHT = H };
  enum { WIDTH  = W };

  /// do not initialize 
  Mat () throw () { ; }

  /// copy matrix
  Mat (const Mat & m) throw()
    : MatExpr<Mat> ()
  {
    (*this) = m;
  }

  /// assign values
  template<typename TB>
  Mat (const Expr<TB> & m)
  {
    MatExpr<Mat>::operator= (m);
  }

  /// fill with scalar
  explicit Mat (TSCAL s) throw()
  {
    (*this) = s;
  }

  /// assign values
  template<typename TB>
  Mat & operator= (const Expr<TB> & m)
  {
    MatExpr<Mat>::operator= (m);
    return *this;
  }

  /// copy matrix
  Mat & operator= (const Mat & m) throw()
  {
    for (int i = 0; i < H*W; i++)
      data[i] = m.data[i];
    return *this;
  }
 
  /// fill values
  Mat & operator= (TSCAL s) throw()
  {
    for (int i = 0; i < H*W; i++)
      data[i] = s;
    return *this;
  }

  ///
  TELEM & operator() (int i) { return data[i]; }
  ///
  TELEM & operator() (int i, int j) { return data[i*W+j]; }
  ///
  const TELEM & operator() (int i) const { return data[i]; }
  ///
  const TELEM & operator() (int i, int j) const { return data[i*W+j]; }

  /// the height
  int Height () const throw() { return H; }
  /// the width
  int Width () const throw() { return W; }
};













/**
   A Matrix with width known at compile time
   No memory allocation/deallocation. User must provide memory.
 */
template <int W, typename T = double>
class FlatMatrixFixWidth : public MatExpr<FlatMatrixFixWidth<W,T> >
{
protected:
  /// the data
  T * data;
  /// the height
  int h;
public:
  /// entry type
  typedef T TELEM;
  /// scalar type of entry
  typedef typename mat_traits<T>::TSCAL TSCAL;

  /// nothing done in default constructor
  FlatMatrixFixWidth () throw()
    : data(0), h(0) { ; }
  
  /// set height and mem
  FlatMatrixFixWidth (int ah, T * adata) throw()
    : data(adata), h(ah) { ; }

  /// allocates at local heap
  FlatMatrixFixWidth (int ah, LocalHeap & lh) throw (LocalHeapOverflow)
    : data((T*)lh.Alloc(ah*W*sizeof(T))), h(ah) { ; }
  
  /// copy constructor. copies pointers, not contents
  FlatMatrixFixWidth (const FlatMatrixFixWidth & m) throw()
    : MatExpr<FlatMatrixFixWidth> (), data(m.data), h(m.h)
  { ; }

  /// copy constructor. copies pointers, not contents
  FlatMatrixFixWidth (FlatMatrix<TSCAL> & m) throw()
    : data(&m(0)), h(m.Height())
  { ; }

  template <int H>
  FlatMatrixFixWidth (const Mat<H,W,TSCAL> & m) throw()
    : data(const_cast<T*>(&m(0,0))), h(H)
  { ; }  

  /// do nothing 
  ~FlatMatrixFixWidth () throw() { ; }

  /// set size, and assign mem
  void AssignMemory (int ah, LocalHeap & lh) throw (LocalHeapOverflow)
  {
    h = ah;
    data = (T*)lh.Alloc(h*W*sizeof(T));
  }
  
  /// set size, and assign mem
  void AssignMemory (int ah, T * mem) throw()
  {
    h = ah;
    data = mem;
  }
  

  /// assign contents
  template<typename TB>
  FlatMatrixFixWidth & operator= (const Expr<TB> & m)
  {
    return MatExpr<FlatMatrixFixWidth>::operator= (m);
  }

  /// copy contents
  FlatMatrixFixWidth & operator= (const FlatMatrixFixWidth & m) throw()
  {
    for (int i = 0; i < h*W; i++)
      data[i] = m(i);
    return *this;
  }

  /// assign constant
  FlatMatrixFixWidth & operator= (TSCAL s) throw()
  {
    for (int i = 0; i < h*W; i++)
      data[i] = s; 
    return *this;
  }

  /// copy size and pointers
  FlatMatrixFixWidth & Assign (const FlatMatrixFixWidth & m) throw()
  {
    h = m.h;
    data = m.data;
    return *this;
  }

  /// access operator, linear access
  TELEM & operator() (int i)
  { 
#ifdef CHECK_RANGE
    CheckMatRange(h,W,i);
#endif
    return data[i]; 
  }

  /// access operator
  TELEM & operator() (int i, int j) 
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,W,i,j);
#endif
    return data[i*W+j]; 
  }


  /// access operator, linear access
  const TELEM & operator() (int i) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,W,i);
#endif
    return data[i]; 
  }

  /// access operator
  const TELEM & operator() (int i, int j) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(h,W,i,j);
#endif
    return data[i*W+j]; 
  }

  /// the height
  int Height () const throw() { return h; }

  /// the width
  int Width () const throw() { return W; }
};


/**
   A Matrix class with memory allocation/deallocation
 */
template <int W, typename T = double>
class MatrixFixWidth : public FlatMatrixFixWidth<W,T>
{
public:
  typedef typename mat_traits<T>::TSCAL TSCAL;

  MatrixFixWidth () { ; }

  /// allocate matrix of size ah * ah
  MatrixFixWidth (int ah) : FlatMatrixFixWidth<W,T> (ah, new T[ah*W]) { ; }

  /// delete memory
  ~MatrixFixWidth() { delete [] this->data; }

  /// sets new size of matrix
  void SetSize(int ah)
  {
    if (this->h == ah) return;
    delete [] this->data;
    this->h = ah;
    this->data = new T[this->h*this->W];
  }

  /// assign matrix, sizes must match
  template<typename TB>
  MatrixFixWidth & operator= (const Expr<TB> & m)
  {
    MatExpr<FlatMatrixFixWidth<W,T> >::operator= (m);
    return *this;
  }

  /// fill matrix with scalar
  MatrixFixWidth & operator= (TSCAL s)
  {
    FlatMatrixFixWidth<W,T>::operator= (s);
    return *this;
  }
};






















/**
   A Matrix which height is known at compile time
   No memory allocation/deallocation. User must provide memory.
   Matrix is stored colum-wise
 */
template <int H, typename T = double>
class FlatMatrixFixHeight : public MatExpr<FlatMatrixFixHeight<H,T> >
{
protected:
  /// the data
  T * data;
  /// the width
  int w;
public:
  ///
  typedef T TELEM;
  ///
  typedef typename mat_traits<T>::TSCAL TSCAL;

  /// 
  FlatMatrixFixHeight () throw()
    : data(0), w(0) { ; }
  
  /// set height and mem
  FlatMatrixFixHeight (int aw, T * adata) throw()
    : data(adata), w(aw) { ; }

  /// allocates at local heap
  FlatMatrixFixHeight (int aw, LocalHeap & lh) throw (LocalHeapOverflow)
    : data((T*)lh.Alloc(aw*H*sizeof(T))), w(aw) { ; }
  
  /// copy constructor. copies pointers, not contents
  FlatMatrixFixHeight (const FlatMatrixFixHeight & m)
    : data(m.data), w(m.w)
  { ; }
  

  /// do nothing 
  ~FlatMatrixFixHeight () throw() { ; }

  /// set size, and assign mem
  void AssignMemory (int aw, LocalHeap & lh) throw (LocalHeapOverflow)
  {
    w = aw;
    data = (T*)lh.Alloc(w*H*sizeof(T));
  }
  
  /// set size, and assign mem
  void AssignMemory (int aw, T * mem) 
  {
    w = aw;
    data = mem;
  }
  

  /// assign contents
  template<typename TB>
  FlatMatrixFixHeight & operator= (const Expr<TB> & m)
  {
    for (int j = 0; j < w; j++)
      for (int i = 0; i < H; i++)
	(*this)(i,j) = m.Spec()(i,j);
    return *this;
  }

  /// copy contents
  FlatMatrixFixHeight & operator= (const FlatMatrixFixHeight & m)
  {
    for (int i = 0; i < w*H; i++)
      data[i] = m(i);
    return *this;
  }

  /// assign constant
  FlatMatrixFixHeight & operator= (TSCAL s)
  {
    for (int i = 0; i < w*H; i++)
      data[i] = s; 
    return *this;
  }

  /// copy size and pointers
  FlatMatrixFixHeight & Assign (const FlatMatrixFixHeight & m)
  {
    w = m.w;
    data = m.data;
    return *this;
  }


  /// access operator, linear access
  TELEM & operator() (int i)
  { 
#ifdef CHECK_RANGE
    CheckMatRange(H,w,i);
#endif
    return data[i]; 
  }

  /// access operator
  TELEM & operator() (int i, int j) 
  {
#ifdef CHECK_RANGE
    CheckMatRange(H,w,i,j);
#endif
    return data[i+j*H]; 
  }


  /// access operator, linear access
  const TELEM & operator() (int i) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(H,w,i);
#endif
    return data[i]; 
  }

  /// access operator
  const TELEM & operator() (int i, int j) const
  {
#ifdef CHECK_RANGE
    CheckMatRange(H,w,i,j);
#endif
    return data[i+j*H]; 
  }


  /// the height
  int Height () const { return H; }

  /// the width
  int Width () const { return w; }

  enum { IS_LINEAR = 0 };
  //  static bool IsLinear() { return false; }
};


/**
   A Matrix class with memory allocation/deallocation
 */
template <int H, typename T = double>
class MatrixFixHeight : public FlatMatrixFixHeight<H,T>
{
public:
  typedef typename mat_traits<T>::TSCAL TSCAL;

  /// allocate matrix of size ah * ah
  MatrixFixHeight (int aw) : FlatMatrixFixHeight<H,T> (aw, new T[aw*H]) { ; }

  /// delete memory
  ~MatrixFixHeight() { delete [] this->data; }

  /// sets new size of matrix
  void SetSize(int aw)
  {
    if (this->w == aw) return;
    delete [] this->data;
    this->w = aw;
    this->data = new T[H*this->w];
  }

  /// assign matrix, sizes must match
  template<typename TB>
  MatrixFixHeight & operator= (const Expr<TB> & m)
  {
    int k = 0;
    for (int j = 0; j < this->w; j++)
      for (int i = 0; i < H; i++)
	{
	  this->data[k] = m.Spec()(i,j);
	  k++;
	}
    return *this;
  }

  /// fill matrix with scalar
  MatrixFixHeight & operator= (TSCAL s)
  {
    FlatMatrixFixHeight<H,T>::operator= (s);
    return *this;
  }
};














/**
   Identity Matrix of fixed size
 */
template <int H>
class Id : public MatExpr<Id<H> >
{
public:
  typedef double TELEM;
  typedef double TSCAL;
  typedef Vec<H, double> TV_COL;
  typedef Vec<H, double> TV_ROW;
  enum { HEIGHT = H };
  enum { WIDTH  = H };
  // static bool IsLinear() { return 0; }
  enum { IS_LINEAR = 0 };

  /// do not initialize 
  Id () { ; }

  ///
  const TELEM operator() (int i) const
  { cerr << "id, linear access" << endl; return 0; }
  ///
  const TELEM operator() (int i, int j) const { return (i == j) ? 1 : 0; }

  /// the height
  int Height () const { return H; }
  /// the width
  int Width () const { return H; }
};







template<int H, int W, typename T>
inline std::ostream & operator<< (std::ostream & s, const Mat<H,W,T> & m)
{
  for (int i = 0; i < H*W; i++)
    s << " " << setw(7) << m(i);
  return s;
}



#endif