"Ben Breen" <news@benbreen.com> wrote in message
news:akds9t$13c$1@paris.btinternet.com...
>
> > I'm wondering if anyone has given any thought to adding support for
> > small vectors and matrices to the C++ standard.
>
[snip]

Small vectors and small matrices are covered pretty good in Intels Small
Matrix Library. You might want to check, it out: it free, and very well
documented: http://www.intel.com/design/PentiumIII/sml/index.htm

Regards
.carsten


---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

Carsten Frigaard wrote:

> "Ben Breen" <news@benbreen.com> wrote in message
> news:akds9t$13c$1@paris.btinternet.com...
>
>>>I'm wondering if anyone has given any thought to adding support for
>>>small vectors and matrices to the C++ standard.
>>>
> [snip]
>
> Small vectors and small matrices are covered pretty good in Intels Small
> Matrix Library. You might want to check, it out: it free, and very well
> documented: http://www.intel.com/design/PentiumIII/sml/index.htm


    It is not, however, written in C++.  That's an encapsulation
of the Intel-only Streaming SIMD Extensions.  The underlying code is
in assembler, and will not run on non-Intel x86 CPUs.  Try
the demo on an AMD CPU, and you'll get an immediate illegal
instruction violation on a "movss" instruction.

    So that's not useful from a standardization standpoint.

    John Nagle
    Animats

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

Thanks for the code. Using this and a suitable class design I have managed
to come up with what I think is a good OO design of a matrix/vector class.

Please take a look at www.benbreen.com for the code, and a document
describing the class design considerations.

Do people think that this sort of design would be the sort of thing that
should be added to the standard. After all, most matrix designs go
completely against OO, so they would never be added. On the other hand, most
class designs are so slow that it would be poinless adding them to the
standard. I think I may have come to a comprimise.

Happy to hear your comments.
Ben Breen



---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

John Nagle wrote:
>
> Valmatrix,h  --  templated matrices, sized at compile time.
>
> ...
>
> valmatrix(const valmatrix& val)                       // copy constructor
> {     for (size_t i=0; i<isize; i++)
> {     for (size_t j=0; j<jsize; j++)
> {     (*this)(i,j) = val(i,j);}       // copy elt by elt
> }
> }
>
> valmatrix& operator=(const valmatrix& val)    // copy assignment
> {     for (size_t i=0; i<isize; i++)
> {     for (size_t j=0; j<jsize; j++)
> {     (*this)(i,j) = val(i,j);}       // copy elt by elt
> }
> return(*this);                                                // returns self
> }

Why are these written out explicitly? The automatically generated ones
will do the right thing, and probably be more efficient (for at least
one of the layout options) unless your compiler is really good at loop
optimisations.

> //    Raw data access
> const T* data() const { return(f_data); }// access as C array
> };

For interfacing with C APIs like OpenGL, you need a non-const version of
this too.

My matrix class is attached below (just the interface, to save space;
the implementation is all pretty standard stuff and probably not of
much interest). It only handles square matrices, but like you I gave it
a column/row-major layout option. I made column-major the default
because I expect most of the code I use it in to be OpenGL stuff.

There's an accompanying vector class that I haven't bothered to include
here; it's basically just a fixed-size array plus some vector
arithmetic operations. The uses of it here should be self-explanatory.

The arithmetic operations return their result by value, rather than
using the supposedly more efficient reference argument that I've seen
in many other matrix libraries, because I expect modern compilers to
fully implement return value optimisations.

(Now that I think about it, the arithmetic operators should be templates
so they can take operands with different value types and layouts, to
avoid using the conversion constructor to create a temporary.)

// Matrix layout options
enum MatrixLayout {
  column_major,
  row_major
};

// Matrix class template
template <typename T, unsigned N, MatrixLayout L = column_major>
class Matrix {
  public:
    // Member types
    typedef unsigned size_type;
    typedef T value_type;
    typedef Vector<T, N> vector_type;
    // Constants
    enum {
      cells = N * N,
      layout = L,
      n = N
    };
    // Default constructor (does no initialisation)
    Matrix();
    // Copy operations between different matrix types
    template <typename T2, MatrixLayout L2>
      Matrix(const Matrix<T2, N, L2>& source);
    template <typename T2, MatrixLayout L2>
      Matrix& operator=(const Matrix<T2, N, L2>& source);
    // Set all cells to the same value
    explicit Matrix(T t);
    // Set the leading diagonal to one value, the rest to another
    Matrix(T t1, T t2);
    // Extract a row or column
    vector_type column(unsigned i) const;
    vector_type row(unsigned i) const;
    // Pointer to the internal data, for use with C APIs
    const T* ptr() const;
    T* ptr();
    // Swap two rows or columns
    void swap_columns(unsigned i, unsigned j);
    void swap_rows(unsigned i, unsigned j);
    // Reference to a cell
    const T& operator()(size_type i, size_type j) const;
    T& operator()(size_type i, size_type j);
    // Matrix arithmetic operations
    Matrix operator-() const;
    Matrix& operator+=(const Matrix& rhs);
    Matrix& operator-=(const Matrix& rhs);
    Matrix& operator*=(const Matrix& rhs);
    friend Matrix operator+(const Matrix& lhs, const Matrix& rhs);
    friend Matrix operator-(const Matrix& lhs, const Matrix& rhs);
    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs);
    // Matrix/vector arithmetic operations
    friend vector_type
      operator*(const Matrix& lhs, const vector_type& rhs);
    friend vector_type
      operator*(const vector_type& lhs, const Matrix& rhs);
    // Matrix/scalar arithmetic operations
    Matrix& operator*=(T rhs);
    Matrix& operator/=(T rhs);
    friend Matrix operator*(const Matrix& lhs, T rhs);
    friend Matrix operator*(T lhs, const Matrix& rhs);
    friend Matrix operator/(const Matrix& lhs, T rhs);
    // Standard matrices
    static Matrix identity();
    static Matrix zero();
};

// Determinant of a matrix
template <typename T, unsigned N, MatrixLayout L>
  T det(const Matrix<T, N, L>& m);

// Inverse of a matrix
template <typename T, unsigned N, MatrixLayout L>
  Matrix<T, N, L> inverse(const Matrix<T, N, L>& m);

// Transpose of a matrix
template <typename T, unsigned N, MatrixLayout L>
  Matrix<T, N, L> transpose(const Matrix<T, N, L>& m);

// Solve a system of linear equations
template <typename T, unsigned N, MatrixLayout L>
  Vector<T, N> solve(const Matrix<T, N, L>& m, const Vector<T, N>& v);

// Comparison operators
template <typename T, unsigned N, MatrixLayout L1, MatrixLayout L2>
  bool operator==(const Matrix<T, N, L1>& lhs,
                  const Matrix<T, N, L2>& rhs);
template <typename T, unsigned N, MatrixLayout L1, MatrixLayout L2>
  bool operator!=(const Matrix<T, N, L1>& lhs,
                  const Matrix<T, N, L2>& rhs);

--
Ross Smith ..................................... Auckland, New Zealand
r-smith@ihug.co.nz ...................................................

        "A specter is haunting Wall Street; it is the specter
        of honest accounting."           -- Arthur D. Hlavaty

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

"Ben Breen" <news@benbreen.com> wrote in message news:<akds9t$13c$1@paris.btinternet.com>...
> > I'm wondering if anyone has given any thought to adding support for
> > small vectors and matrices to the C++ standard.
>
> I was about to start writing highly optimised matrix classes to build my
> dissertation on 3D real-time graphics engines upon.
>
> Now that I have seen so much interest here I will gladly release what I
> manage to come up with.
>
> If people can give me a few specific pointers as to what they want, I will
> try as much as possible to cater for their needs.

Be sure to read Barton & Nackman's "Scientific and Engineering C++".
Most of book is devoted to efficient implementation
of vectors and matrices in C++.

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

> I'm wondering if anyone has given any thought to adding support for
> small vectors and matrices to the C++ standard.

I was about to start writing highly optimised matrix classes to build my
dissertation on 3D real-time graphics engines upon.

Now that I have seen so much interest here I will gladly release what I
manage to come up with.

If people can give me a few specific pointers as to what they want, I will
try as much as possible to cater for their needs.

Would people rather have separate vector and matrix classes (I think this
will help optimisation somewhat)?

Also, have the dimensions of the matrices determined at compile time?

If people have done optimizations for particular platforms, what was
optimized and how was it implemented?




---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

Ben Breen wrote:

> I was about to start writing highly optimised matrix classes to build my
> dissertation on 3D real-time graphics engines upon.


    OK.  Here's something you can start with.  It's
a matrix class that's sized at compile time.  You
also get to decide whether you want a row-major or
column-major internal representation.  (Graphics
programmers known this as "are you using Direct-X,
or OpenGL?")

    I wrote this to see if VC++ 6 could optimize it properly.
It can.  All the encapsulation goes away in the code.  In particular,
with all optimizations on, matrix element access generates
the same code as does access to a C array.

    This needs to be fleshed out with operations on
the matrices. See

   http://www.animats.com/topics/developers.html

for what we've found useful building a physics
engine.  These are mostly
cleaned up versions of code from Graphics Gems,
written for 2D, 3D, and 4D vectors and matrices.
But the code isn't general as regards size.

     John Nagle
     Animats

   =====

//
//
Valmatrix,h  --  templated matrices, sized at compile time.
//
//
John Nagle
//
Animats
//
nagle@animats.com
//
August, 2002
//
//
License: LGPL
//
#ifndef VALMATRIX_H
#define VALMATRIX_H

//
//
class valmatrix  --  2D matrices of fixed size
//
//
Can be represented as either column-major or row-major
//
template <class T, size_t isize, size_t jsize, bool colmajor> class
valmatrix
{
private:
 T f_data[isize*jsize];   // the data
private:
 size_t pos(size_t i, size_t j) const// subscript calculation
 { return( colmajor ? (i * jsize + j) : (j * isize + i)); }
public:
 // Constructors
 valmatrix(const T& val)    // fill with single value
 { for (size_t i=0; i<isize; i++)
  { for (size_t j=0; j<jsize; j++)
   { this->(i,j) = val;}
  }
 }
 valmatrix(const valmatrix& val)   // copy constructor
 { for (size_t i=0; i<isize; i++)
  { for (size_t j=0; j<jsize; j++)
   { (*this)(i,j) = val(i,j);} // copy elt by elt
  }
 }
 valmatrix() {}       // default constructor does nothing
 // Assignment operators
 valmatrix& operator=(const valmatrix& val) // copy assignment
 { for (size_t i=0; i<isize; i++)
  { for (size_t j=0; j<jsize; j++)
   { (*this)(i,j) = val(i,j);} // copy elt by elt
  }
  return(*this);      // returns self
 }
 // Accessors.  () is the subscripting operator.
 T& operator()(size_t i, size_t j)
 { return(f_data[pos(i,j)]); }
 const T& operator()(size_t i, size_t j) const
 { return(f_data[pos(i,j)]); }
 // Info access
 size_t size(unsigned int i) const  // return array dimensions
 { switch (i) {
  case 0: return(isize);
  case 1: return(jsize);
  default: return(-1);
  }
 }
 // Raw data access
 const T* data() const { return(f_data); }// access as C array
};


#endif // VALMATRIX_H

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

===================================== MODERATOR'S COMMENT:
 Standard-related content flatlining, please take followups elsewhere.


===================================== END OF MODERATOR'S COMMENT
John Nagle <nagle@animats.com> wrote

> //
> Valmatrix,h  --  templated matrices, sized at compile time.

This same thing kept happening throughout your source code.
I suspect it happens wherever you used a TAB character --
one more reason to prefer using only the SPACE character.

[snip]
> template <class T, size_t isize, size_t jsize, bool colmajor> class
> valmatrix
> {
> private:
>  T f_data[isize*jsize];   // the data
> private:
>  size_t pos(size_t i, size_t j) const// subscript calculation
>  { return( colmajor ? (i * jsize + j) : (j * isize + i)); }

I suspect that this will optimize correctly on almost all compilers.
But I would recommend using partial specialization on the last
template parameter, whenever possible. That way you could write two
different versions of pos(), and not rely on the compiler doing the
optimization for you.

> public:
>  // Constructors
>  valmatrix(const T& val) // fill with single value
>  { for (size_t i=0; i<isize; i++)
>   { for (size_t j=0; j<jsize; j++)
>    { this->(i,j) = val;}
>   }
>  }

Probably better from an optimization standpoint to use:
    valmatrix(const T&val) {
       for (register size_t i=0; i<sizeof(f_data)/sizeof(f_data[0]); ++i)
          f_data[i]=val;
    }
You lose very little in clarity, and you rely less on the optimizer.

Could your compiler's optimizer reduce the nested loops into one
loop? Maybe. It might even reduce the two variables to one variable,
thereby using less stack space. But your version calls operator(),
which in turn calls pos(). Even if the compiler is smart enough to
make this 100% inline, it still needs to do some multiplication to
come up with an index into f_data. You're going to initialize the
entire array anyway; does it really matter what order you do this?
Avoid the multiplications by accessing the internal array directly.

If you know that T is a POD, you could do even better by using memcpy.

>  valmatrix(const valmatrix& val)   // copy constructor
>  { for (size_t i=0; i<isize; i++)
>   { for (size_t j=0; j<jsize; j++)
>    { (*this)(i,j) = val(i,j);} // copy elt by elt
>   }
>  }

Same concept here, except that this time the assignment calls operator()
twice for every element, thereby performing two multiplications.

Again, I suppose that some optimizing compilers could do almost as
good with the code above -- but why force the optimizer to figure
out something that's so obvious to humans?

>  valmatrix() {}       // default constructor does nothing
>  // Assignment operators
>  valmatrix& operator=(const valmatrix& val) // copy assignment
>  { for (size_t i=0; i<isize; i++)
>   { for (size_t j=0; j<jsize; j++)
>    { (*this)(i,j) = val(i,j);} // copy elt by elt
>   }
>  return(*this);      // returns self
>  }

Same comments again.

>  // Accessors.  () is the subscripting operator.
>  T& operator()(size_t i, size_t j)
>  { return(f_data[pos(i,j)]); }
>  const T& operator()(size_t i, size_t j) const
>  { return(f_data[pos(i,j)]); }
>  // Info access
>  size_t size(unsigned int i) const // return array dimensions
>  { switch (i) {
>   case 0: return(isize);
>   case 1: return(jsize);
>   default: return(-1);
>   }
>  }
>  // Raw data access
>  const T* data() const { return(f_data); }// access as C array
> };

Other than these minor criticisms, your class seems simple and
obvious -- usually the mark of a good class, especially when it
comes to reusable components.

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]

Allan W wrote:

> ===================================== MODERATOR'S COMMENT:
>  Standard-related content flatlining, please take followups elsewhere.
> ===================================== END OF MODERATOR'S COMMENT
> John Nagle <nagle@animats.com> wrote

>>Valmatrix,h  --  templated matrices, sized at compile time.


    Agreed.  To some extent, though, this is motivation
for the "small vectors" discussion and the "operator[]"
discussion.  It shows that the subject isn't entirely
abstract, but leads to useful, efficient code.

    John Nagle
    Animats

---
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@ncar.ucar.edu    ]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.jamesd.demon.co.uk/csc/faq.html                       ]