Added constructor that constructs Matrix from row col val format filename

SpMat.h

@@ -27,6 +27,7 @@
 #include "newmat.h"
 #include "cg.h"
 #include "bicg.h"
+#include "miscmaths.h"
 
 namespace MISCMATHS {
 
@@ -92,6 +93,7 @@ public:
   SpMat(unsigned int m, unsigned int n) : _m(m), _n(n), _nz(0), _ri(n), _val(n) {}
   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);
   ~SpMat() {}
 
   unsigned int Nrows() const {return(_m);}
@@ -360,6 +362,64 @@ SpMat<T>::SpMat(const NEWMAT::GeneralMatrix& M)
   }  
 }
 
+/////////////////////////////////////////////////////////////////////
+//
+// Constructs matrix from row col val format produced by 
+// Save/Print below.
+//
+/////////////////////////////////////////////////////////////////////
+
+template<class T>
+SpMat<T>::SpMat(const std::string&  fname)
+: _m(0), _n(0), _nz(0), _ri(0), _val(0)
+{
+  // First read data into (nz+1)x3 NEWMAT matrix
+  NEWMAT::Matrix rcv;
+  try {
+    rcv = read_ascii_matrix(fname);
+  }
+  catch(...) {
+    throw SpMatException("SpMat::SpMat(string& fname): cannot read file given by fname");
+  }
+  // Then interpret it
+  if (rcv(rcv.Nrows(),3)) throw SpMatException("SpMat::SpMat(string& fname): Last row must have zero value and indicate matrix size");
+  _m = static_cast<unsigned int>(rcv(rcv.Nrows(),1)+0.5);
+  _n = static_cast<unsigned int>(rcv(rcv.Nrows(),2)+0.5);
+  // cout << "rcv = " << endl << rcv << endl << "_n = " << _n << endl;
+  _ri.resize(_n);
+  _val.resize(_n);
+  // First pass to see how many elements in each colum
+  std::vector<unsigned int> col_count(_n,0);
+  unsigned int col = static_cast<unsigned int>(rcv(1,2)+0.5);
+  for (unsigned int indx=1; indx<static_cast<unsigned int>(rcv.Nrows()); indx++) {
+    if (static_cast<unsigned int>(rcv(indx,2)+0.5) != col) {
+      if (static_cast<unsigned int>(rcv(indx,2)+0.5) < col) throw SpMatException("SpMat::SpMat(string& fname): Column index must be monotonously increasing");
+      else col = static_cast<unsigned int>(rcv(indx,2)+0.5);
+      if (col > _n) throw SpMatException("SpMat::SpMat(string& fname): File internally inconsistent");
+    }
+    // cout << "col = " << col << endl;
+    col_count[col-1]++;
+  }
+  // Second pass to allocate and fill vectors
+  unsigned int indx=1;
+  for (col=0; col<_n; col++) {
+    std::vector<unsigned int>& ri = _ri[col];
+    std::vector<T>&            val = _val[col];
+    ri.resize(col_count[col]);
+    val.resize(col_count[col]);
+    for (unsigned int i=0; i<col_count[col]; i++, indx++) {
+      if (i && ri[i] >= static_cast<unsigned int>(rcv(indx,1)+0.5)) throw SpMatException("SpMat::SpMat(string& fname): Row index must be monotonously increasing");
+      if (static_cast<unsigned int>(rcv(indx,1)+0.5) < 1 || static_cast<unsigned int>(rcv(indx,1)+0.5) > _m) {
+        throw SpMatException("SpMat::SpMat(string& fname): Row index outside 1 -- -m range");
+      } 
+      ri[i] = static_cast<unsigned int>(rcv(indx,1)+0.5) - 1;
+      val[i] = rcv(indx,3);
+      _nz++;
+    }
+  }
+}
+
+
 /////////////////////////////////////////////////////////////////////
 //
 // Returns matrix in NEWMAT matrix format. Useful for debugging