Added SetColumn() to SpMat

SpMat.h

@@ -101,8 +101,9 @@ public:
   T Peek(unsigned int r, unsigned int c) const;
   T operator()(unsigned int r, unsigned int c) const {return(Peek(r,c));}    // Read-only
 
-  void Set(unsigned int r, unsigned int c, const T& v) {here(r,c) = v;}
-  void AddTo(unsigned int r, unsigned int c, const T& v) {here(r,c) += v;}
+  void Set(unsigned int r, unsigned int c, const T& v) {here(r,c) = v;}             // Set a single value
+  void SetColumn(unsigned int c, const NEWMAT::ColumnVector& col, double eps=0.0);  // Set a whole column (obliterating what was there before) 
+  void AddTo(unsigned int r, unsigned int c, const T& v) {here(r,c) += v;}          // Add value to a single (possibly existing) value
 
   SpMat<T>& operator+=(const SpMat& M) 
   {
@@ -243,9 +244,9 @@ private:
 //
 // The concept of an accumulator was "borrowed" from Gilbert et al. 
 // 92. It is intended as a helper class for SpMat and is used to
-// hold the content of one column of a matrix C where C is a 
-// sparse matrix resulting from C=A*B where A and B are sparse
-// matrices.
+// hold the content of one column of a matrix. This column can then
+// be accessed both by indexing a certain element, and also by indexing
+// only non-zero elements.
 //
 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
@@ -488,6 +489,45 @@ const SpMat<T> SpMat<T>::t() const
   return(t_mat);  
 }
 
+/////////////////////////////////////////////////////////////////////
+//
+// Sets the values of an entire column, destroying any previous content.
+//
+/////////////////////////////////////////////////////////////////////
+
+template<class T>
+void SpMat<T>::SetColumn(unsigned int                 c,      // Column #
+			 const NEWMAT::ColumnVector&  col,    // The values in that column
+			 double                       eps)    // Any value <= eps is treated as a zero
+{
+  if (c < 1 || c > _n) throw SpMatException("SetColumn: column index out of range");
+  if (static_cast<unsigned int>(col.Nrows()) != _m) throw SpMatException("SetColumn: column size mismatch");
+
+  Accumulator<T>    acc(_m);
+  double            *colp = col.Store();
+  
+  for (unsigned int i=0; i<_m; i++) {
+    if (colp[i] > eps) acc(i) = static_cast<T>(colp[i]);
+  }
+  
+  std::vector<unsigned int>&   ri = _ri[c-1];
+  std::vector<T>&              val = _val[c-1];
+  unsigned int                 old_sz = ri.size();
+  if (old_sz) {
+    ri = std::vector<unsigned int>(acc.NO());
+    val = std::vector<T>(acc.NO());
+  }
+  else {
+    ri.resize(acc.NO()); 
+    val.resize(acc.NO());
+  }
+  for (unsigned int i=0; i<acc.NO(); i++) {
+    ri[i] = acc.ri(i);
+    val[i] = acc.val(i);
+  }
+  _nz += (acc.NO() - old_sz);
+}
+
 /////////////////////////////////////////////////////////////////////
 //
 // Returns value at position i,j (one offset)