From 9b6b65fe0494717dc956c8e4ad63681bb12f0512 Mon Sep 17 00:00:00 2001
From: Victor Brunini <victorbrunini@gmail.com>
Date: Mon, 7 Jan 2013 18:18:10 +0000
Subject: [PATCH] Add function to DenseMatrix that multiplies two square
 matrices.

---
 include/cantera/numerics/DenseMatrix.h |  8 ++++++++
 src/numerics/DenseMatrix.cpp           | 16 ++++++++++++++++
 2 files changed, 24 insertions(+)

diff --git a/include/cantera/numerics/DenseMatrix.h b/include/cantera/numerics/DenseMatrix.h
index ce9030478..8276b9e24 100644
--- a/include/cantera/numerics/DenseMatrix.h
+++ b/include/cantera/numerics/DenseMatrix.h
@@ -138,6 +138,14 @@ public:
      */
     virtual void mult(const double* b, double* prod) const;
 
+    //! Multiply A*B and write result to \c prod.
+    /*!
+     *
+     *  @param    b     input      DenseMatrix B of size NxN
+     *  @param    prod  output     output DenseMatrix prod size NxN
+     */
+    virtual void mult(const DenseMatrix &b, DenseMatrix &prod) const;
+
     //! Left-multiply the matrix by transpose(b), and write the result to prod.
     /*!
      *   @param b    left multiply by this vector. The length must be equal to n
diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp
index 9e09da9ea..85b1ccf9b 100644
--- a/src/numerics/DenseMatrix.cpp
+++ b/src/numerics/DenseMatrix.cpp
@@ -113,6 +113,22 @@ void DenseMatrix::mult(const double* b, double* prod) const
              static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
 }
 //====================================================================================================================
+void DenseMatrix::mult(const DenseMatrix &B, DenseMatrix &prod) const
+{
+  if(m_ncols != B.nColumns() || m_nrows != B.nRows() || m_ncols != m_nrows || m_ncols != prod.nColumns() || m_nrows != prod.nColumns() )
+  {
+    throw CanteraError("mult(const DenseMatrix &B, DenseMatrix &prod)",
+        "Cannot multiply matrices that are not square and/or not the same size.");
+  }
+  const doublereal * const *bcols = B.const_colPts();
+  doublereal * const *prodcols = prod.colPts();
+  for(size_t col=0; col < m_ncols; ++col)
+  {
+    // Loop over ncols multiplying A*column of B and storing in corresponding prod column
+    mult(bcols[col], prodcols[col]);
+  }
+}
+//====================================================================================================================
 void DenseMatrix::leftMult(const double* const b, double* const prod) const
 {
     size_t nc = nColumns();