> On Tue, 23 Mar 1999, E. Robert Tisdale wrote:
>
> > 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.
>
> A better well know strategy is to split M and N in block matrices
> with each block having the memory size
> of just less than one third of the cache,
> and to compute the multiplication in a block wise fashion.
No. This does not solve the striding problem. Striding is a problem
because the memory system fetches two or more elements
to fill the typical cache line but only one of them is used.
The block decomposition strategy is effective
when the footprint of a row of matrix M together with matrix N
exceeds the size of the data cache.
I believe that the block decomposition strategy should be
even more effective with the dot product of two large matrices.
> > M.W. van der Molen wrote:
> > > We have chosen to use M*N for that and to overload the & operator
> > > (arbitrary choice, I admit) as the element-wise multiplication
> > > operator.
>
> > 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)
>
> M.W. van der Molen wrote:
> > Not true. [...]
>
> Please don't discuss about that sort of choice here.
> You cannot decide a consensus in a thread.
> Do it as you want, a simple and short public inheritance class
> lets the user swap the meanings of & and *.
Please refer to an earlier thread on "Scenario of use".
The purpose of this mailing list is to discuss choices just like this.
That's what standardization is all about. Amateur programmers
have different priorities than professional programmers
and perhaps there should be a separate standard for them.
Is anyone passionate enough about this to take on the task?
> > v%M%L
> >
> > instead of
> >
> > L%(M%v)
>
> It is possible to make L*M resulting in a pair class : Product_pair(L,M)
> and that L*M*v results in : Product_pair(Product_pair(L,M),v)
> and to type special members in class Vector, like :
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
> template<class M>
> Vector & Vector::operator+=(
> Product_pair<Product_pair<M,Matrix>,Vector> a_pair) {
> operator+=(Product_pair(
> a_pair.first.first,
> compute_matrix_vector_product(
> a_pair.first.second,
> a_pair.second)
> ));
> return * this;
> }
> /* and other members... */
> _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
>
> to make a correct handling of u += L*M*v.
> Beware, this adds lots of hours understanding TMP stuff
> and lots of additional megabytes at compilation time.
> I tried, for fun. And succeeded.
Bjarne Stroustrup, "The C++ Programming Language: Third Edition",
Chapter 22, "Numerics", Section 4, "Vector Arithmetic",
Subsection 7, "Temporaries, Copying, and Loops", pages 675-6.
But you are evaluating an expression containing left associative operators
right to left which will usually yield a slightly different result
and you will need to explain this to the programmer.
> Just my first 2 pence in this list.
>
> I have done an implementation of matrices and vectors
> relative to algebra basis, to make easier
> the implementation of wavelet transform
> in Partial Differential Equation resolution.
> Not yet stable. Not yet commented. Under GPL.
Please tell us how to get it?
Thanks in advance, Bob Tisdale <edwin@netwood.net>