search

Pressio / pressio

core C++ library
Other
45 stars 7 forks source link

ops: revise level2 kernels #446

Closed fnrizzi closed 1 year ago

fnrizzi commented 1 year ago
image
mzuzek commented 1 year ago

@fnrizzi I caught up with all_have_traits_and_same_scalar semantics and it looks all good - no changes required (for me, it could be simply called all_have_same_scalar - same like all_have_same_rank).

mzuzek commented 1 year ago

@fnrizzi

Eigen implementation notes

Considering y = beta * y + alpha * A * x that triggers Eigen matrix-vector product, scaling, vector summation and vector assignment:

  1. same scalar type for matrix A and vectors x and y is required - Eigen will protest with assertion from Eigen/src/Core/util/XprHelper.h:830-836: 
    // We require Lhs and Rhs to have "compatible" scalar types.
// It is tempting to always allow mixing different types but remember that this is often impossible in the vectorized paths.
// So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to
// add together a float matrix and a double matrix.
#define EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) \
 EIGEN_STATIC_ASSERT((Eigen::internal::has_ReturnType<ScalarBinaryOpTraits<LHS, RHS,BINOP> >::value), \
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)

    (this assertion is not always triggered properly so eventually a compilation may fail with error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<double, float, Eigen::internal::scalar_product_op<double, float> >’)

  2. custom scalar type requirements include that it can be constructed with numeric literal for the calls like Scalar(1) in Eigen/src/Core/ProductEvaluators.h:343-365 to succeed - for example, consider this code that aims to present scalar satisfying minimal Eigen requirements (and compiles and runs fine):

    struct scalar {
 scalar() = default;
 scalar(double v) {} // REQUIRED
 scalar& operator+=(const scalar &o) { return *this; }
};

scalar operator*(const scalar &a, const scalar &b) { return scalar{}; }

// not actually called at runtime but needed for build
scalar operator+(const scalar &a, const scalar &b) { assert("NOT ALLOWED"); return scalar{}; }
bool operator==(const scalar &a, const scalar &b)  { assert("NOT ALLOWED"); return true; }

void test() {
 using matrix = Eigen::Matrix<scalar, Eigen::Dynamic, Eigen::Dynamic>;
 using vector = Eigen::Matrix<scalar, Eigen::Dynamic, 1>;

 const int m = 4, n = 2;
 const auto alpha = scalar{};
 const auto beta = scalar{};
 const matrix A(m, n);
 const vector x(n);
 vector y(m);

 y = beta * y + alpha * A * x;
}

    as of v3.4.90, Eigen docs is missing on scalar types support, so it's hard to tell for sure whether that's a fundamental design rule or implementation detail - but I'd bet the former (low chances for support for scalar types without zero/one in near future).

fnrizzi commented 1 year ago

@mzuzek yes for now we should just constraint on same scalar types

mzuzek commented 1 year ago

@fnrizzi

Kokkos with custom scalar types

Kokkos implementation of matrix-vector product can be called like ::KokkosBlas::gemv("N", alpha, A, x, beta, y) on views with custom scalar types AScalar, XScalar, YScalar under following conditions:

Note: these conditions could be altered, but these alterations would be somewhat artificial (e.g. operator+ returning XType for YType arguments, instead of natural YType).

Here's the full code: 446_ops_level2_kokkos_scalar_types.zip

Important: these constraints can be changed in future Kokkos implementations

mzuzek commented 1 year ago

@fnrizzi It seems that for Teuchos matrix-vector product we have no direct unit tests at all and indirectly, only ::pressio::transpose overload (_ops/teuchos/opslevel2.hpp:89-118) is covered by:

mzuzek commented 1 year ago

Hi @fnrizzi! Would there be anything else to cover under this issue or can we close it ?