Re: [boost] 2D Matrix

From: Oliver Koenig (koenig_at_[hidden])
Date: 2006-11-08 06:38:40

Hi,
we recently investigated available C++ matrix libraries in the context
of a new implementation of our Finite Element library eXperiment FeLyX.
Within the element formulations (e.g. evaluation of stiffness matrices)
of such FE libraries there are lots of fixed-size matrix operations
where the sizes are known at compile time. Usually the matrix sizes are
given through the actual element type used.

We looked at MTL, UBLAS and TVMET (Tiny Vector Matrix library using
Expression Templates, http://tvmet.sourceforge.net/), and finally
decided to use TVMET.
The reasons were:
- As far as I understand your 2D-matrix template, TVMET offers exactly
this type of fixed-size matrix library you are describing, providing
compile-time checks for matrix sizes etc.
- TVMET is fast. For matrices from size 3x3 up to app. 60x60 we did a
very basic performance comparison (Matrix-Matrix / Matrix-Vector
products) with MTL, UBLAS and TVMET (Tiny Vector Matrix library using
Expression Templates, http://tvmet.sourceforge.net/). We were really
impressed by the performance of TVMET, in our basic test it was
significantly faster than the other two.
- The operator overlads of TVMET  (implemented using expression
templates) provide a very convenient and easy-to-use interface
(especially compared to MTL which we used before).

I think TVMET could represent an interesting base for such a library
within BOOST?

Sincerely
ok

-----------
Oliver Koenig
EVEN - Evolutionary Engineering AG
www.even-ag.ch

Jason Kraftcheck wrote:
> Hi,
> 
> I have written a simple 2D matrix template where the extents of the matrix are
> template parameters.  The multi_array container and the uBLAS matrix class both
> provide matrices of fixed dimension (template parameter for multi_array and 2
> for uBLAS) but variable extents.  I think my matrix definition is complementary
> to these. There are several advantages to having an array of fixed extents; but
> the implementation is, of course, unusable when the extents are not known at
> compile-time.
> 
> If the extents of the matrix are part of the type, the compiler can enforce
> correct matrix extents for operations, avoiding the cost of runtime checks.  For
> example:
>   double b[3][2] = { ... };
>   matrix<3,3> A(2.0); // 2 * identity
>   matrix<3,2> B(b);
>   matrix<2,3> C;
>   matrix<3,2> D;
>   D = A * B; // OK
>   C = A * B; // compile error - mismatched types
>   C = transpose(B) * A; // OK
> Optimizing compilers can trivially unroll the loops used to define the
> operations.  Use of this class does not require heap allocation.  The storage
> required for the matrix is the minimal possible (assuming no assumptions about
> the values of the elements.)  An array of matrices can be "serialized" to an
> array of doubles with a single reinterpret_cast.
> 
> The proposed implementation is unsuitable for very large matrices: defining the
> element storage as a contiguous array may be problematic, most of the
> performance gains become insignificant, and the implicit copying will become
> significant.
> 
> I've attached my current implementation.  In it's current form, it is unsuitable
> for inclusion in Boost.  I will fix any required issues if there is any interest
> in including it in Boost.  The things that need to be changed (that I'm aware
> of) are:
>  o license - GPL is not acceptable for Boost
>  o naming (namespace and class name)
>  o change explict 'double' for element type to a template parameter
> Some other changes that should probably be done before inclusion:
>  o implement inverse for other than 2x2 and 3x3 matrices
>  o more efficient determinant implementation for larger matrices
>  o LU factorization
>  o QR factorization
> 
> In timing tests using G++ 4.1 on a Linux/Pentium4 platform, the performance of
> the arithmetic operations was equivalent to those of a special 3x3 matrix
> implementation with all operations explicitly unrolled (the 'old' implementation
> this one was to replace.)
> 
> Is there any interest in including such a matrix definition in Boost?
> 
> Is there already an effort to include something equivalent?  I searched back a
> few years in the mailing list archive and didn't notice anything.
> 
> Any suggested changes (other than those I listed above) for the attached
> implementation?
>