diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C
index 87a39386fc..908c364d9c 100644
--- a/applications/test/Matrix/Test-Matrix.C
+++ b/applications/test/Matrix/Test-Matrix.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "scalarMatrices.H"
+#include "LUscalarMatrix.H"
 #include "LLTMatrix.H"
 #include "QRMatrix.H"
 #include "vector.H"
@@ -113,70 +114,53 @@ int main(int argc, char *argv[])
         Info<< "Solution = " << rhs << endl;
     }
 
+
+    scalarSquareMatrix squareMatrix(3, Zero);
+
+    squareMatrix(0, 0) = 4;
+    squareMatrix(0, 1) = 12;
+    squareMatrix(0, 2) = -16;
+    squareMatrix(1, 0) = 12;
+    squareMatrix(1, 1) = 37;
+    squareMatrix(1, 2) = -43;
+    squareMatrix(2, 0) = -16;
+    squareMatrix(2, 1) = -43;
+    squareMatrix(2, 2) = 98;
+
+    Info<< nl << "Square Matrix = " << squareMatrix << endl;
+
+    const scalarField source(3, 1);
+
     {
-        scalarSquareMatrix squareMatrix(3, Zero);
-
-        squareMatrix(0, 0) = 4;
-        squareMatrix(0, 1) = 12;
-        squareMatrix(0, 2) = -16;
-        squareMatrix(1, 0) = 12;
-        squareMatrix(1, 1) = 37;
-        squareMatrix(1, 2) = -43;
-        squareMatrix(2, 0) = -16;
-        squareMatrix(2, 1) = -43;
-        squareMatrix(2, 2) = 98;
-
-        const scalarSquareMatrix squareMatrixCopy = squareMatrix;
-        Info<< nl << "Square Matrix = " << squareMatrix << endl;
-
-        Info<< "det = " << det(squareMatrixCopy) << endl;
+        {
+            scalarSquareMatrix sm(squareMatrix);
+            Info<< "det = " << det(sm) << endl;
+        }
 
+        scalarSquareMatrix sm(squareMatrix);
         labelList rhs(3, 0);
         label sign;
-        LUDecompose(squareMatrix, rhs, sign);
+        LUDecompose(sm, rhs, sign);
 
-        Info<< "Decomposition = " << squareMatrix << endl;
+        Info<< "Decomposition = " << sm << endl;
         Info<< "Pivots = " << rhs << endl;
         Info<< "Sign = " << sign << endl;
-        Info<< "det = " << detDecomposed(squareMatrix, sign) << endl;
+        Info<< "det = " << detDecomposed(sm, sign) << endl;
     }
 
     {
-        scalarSquareMatrix squareMatrix(3, Zero);
-
-        squareMatrix(0, 0) = 4;
-        squareMatrix(0, 1) = 12;
-        squareMatrix(0, 2) = -16;
-        squareMatrix(1, 0) = 12;
-        squareMatrix(1, 1) = 37;
-        squareMatrix(1, 2) = -43;
-        squareMatrix(2, 0) = -16;
-        squareMatrix(2, 1) = -43;
-        squareMatrix(2, 2) = 98;
-
-        scalarField source(3, 1);
+        LUscalarMatrix LU(squareMatrix);
+        scalarField x((LU.solve(source));
+        Info<< "LU solve residual " << (squareMatrix*x - source) << endl;
+    }
 
+    {
         LLTMatrix<scalar> LLT(squareMatrix);
         scalarField x(LLT.solve(source));
-
         Info<< "LLT solve residual " << (squareMatrix*x - source) << endl;
     }
 
     {
-        scalarSquareMatrix squareMatrix(3, Zero);
-
-        squareMatrix(0, 0) = 4;
-        squareMatrix(0, 1) = 12;
-        squareMatrix(0, 2) = -16;
-        squareMatrix(1, 0) = 12;
-        squareMatrix(1, 1) = 37;
-        squareMatrix(1, 2) = -43;
-        squareMatrix(2, 0) = -16;
-        squareMatrix(2, 1) = -43;
-        squareMatrix(2, 2) = 98;
-
-        scalarField source(3, 1);
-
         QRMatrix<scalarSquareMatrix> QR(squareMatrix);
         scalarField x(QR.solve(source));
 
@@ -184,8 +168,7 @@ int main(int argc, char *argv[])
             << (squareMatrix*x - source) << endl;
 
         Info<< "QR inverse solve residual "
-            << (x - QR.inverse()*source) << endl;
-
+            << (x - QR.inv()*source) << endl;
     }
 
     Info<< "\nEnd\n" << endl;
diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
index 2e8f7421d9..d0f828d4f3 100644
--- a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
+++ b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
@@ -95,6 +95,12 @@ void Foam::LLTMatrix<Type>::solve
     const Field<Type>& source
 ) const
 {
+    // If x and source are different initialize x = source
+    if (&x != &source)
+    {
+        x = source;
+    }
+
     const SquareMatrix<Type>& LLT = *this;
     const label m = LLT.m();
 
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
index c8a83e5e5a..0e0bc01013 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.H
@@ -116,10 +116,16 @@ public:
         //- Perform the LU decomposition of the matrix M
         void decompose(const scalarSquareMatrix& M);
 
-        //- Solve the matrix using the LU decomposition with pivoting
-        //  returning the solution in the source
-        template<class T>
-        void solve(Field<T>& source) const;
+        //- Solve the linear system with the given source
+        //  and returning the solution in the Field argument x.
+        //  This function may be called with the same field for x and source.
+        template<class Type>
+        void solve(Field<Type>& x, const Field<Type>& source) const;
+
+        //- Solve the linear system with the given source
+        //  returning the solution
+        template<class Type>
+        tmp<Field<Type>> solve(const Field<Type>& source) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
index 4d327faeee..605babb852 100644
--- a/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
+++ b/src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
@@ -24,23 +24,34 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LUscalarMatrix.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
+void Foam::LUscalarMatrix::solve
+(
+    Field<Type>& x,
+    const Field<Type>& source
+) const
 {
+    // If x and source are different initialize x = source
+    if (&x != &source)
+    {
+        x = source;
+    }
+
     if (Pstream::parRun())
     {
-        Field<Type> completeSourceSol(m());
+        Field<Type> X(m());
 
         if (Pstream::master(comm_))
         {
             typename Field<Type>::subField
             (
-                completeSourceSol,
-                sourceSol.size()
-            ).assign(sourceSol);
+                X,
+                x.size()
+            ).assign(x);
 
             for
             (
@@ -55,7 +66,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
                     slave,
                     reinterpret_cast<char*>
                     (
-                        &(completeSourceSol[procOffsets_[slave]])
+                        &(X[procOffsets_[slave]])
                     ),
                     (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
                     Pstream::msgType(),
@@ -69,8 +80,8 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
             (
                 Pstream::scheduled,
                 Pstream::masterNo(),
-                reinterpret_cast<const char*>(sourceSol.begin()),
-                sourceSol.byteSize(),
+                reinterpret_cast<const char*>(x.begin()),
+                x.byteSize(),
                 Pstream::msgType(),
                 comm_
             );
@@ -78,12 +89,12 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
 
         if (Pstream::master(comm_))
         {
-            LUBacksubstitute(*this, pivotIndices_, completeSourceSol);
+            LUBacksubstitute(*this, pivotIndices_, X);
 
-            sourceSol = typename Field<Type>::subField
+            x = typename Field<Type>::subField
             (
-                completeSourceSol,
-                sourceSol.size()
+                X,
+                x.size()
             );
 
             for
@@ -99,7 +110,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
                     slave,
                     reinterpret_cast<const char*>
                     (
-                        &(completeSourceSol[procOffsets_[slave]])
+                        &(X[procOffsets_[slave]])
                     ),
                     (procOffsets_[slave + 1]-procOffsets_[slave])*sizeof(Type),
                     Pstream::msgType(),
@@ -113,8 +124,8 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
             (
                 Pstream::scheduled,
                 Pstream::masterNo(),
-                reinterpret_cast<char*>(sourceSol.begin()),
-                sourceSol.byteSize(),
+                reinterpret_cast<char*>(x.begin()),
+                x.byteSize(),
                 Pstream::msgType(),
                 comm_
             );
@@ -122,9 +133,24 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
     }
     else
     {
-        LUBacksubstitute(*this, pivotIndices_, sourceSol);
+        LUBacksubstitute(*this, pivotIndices_, x);
     }
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::LUscalarMatrix::solve
+(
+    const Field<Type>& source
+) const
+{
+    tmp<Field<Type>> tx(new Field<Type>(m()));
+    Field<Type>& x = tx.ref();
+
+    solve(x, source);
+
+    return tx;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index 80e33b0fb2..165056d68e 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -225,7 +225,7 @@ Foam::QRMatrix<MatrixType>::solve
 
 template<class MatrixType>
 typename Foam::QRMatrix<MatrixType>::QMatrixType
-Foam::QRMatrix<MatrixType>::inverse() const
+Foam::QRMatrix<MatrixType>::inv() const
 {
     const label m = Q_.m();
 
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index bbad659be2..e3dab46b4f 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -108,7 +108,7 @@ public:
         tmp<Field<cmptType>> solve(const Field<cmptType>& source) const;
 
         //- Return the inverse of a square matrix
-        QMatrixType inverse() const;
+        QMatrixType inv() const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index 1a547204a7..3e357cdbe7 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -533,8 +533,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
 
     if (directSolveCoarsest_)
     {
-        coarsestCorrField = coarsestSource;
-        coarsestLUMatrixPtr_->solve(coarsestCorrField);
+        coarsestLUMatrixPtr_->solve(coarsestCorrField, coarsestSource);
     }
     //else if
     //(