Combining a linear algebra library with Boost::Units

Question

I'm doing a good amount of scientific programming and made very good experiences with both Boost.Units, which provides compile-time dimensional analysis for quantities (i.e. tags quantities with units and thereby catches many errors with classical physical dimension analysis) and using Eigen 2 for linear algebra.

However, Eigen has no concept of units, and while you can set the scalar quantities in matrices for Eigen, it expects that the multiplication of two quantities yields the same type, which is obviously untrue for units. For example, code like:

using boost::units::quantity;
namespace si = boost::units::si;
Eigen::Matrix< quantity< si::length >, 2, 1 > meter_vector;
quantity< si::area > norm = meter_vector.squaredNorm();

does not work, even though it is logically correct.

Is there any matrix library that supports units? I know that this would have been notoriously hard to implement in the past, and C++11 and decltype will make that much easier, but it surely was possible with C++03 and template specializations.

Answer 1

I believe Blitz++ supports much of Boost.Units functionality.

Edit by the OP: For the reference here is the full test code with which I tested the Blitz matrix multiplication functionality:

#include <blitz/array.h>
#include <boost/units/systems/si/area.hpp>
#include <boost/units/systems/si/length.hpp>
#include <boost/units/quantity.hpp>

using boost::units::quantity;
namespace si = boost::units::si;

namespace blitz {
template< typename U1, typename T1, typename U2, typename T2>
struct Multiply< quantity<U1,T1>, quantity<U2,T2> >
{
    typedef typename boost::units::multiply_typeof_helper< quantity<U1,T1>, quantity<U2,T2> >::type T_numtype;

    static inline T_numtype apply( quantity<U1,T1> a, quantity<U2,T2> b ) { return a*b; }
};

}

using namespace blitz;

int main() {
    Array< quantity<si::length>, 1 > matrix;
    Array< quantity<si::area>, 1 > area;
    area = matrix * matrix;
    return 0;
}

Answer 2

You should check this Wiki page: http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html

Eigen requires some work to use other than primitive data types, but it's generally possible.

Answer 3

The difficulty of using the standard Eigen library plugin option, is that the existing operators +, -, *, etc, needs to be replaced for Boost Units quantities to be used.

For example, for a Boost units custom type to work with the * multiply operator, for an arbitrary CUSTOM_TYPE, it needs to look like this:

template<class X,class Y>
CUSTOM_TYPE<typename boost::units::multiply_typeof_helper<X,Y>::type>
operator*(const CUSTOM_TYPE<X>& x,const CUSTOM_TYPE<Y>& y)
{
    typedef typename boost::units::multiply_typeof_helper<X,Y>::type    type;

    return CUSTOM_TYPE<type>( ... );
}

Notice how the return type is not the same as the input types. Here you use the template helper multiply_typeof_helper to create the return type. This is because multiplying meters with seconds will not give you a quantity of either unit. However, the default Eigen * operator will return the same "type" as that of the inputs - this is the problem.

The other option is to embed the Eigen matrix inside the quantity, rather than embedding the quantity inside the matrix.