diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 6ccd1c79e5..808e00b86c 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -37,19 +37,7 @@ void Foam::Matrix<Form, Type>::allocate()
 }
 
 
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
-
-template<class Form, class Type>
-Foam::Matrix<Form, Type>::~Matrix()
-{
-    if (v_)
-    {
-        delete[] v_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Form, class Type>
 Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
@@ -69,6 +57,33 @@ 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),
+    nCols_(n),
+    v_(NULL)
+{
+    if (nRows_ < 0 || nCols_ < 0)
+    {
+        FatalErrorInFunction
+            << "Bad m, n " << nRows_ << ", " << nCols_
+            << abort(FatalError);
+    }
+
+    allocate();
+
+    if (v_)
+    {
+        const label mn = size();
+        for (label i=0; i<mn; i++)
+        {
+            v_[i] = Zero;
+        }
+    }
+}
+
+
 template<class Form, class Type>
 Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
 :
@@ -79,7 +94,7 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
     if (nRows_ < 0 || nCols_ < 0)
     {
         FatalErrorInFunction
-            << "bad m, n " << nRows_ << ", " << nCols_
+            << "Bad m, n " << nRows_ << ", " << nCols_
             << abort(FatalError);
     }
 
@@ -96,6 +111,24 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
 }
 
 
+template<class Form, class Type>
+inline Foam::Matrix<Form, Type>::Matrix(const mType::Block& block)
+:
+    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)
 :
@@ -116,6 +149,20 @@ Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
 }
 
 
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+template<class Form, class Type>
+Foam::Matrix<Form, Type>::~Matrix()
+{
+    if (v_)
+    {
+        delete[] v_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 template<class Form, class Type>
 void Foam::Matrix<Form, Type>::clear()
 {
@@ -146,6 +193,26 @@ void Foam::Matrix<Form, Type>::transfer(Matrix<Form, Type>& a)
 }
 
 
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::setSize(const label m, const label n)
+{
+    mType newMatrix(m, n, Zero);
+
+    label minM = min(m, nRows_);
+    label minN = min(n, nCols_);
+
+    for (label i=0; i<minM; i++)
+    {
+        for (label j=0; j<minN; j++)
+        {
+            newMatrix(i, j) = (*this)(i, j);
+        }
+    }
+
+    transfer(newMatrix);
+}
+
+
 template<class Form, class Type>
 Form Foam::Matrix<Form, Type>::T() const
 {
@@ -164,29 +231,57 @@ Form Foam::Matrix<Form, Type>::T() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
 template<class Form, class Type>
-void Foam::Matrix<Form, Type>::operator=(const Type& t)
+Foam::Matrix<Form, Type>::ConstBlock::operator Field<cmptType>() const
 {
-    if (v_)
+    if (nCols_ != 1)
     {
-        const label mn = size();
-        for (label i=0; i<mn; i++)
-        {
-            v_[i] = t;
-        }
+        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)
 {
     if (this == &a)
     {
         FatalErrorInFunction
-            << "attempted assignment to self"
+            << "Attempted assignment to self"
             << abort(FatalError);
     }
 
@@ -209,6 +304,247 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 }
 
 
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const mType::Block& block)
+{
+    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 Block& 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 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 j=0; j<nCols_; j++)
+        {
+            (*this)(i, j) = block(i, j);
+        }
+    }
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const Type& t)
+{
+    if (v_)
+    {
+        const label mn = size();
+        for (label i=0; i<mn; i++)
+        {
+            v_[i] = t;
+        }
+    }
+}
+
+
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::operator=(const zero)
+{
+    if (v_)
+    {
+        const label mn = size();
+        for (label i=0; i<mn; i++)
+        {
+            v_[i] = 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>
@@ -234,7 +570,7 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
     else
     {
         FatalErrorInFunction
-            << "matrix is empty"
+            << "Matrix is empty"
             << abort(FatalError);
 
         // Return in error to keep compiler happy
@@ -266,7 +602,7 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
     else
     {
         FatalErrorInFunction
-            << "matrix is empty"
+            << "Matrix is empty"
             << abort(FatalError);
 
         // Return in error to keep compiler happy
@@ -304,7 +640,7 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     if (a.m() != b.m())
     {
         FatalErrorInFunction
-            << "attempted add matrices with different number of rows: "
+            << "Attempt to add matrices with different number of rows: "
             << a.m() << ", " << b.m()
             << abort(FatalError);
     }
@@ -312,7 +648,7 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     if (a.n() != b.n())
     {
         FatalErrorInFunction
-            << "attempted add matrices with different number of columns: "
+            << "Attempt to add matrices with different number of columns: "
             << a.n() << ", " << b.n()
             << abort(FatalError);
     }
@@ -339,7 +675,7 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     if (a.m() != b.m())
     {
         FatalErrorInFunction
-            << "attempted add matrices with different number of rows: "
+            << "Attempt to add matrices with different number of rows: "
             << a.m() << ", " << b.m()
             << abort(FatalError);
     }
@@ -347,7 +683,7 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
     if (a.n() != b.n())
     {
         FatalErrorInFunction
-            << "attempted add matrices with different number of columns: "
+            << "Attempt to add matrices with different number of columns: "
             << a.n() << ", " << b.n()
             << abort(FatalError);
     }
@@ -389,13 +725,55 @@ Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 }
 
 
+template<class Form, class Type>
+Form Foam::operator*(const Matrix<Form, Type>& a, const scalar s)
+{
+    Form sa(a.m(), a.n());
+
+    if (a.m() && a.n())
+    {
+        Type* sav = sa.v();
+        const Type* av = a.v();
+
+        const label mn = a.size();
+        for (label i=0; i<mn; i++)
+        {
+            sav[i] = av[i]*s;
+        }
+    }
+
+    return sa;
+}
+
+
+template<class Form, class Type>
+Form Foam::operator/(const Matrix<Form, Type>& a, const scalar s)
+{
+    Form sa(a.m(), a.n());
+
+    if (a.m() && a.n())
+    {
+        Type* sav = sa.v();
+        const Type* av = a.v();
+
+        const label mn = a.size();
+        for (label i=0; i<mn; i++)
+        {
+            sav[i] = av[i]/s;
+        }
+    }
+
+    return sa;
+}
+
+
 template<class Form, class Type>
 Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
     if (a.n() != b.m())
     {
         FatalErrorInFunction
-            << "attempted to multiply incompatible matrices:" << nl
+            << "Attempt to multiply incompatible matrices:" << nl
             << "Matrix A : " << a.m() << " rows, " << a.n() << " columns" << nl
             << "Matrix B : " << b.m() << " rows, " << b.n() << " columns" << nl
             << "In order to multiply matrices, columns of A must equal "
@@ -405,13 +783,13 @@ Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 
     Form ab(a.m(), b.n(), scalar(0));
 
-    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 l = 0; l < b.m(); l++)
+            for (label k=0; k<b.m(); k++)
             {
-                ab(i, j) += a(i, l)*b(l, j);
+                ab(i, j) += a(i, k)*b(k, j);
             }
         }
     }
@@ -420,6 +798,39 @@ Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 }
 
 
+template<class Form, class Type>
+inline Foam::tmp<Foam::Field<Type>> Foam::operator*
+(
+    const Matrix<Form, Type>& M,
+    const Field<Type>& f
+)
+{
+    if (M.n() != f.size())
+    {
+        FatalErrorInFunction
+            << "Attempt to multiply incompatible matrix and field:" << nl
+            << "Matrix : " << M.m() << " rows, " << M.n() << " columns" << nl
+            << "Field : " << f.size() << " rows" << nl
+            << "In order to multiply a matrix M and field f, "
+               "columns of M must equal rows of f"
+            << abort(FatalError);
+    }
+
+    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 j=0; j<M.n(); ++j)
+        {
+            Mf[i] += M(i, j)*f[j];
+        }
+    }
+
+    return tMf;
+}
+
+
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
 #include "MatrixIO.C"
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 043b9561aa..5f905b465b 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -41,7 +41,9 @@ SourceFiles
 #include "bool.H"
 #include "label.H"
 #include "uLabel.H"
-#include "List.H"
+#include "Field.H"
+#include "MatrixSpace.H"
+#include "zero.H"
 #include "autoPtr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -87,10 +89,168 @@ class Matrix
 
 public:
 
+    //- Matrix type
+    typedef Matrix<Form, Type> mType;
+
+    //- Component type
+    typedef Type cmptType;
+
+
     // Static Member Functions
 
         //- Return a null Matrix
-        inline static const Matrix<Form, Type>& null();
+        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
@@ -102,17 +262,24 @@ public:
         Matrix(const label m, const label n);
 
         //- Construct with given number of rows and columns
-        //  and value for all elements.
+        //  initializing all elements to zero
+        Matrix(const label m, const label n, const zero);
+
+        //- Construct with given number of rows and columns
+        //  initializing all elements to the given value
         Matrix(const label m, const label n, const Type&);
 
         //- Copy constructor.
-        Matrix(const Matrix<Form, Type>&);
+        Matrix(const mType&);
+
+        //- Construct from a block of another matrix
+        Matrix(const mType::Block&);
 
         //- Construct from Istream.
         Matrix(Istream&);
 
         //- Clone
-        inline autoPtr<Matrix<Form, Type>> clone() const;
+        inline autoPtr<mType> clone() const;
 
 
     //- Destructor
@@ -132,6 +299,64 @@ public:
             //- Return the number of elements in matrix (m*n)
             inline label size() const;
 
+            //- Return element vector of the constant Matrix
+            inline const Type* v() const;
+
+            //- Return element vector of the Matrix
+            inline Type* v();
+
+
+        // Block access
+
+            inline ConstBlock block
+            (
+                const label m,
+                const label n,
+                const label mStart,
+                const label nStart
+            ) const;
+
+            template<class VectorSpace>
+            inline ConstBlock block
+            (
+                const label mStart,
+                const label nStart
+            ) const;
+
+            inline ConstBlock col
+            (
+                const label m,
+                const label rowStart
+            ) const;
+
+            inline ConstBlock col
+            (
+                const label m,
+                const label mStart,
+                const label nStart
+            ) const;
+
+
+            inline Block block
+            (
+                const label m,
+                const label n,
+                const label mStart,
+                const label nStart
+            );
+
+            template<class VectorSpace>
+            inline Block block(const label mStart, const label nStart);
+
+            inline Block col(const label m, const label rowStart);
+
+            inline Block col
+            (
+                const label m,
+                const label mStart,
+                const label nStart
+            );
+
 
         // Check
 
@@ -149,7 +374,10 @@ public:
 
             //- Transfer the contents of the argument Matrix into this Matrix
             //  and annul the argument Matrix.
-            void transfer(Matrix<Form, Type>&);
+            void transfer(mType&);
+
+            //- Resize the matrix preserving the elements
+            void setSize(const label m, const label nCols);
 
 
         //- Return the transpose of the matrix
@@ -158,12 +386,6 @@ public:
 
     // Member operators
 
-        //- Return element vector of the constant Matrix
-        inline const Type* v() const;
-
-        //- Return element vector of the Matrix
-        inline Type* v();
-
         //- Return subscript-checked row of Matrix.
         inline Type* operator[](const label);
 
@@ -177,7 +399,13 @@ public:
         inline Type& operator()(const label i, const label j);
 
         //- Assignment operator. Takes linear time.
-        void operator=(const Matrix<Form, Type>&);
+        void operator=(const mType&);
+
+        //- Assignment to a block of another matrix
+        void operator=(const mType::Block&);
+
+        //- Assignment of all entries to zero
+        void operator=(const zero);
 
         //- Assignment of all entries to the given value
         void operator=(const Type&);
@@ -189,49 +417,78 @@ public:
         friend Istream& operator>> <Form, Type>
         (
             Istream&,
-            Matrix<Form, Type>&
+            mType&
         );
 
         // Write Matrix to Ostream.
         friend Ostream& operator<< <Form, Type>
         (
             Ostream&,
-            const Matrix<Form, Type>&
+            const mType&
         );
 };
 
 
 // Global functions and operators
 
-template<class Form, class Type> const Type& max(const Matrix<Form, Type>&);
-template<class Form, class Type> const Type& min(const Matrix<Form, Type>&);
+template<class Form, class Type>
+const Type& max(const Matrix<Form, Type>&);
 
-template<class Form, class Type> Form operator-(const Matrix<Form, Type>&);
+template<class Form, class Type>
+const Type& min(const Matrix<Form, Type>&);
 
-template<class Form, class Type> Form operator+
+template<class Form, class Type>
+Form operator-(const Matrix<Form, Type>&);
+
+template<class Form, class Type>
+Form operator+
 (
     const Matrix<Form, Type>&,
     const Matrix<Form, Type>&
 );
 
-template<class Form, class Type> Form operator-
+template<class Form, class Type>
+Form operator-
 (
     const Matrix<Form, Type>&,
     const Matrix<Form, Type>&
 );
 
-template<class Form, class Type> Form operator*
+template<class Form, class Type>
+Form operator*
 (
     const scalar,
     const Matrix<Form, Type>&
 );
 
-template<class Form, class Type> Form operator*
+template<class Form, class Type>
+Form operator*
+(
+    const Matrix<Form, Type>&,
+    const scalar
+);
+
+template<class Form, class Type>
+Form operator/
+(
+    const Matrix<Form, Type>&,
+    const scalar
+);
+
+template<class Form, class Type>
+Form operator*
 (
     const Matrix<Form, Type>&,
     const Matrix<Form, Type>&
 );
 
+template<class Form, class Type>
+tmp<Field<Type>> operator*
+(
+    const Matrix<Form, Type>&,
+    const Field<Type>&
+);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 532e9df3f8..73992f3285 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -42,6 +42,68 @@ 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>
@@ -76,16 +138,16 @@ template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checki(const label i) const
 {
     #ifdef FULLDEBUG
-    if (!nRows_)
+    if (!nRows_ || !nCols_)
     {
         FatalErrorInFunction
-            << "attempt to access element from zero sized row"
+            << "Attempt to access element from empty matrix"
             << abort(FatalError);
     }
     else if (i<0 || i>=nRows_)
     {
         FatalErrorInFunction
-            << "index " << i << " out of range 0 ... " << nRows_-1
+            << "Index " << i << " out of range 0 ... " << nRows_-1
             << abort(FatalError);
     }
     #endif
@@ -96,10 +158,10 @@ template<class Form, class Type>
 inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 {
     #ifdef FULLDEBUG
-    if (!nCols_)
+    if (!nRows_ || !nCols_)
     {
         FatalErrorInFunction
-            << "attempt to access element from zero sized column"
+            << "Attempt to access element from empty matrix"
             << abort(FatalError);
     }
     else if (j<0 || j>=nCols_)
@@ -112,7 +174,33 @@ inline void Foam::Matrix<Form, Type>::checkj(const label j) const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+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
@@ -128,6 +216,238 @@ inline Type* Foam::Matrix<Form, Type>::v()
 }
 
 
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::ConstBlock
+Foam::Matrix<Form, Type>::block
+(
+    const label m,
+    const label n,
+    const label mStart,
+    const label nStart
+) const
+{
+    return ConstBlock
+    (
+        *this,
+        m,
+        n,
+        mStart,
+        nStart
+    );
+}
+
+
+template<class Form, class Type>
+template<class VectorSpace>
+inline typename Foam::Matrix<Form, Type>::ConstBlock
+Foam::Matrix<Form, Type>::block
+(
+    const label mStart,
+    const label nStart
+) const
+{
+    return ConstBlock
+    (
+        *this,
+        VectorSpace::nRows,
+        VectorSpace::nCols,
+        mStart,
+        nStart
+    );
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::ConstBlock
+Foam::Matrix<Form, Type>::col
+(
+    const label m,
+    const label mStart
+) const
+{
+    return ConstBlock
+    (
+        *this,
+        m,
+        1,
+        mStart,
+        0
+    );
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::ConstBlock
+Foam::Matrix<Form, Type>::col
+(
+    const label m,
+    const label mStart,
+    const label nStart
+) const
+{
+    return ConstBlock
+    (
+        *this,
+        m,
+        1,
+        mStart,
+        nStart
+    );
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::Block
+Foam::Matrix<Form, Type>::block
+(
+    const label m,
+    const label n,
+    const label mStart,
+    const label nStart
+)
+{
+    return Block
+    (
+        *this,
+        m,
+        n,
+        mStart,
+        nStart
+    );
+}
+
+
+template<class Form, class Type>
+template<class VectorSpace>
+inline typename Foam::Matrix<Form, Type>::Block
+Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
+{
+    return Block
+    (
+        *this,
+        VectorSpace::nRows,
+        VectorSpace::nCols,
+        mStart,
+        nStart
+    );
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::Block
+Foam::Matrix<Form, Type>::col(const label m, const label mStart)
+{
+    return Block
+    (
+        *this,
+        m,
+        1,
+        mStart,
+        0
+    );
+}
+
+
+template<class Form, class Type>
+inline typename Foam::Matrix<Form, Type>::Block
+Foam::Matrix<Form, Type>::col
+(
+    const label m,
+    const label mStart,
+    const label nStart
+)
+{
+    return Block
+    (
+        *this,
+        m,
+        1,
+        mStart,
+        nStart
+    );
+}
+
+
+// * * * * * * * * * * * * * * * 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()
 (