Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ARM gemm #41

Merged
merged 4 commits into from
Oct 2, 2024
Merged

ARM gemm #41

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4ffbbf6
Added RegisterMatrixTest.testAxpy
Oct 1, 2024
94d903c
Fixed bug in RegisterMatrix::axpy()
mkatliar Oct 1, 2024
eb05f23
gemm() tests passing on ARM
Oct 1, 2024
0822b47
Removed superfluous comments
Oct 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion include/blast/math/algorithm/Gemm.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@ namespace blast
if (rows(D) != M || columns(D) != N)
BLAST_THROW_EXCEPTION(std::invalid_argument {"Matrix sizes do not match"});

gemm(M, N, K, alpha, ptr(*A), ptr(*B), beta, ptr(*C), ptr(*D));
gemm(M, N, K, alpha, ptr(A), ptr(B), beta, ptr(C), ptr(D));
}


Expand Down
102 changes: 102 additions & 0 deletions include/blast/math/reference/Axpy.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
// 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/TypeTraits.hpp>
#include <blast/math/algorithm/Tile.hpp>
#include <blast/math/register_matrix/Gemm.hpp>
#include <blast/math/Matrix.hpp>
#include <blast/util/Exception.hpp>


namespace blast :: reference
{
/**
* @brief Constant times matrix plus matrix, matrix pointer arguments
*
* C := alpha*A + B
*
* where alpha is a scalar, and A, B and C are M by N matrices
*
* @tparam ST scalar type for @a alpha
* @tparam MPA matrix pointer type for @a A
* @tparam MPB matrix pointer type for @a B
* @tparam MPC matrix pointer type for @a C
*
* @param M the number of rows of the matrices A, B, and C.
* @param N the number of columns of the matrices A, B and C.
* @param alpha the scalar alpha
* @param A pointer to the top left element of matrix A
* @param B pointer to the top left element of matrix B
* @param C pointer to the top left element of matrix C
*/
template <typename ST, MatrixPointer MPA, MatrixPointer MPB, MatrixPointer MPC>
inline void axpy(size_t M, size_t N, ST alpha, MPA A, MPB B, MPC C)
{
for (size_t i = 0; i < M; ++i)
for (size_t j = 0; j < N; ++j)
*(~C)(i, j) = alpha * *(~A)(i, j) + *(~B)(i, j);
}


/**
* @brief Constant times matrix plus matrix, matrix arguments
*
* C := alpha*A + B
*
* where alpha is a scalar, and A, B and C are M by N matrices
*
* @tparam ST scalar type for @a alpha
* @tparam MTA matrix type for @a A
* @tparam MTB matrix type for @a B
* @tparam MTC matrix type for @a C
*
* @param alpha the scalar alpha
* @param A matrix A
* @param B matrix B
* @param C matrix C
*
* @throw @a std::invalid_argument if matrix sizes are not consistent
*/
template <typename ST, Matrix MTA, Matrix MTB, Matrix MTC>
inline void axpy(ST alpha, MTA const& A, MTB const& B, MTC& C)
{
size_t const M = rows(C);
size_t const N = columns(C);

if (rows(A) != M || columns(A) != N ||
rows(B) != M || columns(B) != N)
BLAST_THROW_EXCEPTION(std::invalid_argument {"Inconsistent matrix sizes"});

reference::axpy(M, N, alpha, ptr(A), ptr(B), ptr(C));
}


/**
* @brief Constant times matrix plus matrix, matrix arguments, rvalue reference output argument
*
* C := alpha*A + B
*
* where alpha is a scalar, and A, B and C are M by N matrices
*
* @tparam ST scalar type for @a alpha
* @tparam MTA matrix type for @a A
* @tparam MTB matrix type for @a B
* @tparam MTC matrix view type for @a C
*
* @param alpha the scalar alpha
* @param A matrix A
* @param B matrix B
* @param C matrix view C
*
* @throw @a std::invalid_argument if matrix sizes are not consistent
*/
template <typename ST, Matrix MTA, Matrix MTB, Matrix MTC>
requires IsView_v<MTC>
inline void axpy(ST alpha, MTA const& A, MTB const& B, MTC&& C)
{
reference::axpy(alpha, A, B, C);
}
}
2 changes: 1 addition & 1 deletion include/blast/math/reference/Gemm.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,6 @@ namespace blast :: reference
if (rows(D) != M || columns(D) != N)
BLAST_THROW_EXCEPTION(std::invalid_argument {"Matrix sizes do not match"});

gemm(M, N, K, alpha, ptr(A), ptr(B), beta, ptr(C), ptr(D));
reference::gemm(M, N, K, alpha, ptr(A), ptr(B), beta, ptr(C), ptr(D));
}
}
2 changes: 1 addition & 1 deletion include/blast/math/register_matrix/RegisterMatrix.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ namespace blast
#pragma unroll
for (size_t j = 0; j < N; ++j) if (j < n)
#pragma unroll
for (size_t i = 0; i < RM; ++i) if (i * RM < m)
for (size_t i = 0; i < RM; ++i) if (SS * i < m)
v_[i][j] = fmadd(beta_simd, a(SS * i, j).load(), v_[i][j]);
}

Expand Down
49 changes: 28 additions & 21 deletions test/blast/math/dense/GemmTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,12 @@
#define BLAST_USER_ASSERTION 1

#include <blast/math/algorithm/Gemm.hpp>

#include <test/Testing.hpp>
#include <blast/math/algorithm/Randomize.hpp>
#include <blast/math/dense/DynamicMatrix.hpp>
#include <blast/math/views/Submatrix.hpp>
#include <blast/math/reference/Gemm.hpp>

#include <blast/blaze/Math.hpp>
#include <test/Testing.hpp>


namespace blast :: testing
Expand All @@ -29,22 +30,25 @@ namespace blast :: testing
for (size_t n = 1; n <= 20; n += 1)
for (size_t k = 1; k <= 20; ++k)
{
// Init Blaze matrices
// Init matrices
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Superfluous comment

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

//
blaze::DynamicMatrix<Real, SOA> A(m, k), C(m, n), D(m, n);
blaze::DynamicMatrix<Real, SOB> B(k, n);
DynamicMatrix<Real, SOA> A(m, k), C(m, n), D(m, n);
DynamicMatrix<Real, SOB> B(k, n);
randomize(A);
randomize(B);
randomize(C);

Real alpha {}, beta {};
blaze::randomize(alpha);
blaze::randomize(beta);
randomize(alpha);
randomize(beta);

/// Do gemm
// Do gemm
gemm(alpha, A, B, beta, C, D);

BLAST_ASSERT_APPROX_EQ(D, evaluate(beta * C + alpha * A * B), 1e-10, 1e-10)
DynamicMatrix<Real, columnMajor> D_ref(m, n);
reference::gemm(alpha, A, B, beta, C, D_ref);

BLAST_ASSERT_APPROX_EQ(D, D_ref, 1e-10, 1e-10)
<< "gemm error at size m,n,k=" << m << "," << n << "," << k;
}
}
Expand All @@ -54,31 +58,34 @@ namespace blast :: testing
void testUnalignedImpl()
{
size_t constexpr S_MAX = 20;
blaze::DynamicMatrix<Real, SOA> AA(S_MAX, S_MAX), CC(S_MAX, S_MAX), DD(S_MAX, S_MAX);
blaze::DynamicMatrix<Real, SOB> BB(S_MAX, S_MAX);
DynamicMatrix<Real, SOA> AA(S_MAX, S_MAX), CC(S_MAX, S_MAX), DD(S_MAX, S_MAX);
DynamicMatrix<Real, SOB> BB(S_MAX, S_MAX);

for (size_t m = 1; m <= S_MAX; m += 1)
for (size_t n = 1; n <= S_MAX; n += 1)
for (size_t k = 1; k <= S_MAX; ++k)
{
// Init Blaze matrices
// Init matrices
//
auto A = submatrix(AA, rows(AA) - m, columns(AA) - k, m, k);
auto C = submatrix(CC, rows(CC) - m, columns(CC) - n, m, n);
auto D = submatrix(DD, rows(DD) - m, columns(DD) - n, m, n);
auto B = submatrix(BB, rows(BB) - k, columns(BB) - n, k, n);
auto A = submatrix<unaligned>(AA, rows(AA) - m, columns(AA) - k, m, k);
auto C = submatrix<unaligned>(CC, rows(CC) - m, columns(CC) - n, m, n);
auto D = submatrix<unaligned>(DD, rows(DD) - m, columns(DD) - n, m, n);
auto B = submatrix<unaligned>(BB, rows(BB) - k, columns(BB) - n, k, n);
randomize(A);
randomize(B);
randomize(C);

Real alpha {}, beta {};
blaze::randomize(alpha);
blaze::randomize(beta);
randomize(alpha);
randomize(beta);

/// Do gemm
// Do gemm
gemm(alpha, A, B, beta, C, D);

BLAST_ASSERT_APPROX_EQ(D, evaluate(beta * C + alpha * A * B), 1e-10, 1e-10)
DynamicMatrix<Real, columnMajor> D_ref(m, n);
reference::gemm(alpha, A, B, beta, C, D_ref);

BLAST_ASSERT_APPROX_EQ(D, D_ref, 1e-10, 1e-10)
<< "gemm error at size m,n,k=" << m << "," << n << "," << k;
}
}
Expand Down
82 changes: 74 additions & 8 deletions test/blast/math/simd/RegisterMatrixTest.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include <blast/math/views/Submatrix.hpp>
#include <blast/math/reference/Ger.hpp>
#include <blast/math/reference/Gemm.hpp>
#include <blast/math/reference/Axpy.hpp>
#include <blast/math/StaticPanelMatrix.hpp>
#include <blast/math/dense/StaticMatrix.hpp>
#include <blast/math/expressions/MatTransExpr.hpp>
Expand All @@ -28,14 +29,18 @@ namespace blast :: testing


using MyTypes = Types<
RegisterMatrix<double, 4, 4, columnMajor>,
RegisterMatrix<double, 4, 2, columnMajor>,
RegisterMatrix<double, 4, 1, columnMajor>,
RegisterMatrix<double, 8, 4, columnMajor>,
RegisterMatrix<double, 12, 4, columnMajor>,
RegisterMatrix<float, 8, 4, columnMajor>,
RegisterMatrix<float, 16, 4, columnMajor>,
RegisterMatrix<float, 24, 4, columnMajor>
RegisterMatrix<double, 1 * SimdSize_v<double>, 4, columnMajor>,
RegisterMatrix<double, 1 * SimdSize_v<double>, 2, columnMajor>,
RegisterMatrix<double, 1 * SimdSize_v<double>, 1, columnMajor>,
RegisterMatrix<double, 2 * SimdSize_v<double>, 4, columnMajor>,
RegisterMatrix<double, 2 * SimdSize_v<double>, 2, columnMajor>,
RegisterMatrix<double, 2 * SimdSize_v<double>, 1, columnMajor>,
RegisterMatrix<double, 3 * SimdSize_v<double>, 4, columnMajor>,
RegisterMatrix<double, 3 * SimdSize_v<double>, 2, columnMajor>,
RegisterMatrix<double, 3 * SimdSize_v<double>, 1, columnMajor>,
RegisterMatrix<float, 1 * SimdSize_v<float>, 4, columnMajor>,
RegisterMatrix<float, 2 * SimdSize_v<float>, 4, columnMajor>,
RegisterMatrix<float, 3 * SimdSize_v<float>, 4, columnMajor>
>;


Expand Down Expand Up @@ -629,4 +634,65 @@ namespace blast :: testing
// TODO: should be strictly equal?
BLAST_ASSERT_APPROX_EQ(ker, alpha * B * A, absTol<ET>(), relTol<ET>());
}


TYPED_TEST(RegisterMatrixTest, testAxpy)
{
using RM = TypeParam;
using ET = ElementType_t<RM>;

StaticMatrix<ET, RM::rows(), RM::columns(), columnMajor> A, B;
randomize(A);
randomize(B);

ET alpha {};
randomize(alpha);

RM ker;
ker.load(ptr(B));
ker.axpy(alpha, ptr(A));

StaticMatrix<ET, RM::rows(), RM::columns(), columnMajor> C;
reference::axpy(alpha, A, B, C);

EXPECT_EQ(ker, C);
}


TYPED_TEST(RegisterMatrixTest, testAxpyWithSize)
{
using RM = TypeParam;
using ET = ElementType_t<RM>;

StaticMatrix<ET, RM::rows(), RM::columns(), columnMajor> A, B;
randomize(A);
randomize(B);

ET alpha {};
randomize(alpha);

StaticMatrix<ET, RM::rows(), RM::columns(), columnMajor> C;

for (size_t m = 0; m < RM::rows(); ++m)
{
for (size_t n = 0; n < RM::columns(); ++n)
{
RM ker;
ker.load(ptr(B));
ker.axpy(alpha, ptr(A), ker.rows(), ker.columns());

reference::axpy(alpha,
submatrix<aligned>(A, 0, 0, m, n),
submatrix<aligned>(B, 0, 0, m, n),
submatrix<aligned>(C, 0, 0, m, n)
);

for (size_t i = 0; i < m; ++i)
for (size_t j = 0; j < n; ++j)
ASSERT_EQ(ker(i, j), C(i, j))
<< "element mismatch at (" << i << ", " << j << "), "
<< "axpy() size = " << m << "x" << n;
}
}
}
}
Loading