From 961dc048da4b778be3fc9ede0b3405c7cbe84275 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 22 Mar 2016 14:09:24 +0000
Subject: [PATCH] Matrix: Added support for extracting and assigning blocks
The blocks may be specified directly in terms of the size and location in the
parent matrix or with the size obtained from a template specified
VectorSpace or MatrixSpace type.
---
src/OpenFOAM/matrices/Matrix/Matrix.C | 467 +++++++++++++++++++++++--
src/OpenFOAM/matrices/Matrix/Matrix.H | 301 ++++++++++++++--
src/OpenFOAM/matrices/Matrix/MatrixI.H | 332 +++++++++++++++++-
3 files changed, 1044 insertions(+), 56 deletions(-)
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 6ccd1c79e52..808e00b86c2 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -37,22 +37,28 @@ void Foam::Matrix<Form, Type>::allocate()
}
-// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Form, class Type>
-Foam::Matrix<Form, Type>::~Matrix()
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
+:
+ nRows_(m),
+ nCols_(n),
+ v_(NULL)
{
- if (v_)
+ if (nRows_ < 0 || nCols_ < 0)
{
- delete[] v_;
+ FatalErrorInFunction
+ << "Bad m, n " << nRows_ << ", " << nCols_
+ << abort(FatalError);
}
-}
+ allocate();
+}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const zero)
:
nRows_(m),
nCols_(n),
@@ -66,6 +72,15 @@ Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
}
allocate();
+
+ if (v_)
+ {
+ const label mn = size();
+ for (label i=0; i<mn; i++)
+ {
+ v_[i] = Zero;
+ }
+ }
}
@@ -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>
+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>
-void Foam::Matrix<Form, Type>::operator=(const Type& t)
+Foam::Matrix<Form, Type>::Block::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;
}
+// * * * * * * * * * * * * * * * 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 043b9561aa7..5f905b465bc 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>
+const Type& min(const Matrix<Form, Type>&);
-template<class Form, class Type> Form operator-(const Matrix<Form, Type>&);
+template<class Form, class Type>
+Form operator-(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 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 532e9df3f8d..73992f3285d 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()
(
--
GitLab