Diagonal matrix times dense matrix

[bencteux]

Hello,
I have to compute many times the product of a diagonal matrix D with a dense matrix A, left or right.
I am looking for a very efficient implementation. I don't find any in LAPACK/BLAS, am I right ?

To compute : A * D, I plan to apply dscal on each column vector of A.
But I hesitate about the computation of : D * A.
* computing each row in a loop : for each row, applying dscal with incx equals to the size of a column ?
* or : using dlagtm, which does : T * A, for T tridiagonal ?

Which would be the most efficient ? Another solution ?
Thank you for help,
Guy Bencteux

[Julien Langou]

right multiplication: DLASCL is fine, try an optimized BLAS, compare with a silly loop, compare with reference BLAS
(unrolling of 4), try various unrolling parameters, read a book on low-level optimization, etc. depends on how much time
you have for it! But an optimized DLASCL should be fine. The performance you want to see is basically your memory bandwidth.

left multiplication: DIY: Do It Yourself.
Do your own loop (the outer loop being the column loop):

Code: Select all

            DO 60 J = 1, N
                  DO 50 I = 1, N
                     B( I, J ) = D(I) * B( I, J )
   50            CONTINUE
   60       CONTINUE


and hope for the compiler to do optimize for you.
Some optimization you can do yourself:
* You can try the (I-J) loop (as opposed to J-I) but that should be way slower (since A is stored by columns)
* You can try to unroll by hand (look at dscal.f to see how to unroll)

I am not sure what you mean by:

[I have to compute many times the product of a diagonal matrix D with a dense matrix A, left or right.

If you need to do D1*D2*A*D3*D4 you really want to do: (D1*D2)*A*(D3*D4).

-julien.

[bencteux]

Thank you all for the answers !
Some precisions :
* what I need to do is : D * A, and later in the computation : A * D, with other D and A.
* In a complete computation, I have to do such products up to 5,000 times
* The size of the matrices : up to 2,000

This is not the critical part of the computation, but I want to be sure to take maximal benefit from D being diagonal.
I will do as Julien have suggested.
Thanks once more !