Why does Compile redo the calculation when overflow happens?

Question

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.

Answer 1

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}}]

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.