# Is Mathematica 40 times slower than Fortran with matrixes?

Posted 10 years ago
7400 Views
|
9 Replies
|
9 Total Likes
|
 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}]; // AbsoluteTimingOut[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]= TrueFor some reason, a is not a packed array, but even if we use the packed verion, it's still slow: In[7]:= DeveloperPackedArrayQ[a] Out[7]= False  In[8]:= apk = DeveloperToPackedArray[a];  In[9]:= HL =    Compile[{{t, _Real}}, I ht/2*apk,     CompilationOptions -> {"InlineCompiledFunctions" -> True,       "InlineExternalDefinitions" -> True}];In[10]:= StringFreeQ[CompilePrint[HL], "MainEvaluate"]Out[10]= TrueIn[11]:= Table[HL[200.], {23497}]; // AbsoluteTimingOut[11]= {1.974200, Null} And it's not the problem of Compile:In[12]:= Table[I ht/2*a, {23497}]; // AbsoluteTimingOut[12]= {1.935287, Null} In[13]:= Table[I ht/2*apk, {23497}]; // AbsoluteTimingOut[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 doend docall cpu_time(T1)do i=1,23497   b(:,:)=(0.,1.)*ht/2.*a(:,:)end docall cpu_time(T2)write(*,*) sum(b)print '("Time = ",f12.9," seconds.")', T2-T1end program main
9 Replies
Sort By:
Posted 10 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}]; // RuntimeToolsProfile
Posted 10 years ago
 @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 10 years ago
 It's not big, see my updates to the post.
Posted 10 years ago
 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.0000050Out[66]= 1.22231394*10^13That is one very fast Fortran.
Posted 10 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 10 years ago
 If someone will want to check this all thoroughly, they will probably need your fortran code for benchmarking. Is it big?
Posted 10 years ago
 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 10 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 10 years ago
 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.