diff --git a/RandBLAS/base.hh b/RandBLAS/base.hh index 76280170..40f264b5 100644 --- a/RandBLAS/base.hh +++ b/RandBLAS/base.hh @@ -34,6 +34,7 @@ #include "RandBLAS/config.h" #include "RandBLAS/random_gen.hh" +#include "RandBLAS/exceptions.hh" #include #include @@ -105,6 +106,53 @@ inline submat_spec_64t offset_and_ldim( } } +template +void symmetrize(blas::Layout layout, blas::Uplo uplo, T* A, int64_t n, int64_t lda) { + + auto [inter_row_stride, inter_col_stride] = layout_to_strides(layout, lda); + #define matA(_i, _j) A[(_i)*inter_row_stride + (_j)*inter_col_stride] + if (uplo == blas::Uplo::Upper) { + // copy to lower + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i+1; j < n; ++j) { + matA(j,i) = matA(i,j); + } + } + } else if (uplo == blas::Uplo::Lower) { + // copy to upper + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i+1; j < n; ++j) { + matA(i,j) = matA(j,i); + } + } + } + #undef matA + return; +} + +template +void require_symmetric(blas::Layout layout, T* A, int64_t n, int64_t lda, T tol) { + + auto [inter_row_stride, inter_col_stride] = layout_to_strides(layout, lda); + #define matA(_i, _j) A[(_i)*inter_row_stride + (_j)*inter_col_stride] + + for (int64_t i = 0; i < n; ++i) { + for (int64_t j = i+1; j < n; ++j) { + T Aij = matA(i,j); + T Aji = matA(j,i); + T viol = abs(Aij - Aji); + T rel_tol = (abs(Aij) + abs(Aji) + 1)*tol; + if (viol > rel_tol) { + randblas_error_if_msg( + viol > rel_tol + "Symmetry check failed. |A(%i,%i) - A(%i,%i)| was %d, which exceeds tolerance of %d", i, j, j, i, viol, rel_tol + ); + } + } + } + #undef matA + return; +} template concept SignedInteger = (std::numeric_limits::is_signed && std::numeric_limits::is_integer); diff --git a/RandBLAS/sksy.hh b/RandBLAS/sksy.hh new file mode 100644 index 00000000..016a6ae3 --- /dev/null +++ b/RandBLAS/sksy.hh @@ -0,0 +1,99 @@ +// Copyright, 2024. See LICENSE for copyright holder information. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// (1) Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// (2) Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// (3) Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived from +// this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// + +#ifndef randblas_sksy_hh +#define randblas_sksy_hh + +#include "RandBLAS/base.hh" +#include "RandBLAS/skge.hh" + +namespace RandBLAS { + +using namespace RandBLAS::dense; +using namespace RandBLAS::sparse; + +/* Intended macro definitions. + + .. |op| mathmacro:: \operatorname{op} + .. |mat| mathmacro:: \operatorname{mat} + .. |submat| mathmacro:: \operatorname{submat} + .. |lda| mathmacro:: \texttt{lda} + .. |ldb| mathmacro:: \texttt{ldb} + .. |opA| mathmacro:: \texttt{opA} + .. |opS| mathmacro:: \texttt{opS} +*/ + +template +inline void sketch_symmetric( + // B = alpha*A*S + beta*B, where A is a symmetric matrix stored in the format of a general matrix. + blas::Layout layout, // layout for (A,B) + int64_t n, // number of rows in B + int64_t d, // number of columns in B + T alpha, + const T* A, + int64_t lda, + SKOP &S, + int64_t ro_s, + int64_t co_s, + T beta, + T* B, + int64_t ldb, + T sym_check_tol = -1 +) { + if (sym_check_tol >= 0) { + require_symmetric(layout, A, n, lda, sym_check_tol); + } + sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, n, d, n, alpha, A, lda, S, ro_s, co_s, beta, B, ldb); +} + +template +inline void sketch_symmetric( + // B = alpha*S*A + beta*B + blas::Layout layout, // layout for (A,B) + int64_t d, // number of rows in B + int64_t n, // number of columns in B + T alpha, + SKOP &S, + int64_t ro_s, + int64_t co_s, + const T* A, + int64_t lda, + T beta, + T* B, + int64_t ldb, + T sym_check_tol = -1 +) { + if (sym_check_tol >= 0) { + require_symmetric(layout, A, n, lda, sym_check_tol); + } + sketch_general(layout, blas::Op::NoTrans, blas::Op::NoTrans, d, n, n, alpha, S, ro_s, co_s, A, beta, B, ldb); +} + +} // end namespace RandBLAS +#endif