Message Boards Message Boards

3
|
8617 Views
|
9 Replies
|
9 Total Likes
View groups...
Share
Share this post:
GROUPS:

Is Mathematica 40 times slower than Fortran with matrixes?

Posted 11 years ago
Consider this code:
In[1]:= ht = 0.01;

In[2]:= a = Table[Sin[i] + I Cos[j], {i, 0., 0.5, 0.01}, {j, 0., 0.5, 0.01}];

In[3]:= HL = Compile[{{t, _Real}}, I ht/2*a, CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}];

In[4]:= Table[HL[200.], {23497}]; // AbsoluteTiming
Out[4]= {2.070978, Null}
The same function writtern in Fortran only takes 0.05 seconds to run, so that is a 40 times speed difference.

We can check that there is no calls to MainEvaluator: 
In[5]:= Needs["CompiledFunctionTools`"]

In[6]:= StringFreeQ[CompilePrint[HL], "MainEvaluate"]
Out[6]= True

For some reason, a is not a packed array, but even if we use the packed verion, it's still slow:
 In[7]:= Developer`PackedArrayQ[a]
 Out[7]= False
 
 In[8]:= apk = Developer`ToPackedArray[a];
 
 In[9]:= HL =
   Compile[{{t, _Real}}, I ht/2*apk,
    CompilationOptions -> {"InlineCompiledFunctions" -> True,
      "InlineExternalDefinitions" -> True}];

In[10]:= StringFreeQ[CompilePrint[HL], "MainEvaluate"]
Out[10]= True

In[11]:= Table[HL[200.], {23497}]; // AbsoluteTiming
Out[11]= {1.974200, Null}

And it's not the problem of Compile:
In[12]:= Table[I ht/2*a, {23497}]; // AbsoluteTiming
Out[12]= {1.935287, Null}
In[13]:= Table[I ht/2*apk, {23497}]; // AbsoluteTiming
Out[13]= {1.933595, Null}
So why there is a such a speed difference, and how to improve that?



Update:

This is the Fortran code I'm using 
 program main
 implicit none
 
 double complex :: a(51,51),b(51,51)
 Integer::i,j
 real(8)::ht=0.01
 
 real(8) T1,T2
 
do i=1,51
   do j=1,51
      a(i,j)=cmplx(Sin(0.01*i),Cos(0.01*j))
   end do
end do

call cpu_time(T1)
do i=1,23497
   b(:,:)=(0.,1.)*ht/2.*a(:,:)
end do

call cpu_time(T2)

write(*,*) sum(b)
print '("Time = ",f12.9," seconds.")', T2-T1

end program main
POSTED BY: xslittlegrass *
9 Replies
Posted 11 years ago
I feel it would be convient if we can talk directly with the real code, so here is the code in the bottom of the post. Basically I'm trying to use LinearSolve itorately solve the Shrodinger's equation, and trying to tune the code to get some speedup.
The structure of this function is:
Do[Ct=LinearSolve[LeftMtx[t], RightMtx[t].Ct,{...}]
where LeftMtx is the compiled function HL and the RightMtx.Ct is the compiled function HRDotCt[t,CC]. My problem size is that energyLs has a length of 51, and diplMtx is a 51*51 matrix. I profiled the function and the part of the results is like this :
 Calls Time Evaluation
 1 11.111 TDSEPropagator[energyLs,diplMtx,testLaserParaLs,C0,ht,{211.613,240.71}];
 1 11.111 TDSEPropagator[energyLs,diplMtx,testLaserParaLs,C0,ht,{211.613,240.71}]
 1 11.1089 TDSEPropagator[energyLs_List,diplMtx_List,LaserParaLs_List,C0_List,dt_,solveRange_List]
 1 11.1088 Module[{H0,EE,lth,Cls,tls,Ct,<<5>>},lth=Length[energyLs];H0=DiagonalMatrix[energyLs];EE=Compile[{<<1>>},Evaluate[<<1>>]];Cls=Table[0.,Evaluate[<<1>>]];tls=Range[Apply[<<2>>]];Ct=C0;<<5>>]
 1 11.1055 lth=Length[energyLs];H0=DiagonalMatrix[energyLs];EE=Compile[{{<<2>>}},Evaluate[Sum[<<2>>]]];Cls=Table[0.,Evaluate[Join[<<2>>]]];tls=Range[Sequence@@Append[<<2>>]];Ct=C0;<<5>>
 1 8.8298 Cls=Table[Ct=LinearSolve[HL[<<1>>],HRDotCt[<<2>>]],{t,tls}]
 1 8.82979 Table[Ct=LinearSolve[HL[t],HRDotCt[t,Ct]],{t,tls}]
 11748 8.77454 Ct=LinearSolve[HL[t],HRDotCt[t,Ct]]
11748 8.65778 LinearSolve[HL[t],HRDotCt[t,Ct]]
1 2.22843 Interpolation/@Transpose[Thread/@Transpose[{<<2>>}]]
11748 1.73315 HL[t]
11748 1.59476 HRDotCt[t,Ct]
1 0.452539 Transpose[Thread/@Transpose[{tls,Prepend[<<2>>]}]]
1 0.429521 Thread/@Transpose[{tls,Prepend[Drop[<<2>>],C0]}]
1 0.016248 IdenMtx=Table[If[i==j,1.,0.],{i,1,lth},{j,1,lth}]
1 0.016243 Table[If[i==j,1.,0.],{i,1,lth},{j,1,lth}]
1 0.014355 HL=Compile[{{t,Blank[<<1>>]}},IdenMtx+<<3>> Plus[<<2>>],CompilationOptions->{Rule[<<2>>],Rule[<<2>>]}]
1 0.013833 HRDotCt=Compile[{{t,Blank[<<1>>]},{CC,Blank[<<1>>],1}},(IdenMtx+Times[<<2>>]).CC,CompilationOptions->{Rule[<<2>>],Rule[<<2>>]}]
you can see that the in the total 8.77454 seconds for 11748 of itoration of
Ct=LinearSolve[HL[t],HRDotCt[t,Ct]]
it takes about 3.2s for 11748 times itoration of 
HL[t] and HRDotCt[t,Ct]
which is about the 1/3 of the total time. This is where I think maybe easy to get speadup, since it has such a big performance difference compared with fortran.

You suggest to compile the whole thing, but I think LinearSolve is not compilable, and the call to the MainEvaluator may slow the compiled function (suggested here and here). 
And I don't quite understand your second point in your last reply. Could you have a look at the code and directly refered to the code about where I would able to get some speedup? Thank you very much for your help!


Here is the code:
     TDSEPropagator[energyLs_List, diplMtx_List, LaserParaLs_List, C0_List,
        dt_, solveRange_List] :=
      Module[{H0, EE, lth, Cls, tls, Ct, t, IdenMtx, HL, HR, HRDotCt},
       lth = Length[energyLs];
       H0 = DiagonalMatrix[energyLs];
       EE = Compile[{{t, _Real}},
         Evaluate[Sum[Block[{?, n, eFunc, tc, E0, ?0},
            {?, n, eFunc, tc, E0, ?0} = LaserParaLs[[m]];
            Piecewise[{{E0*eFunc[(? (t - tc))/(2 n)]*
               Cos[? (t - tc) + ?0],
              tc - n ?/? <= t <= tc + n ?/?}}, 0.]
           ], {m, 1, Length[LaserParaLs]}
          ]]];(*compiled piecewise function*)
      
      Cls = Table[0., Evaluate[{t}~Join~solveRange~Join~{dt}]];
      tls = Range[Sequence @@ Append[solveRange + {0., dt}, dt]];
      Ct = C0;
      IdenMtx =
       Table[If[i == j, 1., 0.], {i, 1, lth}, {j, 1,
         lth}];(*identity matrix*)
      
      HL = Compile[{{t, _Real}},
        IdenMtx + (I dt)/2 (H0 + diplMtx*EE[t + dt/2]),
        CompilationOptions -> {"InlineCompiledFunctions" -> True,
          "InlineExternalDefinitions" ->
           True}];(*if possible I would like to speedup this part*)
      
      HRDotCt =
       Compile[{{t, _Real}, {CC, _Complex,
          1}}, (IdenMtx - (I dt)/2 (H0 + diplMtx*EE[t + dt/2])).CC,
        CompilationOptions -> {"InlineCompiledFunctions" -> True,
          "InlineExternalDefinitions" -> True}];(*also speedup this part*)
   
        Cls = Table[
        Ct = LinearSolve[HL[t], HRDotCt[t, Ct]]
        , {t, tls}];
      Interpolation /@
       Transpose[Thread /@ Transpose[{tls, Prepend[Drop[Cls, -1], C0]}]]
      ]

The input parameters for this profile are:
energyLs={0., 30.6061, 31.7552, 33.0657, 33.5432, 33.9684, 35.3902, 29.9259,30.3731, 30.5524, 30.6219, 32.8307, 32.9402, 33.0805, 33.1013, \
33.7457, 33.797, 33.8399, 33.9972, 34.0047, 34.2123, 27.4233,27.8641, 32.1649, 32.3356, 32.7108, 32.7386, 32.9177, 33.5231, \
33.6909, 33.7337, 33.7467, 33.9231, 34.2059, 34.2128, 34.3893,30.3114, 30.4774, 30.5899, 32.9143, 32.9527, 33.1013, 33.7457, \
33.7494, 33.8283, 33.8444, 33.9259, 34.0069, 34.2123, 34.2142,34.3909};
C0 = Prepend[Table[0., {Length[energyLs] - 1}], 1.];
testLaserParaLs = {{29.6265`, 5.`, Cos[#1]^2 &, 14.548390422248389`, 2.5982114443438456`*^27, -(\[Pi]/2)}, {2.37535`, 11.`,
    Cos[#1]^2 &, 226.16134201858867`, 4.504198671763965`*^28, -(\[Pi]/2)}};
ht=00248;
diplMtx=Uncompress["1:eJztm3lUzlvY/pGphESmlKFShiQy7IeSMiUnQkiJSMYMN5WQIVMqaSDCcRKlMoRU9iMNGuiUOkWkJEOK4hhKIfTus9bvXGf9/N271vO28of1dRrOevb93fdwXZ+7/9L1sx1HtmjRYlNr8deMVZs2////ohb/78//yYedo7r+yCkP59TL3alDTX02o6XFn7enq17i5Hz79TIT9WuMbIcP/9Ms4CQjY58tyQobohntTezh9zM1nNGV6Cpf751hnHZxtfH99EMZbTexm9pyTgCjgUn242JzxPc4V89dvCn1DKNJr6qSzmb5cDrzulx+zrIQRq63v4QVRJ5gMnIaTSmwIQN/M3zQ+z2jut7Fa9VTboszruj/2qmkktHM5Fkjg/SfM3LRiHdpHfuK0VP3bfuX7/ubUdv194OWWdRxmpZy1t7cLoVRyO7hW1LkxY+Pr4/8a/KsZEZ3Pr+6kp2Sz+jEwuc1ixWLOHU/q3WwYehFTttV1W487JrGaM9vZ8P72GVxGTmNphTYuSduFU1aEskpUc/eZfeXF4w86v2NnaOKOR2kc0FnR1cxmsWMpijQJ04j1xTc/2O1nJTmJChpZlyq5/T+uWf3SbfSGS0I754UrxPLKcjWr5v9j0BGz8bNHDPhvLjD9WNVq2bIX+eUqrZqgZJVCCezYwdPdV8dw2iX7cb5tyW2zTe28R+8lZ6WGdXGMZrY7XGeqvoRTvk9q1/27v9DxDO2RdTStM+M7g2pKH/8WkTYz7pkyO5bdYwWDcpIqc1rIaWVrT8sb58hvnltRqd5RodqGTkUHlGfeeIToxvnld0k1S0lZJpkFuWzqZWUVKaFfJ6gnsvpbxWHlWsNXjJ6MP3+jsn7Xzff2MZ/6DFHweVsz82M0k5uWjekLInTk6+Xw0YYFDF6NF/6KaLDB06tN6hUdR30jFGvNY8qhn3/yKjg+4Gh46pF9Oz3HuiutkpLSodPxel7V//g9LPXSIfXsXIS0lQvMxsb0VVCj9rEVP10U5DQpNE3yqWaBYw6dOh18ENsBaeZMbNMHLulN9/Yxn9IvP06u3i26GP0+pa2NutzgpPeaOPi7xvyGH2tSns95pY/owuhvXnS6FhRSDN22HZcKGrjwyvVk8ccz2KkMNHoe1yIuJZtF9lZvD/TQULu7mqLLNt9YLR84q273uqtJaTxI3PqGutuUhqjne44aGxLKR1dzSt+1ChJyP/EU8+8QkWpjJxGUwqs7jvt0nC3M5x6rPMPMY0TXbHalzgv+98ec7pdnrB3sJOouk8PFbWNmf2AU2lhVB+tja9Ewxxk095f4Tmn8XONexSffsPpqVHu1AmhIl2nJdYNtKx6wUn+msVfzr6iMk9S3m93LO4jp8LsthaGMeWM0sOnf+447gunkpYt9SbwL82puPEfxmcX6mgFiOzKc8YdttgqrqXx7ADl0+9EZLoH5vVV+SHi+Uguedbri18YLSsodLMIFf8l/8HvW9UqCjk9eDjP6b5iquiHEnx2/3wv5bSmHW087yjytt7x/ttSW4p2am3il/iPLdM4JelP0zD7LY/Thk/r2h1IEG/T6+27El66xjUHtvEfAgsk/SZ8+Sx64EUNgYmhj8SUU/9KctqimlFF7LXAgHKRk+fEmFnOOHiK04ERtaM+mYtLqHohp8ryiBiASu60WL3QUHzPGxOtEp2eiaKpdk34k25O4XR0IjMpMHjI6fdsr2MTt2Uy2rl7r/VNtd2cMm6MU1XWucvpk/KorWUuN5trbOM/TN7xOLXe6y9GO2bqZte2E6Npm0GHPEJTKjnl+g/M9pgqplYHtwND10REc3Kcr+x4w7iGkUn5W41x+iJLn7pfELE9LIeTecu7pVcNRPTabIiaE/Nc1OpsM5UuSvNEqxThmqzXxeEOo295c/qr79kr0vXIaWlrtUXPpKtSOmNV/I3mwDb+Q1mv/CwDdzF5BKuMeCIXLC4q0x6S6m9+mtMsRbMG44hqTp6v12lvvfI3p5DHVhOWl4j584lmxbfOo8WXiuXuHbGbcU0k8GA5pbFJ+ZxanR/+6PRSEWrtHn0nxBw6xSh6vsfXmeEFnL5Khz+5tDqZkVaU4bol1UGMDg5uE3LoXEZzKm78h8X2X9d4Hl3FaLVHzgTd9+KmbUyMVTKY08Doo8Pj9IwOz8XUOvmizoF1baTklmnuvOpaCaMj7ZZ3dHoiXofZNfYWj3fWcBrabsnIiloR2KT0LdWn3NtLKNlKwcnMs4RTlEvhIu2Zjzg9v7S3umKamJ/WW/7+0vtpMic1S7/Nw1VSmwPb+A9Jni9sikLjGJlP/t3Tc9pVEdjlE5epVnSU0E2TFkEa+aLj/Wod1zZcN4FTjfIiPxOH1lLK1/HZ+dy1kpG+r6Lm+LAOUnpj3OL7402iDfKo+V0uYqyYh12dCuqUz7WTUMLoN05WPe5zSlE6tfnGqGhGkxQntUrzesdo1OIWnXumJzQHtvEfzB+pLh31JlBcuXka3voqopCqjKtM6Wblx2ib/4Qrg86LMIYec8kYnp3A6Ja14ZZlBmWc3A9EttOVr+V0+eN8j9HLRQttn+82vbiVGGjr+y3qsKlONEZ3XuwKSHAUwX947PJB70LRVC/4GX8zXnJAvA4Bn1XazRU11m3W33ZLtJ4119jGfxjPFi1dUhvFaNyIEseRUaJ+qi7xC67s8JjRH+v+HKs0WEFCTro3lbv4i/nT8E52sdzUc4yeD7Kui+rRQkL7xj66NrX0HqPyWby7Y3k7KSUYdMs6Wyum1iXGOim6BocZfTrVJffa0ZZSetnFY+83+fuMIhesnzD36XlGd5wy+viblDff2MZ/sB38wWSAMzH684yXZis2S9yi+Ll3I94HMuq7ODlhzoZgRo4G4Vani1pLKCjNXzvTSXTOu6NXWaudER2vxoEbn9pYXeH0YYyX/NVVYkY1jgjucqmui5TiNpkvPbVLZPLUEAfPhm3pnHqczplntOsZJ0uTGedyO2eLgfZ6Rlqo04HmwDb+g9o4n77eQ/ZxylE46Vd5z4eTfMDg+6WqF0Xre48uXlt0i5Pv0UuFnaaIaym3Q0F7VMVFTuE5HX/OfRXLySpsfMLW9/ISyjl4e+iGAdmMfls6LfrT215SWmFs8mBRkkjFDiphlpbfxWxkPeiDdbV9VykNPJUUYPtTJPDiW3+vdh3d3Dz9LzwEzijp1DlWXMtllSsnzbA4zmlAsIXRCmdxG7t9C69vs4RzGqI5akH5vVGc0oJfvtmVc4WRgvK0nfoakZzWJm0Kn+6jKqFJH6rWLQwSHa/5s8NeMR89GaWafH+6/2x3KRU8rDFN2i9+4WSjj8+Hloskv/BSQP3Xr8oSKvywwjvt0u3mwDb+w55wt8Xpius4lTON8LcffDl9ytp+0SXdnNFpT1fDIWV3OAXk+1bZfbRhtFJRUXnGsGBOQZZ7XDLzRUHWm+YQ4qaZxmjY8053c/eIDFzhrbZmz/x/3gvdh2edG15yMjSa77dvp7yUxr96N9DnUCqj9wfrpS9j3jC6+2XxMaYjJ5GR02hKge3kHXd9d6o3o3s35393LQhiNLhoscuTh6GM2qrknmv/Oo/Ry/lDVbQ2XWWU49Rq7E1Tb06RJvOHalVf5/R61moNTysvRp2Kv/W6ZdBTQvFXi/o5nCli9HbqHd1yH5HAjx0bWn7uSRcJqSnFrtEq/c7ox5b2Orb/dMVxK1Z+t0np0OzuNP7D/NzqUGv39YzGns6g3I3ioprpHp64omQ/o6zwjQf1VEXQnAfuWjw1s0o0sYZdXT9PeMqpd+7mQ21GxDLa9+XvtQklxxll7u00rF1qFiOVvVOiXmQqSmh4ZVeVdVo/OXnONk3QrhPN8MjQ8wU7dftLaEW21vh1nftKKbAhcrT5hiQZTcXgwcAPgTcBnwA/G/4n/DL4K9Djod9C74M+9J+e8O/8iXkF/S36IdRP5FvcT8RTNk6UnuyJcTPx/MHpi3PMELcfbzmV+urv+rpdjE2/aUSNOJcSwyn4Sdnlqvo4Tj1tys0VXW9wkjRcLnPe7spoiefupGe7nThtG+T86fAJL07LEma+Tdn1h/iFXu/fziwXn31xj75ationOb0s7x6t2ErE4uPlEOUTP1dycot5UpJrtJr/EljwfeDBwA+BNwGfAD8b/if8Mvgr0OOh30Lvgz4EPQHzJ+YV9Lfoh1A/kW9xP2UksB11VNqfGlzMqMNblRfDXb8y8lx4PdRthLjwgcPje/e/EM0o/Lr7zABbURHo3ZKqYYMTOR3clux1Q3cNo1j7lLhtLWcyWu7TdsfJLuJdVoxaPuRn3hHxSTPsjDLzLEXeqUiM33JK3KRdJnzfsJgVjG4MvZOyef0MUX0crCvnfbBgvwQWvCb4PvBg4IfAm4BPgJ8N/xN+GfwV6PHQb6H3QR+CnoD5E/MK+lv0Q6ifyLcyEli5jrM8yw994RStbbLGIqWW01Xbhx2vx6SJd/nmEJdaPyUppfrNUrm4UBTkEdcOZ2VWiglZ13vHNNuVEYyqjXVqVXc7MbrC25t6vM3/ZyBLMNb5UMApYVjH8zYNi8StPh7/x7YpYkSb+jxSboZUlI+jQdL8yHmTOL2pyH866fLUX28s+Fvwmv/xff/yYOCHwJuAT4CfDf8Tfhn8Fejx0G+h90Efgp6A+RPzCvpb9EOonzIS2OlHl39+P6BQpLbvf9Y6ZrxmdNZrbvqhbg2cTpa7+Bp1E3fCOWpUwZ39tYxCgsdPyZ+tIiXDzYbpZUZ7OJWk+b7T9LrAaVZDz/FHJ6aIU935ZLPjJhGUvBdzBrSoFhe1zLpsamntE06OO29HzdvejdOVb4FbNHVFRjtW1LtPvJX/r4EFTw3+Frwm+D7wYOCHwJuAT4CfDf8Tfhn8Fejx0G+h90Efgp6A+RPzCvpb9EMyEtiaorXpH4eIN7f/06ixWQ1VnIzczwQb1Iusd6t37MQL3T4w2lr6h9pVyxrxetqrvblfL6pufUjJ2D0DRL+993GLFL0O4kuL5UweeKuIHDdOXvL8ooqvaLzlT6r2XnCP04vfh5TLRXtwmqNxp7LjK/GC1NgsibS7VcTpkEn2i9Nu8b8GFnw8eGrwt+A1wfeBBwM/BN4EfAL8bPif8Mvgr0CPh34LvQ/6EPQEzJ+YV9Dfykhg8zJje7zrLOqRbnXWlAtrxZULvnYnoLL9Q5Gk4ut/eAWJt/tYXJ926pmisgSG5gRq6t1k5K63f5HWoFYS+qOPUm3CEGUJjXAcuO9cUpg4Xpe2w73Hb+e0L2/+3a+54qcGDdTKPZwrktTifd9yZlpXiM5lc5B253siOvblJUP2RST/WmOx7wA+Hjw1+FvwmuD7wIOBHwJvAj4Bfjb8T/hl8Fegx0O/hd4HfQh6AuZPzCsyEljTuUeWHewpjiXm64cOF3JFlxm7y5QrtxdD0uFlVtYLcr5xUUjba/QuyRGFqW/v35Yef8JI6VXgpcdHL4ufWrGofoF8CSfTllf2H7cJE2l20/n1FwrsGd2VrjuxV11FQvNsVM8vHiDe7h63ujzTeneQkf+g5AdLbooy/tfRiW3zvOt/DSz2V7DvAD4ePDX4W/Ca4PvAg4EfAm8CPgF+NvxP+GXwV6DHQ7+F3gd9CHoC5k8ZCWxt6Il9GfMecrpY/tXghPojMSmyE3I2MdcZHdFT33V6l6KU7uv8zL2sLj67a7CfWqL2XdHWzjFKVtW/zanhbX50zQkjRsN3mBs88FCTUgP/cd2gZ29xUDq+zmXTxVyoYFM+KXV6OqekPknxDYqXGPnu2JmfayGCPyd9Wf7Pdgq/Bhb7SNhfwb4D+Hjw1OBvwWuC7wMPBn4IvAn4BPjZ8D/hl8FfgR4P/RZ6H/Qh6AkyEtjBYxSDnWtSRCXU1DU1LhKXY7f/IhePUeIS3te0q+M1IiB9wqKNsv3EOBieedgyb2NXKfUc7RdWnSYCG+4XZLmqVHxpaVBJ4i0bcaqRg4MdvuuLzuXGC8PUdfqiGdGb976yVKIhJe/SxK2F78T/InL8hej+3cSpunYKeGJfde7XGov9MuwjYX8F+w7g48FTg78Frwm+DzwY+CHwJuAT4GfD/4RfBn8Fejz0W+h90IdkJLBjktdOr+19h9Gezv0GDHQUH7nzg/w+7ldF1VhWtty4lY9oT4It5WuuTxUFbmnY5VG0UeSm4pzvFqcDxLU0s7Y+/EepeM3Z5YiTbp7ioda+YoNzB9Fu1lgFWBmQyGiprWIOeDmKzz54Kr1MTO8uobL8sS21V4tR+eZ+BXMrF+mvgcW+IPbLsI+E/RXsO4CPB08N/ha8Jvg+8GDgh8CbgE+Anw3/E34Z/BXo8dBvoffJSGDbzjm0Kt4wg5OZuY97kfMZUcXuWpy9Wy06vzbuoydHDxNv7thpppMH5z5glFKvlrfmuZSR+iYbt5/TxKi3PzIry2ZZDylpxWcZPLQVn8tHs2d6r29iWHDeZnlVc4copFc/FtaN7yMy2n3n+V2tzoqf6hLSp9u8ZJG/Ktzz+32dd/3XVIz9T+wLYr8M+0jYX8G+A/h48NTgb8Frgu8DDwZ+CLwJ+AT42fA/4ZfBX4EeD/1WRgL7yK5f//ScKEadOm10m/Y4hNMYLZ02t3zjxETwYv8lm/tJjKx/PtRbGXKZ0/Dyqlam3mXiXe7a0eFaiuhB1o+syffXF4N6needb+02/MkprG/h3Tx/8eNmRxJvjxjQT0p1d1NvXs98J3qmLwZxRd1SGbU/bjXlca043vzfC1c++NZH8ktgsc+L/U/sC2K/DPtI2F/BvgP4ePDU4G/Ba4LvAw8Gfgi8CfgE+NnwP+GXwV+BHi8jgZV/tKtnZEU8p7Iia7O34aLRMFkdbRve6zyjY5bBHrqbxfFv11z7RieiVJxYa9uQ3mtEQjS9ob/4Qfs3Ykawm2vWT16871oRKyqLGkQLvf2F87b7iW0kZJA1bvuCS+I3jy1xGhOoLi68pZ2zx9Gr2hJquX31SI9JSlJaluVw5uGEB7+mYuxnY58X+5/YF8R+GfaRsL+CfQfw8eCpwd+C1wTfBx4M/BB4E/AJ8LPhf8Ivg78iI4G98u5k0jNTcfxhhqYaJ/X2iXiuDLz5qr2YSI0KwyRTWosKJUlZHjstUPTJNW/Kn7mFRzKauGpQa9+/H4tB4NrOs58qxfEW72k9P+eLCOzl5b6+hltFRV3Yc+/ewxXinu/xrFMdGS8Ce1hOV//uWXGlSmvS181uEO3UcROlGnnjv34NLPbtsZ+NfV7sf2JfEPtl2EfC/gr2HcDHg6cGfwteE3wfeDDwQ+BNwCfAz4b/KWN+Ga2P1llv/I96N3nc5rpHCkdEHKoifWPjIxgd7epzbPM18ea+udFxU1xP0VPm25nIPcv8k9GQTHXTt4ZuomX4oaKSlnpKjA/1Nq2PvBNXik49jc6tO89p0MVhxe6hbzlpqrXf3GKCaMIq+yR0k5cT1Xtp2RZV23RxPkuXBl4OTNeTNiU/Fn4Z/BXo8dBvofdBH4KegPkT8wr6W/RDqJ/It7ifiKdsnEZTCiz8T/hl8Fegx0O/hd4HfQh6AuZPzCvob9EPoX4i3+J+ysZpNKXAws+G/wm/DP4K9Hjot9D7oA9BT8D8iXkF/e1//dC/9RP5VjZOoykFFnwC/Gz4n/DL4K9Aj4d++5/e968+BD0B8yfmFfS36IdQP2XjNJpSYMGbgE+Anw3/E34Z/BXo8dBvofdBH4KegPkT8wr6W/RDsnEaTSmw4IfAm4BPgJ8N/xN+GfwV6PHQb6H3QR+CnoD5E/MK+lvZOI2mFFjwYOCHwJuAT4CfDf8Tfhn8Fejx0G+h90Efgp6A+RPzimycRlMKLPg+8GDgh8CbgE+Anw3/E34Z/BXo8dBvofdBH4KegPlTNk6jKQUWvCb4PvBg4IfAm4BPgJ8N/xN+GfwV6PHQb6H3QR+CniAbp9GUAgv+Frwm+D7wYOCHwJuAT4CfDf8Tfhn8Fejx0G+h90Efko3TaEqBBU8N/ha8Jvg+8GDgh8CbgE+Anw3/E34Z/BXo8dBvoffJxmk0pcCCjwdPDf4WvCb4PvBg4IfAm4BPgJ8N/xN+GfwV6PHQb2XjNJpSYLHvAD4ePDX4W/Ca4PvAg4EfAm8CPgF+NvxP+GXwV6DHy8ZpNKXAYn8F+w7g48FTg78Frwm+DzwY+CHwJuAT4GfD/4RfBn9FNk6jKQX2330k7K9g3wF8PHhq8LfgNcH3gQcDPwTeBHwC/Gz4n/DLZOQ0/gc9TmqF"];

and the profile function I use is (in debugger mode):
TDSEPropagator[energyLs, diplMtx, testLaserParaLs, C0, ht, {211.61295159634025`, 240.7097324408371`}]; // RuntimeTools`Profile
POSTED BY: xslittlegrass *
@xslittlegrass

did you ever figure out whether or not your MMA code can be optimized so that it can compete with the FORTRAN version?
Posted 11 years ago
It's not big, see my updates to the post.
POSTED BY: xslittlegrass *
From the computation below, I infer that you are claim Fortran is getting upwards of 10^13 flops/sec.

In[66]:= 23497*(Times @@ Dimensions)/0.0000050

Out[66]= 1.22231394*10^13

That is one very fast Fortran.
POSTED BY: Daniel Lichtblau
Posted 11 years ago
Hi Daniel Lichtblau,

Thanks for the answer. After a careful check, I did make a mistake in timing the fortran program. It seems that the loop part in the Fortran code is not excuted when there is no results dependen on it (explained here). The correct time is 0.05 second, so that means there is about 40 times speed difference. Is it possible to speedup the Mathematica code to be comparable to the Fortran speed or this is the most what we can get from Mathematica? Thanks again for your help and sorry for the mistake. I'll correct the mistake in the post.
POSTED BY: xslittlegrass *
If someone will want to check this all thoroughly, they will probably need your fortran code for benchmarking. Is it big?
POSTED BY: Vitaliy Kaurov
I can make some improvements so that this becomes closer to the Fortran code, both in speed and in methodology. First note that your iteration is not in the Compile code. So we can easily change that. Also it is best to evaluate the matrix argument  prior to those multiplications. So here is my variant.
ht = 0.01;
a = Table[
Sin[i] + N[i]* Cos[j], {i, 0., 0.5, 0.01}, {j, 0., 0.5, 0.01}];
HL = With[{ea = a},
Compile[{{t, _Complex}, {ii, _Integer}}, Do[t *ea, {ii}],
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}]];[/i][/i]

This gives a factor of 6 or so speed improvement on my machine, vs. the earlier code.

[i][i]In[103]:= HL[ht/2.*I, 23497]; // AbsoluteTiming

(* Out[103]= {0.289352, Null} *)[/i][/i]
POSTED BY: Daniel Lichtblau
Posted 11 years ago
Hi  Daniel Lichtblau,

Thanks for the reply. In my real code the matrix a is different in each iteration so it would be difficult to put the iteration in side compile. And if I put the limitation that the iterator should not be inside the compile function, is it still possible to speedup the code? Thanks. 
POSTED BY: xslittlegrass *
I think the real question is whether your real need is something that Mathematica's `Compile` can handle in entirety. If, say, you make modifications to the input matrix in each iteration, and if those modifications are within the realm of `Compile` capabilities (and I suspect they are), then there is no reason not to retain the iteration inside the `Compile`.

If that really is not feasible, then the next suggestion would be to have an outer Compile'd function that calls an inner one. The outer one would do the iterations and modify the input matrix, and the inner one would do the work to obtain the new resulting matrix. Or something along those lines.
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract