Index: /trunk/SophyaLib/TArray/triangmtx.h
===================================================================
--- /trunk/SophyaLib/TArray/triangmtx.h	(revision 2956)
+++ /trunk/SophyaLib/TArray/triangmtx.h	(revision 2957)
@@ -11,5 +11,15 @@
   \class SOPHYA::TriangularMatrix
   \ingroup TArray
-  Class for inferior triangular matrix (base class for the class Alm)
+  \brief Class for inferior triangular matrix (base class for the class Alm)
+  The inferior triangular matrix is represented in memory as column packed,
+  as illustrated below for a 5x5 triangular matrix.
+  \verbatim
+  5x5 Inf.Triang.Matrix, Size= 15 elements (0 ... 14) 
+  | 0                  |
+  | 1   5              |
+  | 2   6   9          |
+  | 3   7   10  12     |
+  | 4   8   11  13  14 |
+  \endverbatim 
 */
 
@@ -24,7 +34,7 @@
 TriangularMatrix()   {;};
 //! instanciate a triangular matrix from the number of rows
-TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) 
+TriangularMatrix(sa_size_t rowSize)  : long_diag_(rowSize) 
   {
-    elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); 
+    elem_.ReSize((rowSize*(rowSize+1)/2) ); 
   }
 //! Copy constructor (possibility of sharing datas)
@@ -32,5 +42,5 @@
 
 //! resize the matrix with a new number of rows
-inline void ReSizeRow(int rowSize) 
+inline void ReSizeRow(sa_size_t rowSize) 
   {
     long_diag_=(uint_4)rowSize;
@@ -47,10 +57,10 @@
 
 //! () operator : access to elements row \b l and column \b m
-inline T& operator()(int l, int m) 
+inline T& operator()(sa_size_t l, sa_size_t m) 
   {
       return  elem_(indexOfElement(l,m));
   }
 
-inline T& operator()(int index) 
+inline T& operator()(sa_size_t index) 
   {
       return  elem_(index);
@@ -59,10 +69,10 @@
 
 //! () operator : access to elements row \b l and column \b m
-inline T const& operator()(int l, int m) const 
+inline T const& operator()(sa_size_t l, sa_size_t m) const 
   {
       return *(elem_.Begin()+ indexOfElement(l,m));
   }
 
-inline T const& operator()(int index) const 
+inline T const& operator()(sa_size_t index) const 
   {
       return *(elem_.Begin()+ index);
@@ -91,5 +101,5 @@
       throw(SzMismatchError("TriangularMatrix<T>::CopyElt(const TriangularMatrix<T>&) SizeMismatch")) ;
     long_diag_ = a.long_diag_;
-    int k;
+    sa_size_t k;
     for (k=0; k< Size(); k++) elem_(k) = a.elem_(k);
     return(*this);
@@ -104,10 +114,10 @@
 
 //! Return number of rows
-inline  int_4  rowNumber() const {return (int_4)long_diag_;}
+inline  sa_size_t rowNumber() const {return (int_4)long_diag_;}
 
 //! Return size of the total array
-  inline int_4 Size() const {return elem_.Size();}
+  inline sa_size_t Size() const {return elem_.Size();}
 
-  inline bool CheckRelativeIndices(int_4 l, int_4 m) const 
+  inline bool CheckRelativeIndices(sa_size_t l, sa_size_t m) const 
    {
      if ( l < m ) 
@@ -117,5 +127,5 @@
      return true;
    }
-  inline bool CheckAbsoluteIndice(int_4 l, int_4 m) const 
+  inline bool CheckAbsoluteIndice(sa_size_t l, sa_size_t m) const 
    {
      if ( indexOfElement(l,m) >= elem_.Size() )
@@ -124,5 +134,5 @@
        }
    }
-  inline bool CheckAbsoluteIndice(int_4 ind) const 
+  inline bool CheckAbsoluteIndice(sa_size_t ind) const 
    {
      if ( ind >= elem_.Size() )
@@ -132,24 +142,33 @@
    }
 
-void Print(int nbLignes=0)
+  //! ASCII dump of the matrix (set nbLignes=-1) for dumping the complete matrix
+void Print(ostream& os, sa_size_t nbLignes=0) const
   {
-    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;
+    os << "TriangularMatrix< " << typeid(T).name() 
+       << " > NRow=" << long_diag_ << " NbElem<>0 : " << Size() << endl;
+    if (nbLignes == 0) return;
+    if (nbLignes < 0 ) nbLignes = long_diag_;
+    if (nbLignes > long_diag_ ) nbLignes = long_diag_;
+    for (sa_size_t k=0; k < nbLignes; k++)  {
+      os << "L[" << k << "]: " ;
+      for (sa_size_t kc = 0; kc <= k ; kc++) 
+	os << " " << elem_(indexOfElement(k,kc));
+      os << endl;
+    }
+    if (nbLignes < long_diag_)  os << " ... ... ... " << endl;
+    return;
   }
 
-  //Return pointer to first element address
-  //inline T* Data()  {return elem_.Begin();}
+inline void Print(sa_size_t nbLignes=0) const  { Print(cout, nbLignes); }
+
+
+//! Return the pointer to the first non zero element in column \b j = &(tmmtx(j,j)) 
+inline const T* columnData(sa_size_t j)  const {return elem_.Begin()+(long_diag_*j-j*(j-1)/2) ;}
+
+//! Return the pointer to the first non zero element in column \b j = &(tmmtx(j,j)) 
+inline T* columnData(sa_size_t j) {return elem_.Begin()+(long_diag_*j-j*(j-1)/2) ;}
 
 //! compute the address of an element in the single array representing the matrix
-inline uint_4 indexOfElement(int i,int j) const 
+inline sa_size_t indexOfElement(sa_size_t i,sa_size_t j) const 
 {
   //  return(i*(i+1)/2+j);
@@ -160,8 +179,12 @@
 private: 
 
-uint_4 long_diag_;    //!< size of the square matrix
+sa_size_t long_diag_;    //!< size of the square matrix
 NDataBlock<T> elem_;  //!< Data block
 
   };
+
+template <class T>
+inline ostream& operator << (ostream& os, const TriangularMatrix<T>& a)
+                            { a.Print(os, 0);    return(os);    }
   
 }   // namespace SOPHYA
