Overhaul of multiplication mechanism

[inakleinbottle]

Currently the multiplication in an algebra is handled by passing the keys to the basis prod member function, which computes the product of keys and returns the result. (This is not actually the case for free_tensor objects. More on this later.) Unfortunately, the product of two basis keys need not be a key itself, and can instead by an arbitrary vector (i.e. a vector). Thus the basis needs to have knowledge of the algebra class that it should return. This, obviously, leads to a few problems:

• The basis class is required, as a template argument, to instantiate the vector (algebra) class, and so this require some careful declarations of classes to make everything work;
• The basis class has to know about coefficients, which it otherwise wouldn't;
• It makes it rather more difficult to separate the logic for densely stored vectors vs sparsely stored vectors.

In the most recent version of libalgebra, this third issue was partially solved for tensors by changing the low-level implementation of multiplication to allow two separate operations to be used on dense and sparse parts of vectors.
For non-tensors, the multiplication remained the same except for a small shim class inserted between which generated a stand-in function object and passed that to the vector mechanism. This was never intended to be a complete solution.

Instead, I think multiplication should be provided by a additional class which is also supplied as a template argument to the algebra class and is stored as a static member of this class. This has several benefits:

• the class implementing multiplication does not need knowledge of the underlying vector type, since it can be derived as a template argument;
• it doesn't even need knowledge of the basis itself, only the key type and the coefficient type that could even be derived from template parameters;
• it gives us complete control over the mechanism used to compute algebra products, we can simply implement different function objects and pass them into the underlying vector mechanics;
• since vector objects are not tied directly to a specific depth, where such a parameter exists, we can more easily implement cross-depth multiplication.

I've coded up tensor multiplication as a template of the kind of thing we need to write for each of the algebra derivatives in libalgebra.

template <typename Coeff>
class free_tensor_multiplication {

    typedef typename Coeff::SCA scalar_t;

    template <typename Transform>
    class index_operator {
        Transform m_transform;

    public:

        index_operator(Transform t) : m_transform(t) {}

        void operator()(scalar_t *result_ptr, scalar_t const *lhs_ptr, scalar_t const *rhs_ptr, DIMN const lhs_target,
                        DIMN const rhs_target, bool assign = false) {
            scalar_t lhs;
            if (assign) {
                for (IDIMN i = 0; i < static_cast<IDIMN>(lhs_target); ++i) {
                    lhs = lhs_ptr[i];
                    for (IDIMN j = 0; j < static_cast<IDIMN>(rhs_target); ++j) {
                        *(result_ptr++) = m_transform(Coeff::mul(lhs, rhs_ptr[j]));
                    }
                }
            } else {
                for (IDIMN i = 0; i < static_cast<IDIMN>(lhs_target); ++i) {
                    lhs = lhs_ptr[i];
                    for (IDIMN j = 0; j < static_cast<IDIMN>(rhs_target); ++j) {
                        *(result_ptr++) += m_transform(Coeff::mul(lhs, rhs_ptr[j]));
                    }
                }
            }
        }
    };

    template <typename Transform>
    class key_operator {
        Transform m_transform;

    public:

        key_operator(Transform t) : m_transform(t) {}

        template <typename Vector>
        void operator()(Vector &result, typename Vector::KEY const &lhs_key, scalar_t const &lhs_val,
                        typename Vector::KEY const &rhs_key, scalar_t const &rhs_val) {
            result.add_scal_prod(lhs_key * rhs_key, m_transform(Coeff::mul(lhs_val, rhs_val)));
        }
    };


public:

    template <typename Algebra, typename Operator>
    Algebra &multiply_and_add(Algebra &result, Algebra const &lhs, Algebra const &rhs, Operator op) const {
        key_operator <Operator> kt(op);
        index_operator<Operator> it(op);
        lhs.buffered_apply_binary_transform(result, rhs, kt, it);
        return result;
    }

    template <typename Algebra, typename Operator>
    Algebra &
    multiply_and_add(Algebra &result, Algebra const &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const {
        key_operator <Operator> kt(op);
        index_operator<Operator> it(op);
        lhs.buffered_apply_binary_transform(result, rhs, kt, it, max_depth);
        return result;
    }

    template <typename Algebra, typename Operator>
    Algebra multiply(Algebra const &lhs, Algebra const &rhs, Operator op) const {
        Algebra result;
        multiply_and_add(result, lhs, rhs, op);
        return result;
    }

    template <typename Algebra, typename Operator>
    Algebra multiply(Algebra const &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const {
        Algebra result;
        multiply_and_add(result, lhs, rhs, op, max_depth);
        return result;
    }

    template <typename Algebra, typename Operator>
    Algebra &multiply_inplace(Algebra &lhs, Algebra const &rhs, Operator op) const {
        key_operator <Operator> kt(op);
        index_operator<Operator> it(op);
        lhs.unbuffered_apply_binary_transform(rhs, kt, it);
        return lhs;
    }

    template <typename Algebra, typename Operator>
    Algebra &multiply_inplace(Algebra &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const {
        key_operator <Operator> kt(op);
        index_operator<Operator> it(op);
        lhs.unbuffered_apply_binary_transform(rhs, kt, it, max_depth);
        return lhs;
    }

};

[terrylyons]

I agree with this. Moreover, I think this is quite elegant. Providing that there is no performance hit, I can see many practical benefits. From: inakleinbottle ***@***.***> Sent: 05 July 2021 17:23 To: terrylyons/libalgebra ***@***.***> Cc: Subscribed ***@***.***> Subject: [terrylyons/libalgebra] Overhaul of multiplication mechanism (#19) Currently the multiplication in an algebra is handled by passing the keys to the basis prod member function, which computes the product of keys and returns the result. (This is not actually the case for free_tensor objects. More on this later.) Unfortunately, the product of two basis keys need not be a key itself, and can instead by an arbitrary vector (i.e. a vector). Thus the basis needs to have knowledge of the algebra class that it should return. This, obviously, leads to a few problems: * The basis class is required, as a template argument, to instantiate the vector (algebra) class, and so this require some careful declarations of classes to make everything work; * The basis class has to know about coefficients, which it otherwise wouldn't; * It makes it rather more difficult to separate the logic for densely stored vectors vs sparsely stored vectors. In the most recent version of libalgebra, this third issue was partially solved for tensors by changing the low-level implementation of multiplication to allow two separate operations to be used on dense and sparse parts of vectors. For non-tensors, the multiplication remained the same except for a small shim class inserted between which generated a stand-in function object and passed that to the vector mechanism. This was never intended to be a complete solution. Instead, I think multiplication should be provided by a additional class which is also supplied as a template argument to the algebra class and is stored as a static member of this class. This has several benefits: * the class implementing multiplication does not need knowledge of the underlying vector type, since it can be derived as a template argument; * it doesn't even need knowledge of the basis itself, only the key type and the coefficient type that could even be derived from template parameters; * it gives us complete control over the mechanism used to compute algebra products, we can simply implement different function objects and pass them into the underlying vector mechanics; * since vector objects are not tied directly to a specific depth, where such a parameter exists, we can more easily implement cross-depth multiplication. I've coded up tensor multiplication as a template of the kind of thing we need to write for each of the algebra derivatives in libalgebra. template <typename Coeff> class free_tensor_multiplication { typedef typename Coeff::SCA scalar_t; template <typename Transform> class index_operator { Transform m_transform; public: index_operator(Transform t) : m_transform(t) {} void operator()(scalar_t *result_ptr, scalar_t const *lhs_ptr, scalar_t const *rhs_ptr, DIMN const lhs_target, DIMN const rhs_target, bool assign = false) { scalar_t lhs; if (assign) { for (IDIMN i = 0; i < static_cast<IDIMN>(lhs_target); ++i) { lhs = lhs_ptr[i]; for (IDIMN j = 0; j < static_cast<IDIMN>(rhs_target); ++j) { *(result_ptr++) = m_transform(Coeff::mul(lhs, rhs_ptr[j])); } } } else { for (IDIMN i = 0; i < static_cast<IDIMN>(lhs_target); ++i) { lhs = lhs_ptr[i]; for (IDIMN j = 0; j < static_cast<IDIMN>(rhs_target); ++j) { *(result_ptr++) += m_transform(Coeff::mul(lhs, rhs_ptr[j])); } } } } }; template <typename Transform> class key_operator { Transform m_transform; public: key_operator(Transform t) : m_transform(t) {} template <typename Vector> void operator()(Vector &result, typename Vector::KEY const &lhs_key, scalar_t const &lhs_val, typename Vector::KEY const &rhs_key, scalar_t const &rhs_val) { result.add_scal_prod(lhs_key * rhs_key, m_transform(Coeff::mul(lhs_val, rhs_val))); } }; public: template <typename Algebra, typename Operator> Algebra &multiply_and_add(Algebra &result, Algebra const &lhs, Algebra const &rhs, Operator op) const { key_operator <Operator> kt(op); index_operator<Operator> it(op); lhs.buffered_apply_binary_transform(result, rhs, kt, it); return result; } template <typename Algebra, typename Operator> Algebra & multiply_and_add(Algebra &result, Algebra const &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const { key_operator <Operator> kt(op); index_operator<Operator> it(op); lhs.buffered_apply_binary_transform(result, rhs, kt, it, max_depth); return result; } template <typename Algebra, typename Operator> Algebra multiply(Algebra const &lhs, Algebra const &rhs, Operator op) const { Algebra result; multiply_and_add(result, lhs, rhs, op); return result; } template <typename Algebra, typename Operator> Algebra multiply(Algebra const &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const { Algebra result; multiply_and_add(result, lhs, rhs, op, max_depth); return result; } template <typename Algebra, typename Operator> Algebra &multiply_inplace(Algebra &lhs, Algebra const &rhs, Operator op) const { key_operator <Operator> kt(op); index_operator<Operator> it(op); lhs.unbuffered_apply_binary_transform(rhs, kt, it); return lhs; } template <typename Algebra, typename Operator> Algebra &multiply_inplace(Algebra &lhs, Algebra const &rhs, Operator op, DEG const max_depth) const { key_operator <Operator> kt(op); index_operator<Operator> it(op); lhs.unbuffered_apply_binary_transform(rhs, kt, it, max_depth); return lhs; } }; — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#19>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABQWN63HZVWR44M6N2HCNF3TWHL57ANCNFSM473AR3TQ>.

[inakleinbottle]

I've merged a first run at this into the develop branch. The tests are passing with these changes but there is some work yet to do.