From 857915dabd27cce56ee43c01fec218862bcb7d2a Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sun, 17 Jul 2016 14:44:50 +0100
Subject: [PATCH] LUscalarMatrix: Added processor-local matrix inverse function

---
 applications/test/Matrix/Test-Matrix.C        |  7 +++
 .../matrices/LUscalarMatrix/LUscalarMatrix.C  | 17 ++++++++
 .../matrices/LUscalarMatrix/LUscalarMatrix.H  |  3 ++
 .../matrices/scalarMatrices/scalarMatrices.C  | 43 ++++++++++++++++++-
 .../matrices/scalarMatrices/scalarMatrices.H  |  8 ++++
 5 files changed, 77 insertions(+), 1 deletion(-)

diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index 523fd7971..8a0e0a3cf 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -152,6 +152,11 @@ int main(int argc, char *argv[])
         LUscalarMatrix LU(squareMatrix);
         scalarField x(LU.solve(source));
         Info<< "LU solve residual " << (squareMatrix*x - source) << endl;
+
+        scalarSquareMatrix inv(3);
+        LU.inv(inv);
+        Info<< "LU inv " << inv << endl;
+        Info<< "LU inv*squareMatrix " << (inv*squareMatrix) << endl;
     }
 
     {
@@ -169,6 +174,8 @@ int main(int argc, char *argv[])
 
         Info<< "QR inverse solve residual "
             << (x - QR.inv()*source) << endl;
+
+        Info<< "QR inv *squareMatrix " << (QR.inv()*squareMatrix) << endl;
     }
 
     Info<< "\nEnd\n" << endl;
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
index ecff7ea54..6c9f0eff7 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
@@ -412,4 +412,21 @@ void Foam::LUscalarMatrix::decompose(const scalarSquareMatrix& M)
 }
 
 
+void Foam::LUscalarMatrix::inv(scalarSquareMatrix& M) const
+{
+    scalarField source(m());
+
+    for (label j=0; j<m(); j++)
+    {
+        source = Zero;
+        source[j] = 1;
+        LUBacksubstitute(*this, pivotIndices_, source);
+        for (label i=0; i<m(); i++)
+        {
+            M(i, j) = source[i];
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index 0e0bc0101..f0948ff14 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
@@ -126,6 +126,9 @@ public:
         //  returning the solution
         template<class Type>
         tmp<Field<Type>> solve(const Field<Type>& source) const;
+
+        //- Set M to the inverse of this square matrix
+        void inv(scalarSquareMatrix& M) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index d96ffb903..0eb52a1be 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -269,7 +269,7 @@ void Foam::multiply
 
     for (label i=0; i<A.m(); i++)
     {
-        for (label g = 0; g < C.n(); g++)
+        for (label g=0; g<C.n(); g++)
         {
             for (label l=0; l<C.m(); l++)
             {
@@ -280,6 +280,47 @@ void Foam::multiply
 }
 
 
+void Foam::multiply
+(
+    scalarSquareMatrix& ans,         // value changed in return
+    const scalarSquareMatrix& A,
+    const DiagonalMatrix<scalar>& B,
+    const scalarSquareMatrix& C
+)
+{
+    if (A.m() != B.size())
+    {
+        FatalErrorInFunction
+            << "A and B must have identical dimensions but A.m = "
+            << A.m() << " and B.m = " << B.size()
+            << abort(FatalError);
+    }
+
+    if (B.size() != C.m())
+    {
+        FatalErrorInFunction
+            << "B and C must have identical dimensions but B.m = "
+            << B.size() << " and C.m = " << C.m()
+            << abort(FatalError);
+    }
+
+    const label size = A.m();
+
+    ans = scalarSquareMatrix(size, Zero);
+
+    for (label i=0; i<size; i++)
+    {
+        for (label g=0; g<size; g++)
+        {
+            for (label l=0; l<size; l++)
+            {
+                ans(i, g) += C(l, g)*A(i, l)*B[l];
+            }
+        }
+    }
+}
+
+
 Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
 (
     const scalarRectangularMatrix& A,
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
index 8cf58611d..214f731c6 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H
@@ -145,6 +145,14 @@ void multiply
     const scalarRectangularMatrix& C
 );
 
+void multiply
+(
+    scalarSquareMatrix& answer,         // value changed in return
+    const scalarSquareMatrix& A,
+    const DiagonalMatrix<scalar>& B,
+    const scalarSquareMatrix& C
+);
+
 //- Return the inverse of matrix A using SVD
 scalarRectangularMatrix SVDinv
 (