diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index c3b701c8f9b6a827a46ba5f693b93ffaa962da47..58b196da091399339655c4e1b20b26a99420a833 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 d3818cfc235bc1d8e0d1d66ad66f11b646ca3f7a..0830b33ae67f49ff0b44d587db8cb66daa4f7061 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 f1ca94018bb68ff4b16c1538030d90f230787481..0489ed012d00f3dce16f9378e98209e4d6e98657 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];
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //