Re: surprising timings for multiplication of diagonalmatrix and full matrix

• To: mathgroup at smc.vnet.net
• Subject: [mg124966] Re: surprising timings for multiplication of diagonalmatrix and full matrix
• From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
• Date: Wed, 15 Feb 2012 04:39:59 -0500 (EST)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com

```Hello Eric,

On Mon, 13 Feb 2012, Eric Michielssen wrote:

> Hi,
>
> I'd like to right multiply an n x n diagonal matrix (specified by a vector)
> by an n x n full matrix.
>
> I tried 3 ways
>
> A. Brute force
>
> B. Using my own routine in which I was hoping to avoid the (supposedly) n^3
> cost of A.
>
> diagtimesmat[diag_, mat_] := MapThread[Times, {diag, mat}];
>
> C. The compiled version of B.
>
> cdiagtimesmat =
>  Compile[{{diag, _Real, 1}, {mat, _Real, 2}},
>   CompilationTarget -> "C"];
>
> I get weird timings.  For n=250:
>

250 is a magic number (100 is also one) where auto compilation for some
functions kicks in.

> n=250;nj=100;
> asd=Table[RandomReal[],{n},{n}];
> bsd=Table[RandomReal[],{n}];
> x=Timing[Do[ed=DiagonalMatrix[bsd].asd,{nj}]][[1]]
> y=Timing[Do[ed=diagtimesmat[bsd,asd],{nj}]][[1]]
> z=Timing[Do[ed=cdiagtimesmat[bsd,asd],{nj}]][[1]]
>
> Out[415]= 0.359
> Out[416]= 6.037
> Out[417]= 0.141
>
> My own uncompiled routine is superslow!!!  There are warnings about arrays
> being unpacked that I did not copy.
>

You may want to use AbsoluteTiming for these kind of tests. (See for
example http://mathematica.stackexchange.com/q/98/21)

Also, asd and bsd are not packed, thus repacking needs to be done in the
loops.

cdiagtimesmat =
Compile[{{diag, _Real, 1}, {mat, _Real, 2}},
CompilationTarget -> "C"];

n = 500; nj = 100;
asd = RandomReal[1, {n, n}];
bsd = RandomReal[1, {n}];
Developer`PackedArrayQ[asd]

x = AbsoluteTiming[Do[ed = DiagonalMatrix[bsd].asd, {nj}]][[1]]

1.847910

Off["Packing"]
y = AbsoluteTiming[Do[ed = diagtimesmat[bsd, asd], {nj}]][[1]]
On["Packing"]

0.577951

z = AbsoluteTiming[Do[ed = cdiagtimesmat[bsd, asd], {nj}]][[1]]

0.354027

> For n=500:
>
> In[418]:= n = 500; nj = 100;
> asd = Table[RandomReal[], {n}, {n}];
> bsd = Table[RandomReal[], {n}];
> x = Timing[Do[ed = DiagonalMatrix[bsd].asd, {nj}]][[1]]
> y = Timing[Do[ed = diagtimesmat[bsd, asd], {nj}]][[1]]
> z = Timing[Do[ed = cdiagtimesmat[bsd, asd], {nj}]][[1]]
>
> Out[421]= 2.777
> Out[422]= 1.31
> Out[423]= 1.014
>
> This is more reasonable. It remains a bit surprising that the routine that
> only touches n^2 numbers is only twice as fast as the one that supposedly
> touches n^3 ones.  Also, the compiled routine still does not achieve 100
> MFlops on my laptop.
>
> How can this behavior be explained?  What is the fastest way of doing this?
> And how about multiplying a full matrix by a diagonal one?

You could use Diagonal[] to extract the diagonal as a vector.

Hope this helps,

Oliver

>
> Thanks!
> Eric
>
> Eric Michielssen
> Department of Electrical Engineering and Computer Science
> University of Michigan
> 1301 Beal Avenue
> Ann Arbor, MI 48109
> Ph: (734) 647 1793 -- Fax: (734) 647 2106
>
>
>
>

```

• Prev by Date: Re: Running MathematicaScript in Mac Terminal
• Next by Date: Re: List Manipulation
• Previous by thread: Re: surprising timings for multiplication of diagonalmatrix and full matrix
• Next by thread: Problem with TableForm