diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 808e00b86c..46508d8e3f 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -30,7 +30,7 @@ License
 template<class Form, class Type>
 void Foam::Matrix<Form, Type>::allocate()
 {
-    if (nRows_ && nCols_)
+    if (mRows_ && nCols_)
     {
         v_ = new Type[size()];
     }
@@ -42,14 +42,14 @@ void Foam::Matrix<Form, Type>::allocate()
 template<class Form, class Type>
 Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
 :
-    nRows_(m),
+    mRows_(m),
     nCols_(n),
     v_(NULL)
 {
-    if (nRows_ < 0 || nCols_ < 0)
+    if (mRows_ < 0 || nCols_ < 0)
     {
         FatalErrorInFunction
-            << "Bad m, n " << nRows_ << ", " << nCols_
+            << "Incorrect m, n " << mRows_ << ", " << nCols_
             << abort(FatalError);
     }
 
@@ -60,14 +60,14 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
 template<class Form, class Type>
 Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const zero)
 :
-    nRows_(m),
+    mRows_(m),
     nCols_(n),
     v_(NULL)
 {
-    if (nRows_ < 0 || nCols_ < 0)
+    if (mRows_ < 0 || nCols_ < 0)
     {
         FatalErrorInFunction
-            << "Bad m, n " << nRows_ << ", " << nCols_
+            << "Incorrect m, n " << mRows_ << ", " << nCols_
             << abort(FatalError);
     }
 
@@ -85,16 +85,16 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const zero)
 
 
 template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& s)
 :
-    nRows_(m),
+    mRows_(m),
     nCols_(n),
     v_(NULL)
 {
-    if (nRows_ < 0 || nCols_ < 0)
+    if (mRows_ < 0 || nCols_ < 0)
     {
         FatalErrorInFunction
-            << "Bad m, n " << nRows_ << ", " << nCols_
+            << "Incorrect m, n " << mRows_ << ", " << nCols_
             << abort(FatalError);
     }
 
@@ -105,45 +105,92 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
         const label mn = size();
         for (label i=0; i<mn; i++)
         {
-            v_[i] = a;
+            v_[i] = s;
         }
     }
 }
 
 
 template<class Form, class Type>
-inline Foam::Matrix<Form, Type>::Matrix(const mType::Block& block)
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& M)
 :
-    nRows_(block.m()),
-    nCols_(block.n())
-{
-    allocate();
-
-    for (label i=0; i<nRows_; i++)
-    {
-        for (label j=0; j<nCols_; j++)
-        {
-            (*this)(i,j) = block(i,j);
-        }
-    }
-}
-
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
-:
-    nRows_(a.nRows_),
-    nCols_(a.nCols_),
+    mRows_(M.mRows_),
+    nCols_(M.nCols_),
     v_(NULL)
 {
-    if (a.v_)
+    if (M.v_)
     {
         allocate();
 
         const label mn = size();
         for (label i=0; i<mn; i++)
         {
-            v_[i] = a.v_[i];
+            v_[i] = M.v_[i];
+        }
+    }
+}
+
+
+template<class Form, class Type>
+template<class Form2>
+Foam::Matrix<Form, Type>::Matrix(const Matrix<Form2, Type>& M)
+:
+    mRows_(M.m()),
+    nCols_(M.n()),
+    v_(NULL)
+{
+    if (M.v())
+    {
+        allocate();
+
+        const label mn = size();
+        for (label i=0; i<mn; i++)
+        {
+            v_[i] = M.v()[i];
+        }
+    }
+}
+
+
+template<class Form, class Type>
+template<class MatrixType>
+inline Foam::Matrix<Form, Type>::Matrix
+(
+    const ConstMatrixBlock<MatrixType>& Mb
+)
+:
+    mRows_(Mb.m()),
+    nCols_(Mb.n())
+{
+    allocate();
+
+    for (label i=0; i<mRows_; i++)
+    {
+        for (label j=0; j<nCols_; j++)
+        {
+            (*this)(i,j) = Mb(i,j);
+        }
+    }
+}
+
+
+template<class Form, class Type>
+template<class MatrixType>
+inline Foam::Matrix<Form, Type>::Matrix
+(
+    const MatrixBlock<MatrixType>& Mb
+)
+:
+    mRows_(Mb.m()),
+    nCols_(Mb.n())
+{
+    allocate();
+
+    for (label i=0; i<mRows_; i++)
+    {
+        for (label j=0; j<nCols_; j++)
+        {
+            (*this)(i,j) = Mb(i,j);
         }
     }
 }
@@ -172,24 +219,24 @@ void Foam::Matrix<Form, Type>::clear()
         v_ = NULL;
     }
 
-    nRows_ = 0;
+    mRows_ = 0;
     nCols_ = 0;
 }
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
+void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& M)
 {
     clear();
 
-    nRows_ = a.nRows_;
-    a.nRows_ = 0;
+    mRows_ = M.mRows_;
+    M.mRows_ = 0;
 
-    nCols_ = a.nCols_;
-    a.nCols_ = 0;
+    nCols_ = M.nCols_;
+    M.nCols_ = 0;
 
-    v_ = a.v_;
-    a.v_ = NULL;
+    v_ = M.v_;
+    M.v_ = NULL;
 }
 
 
@@ -198,7 +245,7 @@ void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
 {
     mType newMatrix(m, n, Zero);
 
-    label minM = min(m, nRows_);
+    label minM = min(m, mRows_);
     label minN = min(n, nCols_);
 
     for (label i=0; i<minM; i++)
@@ -231,65 +278,23 @@ Form Foam::Matrix<Form, Type>::T() const
 }
 
 
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::ConstBlock::operator Field<cmptType>() const
-{
-    if (nCols_ != 1)
-    {
-        FatalErrorInFunction
-            << "Number of columns " << nCols_ << " != 1"
-            << abort(FatalError);
-    }
-
-    Field<cmptType> f(nRows_);
-
-    forAll(f, i)
-    {
-        f[i] = operator()(i, 0);
-    }
-
-    return f;
-}
-
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::Block::operator Field<cmptType>() const
-{
-    if (nCols_ != 1)
-    {
-        FatalErrorInFunction
-            << "Number of columns " << nCols_ << " != 1"
-            << abort(FatalError);
-    }
-
-    Field<cmptType> f(nRows_);
-
-    forAll(f, i)
-    {
-        f[i] = operator()(i, 0);
-    }
-
-    return f;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
+void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& M)
 {
-    if (this == &a)
+    if (this == &M)
     {
         FatalErrorInFunction
             << "Attempted assignment to self"
             << abort(FatalError);
     }
 
-    if (nRows_ != a.nRows_ || nCols_ != a.nCols_)
+    if (mRows_ != M.mRows_ || nCols_ != M.nCols_)
     {
         clear();
-        nRows_ = a.nRows_;
-        nCols_ = a.nCols_;
+        mRows_ = M.mRows_;
+        nCols_ = M.nCols_;
         allocate();
     }
 
@@ -298,106 +303,55 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
         const label mn = size();
         for (label i=0; i<mn; i++)
         {
-            v_[i] = a.v_[i];
+            v_[i] = M.v_[i];
         }
     }
 }
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const mType::Block& block)
+template<class MatrixType>
+void Foam::Matrix<Form, Type>::operator=
+(
+    const ConstMatrixBlock<MatrixType>& Mb
+)
 {
-    for (label i=0; i<nRows_; i++)
+    for (label i=0; i<mRows_; i++)
     {
         for (label j=0; j<nCols_; j++)
         {
-            (*this)(i,j) = block(i,j);
+            (*this)(i,j) = Mb(i,j);
         }
     }
 }
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const Block& block)
+template<class MatrixType>
+void Foam::Matrix<Form, Type>::operator=
+(
+    const MatrixBlock<MatrixType>& Mb
+)
 {
-    if (this != &block)
-    {
-        if (nRows_ != block.m() || nCols_ != block.n())
-        {
-            FatalErrorInFunction
-                << "Attempt to assign blocks of different sizes: "
-                << nRows_ << "x" << nCols_ << " != "
-                << block.m() << "x" << block.n()
-                << abort(FatalError);
-        }
-
-        for (label i=0; i<nRows_; i++)
-        {
-            for (label j=0; j<nCols_; j++)
-            {
-                (*this)(i, j) = block(i, j);
-            }
-        }
-    }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const ConstBlock& block)
-{
-    if (this != &block)
-    {
-        if (nRows_ != block.m() || nCols_ != block.n())
-        {
-            FatalErrorInFunction
-                << "Attempt to assign blocks of different sizes: "
-                << nRows_ << "x" << nCols_ << " != "
-                << block.m() << "x" << block.n()
-                << abort(FatalError);
-        }
-
-        for (label i=0; i<nRows_; i++)
-        {
-            for (label j=0; j<nCols_; j++)
-            {
-                (*this)(i, j) = block(i, j);
-            }
-        }
-    }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const mType& block)
-{
-    if (nRows_ != block.m() || nCols_ != block.n())
-    {
-        FatalErrorInFunction
-            << "Attempt to assign blocks of different sizes: "
-            << nRows_ << "x" << nCols_ << " != "
-            << block.m() << "x" << block.n()
-            << abort(FatalError);
-    }
-
-    for (label i=0; i<nRows_; i++)
+    for (label i=0; i<mRows_; i++)
     {
         for (label j=0; j<nCols_; j++)
         {
-            (*this)(i, j) = block(i, j);
+            (*this)(i,j) = Mb(i,j);
         }
     }
 }
 
 
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Type& t)
+void Foam::Matrix<Form, Type>::operator=(const Type& s)
 {
     if (v_)
     {
         const label mn = size();
         for (label i=0; i<mn; i++)
         {
-            v_[i] = t;
+            v_[i] = s;
         }
     }
 }
@@ -417,155 +371,27 @@ void Foam::Matrix<Form, Type>::operator=(const zero)
 }
 
 
-template<class Form, class Type>
-template
-<
-    template<class, Foam::direction, Foam::direction> class MSBlock,
-    class SubTensor,
-    Foam::direction BRowStart,
-    Foam::direction BColStart
->
-void Foam::Matrix<Form, Type>::Block::operator=
-(
-    const MSBlock<SubTensor, BRowStart, BColStart>& block
-)
-{
-    if (nRows_ != block.nRows || nCols_ != block.nCols)
-    {
-        FatalErrorInFunction
-            << "Attempt to assign blocks of different sizes: "
-            << nRows_ << "x" << nCols_ << " != "
-            << block.nRows << "x" << block.nCols
-            << abort(FatalError);
-    }
-
-    for (direction i=0; i<nRows_; ++i)
-    {
-        for (direction j=0; j<nCols_; ++j)
-        {
-            operator()(i, j) = block(i, j);
-        }
-    }
-}
-
-
-template<class Form, class Type>
-template
-<
-    template<class, Foam::direction> class VSBlock,
-    class SubVector,
-    Foam::direction BStart
->
-void Foam::Matrix<Form, Type>::Block::operator=
-(
-    const VSBlock<SubVector, BStart>& block
-)
-{
-    if (nRows_ != block.nComponents || nCols_ != 1)
-    {
-        FatalErrorInFunction
-            << "Attempt to assign blocks of different sizes: "
-            << nRows_ << "x" << nCols_ << " != "
-            << block.nComponents << "x" << 1
-            << abort(FatalError);
-    }
-
-    for (direction i=0; i<nRows_; ++i)
-    {
-        operator()(i, 0) = block[i];
-    }
-}
-
-
-template<class Form, class Type>
-template<class MSForm, Foam::direction Nrows, Foam::direction Ncols>
-void Foam::Matrix<Form, Type>::Block::operator=
-(
-    const MatrixSpace<MSForm, cmptType, Nrows, Ncols>& ms
-)
-{
-    if (nRows_ != Nrows || nCols_ != Ncols)
-    {
-        FatalErrorInFunction
-            << "Attempt to assign blocks of different sizes: "
-            << nRows_ << "x" << nCols_ << " != "
-            << Nrows << "x" << Ncols
-            << abort(FatalError);
-    }
-
-    for (label i=0; i<nRows_; i++)
-    {
-        for (label j=0; j<nCols_; j++)
-        {
-            (*this)(i, j) = ms(i, j);
-        }
-    }
-}
-
-
-template<class Form, class Type>
-template<class VSForm, Foam::direction Ncmpts>
-void Foam::Matrix<Form, Type>::Block::operator=
-(
-    const VectorSpace<VSForm, cmptType, Ncmpts>& ms
-)
-{
-    if (nRows_ != Ncmpts || nCols_ != 1)
-    {
-        FatalErrorInFunction
-            << "Attempt to assign blocks of different sizes: "
-            << nRows_ << "x" << nCols_ << " != "
-            << Ncmpts << "x" << 1
-            << abort(FatalError);
-    }
-
-    for (direction i=0; i<Ncmpts; ++i)
-    {
-        operator()(i, 0) = ms[i];
-    }
-}
-
-
-template<class Form, class Type>
-void Foam::Matrix<Form, Type>::Block::operator=(const Field<cmptType>& f)
-{
-    if (nRows_ != f.size() || nCols_ != 1)
-    {
-        FatalErrorInFunction
-            << "Error: cannot assign blocks of different size (left is "
-            << nRows_ << "x" << nCols_ << " != "
-            << f.size() << "x" << 1
-            << abort(FatalError);
-    }
-
-    forAll(f, i)
-    {
-        operator()(i, 0) = f[i];
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
-const Type& Foam::max(const Matrix<Form, Type>& a)
+const Type& Foam::max(const Matrix<Form, Type>& M)
 {
-    const label mn = a.size();
+    const label mn = M.size();
 
     if (mn)
     {
         label curMaxI = 0;
-        const Type* v = a.v();
+        const Type* Mv = M.v();
 
         for (label i=1; i<mn; i++)
         {
-            if (v[i] > v[curMaxI])
+            if (Mv[i] > Mv[curMaxI])
             {
                 curMaxI = i;
             }
         }
 
-        return v[curMaxI];
+        return Mv[curMaxI];
     }
     else
     {
@@ -574,30 +400,30 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
             << abort(FatalError);
 
         // Return in error to keep compiler happy
-        return a(0, 0);
+        return M(0, 0);
     }
 }
 
 
 template<class Form, class Type>
-const Type& Foam::min(const Matrix<Form, Type>& a)
+const Type& Foam::min(const Matrix<Form, Type>& M)
 {
-    const label mn = a.size();
+    const label mn = M.size();
 
     if (mn)
     {
         label curMinI = 0;
-        const Type* v = a.v();
+        const Type* Mv = M.v();
 
         for (label i=1; i<mn; i++)
         {
-            if (v[i] < v[curMinI])
+            if (Mv[i] < Mv[curMinI])
             {
                 curMinI = i;
             }
         }
 
-        return v[curMinI];
+        return Mv[curMinI];
     }
     else
     {
@@ -606,7 +432,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
             << abort(FatalError);
 
         // Return in error to keep compiler happy
-        return a(0, 0);
+        return M(0, 0);
     }
 }
 
@@ -614,187 +440,197 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& a)
+Form Foam::operator-(const Matrix<Form, Type>& M)
 {
-    Form na(a.m(), a.n());
+    Form nM(M.m(), M.n());
 
-    if (a.m() && a.n())
+    if (M.m() && M.n())
     {
-        Type* nav = na.v();
-        const Type* av = a.v();
+        Type* nMv = nM.v();
+        const Type* Mv = M.v();
 
-        const label mn = a.size();
+        const label mn = M.size();
         for (label i=0; i<mn; i++)
         {
-            nav[i] = -av[i];
+            nMv[i] = -Mv[i];
         }
     }
 
-    return na;
+    return nM;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+Form Foam::operator+(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
 {
-    if (a.m() != b.m())
+    if (A.m() != B.m())
     {
         FatalErrorInFunction
-            << "Attempt to add matrices with different number of rows: "
-            << a.m() << ", " << b.m()
+            << "Attempt to add matrices with different numbers of rows: "
+            << A.m() << ", " << B.m()
             << abort(FatalError);
     }
 
-    if (a.n() != b.n())
+    if (A.n() != B.n())
     {
         FatalErrorInFunction
-            << "Attempt to add matrices with different number of columns: "
-            << a.n() << ", " << b.n()
+            << "Attempt to add matrices with different numbers of columns: "
+            << A.n() << ", " << B.n()
             << abort(FatalError);
     }
 
-    Form ab(a.m(), a.n());
+    Form AB(A.m(), A.n());
 
-    Type* abv = ab.v();
-    const Type* av = a.v();
-    const Type* bv = b.v();
+    Type* ABv = AB.v();
+    const Type* Av = A.v();
+    const Type* Bv = B.v();
 
-    const label mn = a.size();
+    const label mn = A.size();
     for (label i=0; i<mn; i++)
     {
-        abv[i] = av[i] + bv[i];
+        ABv[i] = Av[i] + Bv[i];
     }
 
-    return ab;
+    return AB;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+Form Foam::operator-(const Matrix<Form, Type>& A, const Matrix<Form, Type>& B)
 {
-    if (a.m() != b.m())
+    if (A.m() != B.m())
     {
         FatalErrorInFunction
-            << "Attempt to add matrices with different number of rows: "
-            << a.m() << ", " << b.m()
+            << "Attempt to add matrices with different numbers of rows: "
+            << A.m() << ", " << B.m()
             << abort(FatalError);
     }
 
-    if (a.n() != b.n())
+    if (A.n() != B.n())
     {
         FatalErrorInFunction
-            << "Attempt to add matrices with different number of columns: "
-            << a.n() << ", " << b.n()
+            << "Attempt to add matrices with different numbers of columns: "
+            << A.n() << ", " << B.n()
             << abort(FatalError);
     }
 
-    Form ab(a.m(), a.n());
+    Form AB(A.m(), A.n());
 
-    Type* abv = ab.v();
-    const Type* av = a.v();
-    const Type* bv = b.v();
+    Type* ABv = AB.v();
+    const Type* Av = A.v();
+    const Type* Bv = B.v();
 
-    const label mn = a.size();
+    const label mn = A.size();
     for (label i=0; i<mn; i++)
     {
-        abv[i] = av[i] - bv[i];
+        ABv[i] = Av[i] - Bv[i];
     }
 
-    return ab;
+    return AB;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
+Form Foam::operator*(const scalar s, const Matrix<Form, Type>& M)
 {
-    Form sa(a.m(), a.n());
+    Form sM(M.m(), M.n());
 
-    if (a.m() && a.n())
+    if (M.m() && M.n())
     {
-        Type* sav = sa.v();
-        const Type* av = a.v();
+        Type* sMv = sM.v();
+        const Type* Mv = M.v();
 
-        const label mn = a.size();
+        const label mn = M.size();
         for (label i=0; i<mn; i++)
         {
-            sav[i] = s*av[i];
+            sMv[i] = s*Mv[i];
         }
     }
 
-    return sa;
+    return sM;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& a, const scalar s)
+Form Foam::operator*(const Matrix<Form, Type>& M, const scalar s)
 {
-    Form sa(a.m(), a.n());
+    Form sM(M.m(), M.n());
 
-    if (a.m() && a.n())
+    if (M.m() && M.n())
     {
-        Type* sav = sa.v();
-        const Type* av = a.v();
+        Type* sMv = sM.v();
+        const Type* Mv = M.v();
 
-        const label mn = a.size();
+        const label mn = M.size();
         for (label i=0; i<mn; i++)
         {
-            sav[i] = av[i]*s;
+            sMv[i] = Mv[i]*s;
         }
     }
 
-    return sa;
+    return sM;
 }
 
 
 template<class Form, class Type>
-Form Foam::operator/(const Matrix<Form, Type>& a, const scalar s)
+Form Foam::operator/(const Matrix<Form, Type>& M, const scalar s)
 {
-    Form sa(a.m(), a.n());
+    Form sM(M.m(), M.n());
 
-    if (a.m() && a.n())
+    if (M.m() && M.n())
     {
-        Type* sav = sa.v();
-        const Type* av = a.v();
+        Type* sMv = sM.v();
+        const Type* Mv = M.v();
 
-        const label mn = a.size();
+        const label mn = M.size();
         for (label i=0; i<mn; i++)
         {
-            sav[i] = av[i]/s;
+            sMv[i] = Mv[i]/s;
         }
     }
 
-    return sa;
+    return sM;
 }
 
 
-template<class Form, class Type>
-Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
+template<class Form1, class Form2, class Type>
+typename Foam::typeOfInnerProduct<Type, Form1, Form2>::type
+Foam::operator*
+(
+    const Matrix<Form1, Type>& A,
+    const Matrix<Form2, Type>& B
+)
 {
-    if (a.n() != b.m())
+    if (A.n() != B.m())
     {
         FatalErrorInFunction
             << "Attempt to multiply incompatible matrices:" << nl
-            << "Matrix A : " << a.m() << " rows, " << a.n() << " columns" << nl
-            << "Matrix B : " << b.m() << " rows, " << b.n() << " columns" << nl
+            << "Matrix A : " << A.m() << " x " << A.n() << nl
+            << "Matrix B : " << B.m() << " x " << B.n() << nl
             << "In order to multiply matrices, columns of A must equal "
             << "rows of B"
             << abort(FatalError);
     }
 
-    Form ab(a.m(), b.n(), scalar(0));
+    typename typeOfInnerProduct<Type, Form1, Form2>::type AB
+    (
+        A.m(),
+        B.n(),
+        Zero
+    );
 
-    for (label i=0; i<ab.m(); i++)
+    for (label i=0; i<AB.m(); i++)
     {
-        for (label j=0; j<ab.n(); j++)
+        for (label j=0; j<AB.n(); j++)
         {
-            for (label k=0; k<b.m(); k++)
+            for (label k=0; k<B.m(); k++)
             {
-                ab(i, j) += a(i, k)*b(k, j);
+                AB(i, j) += A(i, k)*B(k, j);
             }
         }
     }
 
-    return ab;
+    return AB;
 }
 
 
@@ -809,7 +645,7 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
     {
         FatalErrorInFunction
             << "Attempt to multiply incompatible matrix and field:" << nl
-            << "Matrix : " << M.m() << " rows, " << M.n() << " columns" << nl
+            << "Matrix : " << M.m() << " x " << M.n() << nl
             << "Field : " << f.size() << " rows" << nl
             << "In order to multiply a matrix M and field f, "
                "columns of M must equal rows of f"
@@ -819,9 +655,9 @@ inline Foam::tmp<Foam::Field<Type>> Foam::operator*
     tmp<Field<Type>> tMf(new Field<Type>(f.size(), Zero));
     Field<Type>& Mf = tMf.ref();
 
-    for (label i=0; i<M.m(); ++i)
+    for (label i=0; i<M.m(); i++)
     {
-        for (label j=0; j<M.n(); ++j)
+        for (label j=0; j<M.n(); j++)
         {
             Mf[i] += M(i, j)*f[j];
         }
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 5f905b465b..9385d42d87 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -25,8 +25,7 @@ Class
     Foam::Matrix
 
 Description
-    A templated 2D matrix of objects of \<T\>, where the n x m matrix
-    dimensions are known and used for subscript bounds checking, etc.
+    A templated (m x n) matrix of objects of \<T\>.
 
 SourceFiles
     Matrix.C
@@ -42,7 +41,6 @@ SourceFiles
 #include "label.H"
 #include "uLabel.H"
 #include "Field.H"
-#include "MatrixSpace.H"
 #include "zero.H"
 #include "autoPtr.H"
 
@@ -67,6 +65,12 @@ template<class Form, class Type> Ostream& operator<<
     const Matrix<Form, Type>&
 );
 
+template<class MatrixType>
+class ConstMatrixBlock;
+
+template<class MatrixType>
+class MatrixBlock;
+
 
 /*---------------------------------------------------------------------------*\
                            Class Matrix Declaration
@@ -78,7 +82,7 @@ class Matrix
     // Private data
 
         //- Number of rows and columns in Matrix.
-        label nRows_, nCols_;
+        label mRows_, nCols_;
 
         //- Row pointers
         Type* __restrict__ v_;
@@ -102,157 +106,6 @@ public:
         inline static const mType& null();
 
 
-    // Sub-Block Classes
-
-        class ConstBlock
-        {
-            // Private data
-
-                //- Const reference to the parent matrix
-                const mType& matrix_;
-
-                // Block size
-                const label nRows_;
-                const label nCols_;
-
-                // Block location in parent matrix
-                const label rowStart_;
-                const label colStart_;
-
-        public:
-
-            typedef typename mType::cmptType cmptType;
-
-            // Constructors
-
-                //- Construct block for matrix, size and location
-                inline ConstBlock
-                (
-                    const mType& matrix,
-                    const label m,
-                    const label n,
-                    const label mStart,
-                    const label nStart
-                );
-
-
-            // Member Functions
-
-                //- Return the number of rows in the block
-                inline label m() const;
-
-                //- Return the number of columns in the block
-                inline label n() const;
-
-                //- (i, j) const element access operator
-                inline const Type& operator()
-                (
-                    const label i,
-                    const label j
-                ) const;
-
-                //- Convert a column of a matrix to a Field
-                operator Field<Type>() const;
-        };
-
-
-        class Block
-        {
-            // Private data
-
-                //- Reference to the parent matrix
-                mType& matrix_;
-
-                // Block size
-                const label nRows_;
-                const label nCols_;
-
-                // Block location in parent matrix
-                const label rowStart_;
-                const label colStart_;
-
-        public:
-
-            typedef typename mType::cmptType cmptType;
-
-            // Constructors
-
-                //- Construct block for matrix, size and location
-                inline Block
-                (
-                    mType& matrix,
-                    const label m,
-                    const label n,
-                    const label mStart,
-                    const label nStart
-                );
-
-
-            // Member Functions
-
-                //- Return the number of rows in the block
-                inline label m() const;
-
-                //- Return the number of columns in the block
-                inline label n() const;
-
-                //- (i, j) const element access operator
-                inline const Type& operator()
-                (
-                    const label i,
-                    const label j
-                ) const;
-
-                //- (i, j) element access operator
-                inline Type& operator()(const label i, const label j);
-
-                //- Convert a column of a matrix to a Field
-                operator Field<Type>() const;
-
-
-            // Member operators
-
-                //- Assignment to a compatible matrix
-                void operator=(const mType&);
-
-                //- Assignment to a compatible block
-                void operator=(const Block&);
-
-                //- Assignment to a compatible const block
-                void operator=(const ConstBlock&);
-
-                //- Assignment to a compatible MatrixSpace
-                template<class MSForm, direction Nrows, direction Ncols>
-                void operator=(const MatrixSpace<MSForm, Type, Nrows, Ncols>&);
-
-                //- Assignment to a compatible MatrixSpace block
-                template
-                <
-                    template<class, direction, direction> class Block,
-                    class SubTensor,
-                    direction BRowStart,
-                    direction BColStart
-                >
-                void operator=(const Block<SubTensor, BRowStart, BColStart>&);
-
-                //- Assignment to a compatible VectorSpace (column-vector)
-                template<class VSForm, direction Ncmpts>
-                void operator=(const VectorSpace<VSForm, Type, Ncmpts>&);
-
-                //- Assignment to a compatible VectorSpace (column-vector) block
-                template
-                <
-                    template<class, direction> class Block,
-                    class SubVector,
-                    direction BStart
-                >
-                void operator=(const Block<SubVector, BStart>&);
-
-                //- Assignment to a Field (column-vector)
-                void operator=(const Field<Type>&);
-        };
-
-
     // Constructors
 
         //- Null constructor.
@@ -272,8 +125,17 @@ public:
         //- Copy constructor.
         Matrix(const mType&);
 
+        //- Copy constructor from matrix of a different form
+        template<class Form2>
+        explicit Matrix(const Matrix<Form2, Type>&);
+
         //- Construct from a block of another matrix
-        Matrix(const mType::Block&);
+        template<class MatrixType>
+        Matrix(const ConstMatrixBlock<MatrixType>&);
+
+        //- Construct from a block of another matrix
+        template<class MatrixType>
+        Matrix(const MatrixBlock<MatrixType>&);
 
         //- Construct from Istream.
         Matrix(Istream&);
@@ -308,7 +170,7 @@ public:
 
         // Block access
 
-            inline ConstBlock block
+            inline ConstMatrixBlock<mType> block
             (
                 const label m,
                 const label n,
@@ -317,19 +179,19 @@ public:
             ) const;
 
             template<class VectorSpace>
-            inline ConstBlock block
+            inline ConstMatrixBlock<mType> block
             (
                 const label mStart,
                 const label nStart
             ) const;
 
-            inline ConstBlock col
+            inline ConstMatrixBlock<mType> col
             (
                 const label m,
                 const label rowStart
             ) const;
 
-            inline ConstBlock col
+            inline ConstMatrixBlock<mType> col
             (
                 const label m,
                 const label mStart,
@@ -337,7 +199,7 @@ public:
             ) const;
 
 
-            inline Block block
+            inline MatrixBlock<mType> block
             (
                 const label m,
                 const label n,
@@ -346,11 +208,19 @@ public:
             );
 
             template<class VectorSpace>
-            inline Block block(const label mStart, const label nStart);
+            inline MatrixBlock<mType> block
+            (
+                const label mStart,
+                const label nStart
+            );
 
-            inline Block col(const label m, const label rowStart);
+            inline MatrixBlock<mType> col
+            (
+                const label m,
+                const label rowStart
+            );
 
-            inline Block col
+            inline MatrixBlock<mType> col
             (
                 const label m,
                 const label mStart,
@@ -402,7 +272,12 @@ public:
         void operator=(const mType&);
 
         //- Assignment to a block of another matrix
-        void operator=(const mType::Block&);
+        template<class MatrixType>
+        void operator=(const ConstMatrixBlock<MatrixType>&);
+
+        //- Assignment to a block of another matrix
+        template<class MatrixType>
+        void operator=(const MatrixBlock<MatrixType>&);
 
         //- Assignment of all entries to zero
         void operator=(const zero);
@@ -475,11 +350,12 @@ Form operator/
     const scalar
 );
 
-template<class Form, class Type>
-Form operator*
+template<class Form1, class Form2, class Type>
+typename typeOfInnerProduct<Type, Form1, Form2>::type
+operator*
 (
-    const Matrix<Form, Type>&,
-    const Matrix<Form, Type>&
+    const Matrix<Form1, Type>& a,
+    const Matrix<Form2, Type>& b
 );
 
 template<class Form, class Type>
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 73992f3285..73992aac77 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -23,12 +23,14 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "MatrixBlock.H"
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Form, class Type>
 inline Foam::Matrix<Form, Type>::Matrix()
 :
-    nRows_(0),
+    mRows_(0),
     nCols_(0),
     v_(NULL)
 {}
@@ -42,68 +44,6 @@ clone() const
 }
 
 
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::ConstBlock::ConstBlock
-(
-    const mType& matrix,
-    const label m,
-    const label n,
-    const label mStart,
-    const label nStart
-)
-:
-    matrix_(matrix),
-    nRows_(m),
-    nCols_(n),
-    rowStart_(mStart),
-    colStart_(nStart)
-{
-    #ifdef FULLDEBUG
-    if
-    (
-        rowStart_ + nRows_ > matrix.m()
-     || colStart_ + nCols_ > matrix.n()
-    )
-    {
-        FatalErrorInFunction
-            << "Block addresses outside matrix"
-            << abort(FatalError);
-    }
-    #endif
-}
-
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::Block::Block
-(
-    mType& matrix,
-    const label m,
-    const label n,
-    const label mStart,
-    const label nStart
-)
-:
-    matrix_(matrix),
-    nRows_(m),
-    nCols_(n),
-    rowStart_(mStart),
-    colStart_(nStart)
-{
-    #ifdef FULLDEBUG
-    if
-    (
-        rowStart_ + nRows_ > matrix.m()
-     || colStart_ + nCols_ > matrix.n()
-    )
-    {
-        FatalErrorInFunction
-            << "Block addresses outside matrix"
-            << abort(FatalError);
-    }
-    #endif
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
@@ -116,7 +56,7 @@ inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
 template<class Form, class Type>
 inline Foam::label Foam::Matrix<Form, Type>::m() const
 {
-    return nRows_;
+    return mRows_;
 }
 
 
@@ -130,7 +70,7 @@ inline Foam::label Foam::Matrix<Form, Type>::n() const
 template<class Form, class Type>
 inline Foam::label Foam::Matrix<Form, Type>::size() const
 {
-    return nRows_*nCols_;
+    return mRows_*nCols_;
 }
 
 
@@ -138,16 +78,16 @@ template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checki(const label i) const
 {
     #ifdef FULLDEBUG
-    if (!nRows_ || !nCols_)
+    if (!mRows_ || !nCols_)
     {
         FatalErrorInFunction
             << "Attempt to access element from empty matrix"
             << abort(FatalError);
     }
-    else if (i<0 || i>=nRows_)
+    else if (i<0 || i>=mRows_)
     {
         FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << nRows_-1
+            << "Index " << i << " out of range 0 ... " << mRows_-1
             << abort(FatalError);
     }
     #endif
@@ -158,7 +98,7 @@ template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 {
     #ifdef FULLDEBUG
-    if (!nRows_ || !nCols_)
+    if (!mRows_ || !nCols_)
     {
         FatalErrorInFunction
             << "Attempt to access element from empty matrix"
@@ -174,34 +114,6 @@ inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 }
 
 
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::ConstBlock::m() const
-{
-    return nRows_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::ConstBlock::n() const
-{
-    return nCols_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::Block::m() const
-{
-    return nRows_;
-}
-
-
-template<class Form, class Type>
-inline Foam::label Foam::Matrix<Form, Type>::Block::n() const
-{
-    return nCols_;
-}
-
-
 template<class Form, class Type>
 inline const Type* Foam::Matrix<Form, Type>::v() const
 {
@@ -217,7 +129,7 @@ inline Type* Foam::Matrix<Form, Type>::v()
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::block
 (
     const label m,
@@ -226,7 +138,7 @@ Foam::Matrix<Form, Type>::block
     const label nStart
 ) const
 {
-    return ConstBlock
+    return ConstMatrixBlock<mType>
     (
         *this,
         m,
@@ -239,17 +151,17 @@ Foam::Matrix<Form, Type>::block
 
 template<class Form, class Type>
 template<class VectorSpace>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::block
 (
     const label mStart,
     const label nStart
 ) const
 {
-    return ConstBlock
+    return ConstMatrixBlock<mType>
     (
         *this,
-        VectorSpace::nRows,
+        VectorSpace::mRows,
         VectorSpace::nCols,
         mStart,
         nStart
@@ -258,14 +170,14 @@ Foam::Matrix<Form, Type>::block
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::col
 (
     const label m,
     const label mStart
 ) const
 {
-    return ConstBlock
+    return ConstMatrixBlock<mType>
     (
         *this,
         m,
@@ -277,7 +189,7 @@ Foam::Matrix<Form, Type>::col
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::ConstBlock
+inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::col
 (
     const label m,
@@ -285,7 +197,7 @@ Foam::Matrix<Form, Type>::col
     const label nStart
 ) const
 {
-    return ConstBlock
+    return ConstMatrixBlock<mType>
     (
         *this,
         m,
@@ -297,7 +209,7 @@ Foam::Matrix<Form, Type>::col
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::block
 (
     const label m,
@@ -306,7 +218,7 @@ Foam::Matrix<Form, Type>::block
     const label nStart
 )
 {
-    return Block
+    return MatrixBlock<mType>
     (
         *this,
         m,
@@ -319,13 +231,13 @@ Foam::Matrix<Form, Type>::block
 
 template<class Form, class Type>
 template<class VectorSpace>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
 {
-    return Block
+    return MatrixBlock<mType>
     (
         *this,
-        VectorSpace::nRows,
+        VectorSpace::mRows,
         VectorSpace::nCols,
         mStart,
         nStart
@@ -334,10 +246,10 @@ Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::col(const label m, const label mStart)
 {
-    return Block
+    return MatrixBlock<mType>
     (
         *this,
         m,
@@ -349,7 +261,7 @@ Foam::Matrix<Form, Type>::col(const label m, const label mStart)
 
 
 template<class Form, class Type>
-inline typename Foam::Matrix<Form, Type>::Block
+inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
 Foam::Matrix<Form, Type>::col
 (
     const label m,
@@ -357,7 +269,7 @@ Foam::Matrix<Form, Type>::col
     const label nStart
 )
 {
-    return Block
+    return MatrixBlock<mType>
     (
         *this,
         m,
@@ -370,84 +282,6 @@ Foam::Matrix<Form, Type>::col
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-template<class Form, class Type>
-inline const Type& Foam::Matrix<Form, Type>::ConstBlock::operator()
-(
-    const label i,
-    const label j
-) const
-{
-    #ifdef FULLDEBUG
-    if (i<0 || i>=nRows_)
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << nRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
-    #endif
-
-    return matrix_(i + rowStart_, j + colStart_);
-}
-
-
-template<class Form, class Type>
-inline const Type& Foam::Matrix<Form, Type>::Block::operator()
-(
-    const label i,
-    const label j
-) const
-{
-    #ifdef FULLDEBUG
-    if (i<0 || i>=nRows_)
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << nRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
-    #endif
-
-    return matrix_(i + rowStart_, j + colStart_);
-}
-
-
-template<class Form, class Type>
-inline Type& Foam::Matrix<Form, Type>::Block::operator()
-(
-    const label i,
-    const label j
-)
-{
-    #ifdef FULLDEBUG
-    if (i<0 || i>=nRows_)
-    {
-        FatalErrorInFunction
-            << "Index " << i << " out of range 0 ... " << nRows_-1
-            << abort(FatalError);
-    }
-    if (j<0 || j>=nCols_)
-    {
-        FatalErrorInFunction
-            << "Index " << j << " out of range 0 ... " << nCols_-1
-            << abort(FatalError);
-    }
-    #endif
-
-    return matrix_(i + rowStart_, j + colStart_);
-}
-
-
 template<class Form, class Type>
 inline const Type& Foam::Matrix<Form, Type>::operator()
 (
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
index 9be178c12f..51943ada59 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -34,7 +34,7 @@ License
 template<class Form, class Type>
 Foam::Matrix<Form, Type>::Matrix(Istream& is)
 :
-    nRows_(0),
+    mRows_(0),
     nCols_(0),
     v_(NULL)
 {
@@ -59,10 +59,10 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
 
     if (firstToken.isLabel())
     {
-        M.nRows_ = firstToken.labelToken();
+        M.mRows_ = firstToken.labelToken();
         M.nCols_ = readLabel(is);
 
-        label mn = M.nRows_*M.nCols_;
+        label mn = M.mRows_*M.nCols_;
 
         // Read list contents depending on data format
         if (is.format() == IOstream::ASCII || !contiguous<Type>())
@@ -151,7 +151,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
 template<class Form, class Type>
 Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
 {
-    label mn = M.nRows_*M.nCols_;
+    label mn = M.mRows_*M.nCols_;
 
     os  << M.m() << token::SPACE << M.n();
 
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
new file mode 100644
index 0000000000..3aecbe8377
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.C
@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MatrixBlock.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::ConstMatrixBlock<MatrixType>::operator Field<cmptType>() const
+{
+    if (nCols_ != 1)
+    {
+        FatalErrorInFunction
+            << "Number of columns " << nCols_ << " != 1"
+            << abort(FatalError);
+    }
+
+    Field<cmptType> f(mRows_);
+
+    forAll(f, i)
+    {
+        f[i] = operator()(i, 0);
+    }
+
+    return f;
+}
+
+
+template<class MatrixType>
+Foam::MatrixBlock<MatrixType>::operator Field<cmptType>() const
+{
+    if (nCols_ != 1)
+    {
+        FatalErrorInFunction
+            << "Number of columns " << nCols_ << " != 1"
+            << abort(FatalError);
+    }
+
+    Field<cmptType> f(mRows_);
+
+    forAll(f, i)
+    {
+        f[i] = operator()(i, 0);
+    }
+
+    return f;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+template<class Form>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const Matrix<Form, cmptType>& Mb
+)
+{
+    if (mRows_ != Mb.m() || nCols_ != Mb.n())
+    {
+        FatalErrorInFunction
+            << "Attempt to assign blocks of different sizes: "
+            << mRows_ << "x" << nCols_ << " != "
+            << Mb.m() << "x" << Mb.n()
+            << abort(FatalError);
+    }
+
+    for (label i=0; i<mRows_; i++)
+    {
+        for (label j=0; j<nCols_; j++)
+        {
+            (*this)(i, j) = Mb(i, j);
+        }
+    }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const ConstMatrixBlock<MatrixType>& Mb
+)
+{
+    if (this != &Mb)
+    {
+        if (mRows_ != Mb.m() || nCols_ != Mb.n())
+        {
+            FatalErrorInFunction
+                << "Attempt to assign blocks of different sizes: "
+                << mRows_ << "x" << nCols_ << " != "
+                << Mb.m() << "x" << Mb.n()
+                << abort(FatalError);
+        }
+
+        for (label i=0; i<mRows_; i++)
+        {
+            for (label j=0; j<nCols_; j++)
+            {
+                (*this)(i, j) = Mb(i, j);
+            }
+        }
+    }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const MatrixBlock<MatrixType>& Mb
+)
+{
+    if (this != &Mb)
+    {
+        if (mRows_ != Mb.m() || nCols_ != Mb.n())
+        {
+            FatalErrorInFunction
+                << "Attempt to assign blocks of different sizes: "
+                << mRows_ << "x" << nCols_ << " != "
+                << Mb.m() << "x" << Mb.n()
+                << abort(FatalError);
+        }
+
+        for (label i=0; i<mRows_; i++)
+        {
+            for (label j=0; j<nCols_; j++)
+            {
+                (*this)(i, j) = Mb(i, j);
+            }
+        }
+    }
+}
+
+
+template<class MatrixType>
+template<class MatrixType2>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const ConstMatrixBlock<MatrixType2>& Mb
+)
+{
+    if (this != &Mb)
+    {
+        if (mRows_ != Mb.m() || nCols_ != Mb.n())
+        {
+            FatalErrorInFunction
+                << "Attempt to assign blocks of different sizes: "
+                << mRows_ << "x" << nCols_ << " != "
+                << Mb.m() << "x" << Mb.n()
+                << abort(FatalError);
+        }
+
+        for (label i=0; i<mRows_; i++)
+        {
+            for (label j=0; j<nCols_; j++)
+            {
+                (*this)(i, j) = Mb(i, j);
+            }
+        }
+    }
+}
+
+
+template<class MatrixType>
+template<class MatrixType2>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const MatrixBlock<MatrixType2>& Mb
+)
+{
+    if (this != &Mb)
+    {
+        if (mRows_ != Mb.m() || nCols_ != Mb.n())
+        {
+            FatalErrorInFunction
+                << "Attempt to assign blocks of different sizes: "
+                << mRows_ << "x" << nCols_ << " != "
+                << Mb.m() << "x" << Mb.n()
+                << abort(FatalError);
+        }
+
+        for (label i=0; i<mRows_; i++)
+        {
+            for (label j=0; j<nCols_; j++)
+            {
+                (*this)(i, j) = Mb(i, j);
+            }
+        }
+    }
+}
+
+
+template<class MatrixType>
+template
+<
+    template<class, Foam::direction, Foam::direction> class MSBlock,
+    class SubTensor,
+    Foam::direction BRowStart,
+    Foam::direction BColStart
+>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const MSBlock<SubTensor, BRowStart, BColStart>& Mb
+)
+{
+    if (mRows_ != Mb.mRows || nCols_ != Mb.nCols)
+    {
+        FatalErrorInFunction
+            << "Attempt to assign blocks of different sizes: "
+            << mRows_ << "x" << nCols_ << " != "
+            << Mb.mRows << "x" << Mb.nCols
+            << abort(FatalError);
+    }
+
+    for (direction i=0; i<mRows_; ++i)
+    {
+        for (direction j=0; j<nCols_; ++j)
+        {
+            operator()(i, j) = Mb(i, j);
+        }
+    }
+}
+
+
+template<class MatrixType>
+template
+<
+    template<class, Foam::direction> class VSBlock,
+    class SubVector,
+    Foam::direction BStart
+>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const VSBlock<SubVector, BStart>& Mb
+)
+{
+    if (mRows_ != Mb.nComponents || nCols_ != 1)
+    {
+        FatalErrorInFunction
+            << "Attempt to assign blocks of different sizes: "
+            << mRows_ << "x" << nCols_ << " != "
+            << Mb.nComponents << "x" << 1
+            << abort(FatalError);
+    }
+
+    for (direction i=0; i<mRows_; ++i)
+    {
+        operator()(i, 0) = Mb[i];
+    }
+}
+
+
+template<class MatrixType>
+template<class MSForm, Foam::direction Nrows, Foam::direction Ncols>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const MatrixSpace<MSForm, cmptType, Nrows, Ncols>& ms
+)
+{
+    if (mRows_ != Nrows || nCols_ != Ncols)
+    {
+        FatalErrorInFunction
+            << "Attempt to assign blocks of different sizes: "
+            << mRows_ << "x" << nCols_ << " != "
+            << Nrows << "x" << Ncols
+            << abort(FatalError);
+    }
+
+    for (label i=0; i<mRows_; i++)
+    {
+        for (label j=0; j<nCols_; j++)
+        {
+            (*this)(i, j) = ms(i, j);
+        }
+    }
+}
+
+
+template<class MatrixType>
+template<class VSForm, Foam::direction Ncmpts>
+void Foam::MatrixBlock<MatrixType>::operator=
+(
+    const VectorSpace<VSForm, cmptType, Ncmpts>& ms
+)
+{
+    if (mRows_ != Ncmpts || nCols_ != 1)
+    {
+        FatalErrorInFunction
+            << "Attempt to assign blocks of different sizes: "
+            << mRows_ << "x" << nCols_ << " != "
+            << Ncmpts << "x" << 1
+            << abort(FatalError);
+    }
+
+    for (direction i=0; i<Ncmpts; ++i)
+    {
+        operator()(i, 0) = ms[i];
+    }
+}
+
+
+template<class MatrixType>
+void Foam::MatrixBlock<MatrixType>::operator=(const Field<cmptType>& f)
+{
+    if (mRows_ != f.size() || nCols_ != 1)
+    {
+        FatalErrorInFunction
+            << "Error: cannot assign blocks of different size (left is "
+            << mRows_ << "x" << nCols_ << " != "
+            << f.size() << "x" << 1
+            << abort(FatalError);
+    }
+
+    forAll(f, i)
+    {
+        operator()(i, 0) = f[i];
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
new file mode 100644
index 0000000000..67ece198ee
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlock.H
@@ -0,0 +1,240 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::MatrixBlock
+
+Description
+    A templated block of an (m x n) matrix of type \<MatrixType\>.
+
+        Foam::ConstMatrixBlock: block of a const matrix
+        Foam::MatrixBlock: block of a non-const matrix
+
+    The block may be assigned to a block of another matrix or to a VectorSpace
+    or MatrixSpace e.g. \c tensor.  Conversion of a column block to a \c
+    Field<T> is also provide.
+
+SourceFiles
+    MatrixBlock.C
+    MatrixBlockI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MatrixBlock_H
+#define MatrixBlock_H
+
+#include "Matrix.H"
+#include "MatrixSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class ConstMatrixBlock Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MatrixType>
+class ConstMatrixBlock
+{
+    // Private data
+
+        //- Const reference to the parent matrix
+        const MatrixType& matrix_;
+
+        // Block size
+        const label mRows_;
+        const label nCols_;
+
+        // Block location in parent matrix
+        const label rowStart_;
+        const label colStart_;
+
+public:
+
+    typedef typename MatrixType::cmptType cmptType;
+
+    // Constructors
+
+        //- Construct block for matrix, size and location
+        inline ConstMatrixBlock
+        (
+            const MatrixType& matrix,
+            const label m,
+            const label n,
+            const label mStart,
+            const label nStart
+        );
+
+
+    // Member Functions
+
+        //- Return the number of rows in the block
+        inline label m() const;
+
+        //- Return the number of columns in the block
+        inline label n() const;
+
+        //- (i, j) const element access operator
+        inline const cmptType& operator()
+        (
+            const label i,
+            const label j
+        ) const;
+
+        //- Convert a column of a matrix to a Field
+        operator Field<cmptType>() const;
+};
+
+
+/*---------------------------------------------------------------------------*\
+                          Class MatrixBlock Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MatrixType>
+class MatrixBlock
+{
+    // Private data
+
+        //- Reference to the parent matrix
+        MatrixType& matrix_;
+
+        // Block size
+        const label mRows_;
+        const label nCols_;
+
+        // Block location in parent matrix
+        const label rowStart_;
+        const label colStart_;
+
+public:
+
+    typedef typename MatrixType::cmptType cmptType;
+
+    // Constructors
+
+        //- Construct block for matrix, size and location
+        inline MatrixBlock
+        (
+            MatrixType& matrix,
+            const label m,
+            const label n,
+            const label mStart,
+            const label nStart
+        );
+
+
+    // Member Functions
+
+        //- Return the number of rows in the block
+        inline label m() const;
+
+        //- Return the number of columns in the block
+        inline label n() const;
+
+        //- (i, j) const element access operator
+        inline const cmptType& operator()
+        (
+            const label i,
+            const label j
+        ) const;
+
+        //- (i, j) element access operator
+        inline cmptType& operator()(const label i, const label j);
+
+        //- Convert a column of a matrix to a Field
+        operator Field<cmptType>() const;
+
+
+    // Member operators
+
+        //- Assignment to a compatible matrix
+        template<class Form>
+        void operator=(const Matrix<Form, cmptType>&);
+
+        //- Assignment to a compatible const block
+        void operator=(const ConstMatrixBlock<MatrixType>&);
+
+        //- Assignment to a compatible block
+        void operator=(const MatrixBlock<MatrixType>&);
+
+        //- Assignment to a compatible const block
+        template<class MatrixType2>
+        void operator=(const ConstMatrixBlock<MatrixType2>&);
+
+        //- Assignment to a compatible block
+        template<class MatrixType2>
+        void operator=(const MatrixBlock<MatrixType2>&);
+
+        //- Assignment to a compatible MatrixSpace
+        template<class MSForm, direction Nrows, direction Ncols>
+        void operator=(const MatrixSpace<MSForm, cmptType, Nrows, Ncols>&);
+
+        //- Assignment to a compatible MatrixSpace block
+        template
+        <
+            template<class, direction, direction> class Block,
+            class SubTensor,
+            direction BRowStart,
+            direction BColStart
+        >
+        void operator=(const Block<SubTensor, BRowStart, BColStart>&);
+
+        //- Assignment to a compatible VectorSpace (column-vector)
+        template<class VSForm, direction Ncmpts>
+        void operator=(const VectorSpace<VSForm, cmptType, Ncmpts>&);
+
+        //- Assignment to a compatible VectorSpace (column-vector) block
+        template
+        <
+            template<class, direction> class Block,
+            class SubVector,
+            direction BStart
+        >
+        void operator=(const Block<SubVector, BStart>&);
+
+        //- Assignment to a Field (column-vector)
+        void operator=(const Field<cmptType>&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MatrixBlockI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "MatrixBlock.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
new file mode 100644
index 0000000000..e764dea032
--- /dev/null
+++ b/src/OpenFOAM/matrices/MatrixBlock/MatrixBlockI.H
@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MatrixType>
+Foam::ConstMatrixBlock<MatrixType>::ConstMatrixBlock
+(
+    const MatrixType& matrix,
+    const label m,
+    const label n,
+    const label mStart,
+    const label nStart
+)
+:
+    matrix_(matrix),
+    mRows_(m),
+    nCols_(n),
+    rowStart_(mStart),
+    colStart_(nStart)
+{
+    #ifdef FULLDEBUG
+    if
+    (
+        rowStart_ + mRows_ > matrix.m()
+     || colStart_ + nCols_ > matrix.n()
+    )
+    {
+        FatalErrorInFunction
+            << "Block addresses outside matrix"
+            << abort(FatalError);
+    }
+    #endif
+}
+
+
+template<class MatrixType>
+Foam::MatrixBlock<MatrixType>::MatrixBlock
+(
+    MatrixType& matrix,
+    const label m,
+    const label n,
+    const label mStart,
+    const label nStart
+)
+:
+    matrix_(matrix),
+    mRows_(m),
+    nCols_(n),
+    rowStart_(mStart),
+    colStart_(nStart)
+{
+    #ifdef FULLDEBUG
+    if
+    (
+        rowStart_ + mRows_ > matrix.m()
+     || colStart_ + nCols_ > matrix.n()
+    )
+    {
+        FatalErrorInFunction
+            << "Block addresses outside matrix"
+            << abort(FatalError);
+    }
+    #endif
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline Foam::label Foam::ConstMatrixBlock<MatrixType>::m() const
+{
+    return mRows_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::ConstMatrixBlock<MatrixType>::n() const
+{
+    return nCols_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::MatrixBlock<MatrixType>::m() const
+{
+    return mRows_;
+}
+
+
+template<class MatrixType>
+inline Foam::label Foam::MatrixBlock<MatrixType>::n() const
+{
+    return nCols_;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class MatrixType>
+inline const typename MatrixType::cmptType&
+Foam::ConstMatrixBlock<MatrixType>::operator()
+(
+    const label i,
+    const label j
+) const
+{
+    #ifdef FULLDEBUG
+    if (i<0 || i>=mRows_)
+    {
+        FatalErrorInFunction
+            << "Index " << i << " out of range 0 ... " << mRows_-1
+            << abort(FatalError);
+    }
+    if (j<0 || j>=nCols_)
+    {
+        FatalErrorInFunction
+            << "Index " << j << " out of range 0 ... " << nCols_-1
+            << abort(FatalError);
+    }
+    #endif
+
+    return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+template<class MatrixType>
+inline const typename MatrixType::cmptType&
+Foam::MatrixBlock<MatrixType>::operator()
+(
+    const label i,
+    const label j
+) const
+{
+    #ifdef FULLDEBUG
+    if (i<0 || i>=mRows_)
+    {
+        FatalErrorInFunction
+            << "Index " << i << " out of range 0 ... " << mRows_-1
+            << abort(FatalError);
+    }
+    if (j<0 || j>=nCols_)
+    {
+        FatalErrorInFunction
+            << "Index " << j << " out of range 0 ... " << nCols_-1
+            << abort(FatalError);
+    }
+    #endif
+
+    return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+template<class MatrixType>
+inline typename MatrixType::cmptType&
+Foam::MatrixBlock<MatrixType>::operator()
+(
+    const label i,
+    const label j
+)
+{
+    #ifdef FULLDEBUG
+    if (i<0 || i>=mRows_)
+    {
+        FatalErrorInFunction
+            << "Index " << i << " out of range 0 ... " << mRows_-1
+            << abort(FatalError);
+    }
+    if (j<0 || j>=nCols_)
+    {
+        FatalErrorInFunction
+            << "Index " << j << " out of range 0 ... " << nCols_-1
+            << abort(FatalError);
+    }
+    #endif
+
+    return matrix_(i + rowStart_, j + colStart_);
+}
+
+
+// ************************************************************************* //