//  Triangular Matrix Class Template
//  --------------------------------

//  A triangular matrix is one in which all elements below the main
//  diagonal are 0. Therefore, they need not be stored.

#include "NumVec.hpp"

template<class T>            //  T must have a 0 value
class TMatrix : public NumVec<T> {
protected:                   //  Member data items (in addition to NumVec's)
  int TM_size;               //  Dimension -- needed only for efficiency
                             //               in subscript mapping
  T   dummy;                 //  Harmless target for any attempt to store
                             //     an element below the diagonal

public:                      //  Constructors:

//  Since the two dimensions are equal, we need only one set of bounds,
//  and, to simplify the example, we implement only the 0-origin version:

  TMatrix(INT dim)      : NumVec<T>(dim * (dim+1) / 2), TM_size(dim) {}
  TMatrix(TMatrix& x)   : NumVec<T>(x),           TM_size(x.TM_size) {}


//  Since the subscript operator [] allows only one parameter, we use
//  ordinary parentheses (the "function call" operator) for indexing:

  T  operator()(INT i, INT j) const   //  Const version
     {return i > j ?   0   : data[dim()+j-i-(TM_size-i)*(TM_size-i+1)/2];}

  T& operator()(INT i, INT j)         //  Lvalue version
     {dummy = 0;                      //  (in case of fetch)
      return i > j ? dummy : data[dim()+j-i-(TM_size-i)*(TM_size-i+1)/2];}

  private: T operator[](INT rs); T& operator[](INT rs) const; public:


//  Assignment operator

  TMatrix& operator=(const TMatrix& rs) //  From a matrix
  {assert(TM_size == rs.TM_size);       //  Verify compatible size
   NumVec<T>::operator=(rs);            //  Assign the data
   return *this;
  }

  TMatrix& operator=(const T  rs)       //  From a scalar
  {NumVec<T>::operator=(rs);
   return *this;
  }
};