From 0ea0848047b6576ac6ae885a91991740416202d1 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sun, 20 Mar 2016 15:00:36 +0000
Subject: [PATCH] Matrix: Replace the row-start pointer array with computed
offsets
The row-start pointer array provided performance benefits on old
computers but now that computation is often cache-miss limited the
benefit of avoiding a integer multiply is more than offset by the
addition memory access into a separately allocated array.
With the new addressing scheme LUsolve is 15% faster.
---
src/OpenFOAM/matrices/Matrix/Matrix.C | 91 ++++++++-----------
src/OpenFOAM/matrices/Matrix/Matrix.H | 30 ++++--
src/OpenFOAM/matrices/Matrix/MatrixI.H | 55 +++++++++--
src/OpenFOAM/matrices/Matrix/MatrixIO.C | 8 +-
.../scalarMatrices/scalarMatricesTemplates.C | 2 -
.../primitives/MatrixSpace/MatrixSpace.H | 6 +-
.../primitives/MatrixSpace/MatrixSpaceI.H | 16 ++--
7 files changed, 120 insertions(+), 88 deletions(-)
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index ec9ae8b1b0f..829156693ca 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -32,13 +32,7 @@ void Foam::Matrix<Form, Type>::allocate()
{
if (nRows_ && nCols_)
{
- v_ = new Type*[nRows_];
- v_[0] = new Type[nRows_*nCols_];
-
- for (label i=1; i<nRows_; i++)
- {
- v_[i] = v_[i-1] + nCols_;
- }
+ v_ = new Type[size()];
}
}
@@ -50,7 +44,6 @@ Foam::Matrix<Form, Type>::~Matrix()
{
if (v_)
{
- delete[] (v_[0]);
delete[] v_;
}
}
@@ -59,16 +52,16 @@ Foam::Matrix<Form, Type>::~Matrix()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label n, const label m)
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n)
:
- nRows_(n),
- nCols_(m),
+ nRows_(m),
+ nCols_(n),
v_(NULL)
{
if (nRows_ < 0 || nCols_ < 0)
{
FatalErrorInFunction
- << "bad n, m " << nRows_ << ", " << nCols_
+ << "Bad m, n " << nRows_ << ", " << nCols_
<< abort(FatalError);
}
@@ -77,16 +70,16 @@ Foam::Matrix<Form, Type>::Matrix(const label n, const label m)
template<class Form, class Type>
-Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
+Foam::Matrix<Form, Type>::Matrix(const label m, const label n, const Type& a)
:
- nRows_(n),
- nCols_(m),
+ nRows_(m),
+ nCols_(n),
v_(NULL)
{
if (nRows_ < 0 || nCols_ < 0)
{
FatalErrorInFunction
- << "bad n, m " << nRows_ << ", " << nCols_
+ << "bad m, n " << nRows_ << ", " << nCols_
<< abort(FatalError);
}
@@ -94,13 +87,10 @@ Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
if (v_)
{
- Type* v = v_[0];
-
- label mn = nRows_*nCols_;
-
+ const label mn = size();
for (label i=0; i<mn; i++)
{
- v[i] = a;
+ v_[i] = a;
}
}
}
@@ -116,13 +106,11 @@ Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
if (a.v_)
{
allocate();
- Type* v = v_[0];
- const Type* av = a.v_[0];
- label mn = nRows_*nCols_;
+ const label mn = size();
for (label i=0; i<mn; i++)
{
- v[i] = av[i];
+ v_[i] = a.v_[i];
}
}
}
@@ -133,12 +121,12 @@ void Foam::Matrix<Form, Type>::clear()
{
if (v_)
{
- delete[] (v_[0]);
delete[] v_;
+ v_ = NULL;
}
+
nRows_ = 0;
nCols_ = 0;
- v_ = NULL;
}
@@ -183,12 +171,10 @@ void Foam::Matrix<Form, Type>::operator=(const Type& t)
{
if (v_)
{
- Type* v = v_[0];
-
- label mn = nRows_*nCols_;
+ const label mn = size();
for (label i=0; i<mn; i++)
{
- v[i] = t;
+ v_[i] = t;
}
}
}
@@ -214,13 +200,10 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
if (v_)
{
- Type* v = v_[0];
- const Type* av = a.v_[0];
-
- label mn = nRows_*nCols_;
+ const label mn = size();
for (label i=0; i<mn; i++)
{
- v[i] = av[i];
+ v_[i] = a.v_[i];
}
}
}
@@ -231,12 +214,12 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
template<class Form, class Type>
const Type& Foam::max(const Matrix<Form, Type>& a)
{
- label mn = a.m()*a.n();
+ const label mn = a.size();
if (mn)
{
label curMaxI = 0;
- const Type* v = a[0];
+ const Type* v = a.v();
for (label i=1; i<mn; i++)
{
@@ -263,12 +246,12 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
template<class Form, class Type>
const Type& Foam::min(const Matrix<Form, Type>& a)
{
- label mn = a.m()*a.n();
+ const label mn = a.size();
if (mn)
{
label curMinI = 0;
- const Type* v = a[0];
+ const Type* v = a.v();
for (label i=1; i<mn; i++)
{
@@ -301,10 +284,10 @@ Form Foam::operator-(const Matrix<Form, Type>& a)
if (a.m() && a.n())
{
- Type* nav = na[0];
- const Type* av = a[0];
+ Type* nav = na.v();
+ const Type* av = a.v();
- label mn = a.m()*a.n();
+ const label mn = a.size();
for (label i=0; i<mn; i++)
{
nav[i] = -av[i];
@@ -336,11 +319,11 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
Form ab(a.m(), a.n());
- Type* abv = ab[0];
- const Type* av = a[0];
- const Type* bv = b[0];
+ Type* abv = ab.v();
+ const Type* av = a.v();
+ const Type* bv = b.v();
- label mn = a.m()*a.n();
+ const label mn = a.size();
for (label i=0; i<mn; i++)
{
abv[i] = av[i] + bv[i];
@@ -371,11 +354,11 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
Form ab(a.m(), a.n());
- Type* abv = ab[0];
- const Type* av = a[0];
- const Type* bv = b[0];
+ Type* abv = ab.v();
+ const Type* av = a.v();
+ const Type* bv = b.v();
- label mn = a.m()*a.n();
+ const label mn = a.size();
for (label i=0; i<mn; i++)
{
abv[i] = av[i] - bv[i];
@@ -392,10 +375,10 @@ Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
if (a.m() && a.n())
{
- Type* sav = sa[0];
- const Type* av = a[0];
+ Type* sav = sa.v();
+ const Type* av = a.v();
- label mn = a.m()*a.n();
+ const label mn = a.size();
for (label i=0; i<mn; i++)
{
sav[i] = s*av[i];
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 49bf41fc6fc..043b9561aa7 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -79,10 +79,9 @@ class Matrix
label nRows_, nCols_;
//- Row pointers
- Type** __restrict__ v_;
+ Type* __restrict__ v_;
- //- Allocate the storage for the row-pointers and the data
- // and set the row pointers
+ //- Allocate the storage for the element vector
void allocate();
@@ -100,11 +99,11 @@ public:
inline Matrix();
//- Construct given number of rows and columns.
- Matrix(const label n, const label m);
+ Matrix(const label m, const label n);
//- Construct with given number of rows and columns
// and value for all elements.
- Matrix(const label n, const label m, const Type&);
+ Matrix(const label m, const label n, const Type&);
//- Copy constructor.
Matrix(const Matrix<Form, Type>&);
@@ -130,16 +129,16 @@ public:
//- Return the number of columns
inline label n() const;
- //- Return the number of elements in matrix (n*m)
+ //- Return the number of elements in matrix (m*n)
inline label size() const;
// Check
- //- Check index i is within valid range (0 ... n-1).
+ //- Check index i is within valid range (0 ... m-1).
inline void checki(const label i) const;
- //- Check index j is within valid range (0 ... m-1).
+ //- Check index j is within valid range (0 ... n-1).
inline void checkj(const label j) const;
@@ -159,12 +158,24 @@ 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);
//- Return subscript-checked row of constant Matrix.
inline const Type* operator[](const label) 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);
+
//- Assignment operator. Takes linear time.
void operator=(const Matrix<Form, Type>&);
@@ -221,13 +232,14 @@ template<class Form, class Type> Form operator*
const Matrix<Form, Type>&
);
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- #include "MatrixI.H"
+#include "MatrixI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index fbfd4462819..b9ecb4c744b 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -51,7 +51,6 @@ inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
}
-//- Return the number of rows
template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::m() const
{
@@ -76,6 +75,7 @@ inline Foam::label Foam::Matrix<Form, Type>::size() const
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checki(const label i) const
{
+ #ifdef FULLDEBUG
if (!nRows_)
{
FatalErrorInFunction
@@ -88,12 +88,14 @@ inline void Foam::Matrix<Form, Type>::checki(const label i) const
<< "index " << i << " out of range 0 ... " << nRows_-1
<< abort(FatalError);
}
+ #endif
}
template<class Form, class Type>
inline void Foam::Matrix<Form, Type>::checkj(const label j) const
{
+ #ifdef FULLDEBUG
if (!nCols_)
{
FatalErrorInFunction
@@ -106,28 +108,65 @@ inline void Foam::Matrix<Form, Type>::checkj(const label j) const
<< "index " << j << " out of range 0 ... " << nCols_-1
<< abort(FatalError);
}
+ #endif
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Form, class Type>
-inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
+inline const Type* Foam::Matrix<Form, Type>::v() const
+{
+ return v_;
+}
+
+
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::v()
+{
+ return v_;
+}
+
+
+template<class Form, class Type>
+inline const Type& Foam::Matrix<Form, Type>::operator()
+(
+ const label i,
+ const label j
+) const
{
- #ifdef FULLDEBUG
checki(i);
- #endif
- return v_[i];
+ checki(j);
+ return v_[i*nCols_ + j];
+}
+
+
+template<class Form, class Type>
+inline Type& Foam::Matrix<Form, Type>::operator()
+(
+ const label i,
+ const label j
+)
+{
+ checki(i);
+ checki(j);
+ return v_[i*nCols_ + j];
}
template<class Form, class Type>
inline const Type* Foam::Matrix<Form, Type>::operator[](const label i) const
{
- #ifdef FULLDEBUG
checki(i);
- #endif
- return v_[i];
+ return v_ + i*nCols_;
+}
+
+
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
+{
+ checki(i);
+ return v_ + i*nCols_;
}
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixIO.C b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
index edf2c54b22d..9be178c12f3 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixIO.C
+++ b/src/OpenFOAM/matrices/Matrix/MatrixIO.C
@@ -73,7 +73,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
if (mn)
{
M.allocate();
- Type* v = M.v_[0];
+ Type* v = M.v_;
if (listDelimiter == token::BEGIN_LIST)
{
@@ -124,7 +124,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
if (mn)
{
M.allocate();
- Type* v = M.v_[0];
+ Type* v = M.v_;
is.read(reinterpret_cast<char*>(v), mn*sizeof(Type));
@@ -162,7 +162,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
{
bool uniform = false;
- const Type* v = M.v_[0];
+ const Type* v = M.v_;
if (mn > 1 && contiguous<Type>())
{
@@ -248,7 +248,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
{
if (mn)
{
- os.write(reinterpret_cast<const char*>(M.v_[0]), mn*sizeof(Type));
+ os.write(reinterpret_cast<const char*>(M.v_), mn*sizeof(Type));
}
}
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
index 59ca2d34e5c..077fa60ced7 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
@@ -56,8 +56,6 @@ void Foam::solve
if (i != iMax)
{
- //Info<< "Pivoted on " << i << " " << iMax << endl;
-
for (label k=i; k<m; k++)
{
Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
diff --git a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
index 7d1b1186f4b..fd0beaca09d 100644
--- a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
+++ b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpace.H
@@ -279,12 +279,12 @@ public:
//- (i, j) const element access operator
inline const Cmpt& operator()
(
- const direction& row,
- const direction& col
+ const direction& i,
+ const direction& j
) const;
//- (i, j) element access operator
- inline Cmpt& operator()(const direction& row, const direction& col);
+ inline Cmpt& operator()(const direction& i, const direction& j);
// Member Operators
diff --git a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H
index a7afa55fdfc..3fca897851c 100644
--- a/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H
+++ b/src/OpenFOAM/primitives/MatrixSpace/MatrixSpaceI.H
@@ -322,12 +322,12 @@ Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::block()
template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
(
- const direction& row,
- const direction& col
+ const direction& i,
+ const direction& j
) const
{
#ifdef FULLDEBUG
- if (row < 0 || row > Nrows-1 || col < 0 || col > Ncols-1)
+ if (i < 0 || i > Nrows-1 || j < 0 || j > Ncols-1)
{
FatalErrorInFunction
<< "indices out of range"
@@ -335,19 +335,19 @@ inline const Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
}
#endif
- return this->v_[row*Ncols + col];
+ return this->v_[i*Ncols + j];
}
template<class Form, class Cmpt, Foam::direction Nrows, Foam::direction Ncols>
inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
(
- const direction& row,
- const direction& col
+ const direction& i,
+ const direction& j
)
{
#ifdef FULLDEBUG
- if (row < 0 || row > Nrows-1 || col < 0 || col > Ncols-1)
+ if (i < 0 || i > Nrows-1 || j < 0 || j > Ncols-1)
{
FatalErrorInFunction
<< "indices out of range"
@@ -355,7 +355,7 @@ inline Cmpt& Foam::MatrixSpace<Form, Cmpt, Nrows, Ncols>::operator()
}
#endif
- return this->v_[row*Ncols + col];
+ return this->v_[i*Ncols + j];
}
--
GitLab