- Added blast::StaticMatrix class

- RegisterMatrixTest.testGer using reference implementation of ger()

include/blast/math/dense/StaticMatrix.hpp

@@ -0,0 +1,201 @@
+// Copyright (c) 2019-2020 Mikhail Katliar All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+#pragma once
+
+#include <blast/math/Forward.hpp>
+#include <blast/math/StorageOrder.hpp>
+#include <blast/math/Simd.hpp>
+#include <blast/math/TypeTraits.hpp>
+#include <blast/system/CacheLine.hpp>
+#include <blast/util/NextMultiple.hpp>
+#include <blast/util/Types.hpp>
+
+#include <initializer_list>
+
+
+namespace blast
+{
+    /// @brief Matrix with statically defined size.
+    ///
+    /// @tparam T element type of the matrix
+    /// @tparam M number of rows
+    /// @tparam N number of columns
+    /// @tparam SO storage order
+    template <typename T, size_t M, size_t N, bool SO>
+    class StaticMatrix
+    {
+    public:
+        using ElementType = T;
+        static bool constexpr storageOrder = SO;
+
+
+        StaticMatrix() noexcept
+        {
+            // Initialize padding elements to 0 to prevent denorms in calculations.
+            // Denorms can significantly impair performance, see https://github.com/giaf/blasfeo/issues/103
+            std::fill_n(v_, capacity_, T {});
+        }
+
+
+        StaticMatrix(T const& v) noexcept
+        {
+            std::fill_n(v_, capacity_, v);
+        }
+
+
+        /**
+         * @brief Construct from an initializer list.
+         *
+         * The initializer list is row-major regardless of the @a RegsterMatrix storage order,
+         * such that the initializer elements are written in a natural way.
+         *
+         * \code
+         * StaticMatrix<double, 2, 3, columnMajor> m {
+         *     {1., 2., 3.},
+         *     {4., 5., 6.},
+         *     {7., 8., 9.}
+         * };
+         * \endcode
+         *
+         * @param list list of lists of matrix row elements
+         *
+         */
+        constexpr StaticMatrix(std::initializer_list<std::initializer_list<T>> list)
+        {
+            std::fill_n(v_, capacity_, T {});
+
+            if (list.size() != M || determineColumns(list) > N)
+                throw std::invalid_argument {"Invalid setup of static matrix"};
+
+            size_t i = 0;
+
+            for (auto const& row : list)
+            {
+                size_t j = 0;
+
+                for (const auto& element : row)
+                {
+                    v_[elementIndex(i, j)] = element;
+                    ++j;
+                }
+
+                ++i;
+            }
+        }
+
+
+        StaticMatrix& operator=(T val) noexcept
+        {
+            for (size_t i = 0; i < M; ++i)
+                for (size_t j = 0; j < N; ++j)
+                    (*this)(i, j) = val;
+
+            return *this;
+        }
+
+
+        constexpr T const& operator()(size_t i, size_t j) const noexcept
+        {
+            return v_[elementIndex(i, j)];
+        }
+
+
+        constexpr T& operator()(size_t i, size_t j)
+        {
+            return v_[elementIndex(i, j)];
+        }
+
+
+        static size_t constexpr rows() noexcept
+        {
+            return M;
+        }
+
+
+        static size_t constexpr columns() noexcept
+        {
+            return N;
+        }
+
+
+        static size_t constexpr spacing() noexcept
+        {
+            return spacing_;
+        }
+
+
+        T * data() noexcept
+        {
+            return v_;
+        }
+
+
+        T const * data() const noexcept
+        {
+            return v_;
+        }
+
+
+    private:
+        static size_t constexpr spacing_ = nextMultiple(SO == columnMajor ? M : N, SimdSize_v<T>);
+        static size_t constexpr capacity_ = spacing_ * (SO == columnMajor ? N : M);
+
+        // Alignment of the data elements.
+        static size_t constexpr alignment_ = CACHE_LINE_SIZE;
+
+        // Aligned element storage.
+        alignas(alignment_) T v_[capacity_];
+
+
+        size_t elementIndex(size_t i, size_t j) const
+        {
+            return SO == columnMajor ? i + spacing_ * j : spacing_ * i + j;
+        }
+    };
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    inline size_t constexpr rows(StaticMatrix<T, M, N, SO> const& m) noexcept
+    {
+        return m.rows();
+    }
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    inline size_t constexpr columns(StaticMatrix<T, M, N, SO> const& m) noexcept
+    {
+        return m.columns();
+    }
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    inline constexpr T * data(StaticMatrix<T, M, N, SO>& m) noexcept
+    {
+        return m.data();
+    }
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    inline constexpr T const * data(StaticMatrix<T, M, N, SO> const& m) noexcept
+    {
+        return m.data();
+    }
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    struct IsDenseMatrix<StaticMatrix<T, M, N, SO>> : std::true_type {};
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    struct IsStatic<StaticMatrix<T, M, N, SO>> : std::true_type {};
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    struct Spacing<StaticMatrix<T, M, N, SO>> : std::integral_constant<size_t, StaticMatrix<T, M, N, SO>::spacing()> {};
+
+
+    template <typename T, size_t M, size_t N, bool SO>
+    struct StorageOrderHelper<StaticMatrix<T, M, N, SO>> : std::integral_constant<StorageOrder, StorageOrder(SO)> {};
+}

include/blast/math/reference/Ger.hpp

@@ -0,0 +1,40 @@
+// Copyright (c) 2019-2024 Mikhail Katliar All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+#pragma once
+
+#include <blast/math/typetraits/MatrixPointer.hpp>
+#include <blast/math/typetraits/VectorPointer.hpp>
+#include <blast/util/Types.hpp>
+
+
+namespace blast :: reference
+{
+    /**
+     * @brief Reference implementation of rank-1 update with multiplier
+     *
+     * a(i, j) += alpha * x(i) * y(j)
+     * for i=0...m-1, j=n-1
+     *
+     * @tparam Real real number type
+     * @tparam VPX vector pointer type for the column vector @a x
+     * @tparam VPY vector pointer type for the row vector @a y
+     * @tparam MPA
+     *
+     * @param m number of rows in the matrix
+     * @param n number of columns in the matrix
+     * @param alpha scalar multiplier
+     * @param x column vector
+     * @param y row vector
+     * @param a matrix to perform update on
+     *
+     */
+    template <typename Real, typename VPX, typename VPY, typename MPA>
+    requires VectorPointer<VPX, Real> && VectorPointer<VPY, Real> && MatrixPointer<MPA, Real>
+    inline void ger(size_t m, size_t n, Real alpha, VPX x, VPY y, MPA a)
+    {
+        for (size_t i = 0; i < m; ++i)
+            for (size_t j = 0; j < n; ++j)
+                *(~a)(i, j) += alpha * *(~x)(i) * *(~y)(j);
+    }
+}

include/blast/util/NextMultiple.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <blast/util/Types.hpp>
+
+
+namespace blast
+{
+    /**
+     * @brief Rounds up an integral value to the next multiple of a given factor.
+    //
+    // @param value The integral value to be rounded up \f$[1..\infty)\f$.
+    // @param factor The factor of the multiple \f$[1..\infty)\f$.
+    // @return The next multiple of the given factor.
+    //
+    // This function rounds up the given integral value to the next multiple of the given integral
+    // factor. In case the integral value is already a multiple of the given factor, the value itself
+    // is returned.
+     */
+    inline size_t constexpr nextMultiple(size_t value, size_t factor) noexcept
+    {
+        return value + (factor - value % factor) % factor;
+    }
+}

test/blast/math/simd/RegisterMatrixTest.cpp

@@ -3,11 +3,11 @@
 // license that can be found in the LICENSE file.
 
 #include <blast/math/RegisterMatrix.hpp>
-#include <blast/math/StaticPanelMatrix.hpp>
 #include <blast/math/Matrix.hpp>
 #include <blast/math/Vector.hpp>
 #include <blast/math/views/submatrix/Panel.hpp>
-#include <blast/blaze/Math.hpp>
+#include <blast/math/reference/Ger.hpp>
+#include <blast/math/panel/StaticPanelMatrix.hpp>
 
 #include <test/Testing.hpp>
 #include <test/Randomize.hpp>
@@ -131,7 +131,7 @@ namespace blast :: testing
         StaticMatrix<ET, Traits::rows, Traits::columns, columnMajor> A;
         for (size_t i = 0; i < rows(A); ++i)
             for (size_t j = 0; j < columns(A); ++j)
-                (*A)(i, j) = 1000 * i + j;
+                A(i, j) = 1000 * i + j;
 
         for (size_t m = 1; m <= rows(A); ++m)
         {
@@ -179,7 +179,7 @@ namespace blast :: testing
         using ET = ElementType_t<RM>;
 
         StaticMatrix<ET, Traits::rows, Traits::columns, columnMajor> A, B(0.);
-        randomize(A);
+        blast::randomize(A);
 
         RM ker;
         ker.load(1., ptr<aligned>(A, 0, 0));
@@ -374,10 +374,12 @@ namespace blast :: testing
         blaze::randomize(alpha);
 
         TypeParam ker;
-        ker.load(1., ptr(*C));
+        ker.load(1., ptr(C));
         ker.ger(alpha, ptr(a), ptr(b));
 
-        BLAST_EXPECT_APPROX_EQ(ker, evaluate(C + alpha * a * b), absTol<ET>(), relTol<ET>());
+        reference::ger(rows(C), columns(C), alpha, ptr(a), ptr(b), ptr(C));
+
+        BLAST_EXPECT_APPROX_EQ(ker, C, absTol<ET>(), relTol<ET>());
     }
 
 
@@ -443,7 +445,7 @@ namespace blast :: testing
             for (size_t n = 0; n <= columns(C); ++n)
             {
                 TypeParam ker;
-                ker.load(1., ptr(*C));
+                ker.load(1., ptr(C));
                 ker.ger(ET(1.), ptr(a), ptr(trans(b)), m, n);
 
                 for (size_t i = 0; i < m; ++i)
@@ -471,9 +473,9 @@ namespace blast :: testing
         randomize(C);
 
         TypeParam ker;
-        ker.load(1., ptr(*C));
+        ker.load(1., ptr(C));
         ker.ger(ET(1.), ptr(a), ptr(trans(b)));
-        ker.store(ptr(*D));
+        ker.store(ptr(D));
 
         BLAST_EXPECT_EQ(D, evaluate(C + a * trans(b)));
     }