> A colleague and I have started recently developing
> a C++ library for Vectors, Matrices and Tensors.
> Then another colleague informed me of your work,
> so I had a look at it. I must say that I am impressed
> by the code (I myself am quite new to C++)
> and by the completeness of your library.
>
> On the other hand I would like to make some remarks:
>
> - It is easy to access a row of a matrix (M[i]) but not
> a column, so I started off by defining M.column(j) =
> M.sub(0,M.length2(),1,j,1,1) which make things easier.
> Couldn't you include that in your lib?
What was wrong with M.t()[j]?
> - Secondly I miss a function to invert a matrix,
> which is quite a basic function.
Make sure that you have the very latest version of SVMT
http://www.netwood.net/~edwin/svmt/
It includes SubSquare and Square matrix classes.
Member functions i(<type> r = 0) and invert(<type> r = 0)
are declared but not defined in the floating-point
SubSquare class definitions. The problem is that
there are several ways to invert a square matrix
so the definition of these member functions
is left up to the application programmer.
> - The most important problem I have with your library however
> is that you define M*N as an _element-wise_ multiplication
> of to matrices !!??? Do I really need to write
> M.dot(N.t()) for a normal matrix multiplication?
You could define
inline
const
<Type><System>Matrix operator %(
const <Type><System>SubMatrix& M,
const <Type><System>SubMatrix& N) {
return M.dot(N.t()); }
for floating-point types but the use of matrix-matrix multiplication
is discouraged because striding across rows of a large matrix N
results in a cache miss on every access to an element of N
which can have a very negative impact on performance.
More than likely, you will find that it makes more sense
to store the transpose of N instead of N in memory
and use the dot product. Store both N and its transpose
if you really need N and have the required memory.
> We have chosen to use M*N for that
> and to overload the & operator (arbitrary choice, I admit)
> as the element-wise multiplication operator.
> (The logic being that & normally is a bit-wise AND
> - also element-wise).
> The | would then be the element-wise division.
I think that you would do well to reconsider your choices.
Matrix-matrix multiplication is a very important operation
but element-by-element multiplication (and division)
will appear far more often in your application source code.
The meaning of operator * (...) and operator /(...) should be
consistent with the meaning of operator +(..) and operator -(...)
so that they are easy for C and C++ programmers to learn
and remember. Matrix-matrix multiplication has
a very special meaning and should have a special operator.
Matrix-matrix multiplication operator %(...) would have
the correct precedence and obviate the extra parentheses
that your choices require in expressions combining
element-by-element addition, subtraction, multiplication
and division with matrix-matrix multiplication.
Using operator &(...) and operator |(...) respectively
for element-by-element multiplication and division
will only serve to confuse application programmers
if (when?) you decide to implement
vector, matrix and tensor classes for integral types.
[In C and C++, operator *(...) and operator |(...)
are normally commutative operations
but operator %(...) and operator /(...)
like matrix-matrix multiplication
and element-by-element division
are not normally commutative.]
> - Finally I get from expressions like:
> v.dot(w) = vwT and M.dot(v) = MvT
> that you define a vector as being a row-vector,
> while the whole scientific community (as far I know at least)
> defines vectors as being column-vectors.
You may have heard people say that,
"Fortran stores two dimensional arrays in column major order
but C stores two dimensional arrays in row major order."
In fact, Fortran and C compilers almost always store
multidimensional arrays in exactly the same way
but the order of the subscripted indices is reversed.
Unfortunately, C and C++ programmers don't interpret M[i][j]
as a reference to the element in row j and column i of matrix M.
They think that the leftmost subscript is a row index
and the rightmost subscript is a column index.
Two dimensional arrays appear to be transposed
when a C program calls a Fortran subroutine
or when a Fortran program calls a C function.
SVMT class libraries are consistent with C and C++.
The subvector M[i] references row i of submatrix M
so it is natural to interpret SVMT vectors as row vectors.
Two dimensional Fortran arrays conform
to the conventional interpretation of matrices
as a collection of column vectors so it is natural to interpret
one dimensional Fortran arrays as column vectors.
I use the conventional notation to derive expressions
which involve vector and matrix objects but then
I transpose the entire expression and replace
each individual vector and matrix object with its transpose
before I code it using SVMT class libraries.
This interpretation of vector and matrix objects has subtle
but practical advantages over the conventional interpretation
for C and C++ programmers.
The elements of vector and matrix objects are stored in memory
in the order in which they appear in the I/O stream.
The elements of a (short) row vector
can be displayed on a single line of output.
The application programmer simply writes
v%M%L
instead of
L%(M%v)
to evaluate the combined matrix-matrix-vector product
efficiently since the operators are left associative.
> PS: We are aware that software published under the GNU-GPL
> cannot be used to write commercial software.
> We only use it internally for development purposes.
You certainly _may_ use the SVMT portable reference library
to write commercial software. You _may_ copyright and distribute
application programs which use the SVMT portable reference library
without disclosing the application source code to anyone.
If you write your own SVMT class library from scratch,
you are permitted and encouraged to copyright it and distribute it for profit
without disclosing any of the source code other than the header files.
You may _not_ modify any part of the existing SVMT portable reference library
and distribute it unless you comply with the terms of the GNU-GPL.
With your permission, I would like to publish this reply
to the Object Oriented Numerics Standard mailing list.
http://monet.uwaterloo.ca/oon/oonstd/
Thanks in advance, E. Robert Tisdale <edwin@netwood.net>
Edwin,
First of all thanks for answering so quickly.
> > - it is easy to access a row of a matrix (M[i]) but not
> > a column, so I started off by defining M.column(j) =
> > M.sub(0,M.length2(),1,j,1,1) which make things easier.
> > Couldn't you include that in your lib?
>
> What was wrong with M.t()[j]?
It returns a row, not a column. Well, see later on...
> <SNAP>
> > - The most important problem I have with your library however
> > is that you define M*N as an _element-wise_ multiplication
> > of to matrices !!??? Do I really need to write
> > M.dot(N.t()) for a normal matrix multiplication?
>
> You could define
>
> inline
> const
> <Type><System>Matrix operator %(
> const <Type><System>SubMatrix& M,
> const <Type><System>SubMatrix& N) {
> return M.dot(N.t()); }
>
> for floating-point types but the use of matrix-matrix multiplication
> is discouraged because striding across rows of a large matrix N
> results in a cache miss on every access to a element of N
> which can have a very negative impact on performance.
> More than likely, you will find that it makes more sense
> to store the transpose of N instead of N in memory
> and use the dot product. Store both N and its transpose
> if you really need N and have the required memory.
That is a way out, but I would much prefer
being able to write M*N than M.dot(N_t). Personally
I prefer readable and aesthetic code over fast code.
But who am I ?
> > We have chosen to use M*N for that
> > and to overload the & operator (arbitrary choice, I admit)
> > as the element-wise multiplication operator.
> > (The logic being that & normally is a bit-wise AND
> > - also element-wise).
> > The | would then be the element-wise division.
>
> I think that you would do well to reconsider your choices.
> Matrix-matrix multiplication is a very important operation
> but element-by-element multiplication (and division)
> will appear far more often in your application source code.
Not true.
It would look like you guessed I am an image processor.
For those applications element-wise op's are indeed most frequent.
But we have special software for that.
My C++ code is about linear algebra and 80-90% of the operations
are Matrix-Matrix or Matrix-Vector multiplications
or additions/sustractions.
> The meaning of operator * (...) and operator /(...) should be
> consistent with the meaning of operator +(..) and operator -(...)
> so that they are easy for C and C++ programmers to learn
> and remember. Matrix-matrix multiplication has
> a very special meaning and should have a special operator.
I don't agree.
You think very much from the point of view of a (C) programmer,
while I take the point of view of a researcher/mathematician.
In linear algebra the only meaning that exists for 'multiplication'
in association with a matrix, is matrix-matrix multiplication.
This is not a 'special meaning' but the common meaning.
Element-wise multiplication is only a programmer's trick
(except for scalar-matrix multiplication).
I think the library shouldn't be easy
for programmers to use and remember, but for researchers.
We use it to develop test code, not so much for production code.
(That also explains why I don't care too much for maximum speed code).
The math of our algorithms is so complex
that if you add to that by changing operators, we will get lost.
> Matrix-matrix multiplication operator %(...) would have
> the correct precedence and obviate the extra parentheses
> that your choices require in expressions combining
> element-by-element addition, subtraction, multiplication
> and division with matrix-matrix multiplication.
> Using operator &(...) and operator |(...) respectively
> for element-by-element multiplication and division
> will only serve to confuse application programmers
> if (when?) you decide to implement
> vector, matrix and tensor classes for integral types.
We use template based classes,
so they are defined for integral data-types.
But we don't see an int as a set of bits, but always as a number,
so the original definition of & and | are meaningless.
> > - Finally I get from expressions like:
> > v.dot(w) = vwT and M.dot(v) = MvT
> > that you define a vector as being a row-vector,
> > while the whole scientific community (as far I know at least)
> > defines vectors as being column-vectors.
>
> You may have heard people say that,
> "Fortran stores two dimensional arrays in column major order
> but C stores two dimensional arrays in row major order."
> In fact, Fortran and C compilers almost always store
> multidimensional arrays in exactly the same way
> but the order of the subscripted indices is reversed.
> Unfortunately, C and C++ programmers don't interpret M[i][j]
> as a reference to the element in row j and column i of matrix M.
> They think that the leftmost subscript is a row index
> and the rightmost subscript is a column index.
> Two dimensional arrays appear to be transposed
> when a C program calls a Fortran subroutine
> or when a Fortran program calls a C function.
> SVMT class libraries are consistent with C and C++.
> The subvector M[i] references row i of submatrix M
> so it is natural to interpret SVMT vectors as row vectors.
>
> Two dimensional Fortran arrays conform
> to the conventional interpretation of matrices
> as a collection of column vectors so it is natural to interpret
> one dimensional Fortran arrays as column vectors.
>
> I use the conventional notation to derive expressions
> which involve vector and matrix objects but then
> I transpose the entire expression and replace
> each individual vector and matrix object with its transpose
> before I code it using SVMT class libraries.
>
> This interpretation of vector and matrix objects has subtle
> but practical advantages over the conventional interpretation
> for C and C++ programmers.
> The elements of vector and matrix objects are stored in memory
> in the order in which they appear in the I/O stream.
> The elements of a (short) row vector
> can be displayed on a single line of output.
> The application programmer simply writes
>
> v%M%L
>
> instead of
>
> L%(M%v)
>
> to evaluate the combined matrix-matrix-vector product
> efficiently since the operators are left associative.
That is exactly what I DON'T want to do.
Again, I'm writing this as a scientist, not a a programmer.
I don't even know Fortran.
My interpretation of a vector as a column does not stem
from Fortran or C/C++ conventions but from linear algebra.
There Mij is the j'th (column) element of the i'th row.
Btw. had you been a mathematician,
you would have known that (L*M)*v = L*(M*v)
so the parentheses are not necessary.
Again, as a mathematician I want my code
to resemble as much as possible what I write on paper
or (as someone else mentioned in the OONSTD newsgroup) Matlab code.
We usually start with a piece of paper, try it out in Matlab
and then make a more elaborate test and/or demo program in C++.
It would be nice
if the conventions would be the same across these 'platforms'.
> > PS: We are aware that software published under the GNU-GPL
> > cannot be used to write commercial software.
> > We only use it internally for development purposes.
>
> You certainly _may_ use the SVMT portable reference library
> to write commercial software. You _may_ copyright and distribute
> application programs which use the SVMT portable reference library
> without disclosing the application source code to anyone.
> If you write your own SVMT class library from scratch,
> you are permitted and encouraged to copyright it and distribute it for profit
> without disclosing any of the source code other than the header files.
> You may _not_ modify any part of the existing SVMT portable reference library
> and distribute it unless you comply with the terms of the GNU-GPL.
Thanks, that is generous of you.
I hope we will be able to use it
and maybe in return contribute some to it.
> With your permission, I would like to publish this reply
> to the Object Oriented Numerics Standard mailing list.
>
> http://monet.uwaterloo.ca/oon/oonstd/
>
> Thanks in advance, E. Robert Tisdale <edwin@netwood.net>
No problem. You can also add this reply to it.
I think I will subscribe,
so we can continue this discussion there with the others...
Kind regards,
Matthijs
-- --------------------------------------------------------------------- | ir. Matthijs W. van der Molen | Volmerlaan 8, L5-300 | | Shell Int'l Exploration & Production | P.O.Box 60, 2280 AB Rijswijk | | Research and Technical Services | the Netherlands | | New Technologies Department | tel. ++31 - 70 - 311 6114 | | m.w.vandermolen@siep.shell.com | fax. ++31 - 70 - 311 3366 | ---------------------------------------------------------------------