research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Numerics] Fix DenseMatrix::mult to work for rectangular matrices

Browse source
This commit is contained in:
Ray Speth 2016-04-09 00:48:45 -04:00
parent b94e272291
commit bd24aeb926
3 changed files with 41 additions and 6 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
include/cantera/numerics/DenseMatrix.h
View file

@ -112,8 +112,9 @@ public:
//! Multiply A*B and write result to \c prod.
/*!
* @param[in] b DenseMatrix B of size NxN
* @param[out] prod DenseMatrix prod size NxN
* Take this matrix to be of size NxM.
* @param[in] b DenseMatrix B of size MxP
* @param[out] prod DenseMatrix prod size NxP
*/
virtual void mult(const DenseMatrix& b, DenseMatrix& prod) const;

13
src/numerics/DenseMatrix.cpp
View file

@ -93,13 +93,18 @@ void DenseMatrix::mult(const double* b, double* prod) const
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.");
if (nColumns() != B.nRows()) {
throw CanteraError("DenseMatrix::mult",
"Inner matrix dimensions do not agree: {} != {}", nColumns(), B.nRows());
}
if (nRows() != prod.nRows() || B.nColumns() != prod.nColumns()) {
throw CanteraError("DenseMatrix::mult",
"Output matrix has wrong dimensions: {}x{} != {}x{}",
prod.nRows(), prod.nColumns(), nRows(), B.nColumns());
}
const doublereal* const* bcols = B.const_colPts();
doublereal* const* prodcols = prod.colPts();
for (size_t col=0; col < m_ncols; ++col) {
for (size_t col=0; col < B.nColumns(); ++col) {
// Loop over ncols multiplying A*column of B and storing in
// corresponding prod column
mult(bcols[col], prodcols[col]);

29
test/general/test_matrices.cpp
View file

@ -129,6 +129,16 @@ public:
A1(3,3) = -9;
}
double special_sum(DenseMatrix M) {
double sum = 0;
for (size_t i = 0; i < M.nRows(); i++) {
for (size_t j = 0; j < M.nColumns(); j++) {
sum += M(i,j) * (i+2*j+1);
}
}
return sum;
}
DenseMatrix A1, A2, A3;
vector_fp x4, x3;
vector_fp b1, b2, b3;
@ -150,3 +160,22 @@ TEST_F(DenseMatrixTest, matrix_times_vector)
EXPECT_DOUBLE_EQ(b3[i], c[i]);
}
}
TEST_F(DenseMatrixTest, matrix_times_matrix)
{
DenseMatrix c(4, 3);
A1.mult(A2, c);
EXPECT_DOUBLE_EQ(3033, special_sum(c));
c.resize(3, 4);
A3.mult(A1, c);
EXPECT_DOUBLE_EQ(3386, special_sum(c));
c.resize(3, 3);
A3.mult(A2, c);
EXPECT_DOUBLE_EQ(2989, special_sum(c));
c.resize(4, 4);
A2.mult(A3, c);
EXPECT_DOUBLE_EQ(4014, special_sum(c));
}