LUscalarMatrix: Added processor-local matrix inverse function

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;

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];
+        }
+    }
+}
+
+
 // ************************************************************************* //

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;
 };
 
 

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,

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
 (