#ifndef __sparsemat
#define __sparsemat

#include <vector> // STL vector
#include <list>   // STL list
#include <utility>   // STL pair
#include <iostream>
#include <iterator> // STL ostream_iterator

using namespace std;

template <class Field> class sparsemat; // forward declare template class,
                                        // which is used in << and >>

template <class Field> // must declare before friend declaration in
                       // template class sparsemat
ostream& operator<<(ostream& os, const sparsemat<Field>& A) {
  ostream_iterator< pair< int, Field > > out(os, " ");

  for (int i = 1; i <= A.m; i++) {
    os << "Row " << i << ": ";
    os.flush();
    copy(A.A[i].begin(), A.A[i].end(), out);
    os << endl;
  } // for (;;)

  return os;
}// end operator<< sparsemat<Field>

template <class Field>
istream& operator>>(istream& is, sparsemat<Field>& A) {
// reads the elements of a sparse matrix by reading in 
// a list of 3 items, first the row and column index and then a field element.
  int i; int j;
  Field a(0);

  while (is >> i) {
    // operator>> returns a reference to an istream
    // then istream::operator void*() is called which
    // returns !basic_ios::fail() [Stroutstrup, p.617, ???]
    if(i < 0) break; // return also if row index is negative
    is >> j >> a;
    A.put_value(make_pair(i,j), a);
  } // while-loop

  return is;
} // operator>> sparsemat<Field>

template <class Field>
class sparsemat {

// specialized templates (in this case, with a template parameter type)
friend ostream& operator<< <> (ostream&, const sparsemat<Field>&);
friend istream& operator>> <> (istream&, sparsemat<Field>&);

vector< list< pair< int, Field > > > A; // sparse matrix data struct
// A[i], the i-th row, is an STL list< pair< int, Field > >.
// Each list entry is a column index/value pair

int m; // rowdimension from 1..m; the actual dimensions
int n; // coldimension from 1..n; the actual dimensions

typedef list< pair< int, Field > > row_type;
typedef typename list< pair< int, Field > >::iterator row_iter;
typedef typename list< pair< int, Field > >::const_iterator const_row_iter;
// shorthands:
// can write template member functions with sparsemat<Field>::row_type
// and sparsemat<Field>::row_iter (in put_value, e.g,),
//     sparsemat<Field>::const_row_iter (in operator[])

class comp_w_col_ind {
public:
   bool operator() (const pair< int, Field >& entry, int col_in)
     {return entry.first < col_in; }
};
// used in lower_bound as function object

public:

sparsemat(int m_init, int n_init) :
// constructs a matrix of the given dimensions with all 0s
   A(m_init+1, list< pair< int, Field > >() ),
   // initialize each list in the vector as empty list
   m(m_init), n(n_init)
   // set the dimensions
   {}
// note: the copy constructor and operator= will work as intended
//       because of STL's container design

int get_rowdim() const {return m;}
int get_coldim() const {return n;}

Field operator[] (const pair<int, int>& ind) const;
// return:     a copy of the element A[ind.first, ind.second]
// example:    A[make_pair(3, 5)]
// exceptions: if the indices are out of range, an error is printed
//             and Field(0) is returned
// note:       unlike STL vector's operator[] this member does
//             not return a reference to the Field element in the matrix
//             it is therefore impossible to write
//             A[make_pair(3, 5)] = Field(1);
//             the reason is that zero entries have no memory where
//             one could place the right side field element

void put_value (const pair<int, int>& ind, const Field& a);
// sets A[ind.first, ind.second] = a
// example:    A.put_value(make_pair(3, 5), Field(-4))
// exceptions: if the indices are out of range, an error is printed
// note:       if a==Field(0) no cell will be inserted
//             and an already existing cell for the entry will be erased


sparsemat& linsolve(sparsemat& b);
// arguments:  b must be an m by 1 matrix
// return:     an n by 1 matrix x such that A * x = b
// exceptions: if the system does not have a unique solution,
//             an error is printed and a zero matrix returned
// See: www.math.ncsu.edu/~kaltofen/courses/LinAlgebra/Maple/refpkg/Src/ref.mpl
//      for a Maple row echelon form algorithm

void swaprow(int i, int j); // exchanges row i and j in A
// auxiliary member function to linsolve
// note:      uses list< <int, Field> >::swap

void addrow(int i, int j, const Field& a);
// auxiliary member function to linsolve
// adds a*(row i) to (row j) in A


// ALTERNATIVE to linsolve: matrix times vector product
sparsemat& lintrafo(sparsemat& b);
// arguments:  b must be an n by 1 matrix
// return:     an m by 1 matrix y such that y = A * b

};// end template class sparsemat

namespace std {// peculiarity of overloading
// << is called from a member function in std, which already
// has operator<< declared, so it will not overload outside the std namespace
template <class Field>
ostream& operator<<(ostream& os, const pair< int, Field >& entry)
// needed in sparsemat<Field>::print
{os << "(" << entry.first << ", " << entry.second << ")";
 return os;
}// end << pair< int, Field >
}// end namespace std
#endif