MET - Matrix Expression Templates

Goals of MET  - Promoting the code Reusability

MET is a C++ matrix class library which promotes the notational convenience of linear algebraic codes but is free of the overhead of superfluous temporary matrix objects. Instead of coding tedious loops for the linear algebraic expression in C or fortran,  C++ overloaded operators enable one to write the simplest code like u = m*v + w. This advantage, however, came with a performance penalty until recently. Due to the overhead of temporaries and copying of matrix objects,  C++ lagged behind fortran's performance by an order of magnitude.

This problem was resolved by expression templates, whose effectiveness was proved in various array libraries, Blitz++,POOMA, PETE. But these libraries are just for array and vector and do not provide an ordinary matrix class. This is why I decided to write such a class using expression templates.

Now there already exists quite a few matrix classes over the internet. Moreover it is very likely that you would make another matrix class with the different data layout or different feature. Then it will be more desirable to attach (or "glom") an expression template facility onto your matrix class than to replace it.

Theoretically, "disambiguated glommable" expression templates can be added to yours.( Computer in Physics,  vol. 11, no. 3, 263p ) Since this newer technique is implemented in MET,  you can easily register your own matrix class with the expression template. All you need to do is to derive your class from that of MET and provide a few specific methods :  int rows() , int cols()  and operator()(int i, int j) .

Thus, in addition to providing the ready-to-use matrix and vector classes, MET serves as a code basis for extending your classes to be more notationally convenient but suppress the superfluous overhead.
 

Features

The Matrix Expression Templates library is,

Latest Release

can be downloaded here.

Installation

Makefile is prepared for GNU make and C++ compiler (gcc 2.95.2). Expand gziped tar file  and just type,

cd met/
make

Getting Started

Let's start with a simple example, exam01.C. Consider u = m*v + w, where u,v and w are vectors and m is matrix.

const int size = 3;
Vec<double> u(size), v(size), w(size);  // the element type is double precision.
Mat<double> m(size);                    // 3x3 square matrix

v(0) = 1;  v(1) = 3;  v(2) = 5;         // access to the vector element
w = v;

m(0,0) = 2;  m(0,1) = 3;  m(0,2) = 4;   // access to the matrix element
// set the rest of elements as you like.

u = m * v + w;

cout << "Answer is (";
for (int i=0; i<size; i++) cout << " " << u(i);
cout << " ) ." << endl;
 

How to attach an expression template facility onto your class?

The expression template is implemented with container classes for holding the data of vector or matrix elements, along with a set of global template functions which provides the closure object construction mechanism. The evaluation of any expression is automatically transfered into an a composite closure object which resolves to simple loops because of inlining. This composite closure is called as an expression template, which is parsed in the compilation. ( Detailed explanation of the technique and the outstanding benchmark results are shown in blitz++ web site. )

Since terminals of the parse tree are container objects, that implementation is easily generalized by using abstract base classes to represent the container of element data. User-defined matrix and vector classes would be derived from these base classes to be eligible to be an terminal in an expression-template parse tree. You can override the indexing operator for accessing an element to conform to your own data layout. The base classes, however, must be argumented with class template, otherwise the speed penalty of virtual functions is not acceptable for accessing vector and matrix elements.

Hence the base class of vector classes in MET is declared in xpr1.h as ;

template <class ElementType, class UserContainer>
class Dim1 {
public:
  // Element access
  ElementType operator() (int n) const;

  // Assignment from an expression, vector, scalar
  template <class E>
    UserContainer& assignFrom(const Xpr1<ElementType,E>& x);
  template <class V>
    UserContainer& assignFrom(const Dim1<ElementType,V>& x);
  UserContainer& assignFrom(ElementType x);

  // operator += with an expression, vector, scalar
  UserContainer&
    Dim1<ElementType,UserContainer>&
      operator+=(const Xpr1<ElementType,E>& x);
  // ....

  // other operator declarations
};

so the user-defined vector class can be derived as ;

class UserVec : public Dim1<double,UserVec> {
  int sz;        // size of vector
  double *data;  // or your own data structure

public:
  explicit UserVec(int sz_) : Dim1<double,UserVec>(),
                              sz(x.sz), data(new double[sz]) {}
  UserVec(const UserVec& x) : Dim1<double,UserVec>(),
                              sz(x.sz), data(new double[sz]) {
    assignFrom(x);
  }
  ~UserVec() { delete [] data; }

  // Assignment from an expression, vector, scalar
  template <class E>
    UserVec& operator=(const Xpr1<ElementType,E>& x) {
      return assignFrom(x);
    }
  template <class V>
    UserVec& operator=(const Dim1<ElementType,V>& x) {
      return assignFrom(x);
    }
  UserVec& assignFrom(double x) {
    return assignFrom(x);
  }

  int size() const { return sz; }
  double operator()(int i) const { return data[i]; }
  double& operator()(int i) { return data[i]; }

  // your own methods ...
};
 

Since the assignment operators cannot be inherited,  you have to define them using  asignFrom(x) function which is implemented in the base class.  asignFrom(x) can also be used to define the copy constructor. For registering your class with the expression template facility,  you have to provide  int size() to indicate the indexing range and operator()(int i) to actually access an element of your vector.

The user-defined matrix class can be derived from Dim2<ElementType, UserMat> class  (defined in xpr2.h) in a similar way.

class UserMat : public Dim2<double,UserMat> {
  int r,c;        // the num. of rows and columns
  double *data, **m;  // or your own data structure

public:
  explicit UserMat(int r_, int c_) : Dim2<double,UserMat>(),
    r(r_), c(c_), data(new double[r*c]), m(new double*[r] {
    for (int i=0; i<c; i++) m[i] = data + i*r;
  }
  UserMat(const UserMat& x) : im2<double,UserMat>(),
    r(r_), c(c_), data(new double[r*c]), m(new double*[r] {
    for (int i=0; i<c; i++) m[i] = data + i*r;
    assignFrom(x);
  }
  ~UserMat() {
    delete [] m;
    delete [] data;
  }

  // Assignment from an expression, matrix, scalar
  template <class E>
    UserMat& operator=(const Xpr2<ElementType,E>& x) {
      return assignFrom(x);
    }
  template <class M>
    UserMat& operator=(const Dim2<ElementType,M>& x) {
      return assignFrom(x);
    }
  UserMat& assignFrom(double x) {
    return assignFrom(x);
  }

  int rows() const { return r; }
  int cols() const { return c; }
  double operator()(int i, int j) const { return m[j][i]; }
  double& operator()(int i, int j) { return m[j][i]; }

  // your own methods ...
};
 
 

Copyright & Disclaimer

To Do

  • Wrap LAPACK library for use in C++.
  • Write documents. ( If you find something wrong, please correct me via e-mail. )
  • Support for cross-compiling through GNU Autotools
  • Tune codes to rival fortran.   :-)   (Currently MET is not so fast as f77.)


    • Link

  • The Object-Oriented Numerics Page
  • The Open Science Project

    •         > Mathematics > Linear Algebra
  • SAL (Scientific Applications on Linux)

    •         > Numerical Analysis > Miscellaneous software
  • freshmeat
  • Blitz++ :

    •         Array and Vector classes which rival Fortran's performance.
     

    Project Summary

    If you would like to help to develop MET, please e-mail to me or visit the project summary page. This page has :

    © 2001, Copyright by Masakatsu Ito; All Rights Reserved. mailto:masakatsu@users.sourceforge.net





    This site is hosted by SourceForge Logo