// This may look like C code, but it is really -*- C++ -*-

#ifndef TRIANGMTX_H_SEEN
#define TRIANGMTX_H_SEEN

#include "ndatablock.h"
#include "pexceptions.h"

// doit etre mis en dehors du namespace
/*!
  \class SOPHYA::TriangularMatrix
  \ingroup TArray
  Class for inferior triangular matrix (base class for the class Alm)
*/

namespace SOPHYA {
  
//! Class for inferior triangular matrix (base class for the class Alm)
template <class T>
class TriangularMatrix {
public :

//! Default constructor
TriangularMatrix()   {};
//! instanciate a triangular matrix from the number of rows
TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) {elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); };
//! Copy constructor (possibility of sharing datas)
TriangularMatrix(const TriangularMatrix<T>& a,  bool share=false)  : elem_(a.elem_, share),  long_diag_(a.long_diag_) {;}

//! resize the matrix with a new number of rows
inline void ReSizeRow(int rowSize) 
  {
    long_diag_=(uint_4)rowSize;
    elem_.ReSize(long_diag_*(long_diag_+1)/2);
  }

TriangularMatrix<T>& SetT(T a)
  {
    if (long_diag_ < 1)
    throw RangeCheckError("TriangularMatrix<T>::SetT(T )  - TriangularMatrix not dimensionned ! ");
    elem_ = a;
    return (*this);
  }

//! () operator : access to elements row \b l and column \b m
inline T& operator()(int l, int m) 
  {
      return  elem_(indexOfElement(l,m));
  }

inline T& operator()(int index) 
  {
      return  elem_(index);
  }


//! () operator : access to elements row \b l and column \b m
inline T const& operator()(int l, int m) const 
  {
      return *(elem_.Begin()+ indexOfElement(l,m));
  }

inline T const& operator()(int index) const 
  {
      return *(elem_.Begin()+ index);
  }


//! Return number of rows
inline  int_4  rowNumber() const {return (int_4)long_diag_;}

void Print(int nbLignes=0)
  {
    if (nbLignes == 0 ) nbLignes = long_diag_;
    cout << " ***** matrice triangulaire : ********* " << endl;
    for (int k=0; k < nbLignes; k++)
      {
	for (int kc = 0; kc <= k ; kc++)
	  {
	    cout << " " << elem_(indexOfElement(k,kc));
	  }
	cout << endl;
      }
    cout << "---------------- fin matrice ------------" << endl;
  }

  //Return pointer to first element address
  //inline T* Data()  {return elem_.Begin();}

//! compute the address of an element in the single array representing the matrix
inline uint_4 indexOfElement(int i,int j) const 
{
  //  return(i*(i+1)/2+j);
  // the (inferior triangular )matrix is stored column by column
  return(i+ long_diag_*j-j*(j+1)/2);
}

private: 

uint_4 long_diag_;    //!< size of the square matrix
NDataBlock<T> elem_;  //!< Data block

  };
  
}   // namespace SOPHYA

#endif