Commit 961dc048 authored by Henry Weller's avatar Henry Weller
Browse files

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.
parent 27ec0178
......@@ -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"
......
......@@ -41,7 +41,9 @@ SourceFiles
#include "bool.H"
#include "label.H"