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.