Open inakleinbottle opened 3 years ago
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:
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:
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
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 GitHubhttps://github.com/terrylyons/libalgebra/issues/19, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABQWN63HZVWR44M6N2HCNF3TWHL57ANCNFSM473AR3TQ.
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.
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 forfree_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: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:
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.