Commit 1e3d4eb0 authored by laurence's avatar laurence
Browse files

ENH: multiply(ans, A, B) for matrices templated on Matrix type. Previously...

ENH: multiply(ans, A, B) for matrices templated on Matrix type. Previously only worked for rectangular matrices.
parent 3b9cdac8
......@@ -136,41 +136,6 @@ void Foam::LUDecompose
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::multiply
(
scalarRectangularMatrix& ans, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B
)
{
if (A.m() != B.n())
{
FatalErrorIn
(
"multiply("
"scalarRectangularMatrix& answer "
"const scalarRectangularMatrix& A, "
"const scalarRectangularMatrix& B)"
) << "A and B must have identical inner dimensions but A.m = "
<< A.m() << " and B.n = " << B.n()
<< abort(FatalError);
}
ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
for (register label i = 0; i < A.n(); i++)
{
for (register label j = 0; j < B.m(); j++)
{
for (register label l = 0; l < B.n(); l++)
{
ans[i][j] += A[i][l]*B[l][j];
}
}
}
}
void Foam::multiply
(
scalarRectangularMatrix& ans, // value changed in return
......
......@@ -88,11 +88,12 @@ void LUBacksubstitute
template<class Type>
void LUsolve(scalarSquareMatrix& matrix, Field<Type>& source);
template<class Form, class Type>
void multiply
(
scalarRectangularMatrix& answer, // value changed in return
const scalarRectangularMatrix& A,
const scalarRectangularMatrix& B
Matrix<Form, Type>& answer, // value changed in return
const Matrix<Form, Type>& A,
const Matrix<Form, Type>& B
);
void multiply
......
......@@ -176,4 +176,40 @@ void Foam::LUsolve
}
template<class Form, class Type>
void Foam::multiply
(
Matrix<Form, Type>& ans, // value changed in return
const Matrix<Form, Type>& A,
const Matrix<Form, Type>& B
)
{
if (A.m() != B.n())
{
FatalErrorIn
(
"multiply("
"Matrix<Form, Type>& answer "
"const Matrix<Form, Type>& A, "
"const Matrix<Form, Type>& B)"
) << "A and B must have identical inner dimensions but A.m = "
<< A.m() << " and B.n = " << B.n()
<< abort(FatalError);
}
ans = Matrix<Form, Type>(A.n(), B.m(), scalar(0));
for (register label i = 0; i < A.n(); i++)
{
for (register label j = 0; j < B.m(); j++)
{
for (register label l = 0; l < B.n(); l++)
{
ans[i][j] += A[i][l]*B[l][j];
}
}
}
}
// ************************************************************************* //
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment