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

@@ -32,13 +32,7 @@ void Foam::Matrix<Form, Type>::allocate()
 {
     if (nRows_ && nCols_)
     {
-        v_ = new Type*[nRows_];
+        v_ = new Type[size()];
-        v_[0] = new Type[nRows_*nCols_];
-
-        for (label i=1; i<nRows_; i++)
-        {
-            v_[i] = v_[i-1] + nCols_;
-        }
     }
 }
 
@@ -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),
+    nRows_(m),
-    nCols_(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),
+    nRows_(m),
-    nCols_(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];
+        const label mn = size();
-
-        label mn = nRows_*nCols_;
-
         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];
+        const label mn = size();
-
-        label mn = nRows_*nCols_;
         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 label mn = size();
-        const Type* av = a.v_[0];
-
-        label mn = nRows_*nCols_;
         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];
+        Type* nav = na.v();
-        const Type* av = a[0];
+        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];
+    Type* abv = ab.v();
-    const Type* av = a[0];
+    const Type* av = a.v();
-    const Type* bv = b[0];
+    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];
+    Type* abv = ab.v();
-    const Type* av = a[0];
+    const Type* av = a.v();
-    const Type* bv = b[0];
+    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];
+        Type* sav = sa.v();
-        const Type* av = a[0];
+        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];

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
+        //- Allocate the storage for the element vector
-        //  and set the row pointers
         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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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
+    checki(j);
-    return v_[i];
+    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*nCols_;
-    return v_[i];
+}
+
+
+template<class Form, class Type>
+inline Type* Foam::Matrix<Form, Type>::operator[](const label i)
+{
+    checki(i);
+    return v_ + i*nCols_;
 }
 
 

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));
         }
     }
 

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]);

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& i,
-            const direction& col
+            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

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& i,
-    const direction& col
+    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& i,
-    const direction& col
+    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];
 }