From b601bbd5bc4a0e3b44be170ebec5021959497850 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miko=C5=82aj=20Zuzek?= Date: Fri, 27 May 2022 20:58:59 +0200 Subject: [PATCH 1/2] kokkos: fix reversed triangles in TRSM and TRSV --- tests/kokkos-based/triangular_matrix_matrix_left_solve.cpp | 6 +++--- tests/kokkos-based/triangular_matrix_matrix_right_solve.cpp | 6 +++--- tests/kokkos-based/triangular_matrix_vector_solve.cpp | 6 +++--- .../kokkos-kernels/blas3_triangular_matrix_matrix_solve.hpp | 3 +-- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/tests/kokkos-based/triangular_matrix_matrix_left_solve.cpp b/tests/kokkos-based/triangular_matrix_matrix_left_solve.cpp index 32802f17..86de4dec 100644 --- a/tests/kokkos-based/triangular_matrix_matrix_left_solve.cpp +++ b/tests/kokkos-based/triangular_matrix_matrix_left_solve.cpp @@ -66,10 +66,10 @@ void triangular_matrix_matrix_left_solve_gold_solution( for (size_type k = 0; k < num_vectors; ++k) { for (size_type ii = 0; ii < A_ext; ++ii) { - const size_type i = lower_triangle ? A_ext - 1 - ii : ii; + const size_type i = lower_triangle ? ii : A_ext - 1 - ii; // A(i, j) has lower triangle in i <= j - const size_type j0 = lower_triangle ? i + 1 : 0; - const size_type j1 = lower_triangle ? A_ext : i; + const size_type j0 = lower_triangle ? 0 : i + 1; + const size_type j1 = lower_triangle ? i : A_ext; for (size_type j = j0; j < j1; ++j) { X(i, k) -= A(i, j) * X(j, k); } diff --git a/tests/kokkos-based/triangular_matrix_matrix_right_solve.cpp b/tests/kokkos-based/triangular_matrix_matrix_right_solve.cpp index 9ec1969c..63f1690a 100644 --- a/tests/kokkos-based/triangular_matrix_matrix_right_solve.cpp +++ b/tests/kokkos-based/triangular_matrix_matrix_right_solve.cpp @@ -66,10 +66,10 @@ void triangular_matrix_matrix_right_solve_gold_solution( for (size_type i = 0; i < num_vectors; ++i) { for (size_type kk = 0; kk < A_ext; ++kk) { - const size_type k = lower_triangle ? kk : A_ext - 1 - kk; + const size_type k = lower_triangle ? A_ext - 1 - kk : kk; // A(j, k) has lower triangle in j <= k - const size_type j0 = lower_triangle ? 0 : k + 1; - const size_type j1 = lower_triangle ? k : A_ext; + const size_type j0 = lower_triangle ? k + 1 : 0; + const size_type j1 = lower_triangle ? A_ext : k; for (size_type j = j0; j < j1; ++j) { X(i, k) -= X(i, j) * A(j, k); } diff --git a/tests/kokkos-based/triangular_matrix_vector_solve.cpp b/tests/kokkos-based/triangular_matrix_vector_solve.cpp index 83def41e..d3667708 100644 --- a/tests/kokkos-based/triangular_matrix_vector_solve.cpp +++ b/tests/kokkos-based/triangular_matrix_vector_solve.cpp @@ -64,9 +64,9 @@ void triangular_matrix_vector_solve_gold_solution(A_t A, Triangle t, DiagonalSto const size_type ext1 = A.extent(1); for (size_type ii = 0; ii < ext0; ++ii) { - const size_type i = lower_triangle ? ext0 - 1 - ii : ii; - const size_type j0 = lower_triangle ? i + 1 : 0; - const size_type j1 = lower_triangle ? ext1 : i; + const size_type i = lower_triangle ? ii : ext0 - 1 - ii; + const size_type j0 = lower_triangle ? 0 : i + 1; + const size_type j1 = lower_triangle ? i : ext1; for (size_type j = j0; j < j1; ++j) { x(i) -= A(i, j) * x(j); } diff --git a/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_triangular_matrix_matrix_solve.hpp b/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_triangular_matrix_matrix_solve.hpp index 1ebf59e4..7ec0a7c0 100644 --- a/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_triangular_matrix_matrix_solve.hpp +++ b/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_triangular_matrix_matrix_solve.hpp @@ -63,9 +63,8 @@ void trsm(Side /*s*/, Triangle /*t*/, DiagonalStorage /*d*/, AViewType A, XViewT { const auto side = std::is_same_v ? "L" : "R"; - // KK and stdBLAS use REVERSED triangle definitions const auto uplo = std::is_same_v ? "U" : "L"; + std::experimental::linalg::lower_triangle_t> ? "L" : "U"; const auto diag = std::is_same_v ? "N" : "U"; From 7784ceed66b9feab122bf908034cbd3f83185b57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Miko=C5=82aj=20Zuzek?= Date: Fri, 27 May 2022 20:57:22 +0200 Subject: [PATCH 2/2] kokkos: fix reversed triangles in TRMM #226 --- .../triangular_matrix_left_product_kokkos.cpp | 6 +++--- .../triangular_matrix_right_product_kokkos.cpp | 6 +++--- .../kokkos-kernels/blas3_matrix_product_kk.hpp | 12 ++++++------ 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/kokkos-based/triangular_matrix_left_product_kokkos.cpp b/tests/kokkos-based/triangular_matrix_left_product_kokkos.cpp index e1179bc4..d5924d77 100644 --- a/tests/kokkos-based/triangular_matrix_left_product_kokkos.cpp +++ b/tests/kokkos-based/triangular_matrix_left_product_kokkos.cpp @@ -65,11 +65,11 @@ void updating_triangular_matrix_left_product_gold_solution(A_t A, // because (i,j) indexing respects C updating order // and parallelism is restricted accordingly. for (size_type ii = 0; ii < C_ext0; ++ii) { - const size_type i = lower ? ii : C_ext0 - 1 - ii; + const size_type i = lower ? C_ext0 - 1 - ii : ii; for (size_type j = 0; j < C_ext1; ++j) { c_element_type t = E(i, j); - const size_type k0 = lower ? (explicit_diag ? i : i + 1) : 0; - const size_type k1 = lower ? C_ext0 :(explicit_diag ? i + 1 : i); + const size_type k0 = lower ? 0 : (explicit_diag ? i : i + 1); + const size_type k1 = lower ? (explicit_diag ? i + 1 : i) : C_ext0; for (size_type k = k0; k < k1; ++k) { t += A(i, k) * B(k, j); } diff --git a/tests/kokkos-based/triangular_matrix_right_product_kokkos.cpp b/tests/kokkos-based/triangular_matrix_right_product_kokkos.cpp index 4a4e14cd..a0600f65 100644 --- a/tests/kokkos-based/triangular_matrix_right_product_kokkos.cpp +++ b/tests/kokkos-based/triangular_matrix_right_product_kokkos.cpp @@ -65,12 +65,12 @@ void updating_triangular_matrix_right_product_gold_solution(A_t A, // because (i,j) indexing respects C updating order // and parallelism is restricted accordingly. for (size_type jj = 0; jj < C_ext1; ++jj) { - const size_type j = lower ? C_ext1 - 1 - jj : jj; + const size_type j = lower ? jj : C_ext1 - 1 - jj; for (size_type i = 0; i < C_ext0; ++i) { c_element_type t = E(i, j); // Note: lower triangle of A(k, j) means k <= j - const auto k0 = lower ? 0 : (explicit_diag ? j : j + 1); - const auto k1 = lower ? (explicit_diag ? j + 1 : j) : C.extent(1); + const auto k0 = lower ? (explicit_diag ? j : j + 1) : 0; + const auto k1 = lower ? C_ext1 : (explicit_diag ? j + 1 : j); for (size_type k = k0; k < k1; ++k) { t += B(i, k) * A(k, j); } diff --git a/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_matrix_product_kk.hpp b/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_matrix_product_kk.hpp index a9ecad0b..2d17c466 100644 --- a/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_matrix_product_kk.hpp +++ b/tpl-implementations/include/experimental/__p1673_bits/kokkos-kernels/blas3_matrix_product_kk.hpp @@ -261,12 +261,12 @@ void trmm_left(KokkosExecSpace &&exec, Triangle t, DiagonalStorage d, // because (i,j) indexing respects C updating order // and parallelism is restricted accordingly for (size_type ii = 0; ii < C_ext0; ++ii) { - const auto i = lower ? ii : C_ext0 - 1 - ii; + const auto i = lower ? C_ext0 - 1 - ii : ii; Kokkos::parallel_for(Kokkos::RangePolicy(exec, 0, C_ext1), KOKKOS_LAMBDA (const size_type j) { c_element_type t{}; - const size_type k0 = lower ? (explicit_diag ? i : i + 1) : 0; - const size_type k1 = lower ? C_ext0 :(explicit_diag ? i + 1 : i); + const size_type k0 = lower ? 0 : (explicit_diag ? i : i + 1); + const size_type k1 = lower ? (explicit_diag ? i + 1 : i) : C_ext0; for (size_type k = k0; k < k1; ++k) { t += A_view(i, k) * B_view(k, j); } @@ -300,13 +300,13 @@ void trmm_right(KokkosExecSpace &&exec, Triangle t, DiagonalStorage d, // because (i,j) indexing respects C updating order // and parallelism is restricted accordingly for (size_type jj = 0; jj < C_ext1; ++jj) { - const size_type j = lower ? C_ext1 - 1 - jj : jj; + const size_type j = lower ? jj : C_ext1 - 1 - jj; Kokkos::parallel_for(Kokkos::RangePolicy(exec, 0, C_ext0), KOKKOS_LAMBDA (const size_type i) { c_element_type t{}; // Note: lower triangle of A(k, j) means k <= j - const auto k0 = lower ? 0 : (explicit_diag ? j : j + 1); - const auto k1 = lower ? (explicit_diag ? j + 1 : j) : C_ext1; + const auto k0 = lower ? (explicit_diag ? j : j + 1) : 0; + const auto k1 = lower ? C_ext1 : (explicit_diag ? j + 1 : j); for (size_type k = k0; k < k1; ++k) { t += B_view(i, k) * A_view(k, j); }