Commit 6f7c7d8c authored by Kutalmis Bercin's avatar Kutalmis Bercin Committed by Mark Olesen

ENH: partial overhaul of Matrix type (#1220)

- additional operators:
  + compound assignment
  + inner product: operator&
  + outer product: operator^

- additional functions:
   - MatrixBlock methods: subColumn, subRow, subMatrix
   - L2 norms for matrix or column
   - trace, diag, round, transpose

- MatrixBlock methods: col(), block() are deprecated since their
  access patterns with (size, offset) are unnatural/unwieldy.

- verifications by test/Matrix/Test-Matrix
parent 165ea4ed
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -128,40 +128,36 @@ inline bool Foam::Matrix<Form, Type>::empty() const noexcept
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checki(const label i) const
{
#ifdef FULLDEBUG
if (!mRows_ || !nCols_)
{
FatalErrorInFunction
<< "Attempt to access element from empty matrix"
<< abort(FatalError);
}
if (i < 0 || i >= mRows_)
if (i < 0 || mRows_ <= i)
{
FatalErrorInFunction
<< "Index " << i << " out of range 0 ... " << mRows_-1
<< abort(FatalError);
}
#endif
}
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checkj(const label j) const
{
#ifdef FULLDEBUG
if (!mRows_ || !nCols_)
{
FatalErrorInFunction
<< "Attempt to access element from empty matrix"
<< abort(FatalError);
}
if (j < 0 || j >= nCols_)
if (j < 0 || nCols_ <= j)
{
FatalErrorInFunction
<< "index " << j << " out of range 0 ... " << nCols_-1
<< abort(FatalError);
}
#endif
}
......@@ -238,7 +234,7 @@ template<class Form, class Type>
inline const Type& Foam::Matrix<Form, Type>::at(const label idx) const
{
#ifdef FULLDEBUG
if (idx < 0 || idx >= this->size())
if (idx < 0 || this->size() <= idx)
{
FatalErrorInFunction
<< "index " << idx << " out of range 0 ... " << this->size()
......@@ -253,7 +249,7 @@ template<class Form, class Type>
inline Type& Foam::Matrix<Form, Type>::at(const label idx)
{
#ifdef FULLDEBUG
if (idx < 0 || idx >= this->size())
if (idx < 0 || this->size() <= idx)
{
FatalErrorInFunction
<< "index " << idx << " out of range 0 ... " << this->size()
......@@ -266,152 +262,188 @@ inline Type& Foam::Matrix<Form, Type>::at(const label idx)
template<class Form, class Type>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
Foam::Matrix<Form, Type>::subColumn
(
const label m,
const label n,
const label mStart,
const label nStart
const label colIndex,
const label rowIndex,
label len
) const
{
if (len < 0)
{
len = mRows_ - rowIndex;
}
return ConstMatrixBlock<mType>
(
*this,
m,
n,
mStart,
nStart
len, // rows
1,
rowIndex,
colIndex
);
}
template<class Form, class Type>
template<class VectorSpace>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
Foam::Matrix<Form, Type>::subRow
(
const label mStart,
const label nStart
const label rowIndex,
const label colIndex,
label len
) const
{
if (len < 0)
{
len = nCols_ - colIndex;
}
return ConstMatrixBlock<mType>
(
*this,
VectorSpace::mRows,
VectorSpace::nCols,
mStart,
nStart
1,
len, // columns
rowIndex,
colIndex
);
}
template<class Form, class Type>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
Foam::Matrix<Form, Type>::subMatrix
(
const label m,
const label mStart
const label rowIndex,
const label colIndex,
label szRows,
label szCols
) const
{
if (szRows < 0) szRows = mRows_ - rowIndex;
if (szCols < 0) szCols = nCols_ - colIndex;
return ConstMatrixBlock<mType>
(
*this,
m,
1,
mStart,
0
szRows,
szCols,
rowIndex,
colIndex
);
}
template<class Form, class Type>
template<class VectorSpace>
inline Foam::ConstMatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
Foam::Matrix<Form, Type>::block
(
const label m,
const label mStart,
const label nStart
const label rowIndex,
const label colIndex
) const
{
return ConstMatrixBlock<mType>
(
*this,
m,
1,
mStart,
nStart
VectorSpace::mRows,
VectorSpace::nCols,
rowIndex,
colIndex
);
}
template<class Form, class Type>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block
Foam::Matrix<Form, Type>::subColumn
(
const label m,
const label n,
const label mStart,
const label nStart
const label colIndex,
const label rowIndex,
label len
)
{
if (len < 0)
{
len = mRows_ - rowIndex;
}
return MatrixBlock<mType>
(
*this,
m,
n,
mStart,
nStart
len, // rows
1,
rowIndex,
colIndex
);
}
template<class Form, class Type>
template<class VectorSpace>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::block(const label mStart, const label nStart)
Foam::Matrix<Form, Type>::subRow
(
const label rowIndex,
const label colIndex,
label len
)
{
if (len < 0)
{
len = nCols_ - colIndex;
}
return MatrixBlock<mType>
(
*this,
VectorSpace::mRows,
VectorSpace::nCols,
mStart,
nStart
1,
len, // columns
rowIndex,
colIndex
);
}
template<class Form, class Type>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col(const label m, const label mStart)
Foam::Matrix<Form, Type>::subMatrix
(
const label rowIndex,
const label colIndex,
label szRows,
label szCols
)
{
if (szRows < 0) szRows = mRows_ - rowIndex;
if (szCols < 0) szCols = nCols_ - colIndex;
return MatrixBlock<mType>
(
*this,
m,
1,
mStart,
0
szRows,
szCols,
rowIndex,
colIndex
);
}
template<class Form, class Type>
template<class VectorSpace>
inline Foam::MatrixBlock<Foam::Matrix<Form, Type>>
Foam::Matrix<Form, Type>::col
Foam::Matrix<Form, Type>::block
(
const label m,
const label mStart,
const label nStart
const label rowIndex,
const label colIndex
)
{
return MatrixBlock<mType>
(
*this,
m,
1,
mStart,
nStart
VectorSpace::mRows,
VectorSpace::nCols,
rowIndex,
colIndex
);
}
......@@ -532,8 +564,10 @@ inline const Type& Foam::Matrix<Form, Type>::operator()
const label jcol
) const
{
#ifdef FULLDEBUG
checki(irow);
checkj(jcol);
#endif
return v_[irow*nCols_ + jcol];
}
......@@ -545,8 +579,10 @@ inline Type& Foam::Matrix<Form, Type>::operator()
const label jcol
)
{
#ifdef FULLDEBUG
checki(irow);
checkj(jcol);
#endif
return v_[irow*nCols_ + jcol];
}
......@@ -554,7 +590,9 @@ inline Type& Foam::Matrix<Form, Type>::operator()
template<class Form, class Type>
inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
{
#ifdef FULLDEBUG
checki(irow);
#endif
return v_ + irow*nCols_;
}
......@@ -562,7 +600,9 @@ inline const Type* Foam::Matrix<Form, Type>::operator[](const label irow) const
template<class Form, class Type>
inline Type* Foam::Matrix<Form, Type>::operator[](const label irow)
{
#ifdef FULLDEBUG
checki(irow);
#endif
return v_ + irow*nCols_;
}
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ 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 "MatrixTools.H"
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * * //
template<class Form1, class Form2, class Type>
bool Foam::MatrixTools::equal
(
const Matrix<Form1, Type>& A,
const Matrix<Form2, Type>& B,
const bool verbose,
const scalar relTol,
const scalar absTol
)
{
const label len = A.size();
if (len != B.size())
{
if (verbose)
{
Info<< "Matrices have different sizes: "
<< len << " vs " << B.size() << nl;
}
return false;
}
auto iter1 = A.cbegin();
auto iter2 = B.cbegin();
for (label i = 0; i < len; ++i)
{
if ((absTol + relTol*mag(*iter2)) < Foam::mag(*iter1 - *iter2))
{
if (verbose)
{
Info<< "Matrix element " << i
<< " differs beyond tolerance: "
<< *iter1 << " vs " << *iter2 << nl;
}
return false;
}
++iter1;
++iter2;
}
if (verbose)
{
Info<< "All elements equal within the tolerances" << nl;
}
return true;
}
template<class Container>
Foam::Ostream& Foam::MatrixTools::printMatrix
(
Ostream& os,
const Container& mat
)
{
os << mat.m() << ' ' << mat.n();
if (mat.m() == 1)
{
// row-vector
os << " (";
for (label j = 0; j < mat.n(); ++j)
{
if (j) os << ' ';
os << mat(0,j);
}
os << ')' << nl;
}
else if (mat.n() == 1)
{
// col-vector
os << " (";
for (label i = 0; i < mat.m(); ++i)
{
if (i) os << ' ';
os << mat(i,0);
}
os << ')' << nl;
}
else
{
// Regular
os << nl << '(' << nl;
for (label i = 0; i < mat.m(); ++i)
{
os << '(';
for (label j = 0; j < mat.n(); ++j)
{
if (j) os << ' ';
os << mat(i,j);
}
os << ')' << nl;
}
os << ')' << nl;
}
return os;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 OpenCFD Ltd.
\\/ 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/>.
Namespace
Foam::MatrixTools
Description
Collection of static functions for matrix-related verifications.
SourceFiles
MatrixTools.C
\*---------------------------------------------------------------------------*/
#ifndef MatrixTools_H
#define MatrixTools_H
#include "Matrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
class Ostream;
/*---------------------------------------------------------------------------*\
Namespace MatrixTools Declaration
\*---------------------------------------------------------------------------*/
namespace MatrixTools
{
//- Compare matrix elements for absolute or relative equality
template<class Form1, class Form2, class Type>
bool equal
(
const Matrix<Form1, Type>& A,
const Matrix<Form2, Type>& B,
const bool verbose = false,
const scalar relTol = 1e-5,
const scalar absTol = 1e-8
);
//- Simple ASCII output of Matrix, MatrixBlock
template<class Container>
Ostream& printMatrix(Ostream& os, const Container& mat);
} // End namespace MatrixTools
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "MatrixTools.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd |
\\ / A nd | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
| Copyright (C) 2016 OpenFOAM Foundation
......@@ -71,6 +71,76 @@ Foam::MatrixBlock<MatrixType>::operator Field<cmptType>() const
}
template<class MatrixType>
Foam::label Foam::ConstMatrixBlock<MatrixType>::disallow
(
const char* what
) const
{
FatalErrorInFunction
<< "Block addresses " << what
<< " outside matrix or invalid matrix components"
<< abort(FatalError);
return 0;
}
template<class MatrixType>
Foam::label Foam::MatrixBlock<MatrixType>::disallow
(
const char* what
) const
{
FatalErrorInFunction
<< "Block addresses " << what
<< " outside matrix or invalid matrix components"
<< abort(FatalError);
return 0;
}
template<class MatrixType> void Foam::ConstMatrixBlock<MatrixType>::checkIndex
(
const label i,
const label j
) const
{
if (i < 0 || i >= mRows_)
{
FatalErrorInFunction
<< "Index " << i << " is out of range 0 ... " << mRows_ - 1
<< abort(FatalError);
}
else if (j < 0 || j >= nCols_)
{
FatalErrorInFunction
<< "Index " << j << " is out of range 0 ... " << nCols_ - 1
<< abort(FatalError);
}
}