diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
index 9208048b28..9f0422ca43 100644
--- a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
+++ b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.C
@@ -91,10 +91,11 @@ void Foam::LLTMatrix<Type>::decompose(const SquareMatrix<Type>& mat)
 
 
 template<class Type>
-void Foam::LLTMatrix<Type>::solve
+template<template<typename> class ListContainer>
+void Foam::LLTMatrix<Type>::solveImpl
 (
     List<Type>& x,
-    const UList<Type>& source
+    const ListContainer<Type>& source
 ) const
 {
     // If x and source are different, copy initialize x = source
@@ -132,6 +133,28 @@ void Foam::LLTMatrix<Type>::solve
 }
 
 
+template<class Type>
+void Foam::LLTMatrix<Type>::solve
+(
+    List<Type>& x,
+    const UList<Type>& source
+) const
+{
+    solveImpl(x, source);
+}
+
+
+template<class Type>
+template<class Addr>
+void Foam::LLTMatrix<Type>::solve
+(
+    List<Type>& x,
+    const IndirectListBase<Type, Addr>& source
+) const
+{
+    solveImpl(x, source);
+}
+
 template<class Type>
 Foam::tmp<Foam::Field<Type>> Foam::LLTMatrix<Type>::solve
 (
@@ -146,4 +169,19 @@ Foam::tmp<Foam::Field<Type>> Foam::LLTMatrix<Type>::solve
 }
 
 
+template<class Type>
+template<class Addr>
+Foam::tmp<Foam::Field<Type>> Foam::LLTMatrix<Type>::solve
+(
+    const IndirectListBase<Type, Addr>& source
+) const
+{
+    auto tresult(tmp<Field<Type>>::New(source.size()));
+
+    solve(tresult.ref(), source);
+
+    return tresult;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H
index a9f853e41c..ac7753ddd6 100644
--- a/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H
+++ b/src/OpenFOAM/matrices/LLTMatrix/LLTMatrix.H
@@ -58,6 +58,20 @@ class LLTMatrix
 :
     public SquareMatrix<Type>
 {
+
+    // Private Member Functions
+
+        //- Solve the linear system with the given source
+        //- and return the solution in x
+        //  This function may be called with the same field for x and source.
+        template<template<typename> class ListContainer>
+        void solveImpl
+        (
+            List<Type>& x,
+            const ListContainer<Type>& source
+        ) const;
+
+
 public:
 
     // Constructors
@@ -75,13 +89,38 @@ public:
         void decompose(const SquareMatrix<Type>& mat);
 
         //- Solve the linear system with the given source
-        //- and returning the solution in the Field argument x.
+        //- and return the solution in the argument x.
         //  This function may be called with the same field for x and source.
-        void solve(List<Type>& x, const UList<Type>& source) const;
+        void solve
+        (
+            List<Type>& x,
+            const UList<Type>& source
+        ) const;
 
         //- Solve the linear system with the given source
-        //- returning the solution
-        tmp<Field<Type>> solve(const UList<Type>& source) const;
+        //- and return the solution in the argument x.
+        //  This function may be called with the same field for x and source.
+        template<class Addr>
+        void solve
+        (
+            List<Type>& x,
+            const IndirectListBase<Type, Addr>& source
+        ) const;
+
+        //- Solve the linear system with the given source
+        //- return the solution
+        tmp<Field<Type>> solve
+        (
+            const UList<Type>& source
+        ) const;
+
+        //- Solve the linear system with the given source
+        //- return the solution
+        template<class Addr>
+        tmp<Field<Type>> solve
+        (
+            const IndirectListBase<Type, Addr>& source
+        ) const;
 };
 
 
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 6024d21c18..bfe1c4e51f 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -326,10 +326,7 @@ public:
         Form T() const;
 
         //- Multiply matrix with vector (A * x)
-        inline tmp<Field<Type>> Amul
-        (
-            const UList<Type>& x
-        ) const;
+        inline tmp<Field<Type>> Amul(const UList<Type>& x) const;
 
         //- Multiply matrix with vector (A * x)
         template<class Addr>
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
index 6c1ad97928..ea3c437a86 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.C
@@ -188,10 +188,11 @@ void Foam::QRMatrix<MatrixType>::solvex
 
 
 template<class MatrixType>
-void Foam::QRMatrix<MatrixType>::solve
+template<template<typename> class ListContainer>
+void Foam::QRMatrix<MatrixType>::solveImpl
 (
     List<cmptType>& x,
-    const UList<cmptType>& source
+    const ListContainer<cmptType>& source
 ) const
 {
     // Assert (&x != &source) ?
@@ -211,6 +212,29 @@ void Foam::QRMatrix<MatrixType>::solve
 }
 
 
+template<class MatrixType>
+void Foam::QRMatrix<MatrixType>::solve
+(
+    List<cmptType>& x,
+    const UList<cmptType>& source
+) const
+{
+    solveImpl(x, source);
+}
+
+
+template<class MatrixType>
+template<class Addr>
+void Foam::QRMatrix<MatrixType>::solve
+(
+    List<cmptType>& x,
+    const IndirectListBase<cmptType, Addr>& source
+) const
+{
+    solveImpl(x, source);
+}
+
+
 template<class MatrixType>
 Foam::tmp<Foam::Field<typename MatrixType::cmptType>>
 Foam::QRMatrix<MatrixType>::solve
@@ -218,9 +242,25 @@ Foam::QRMatrix<MatrixType>::solve
     const UList<cmptType>& source
 ) const
 {
-    auto tresult(tmp<Field<cmptType>>::New(Q_.m()));
+    auto tresult(Q_.Tmul(source));
 
-    solve(tresult.ref(), source);
+    solvex(tresult.ref());
+
+    return tresult;
+}
+
+
+template<class MatrixType>
+template<class Addr>
+Foam::tmp<Foam::Field<typename MatrixType::cmptType>>
+Foam::QRMatrix<MatrixType>::solve
+(
+    const IndirectListBase<cmptType, Addr>& source
+) const
+{
+    auto tresult(Q_.Tmul(source));
+
+    solvex(tresult.ref());
 
     return tresult;
 }
diff --git a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
index b716cd1954..0dc739ce93 100644
--- a/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
+++ b/src/OpenFOAM/matrices/QRMatrix/QRMatrix.H
@@ -79,6 +79,15 @@ private:
         template<template<typename> class ListContainer>
         void solvex(ListContainer<cmptType>& x) const;
 
+        //- Solve the linear system with the given source
+        //- and return the solution in x
+        template<template<typename> class ListContainer>
+        void solveImpl
+        (
+            List<cmptType>& x,
+            const ListContainer<cmptType>& source
+        ) const;
+
 
 public:
 
@@ -103,12 +112,36 @@ public:
         void decompose(const MatrixType& M);
 
         //- Solve the linear system with the given source
-        //- and return the solution in the Field argument x
-        void solve(List<cmptType>& x, const UList<cmptType>& source) const;
+        //- and return the solution in the argument x
+        void solve
+        (
+            List<cmptType>& x,
+            const UList<cmptType>& source
+        ) const;
+
+        //- Solve the linear system with the given source
+        //- and return the solution in the argument x
+        template<class Addr>
+        void solve
+        (
+            List<cmptType>& x,
+            const IndirectListBase<cmptType, Addr>& source
+        ) const;
 
         //- Solve the linear system with the given source
         //- and return the solution
-        tmp<Field<cmptType>> solve(const UList<cmptType>& source) const;
+        tmp<Field<cmptType>> solve
+        (
+            const UList<cmptType>& source
+        ) const;
+
+        //- Solve the linear system with the given source
+        //- and return the solution
+        template<class Addr>
+        tmp<Field<cmptType>> solve
+        (
+            const IndirectListBase<cmptType, Addr>& source
+        ) const;
 
         //- Return the inverse of (Q*R), so that solving x = (Q*R).inv()*source
         QMatrixType inv() const;