diff --git a/SpMat.h b/SpMat.h
index 625dc916105af0994e9285f607f7edcf1234e68d..ded38742a83bd4ad988ba7532a4064538010a3d3 100644
--- a/SpMat.h
+++ b/SpMat.h
@@ -89,8 +89,41 @@ template<class T>
 class SpMat
 {
 public:
-  SpMat() : _m(0), _n(0), _nz(0), _ri(0), _val(0), _pw(false) {}
-  SpMat(unsigned int m, unsigned int n) : _m(m), _n(n), _nz(0), _ri(n), _val(n), _pw(false) {}
+
+    class Iterator {
+    public:
+    Iterator(SpMat<T>& mat, bool oob=false) : _mat(mat), _i(0), _oob(oob)
+      {
+	_j = 0;
+	while (_j < _mat._n && !_mat._ri[_j].size()) _j++;
+	if (_j == _mat._n) _oob = true;
+      }
+      ~Iterator() {}
+      bool operator==(const Iterator& other) { return(&_mat==&other._mat && ((_oob && other._oob) || (_i==other._i && _j==other._j))); }
+      bool operator!=(const Iterator& other) { return(!(*this==other)); }
+      T& operator*() { return((_mat._val[_j])[_i]); }
+      Iterator& operator++() 
+      {
+	if (++_i < _mat._ri[_j].size()) return(*this);
+        else {
+          while (++_j < _mat._n && !_mat._ri[_j].size()) ;
+        }
+        if (_j == _mat._n) _oob = true;
+        else _i = 0;
+        return(*this); 
+      }
+      unsigned int Row() { return((_mat._ri[_j])[_i]+1); }
+      unsigned int Col() { return(_j+1); }
+    private:
+      SpMat<T>&     _mat;
+      unsigned int  _i;
+      unsigned int  _j;
+      bool          _oob;
+    };
+
+
+  SpMat() : _m(0), _n(0), _nz(0), _ri(0), _val(0), _pw(false), _ei(*this,true) {}
+  SpMat(unsigned int m, unsigned int n) : _m(m), _n(n), _nz(0), _ri(n), _val(n), _pw(false), _ei(*this,true) {}
   SpMat(unsigned int m, unsigned int n, const unsigned int *irp, const unsigned int *jcp, const double *sp);
   SpMat(const NEWMAT::GeneralMatrix& M);
   SpMat(const std::string& fname);
@@ -100,6 +133,9 @@ public:
   unsigned int Ncols() const {return(_n);}
   unsigned int NZ() const {return(_nz);}
 
+  Iterator         begin() { return(Iterator(*this)); }
+  const Iterator&  end() { return(_ei); }
+
   NEWMAT::ReturnMatrix AsNEWMAT() const;
   void Save(const std::string&   fname,
             unsigned int         precision) const;
@@ -112,7 +148,6 @@ public:
   void WarningsOn() {_pw=true;}
   void WarningsOff() {_pw=false;}
 
-
   T Peek(unsigned int r, unsigned int c) const;
   T operator()(unsigned int r, unsigned int c) const {return(Peek(r,c));}    // Read-only
 
@@ -177,6 +212,7 @@ private:
   std::vector<std::vector<unsigned int> >   _ri;
   std::vector<std::vector<T> >              _val;
   bool                                      _pw;   // Print Warnings
+  Iterator                                  _ei;
 
   bool found(const std::vector<unsigned int>&  ri, unsigned int key, int& pos) const;
   T& here(unsigned int r, unsigned int c);
@@ -309,7 +345,7 @@ private:
 
 template<class T>
 SpMat<T>::SpMat(unsigned int m, unsigned int n, const unsigned int *irp, const unsigned int *jcp, const double *sp)
-  : _m(m), _n(n), _nz(0), _ri(n), _val(n), _pw(false)
+  : _m(m), _n(n), _nz(0), _ri(n), _val(n), _pw(false), _ei(*this,true)
 {
   _nz = jcp[n];
   unsigned long nz = 0;
@@ -339,7 +375,7 @@ SpMat<T>::SpMat(unsigned int m, unsigned int n, const unsigned int *irp, const u
 
 template<class T>
 SpMat<T>::SpMat(const NEWMAT::GeneralMatrix& M)
-  : _m(M.Nrows()), _n(M.Ncols()), _nz(0), _ri(M.Ncols()), _val(M.Ncols()), _pw(false)
+  : _m(M.Nrows()), _n(M.Ncols()), _nz(0), _ri(M.Ncols()), _val(M.Ncols()), _pw(false), _ei(*this,true)
 {
   double *m = static_cast<double *>(M.Store());
 
@@ -375,7 +411,7 @@ SpMat<T>::SpMat(const NEWMAT::GeneralMatrix& M)
 
 template<class T>
 SpMat<T>::SpMat(const std::string&  fname)
-: _m(0), _n(0), _nz(0), _ri(0), _val(0), _pw(false)
+  : _m(0), _n(0), _nz(0), _ri(0), _val(0), _pw(false), _ei(*this,true)
 {
   // First read data into (nz+1)x3 NEWMAT matrix
   NEWMAT::Matrix rcv;