diff --git a/SpMat.h b/SpMat.h
index 2bfec92a6a8684e8a3cdc13c3544779e9bf764a9..569106a8d78ac732173d73c27abb408406ab6ff6 100644
--- a/SpMat.h
+++ b/SpMat.h
@@ -112,12 +112,12 @@ public:
   void WarningsOff() {_pw=false;}
   bool IsSorted() const { return(is_sorted()); }  // Returns true if all _ri arrays are sorted. For debugging.
 
-
   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;}             // 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 SetColumn(unsigned int c, const NEWMAT::ColumnVector& col, double eps=0.0);  // Set a whole column (obliterating what was there before)
+  SpMat<T>& MultiplyColumns(const NEWMAT::Matrix& vals);                             // Multiply each column of a matrix by a value
   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) 
@@ -657,6 +657,41 @@ const SpMat<T> SpMat<T>::t() const
   return(t_mat);
 }
 
+/////////////////////////////////////////////////////////////////////
+//
+// Takes a row- or column-vector with as many elements as there are
+// columns in the SpMat and multiplies each column by those values.
+// It is the same as post-multiplication of SpMat matrix by a diagonal
+// matrix with the vals values on the diagonal.
+//
+/////////////////////////////////////////////////////////////////////
+
+template<class T>
+SpMat<T>& SpMat<T>::MultiplyColumns(const NEWMAT::Matrix& vals)
+{
+  if (vals.Nrows() != 1 && vals.Ncols() != 1) { // vals must be a row- or column-vector
+    throw SpMatException("MultiplyColumns: vals must be a row- or column-vector");
+  }
+  if (vals.Nrows() != _n && vals.Ncols() != _n) { // vals must have as many elements as there are columns
+    throw SpMatException("MultiplyColumns: vals must have as many elements as there are columns");
+  }
+  if (vals.Nrows() == 1) { // If it is a row-vector
+    for (unsigned int ci=0; ci<_n; ci++) {
+      for (unsigned int i=0; i<_val[ci].size(); i++) {
+	_val[ci][i] *= vals(1,ci+1);
+      }
+    }
+  }
+  else if (vals.Ncols() == 1) { // If it is a column-vector
+    for (unsigned int ci=0; ci<_n; ci++) {
+      for (unsigned int i=0; i<_val[ci].size(); i++) {
+	_val[ci][i] *= vals(ci+1,1);
+      }
+    }
+  } 
+  return(*this);
+}                                  
+
 /////////////////////////////////////////////////////////////////////
 //
 // Sets the values of an entire column, destroying any previous content.
@@ -1212,7 +1247,7 @@ SpMat<T>& SpMat<T>::add_diff_sparsity_mat_to_me(const SpMat<T>& M, double s)
       std::vector<unsigned int>&        ri = _ri[c];
       std::vector<T>&                   val = _val[c];
       for (unsigned int i=0; i<ri.size(); i++) {
-        acc(ri[i]) += s*val[i];
+        acc(ri[i]) += val[i];
       }
       ri.resize(acc.NO());
       val.resize(acc.NO());