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}}, > Module[{}, MapThread[Times, {diag, mat}]], > 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}}, Block[{}, MapThread[Times, {diag, mat}]], 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 > Radiation Laboratory > 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 > > > >