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 blockwise fashion.
> 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 let the user swap the meanings of & and *.
> 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 additionnal megabytes at compilation
time. I tried, for fun. And succeeded.
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.
Daniel Goujot, from university Paris-6. http://www.ann.jussieu.fr