#include "spareM.h"


#include <iostream.h>


#include "../Basics/exception.h"



#define _CHECK_INDEX



// constructors

tcT /*inline*/ spareM<T>::spareM () : 
   coefPtr(0), rowPtr(0), N(0), val(), ind(), beg(), cnt(), finger() {}


tcT /*inline*/ spareM<T>::spareM (const int m, const int nbMaxCoef) :
   val(nbMaxCoef), ind(nbMaxCoef), beg(1,m+1,0), cnt(1,m,0), finger(1,m,0), 
   coefPtr(0), rowPtr(0), N(0) {}


tcT void spareM<T>::newCoef (const int i, const int j, const T v)
{
#ifdef _CHECK_INDEX
   if ( i<0 || i>m() )
      throw whereOutOfRange(
	 "void spareM<T>::newCoef (const int i, const int j, const T v)",
	 "i",i,0,m()-1);
   if ( j<0 )
      throw whereOutOfRange(
	 "void spareM<T>::newCoef (const int i, const int j, const T v)",
	 "j",j,0,val.n()-1);
   if ( coefPtr>=val.n() )
      throw whereWhich(
	 "void spareM<T>::newCoef (const int i, const int j, const T v)",
	 "More non zero coefficient than expected");
   if ( i<rowPtr )
      throw whereWhich(
	 "void spareM<T>::newCoef (const int i, const int j, const T v)",
	 "Rows must be introduced in incremental order");
#endif
   rowPtr=i;
   beg.s(0,1,rowPtr+1) = beg(rowPtr+1)+1;
   cnt(rowPtr)++;
   ind(coefPtr)=j;
   val(coefPtr)=v;
   coefPtr++;
   if ( j >= N ) N=j+1;
}


// selectors

tcT /*inline*/ int spareM<T>::m() const { return cnt.n(); }

tcT /*inline*/ int spareM<T>::n() const { return N; }


// indexation

tcT T spareM<T>::operator () (const int i, const int j) const
{
   int mini = beg(i);
   int maxi = beg(i+1)-1;
   if ( maxi<mini || ind(mini)>j || ind(maxi)<j ) return T(0);
   if ( ind(mini)==j ) return val(mini);
   if ( ind(maxi)==j ) return val(maxi);
   int midi;
   while ( mini+1<maxi )
   {
      midi=(mini+maxi)/2;
      if ( ind(midi)==j ) return val(midi);
      if ( ind(midi)<j ) mini=midi;
      else               maxi=midi;
   }
   if ( ind(mini)==j ) return val(mini);
   if ( ind(maxi)==j ) return val(maxi);
   return T(0);
}


// rows traversal

tcT boolean spareM<T>::first (const int i, int& j, T& v)
{
   if ( !cnt(i) ) return false;
   finger(i)=beg(i);
   j=ind(finger(i));
   v=val(finger(i));
   return true;
}


tcT boolean spareM<T>::next (const int i, int& j, T& v)
{
   if ( (++finger(i))>=beg(i+1) ) return false;
   j=ind(finger(i)); // 
   v=val(finger(i));
   return true;
}
   

// output of a spare Matrix

tcT ostream& operator << (ostream& s, const spareM<T>& a)
{
   s << "beg " << a.beg << "cnt " << a.cnt;
   for (int ii=0; ii<a.m(); ii++)
   { 
      for (int jj=0; jj<a.n(); jj++)
	 s << setw(6) << int(a(ii,jj));
      s<<endl<<flush;
   }
   return s;
}