0
$\begingroup$

When implementing the Numerov method solving Schrodinger equation, I encountered this problem. In order to show the order of the evaluation, Print is added.

calU=Compile[{{x,_Real,1},{energy,_Real},{m,_Real},{a,_Real}},
Module[{i,node,xn,nn,phi,V,h,f,temp},
  Print["begin ini"];

  h=x[[3]];
  xn=Range[x[[1]],x[[2]],h];
  nn=Length@xn;
  phi=Table[0.0,{nn}];
  V=1/2 m a^2 xn^2;
  f=2m(energy-V);
  phi[[-2]]=h;
  node=0;

  Print["end ini"];

  Print["begin cal"];
  For[i=nn-2,i>1,i--,
    phi[[i]]=1/(1+h^2/12 f[[i]]) (2(1-(5h^2)/12 f[[i+1]])phi[[i+1]]-(1+h^2/12 f[[i+2]])phi[[i+2]]);
    If[phi[[i]]*phi[[i+1]]<0,node++];
  ];
  Print["end cal"];

  Print["begin renormalize"];
  temp=Total[phi^2];
  Print[temp];
  phi=phi/Sqrt[temp h];
  Print["end renormalize"];
{node-1,Transpose[{xn,phi}]}
]];

Evaluating Timing[calU[{0, 14, 0.0001}, 1.5001128419302405, 3.707412760373, 1.0];] gives

begin ini
end ini
begin cal
end cal
begin renormalize (* <==weird happens here *)
begin ini
end ini
begin cal
end cal
begin renormalize
1.124558754621281*10^312
end renormalize
{2.468750, Null}

Please be noted that when comes to the "normalize" part, the process jump to the initialization without warning of overflow (it actually happens, since the maximum number in double precission is about 1.8*10^308) or using uncompiled one to preceed.

The output becomes weirder when evaluate Timing[calU[{0, 10, 0.0001}, 1.5001128419302405, 3.707412760374, 1.0];] three times, the first two give

begin ini
end ini
begin cal
end cal
begin renormalize
1.16629*10^158
end renormalize
CompiledFunction::cfex: Could not complete external evaluation at instruction 114; proceeding with uncompiled evaluation. >>
begin ini
end ini
begin cal
end cal
begin renormalize
1.16629*10^158
end renormalize
{1.765625, Null}

even the "normalize" part is worked out they still jump to the "initialization" again; the third one gives the normal output as I expected (not the evaluation time)

So how to understand these?

Update

As blochwave point it out, the time cost is due to returning multiple results in compile and overflow.

Here comes two more related questions:

It seems that if overflow happens (it's turned off in default), Mathematica will drop the data already get and start over, why not taking over the data and using the main evaluation routine to do it?

I just find out that even if I set the RuntimeOptions -> "Quality" or turn on the "CatchMachineOverflow" , still there aren't any warnings. While "RuntimeErrorHandler" -> Function[Throw[\$Failed]] and catch it will get \$Failed. This inconsistence beats me again.

$\endgroup$
  • $\begingroup$ Does it happen without the Print? It's worth noting that Print is not a compilable function in Mathematica, and so the function calls MainEvaluate each time it encounters Print. You can check this by Needs["CompiledFunctionTools`"]; CompilePrint@calU. $\endgroup$ – dr.blochwave Apr 14 '15 at 7:00
  • $\begingroup$ @blochwave Actually I'm find this problem because of the long time it cost. If I remove the Print , I can only tell that the evaluation time is very close to the above. So my guess is yes. BTW, thank you for the link :) $\endgroup$ – luyuwuli Apr 14 '15 at 7:18
5
$\begingroup$

Looking at your code, the problem occurs when you're trying to return {node-1, Transpose[{xn,phi}]}.

If instead you run the following code, which only returns Transpose[{xn,phi}],

calU=Compile[{{x,_Real,1},{energy,_Real},{m,_Real},{a,_Real}},
Module[{i,node,xn,nn,phi,V,h,f,temp},
  h=x[[3]];
  xn=Range[x[[1]],x[[2]],h];
  nn=Length@xn;
  phi=Table[0.0,{nn}];
  V=1/2 m a^2 xn^2;
  f=2m(energy-V);
  phi[[-2]]=h;
  node=0;
  For[i=nn-2,i>1,i--,
    phi[[i]]=1/(1+h^2/12 f[[i]]) (2(1-(5h^2)/12 f[[i+1]])phi[[i+1]]-(1+h^2/12 f[[i+2]])phi[[i+2]]);
    If[phi[[i]]*phi[[i+1]]<0,node++];
  ];
  temp=Total[phi^2];
  phi=phi/Sqrt[temp h];
  Transpose[{xn,phi}]
]];

calU[{0, 10, 0.0001}, 1.5001128419302405, 3.707412760374, 1.0]; // AbsoluteTiming

(* 0.0657 seconds *)

It is much quicker.

The problem is in your last line, where {node-1,Transpose[{xn,phi}]} doesn't have consistent dimensions (node is an integer, and {xn,phi} a list), and so the compiled function instead calls MainEvaluate to return your data - thus losing all the benefits of compilation.

You can check this via:

Needs["CompiledFunctionTools`"]
CompilePrint@calU

(* [line] 115   T(R2)6 = MainEvaluate[ Hold[List][ I12, T(R2)8]] *)

If you do indeed want your function to return multiple results, have a look at How to force Compile to return multiple results?.


Finally, regarding the point about machine overflow, you can use the RuntimeOptions of Compile to catch it or not depending on preference.


Update

"The problem is if I change the upper bound of x to 14, i.e.Timing[calU[{0, 14, 0.0001}, 1.5001128419302405, 3.707412760373, 1.0];] , it still takes the same time on my machine. In my code posted, since it actually work it out, why bother to start over?"

Let's take a look at the timing behaviour as the upper bound of x increases:

ListLogPlot[
 First@AbsoluteTiming[
     calU[{0, #, 0.0001}, 1.5001128419302405, 3.707412760374, 
      1.0]] & /@ Range[15], PlotRange -> {{0, 15}, {10^-3, 10}}]

enter image description here

There's a significant jump for x = 14. Why? Because Total[phi^2] becomes larger than the maximum machine number at that point. If you want to use an upper bound of this or larger, then compilation (and thus machine precision) probably isn't of use to you.

$\endgroup$
  • $\begingroup$ Thanks a lot (and sorry for my late response). The problem is if I change the upper bound of x to 14, i.e.Timing[calU[{0, 14, 0.0001}, 1.5001128419302405, 3.707412760373, 1.0];] , it still takes the same time on my machine. In my code posted, since it actually work it out, why bother to start over? $\endgroup$ – luyuwuli Apr 14 '15 at 8:35
  • $\begingroup$ @luyuwuli in that case don't bother using compilation - it won't speed anything up for you because you want to use numbers larger than $MaxMachineNumber. $\endgroup$ – dr.blochwave Apr 14 '15 at 8:55
  • $\begingroup$ Did you mean that if overflow happens (it's turned off in default), Mathematica will drop the data already get and start over? If so, why not taking over the data and using the main evaluation routine to do it? Are there some technical considerations for this? $\endgroup$ – luyuwuli Apr 14 '15 at 9:08
  • $\begingroup$ I'm afraid I don't have an answer to that. $\endgroup$ – dr.blochwave Apr 14 '15 at 9:10
  • $\begingroup$ I just find out that even if I set the RuntimeOptions -> "Quality" or turn on the "CatchMachineOverflow" , still there aren't any warnings. While "RuntimeErrorHandler" -> Function[Throw[$Failed]] and catch it will get $Failed. This inconsistence beats me again. $\endgroup$ – luyuwuli Apr 14 '15 at 9:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.