From 1e3d4eb0676d7e58682e838559b6eb21fe133f0d Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Thu, 20 Dec 2012 12:00:14 +0000
Subject: [PATCH] ENH: multiply(ans, A, B) for matrices templated on Matrix
 type. Previously only worked for rectangular matrices.

---
 .../matrices/scalarMatrices/scalarMatrices.C  | 35 ------------------
 .../matrices/scalarMatrices/scalarMatrices.H  |  7 ++--
 .../scalarMatrices/scalarMatricesTemplates.C  | 36 +++++++++++++++++++
 3 files changed, 40 insertions(+), 38 deletions(-)

diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index c3b701c8f9b..58b196da091 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -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
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
index d3818cfc235..0830b33ae67 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
@@ -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
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index f1ca94018bb..0489ed012d0 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -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];
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
-- 
GitLab