MathGroup Archive: August 2011 [00478]

[Date Index] [Thread Index] [Author Index]
Problem with NIntegrating a compiled function
To: mathgroup at smc.vnet.net
Subject: [mg121077] Problem with NIntegrating a compiled function
From: "Eric Michielssen" <emichiel at eecs.umich.edu>
Date: Fri, 26 Aug 2011 05:24:26 -0400 (EDT)
Delivered-to: l-mathgroup@mail-archive0.wolfram.com
Hi,

I am running into unexpected problems using NIntegrate on a compiled
function.  In the code below, I NIntegrate the functions "myint" and
"cmyint"; the latter is a compiled version of the former -- both routines
are listed at the bottom of this email.

ClearAll[dt,r,c,i];
dt=0.1;
c=1.;
r=0.95;
i=11;
Print["Integrate using noncompiled modules"];
NIntegrate[myint[r,c,dt,i],{r,0.95,1.05},EvaluationMonitor:>Print[r,"
",myint[r,c,dt,i]]]
Print["Integrate using compiled modules"];
NIntegrate[cmyint[r,c,dt,i],{r,0.95,1.05},EvaluationMonitor:>Print[r "
",cmyint[r,c,dt,i],"  ",myint[r,c,dt,i] ,"  ",myint[r,c,dt,i]
-cmyint[r,c,dt,i]]]

The noncompiled NIntegrate yields a nonzero answer. The compiled one yields
zero, even though EvaluationMonitor indicates the integrands evaluate to the
same numerical value for all values of the integration variable. Changing
the NIntegrate options (PrecisionGoal, MinRecursion etc) does not change the
behavior.

During evaluation of In[2288]:= Integrate using noncompiled modules [[I
print integration variable and myint]]

During evaluation of In[2288]:= 0.950796 0.000159148
During evaluation of In[2288]:= 0.954691 0.000169447
During evaluation of In[2288]:= 0.962292 0.000191239
During evaluation of In[2288]:= 0.973077 0.000225252
During evaluation of In[2288]:= 0.986018 0.000269762
During evaluation of In[2288]:= 1. 0.00032485
During evaluation of In[2288]:= 1.01398 0.000381631
During evaluation of In[2288]:= 1.02692 0.000407476
During evaluation of In[2288]:= 1.03771 0.000412319
During evaluation of In[2288]:= 1.04531 0.00040786
During evaluation of In[2288]:= 1.0492 0.000403319
During evaluation of In[2288]:= 0.950398 0.000158134
During evaluation of In[2288]:= 0.952346 0.000163168
During evaluation of In[2288]:= 0.956146 0.000173454
During evaluation of In[2288]:= 0.961538 0.000188989
During evaluation of In[2288]:= 0.968009 0.000208878
During evaluation of In[2288]:= 0.975 0.000231631
During evaluation of In[2288]:= 0.981991 0.000255507
During evaluation of In[2288]:= 0.988462 0.000278613
During evaluation of In[2288]:= 0.993854 0.000298889
During evaluation of In[2288]:= 0.997654 0.000314236
During evaluation of In[2288]:= 0.999602 0.000322896
During evaluation of In[2288]:= 1.0004 0.000326926
During evaluation of In[2288]:= 1.00235 0.000336681
During evaluation of In[2288]:= 1.00615 0.000353817
During evaluation of In[2288]:= 1.01154 0.000374006
During evaluation of In[2288]:= 1.01801 0.000392226
During evaluation of In[2288]:= 1.025 0.000405096
During evaluation of In[2288]:= 1.03199 0.00041149
During evaluation of In[2288]:= 1.03846 0.000412151
During evaluation of In[2288]:= 1.04385 0.000409176
During evaluation of In[2288]:= 1.04765 0.000405299
During evaluation of In[2288]:= 1.0496 0.000402774
During evaluation of In[2288]:= 0.950199 0.000157629
During evaluation of In[2288]:= 0.951173 0.000160117
During evaluation of In[2288]:= 0.953073 0.000165091
During evaluation of In[2288]:= 0.955769 0.000172408
During evaluation of In[2288]:= 0.959005 0.000181561
During evaluation of In[2288]:= 0.9625 0.000191865
During evaluation of In[2288]:= 0.965995 0.000202557
During evaluation of In[2288]:= 0.969231 0.000212766
During evaluation of In[2288]:= 0.971927 0.000221482
During evaluation of In[2288]:= 0.973827 0.000227732
During evaluation of In[2288]:= 0.974801 0.000230967
During evaluation of In[2288]:= 0.975199 0.000232295
During evaluation of In[2288]:= 0.976173 0.000235561
During evaluation of In[2288]:= 0.978073 0.000241993
During evaluation of In[2288]:= 0.980769 0.000251257
During evaluation of In[2288]:= 0.984005 0.000262586
During evaluation of In[2288]:= 0.9875 0.000275109
During evaluation of In[2288]:= 0.990995 0.00028799
During evaluation of In[2288]:= 0.994231 0.000300358
During evaluation of In[2288]:= 0.996927 0.000311187
During evaluation of In[2288]:= 0.998827 0.000319334
During evaluation of In[2288]:= 0.999801 0.000323855
During evaluation of In[2288]:= 0.975099 0.000231963
During evaluation of In[2288]:= 0.975586 0.000233592
During evaluation of In[2288]:= 0.976536 0.000236786
During evaluation of In[2288]:= 0.977885 0.000241352
During evaluation of In[2288]:= 0.979502 0.000246885
During evaluation of In[2288]:= 0.98125 0.000252926
During evaluation of In[2288]:= 0.982998 0.000259035
During evaluation of In[2288]:= 0.984615 0.000264752
During evaluation of In[2288]:= 0.985964 0.000269565
During evaluation of In[2288]:= 0.986914 0.000272986
During evaluation of In[2288]:= 0.987401 0.000274748
During evaluation of In[2288]:= 0.987599 0.00027547
During evaluation of In[2288]:= 0.988086 0.000277242
During evaluation of In[2288]:= 0.989036 0.000280721
During evaluation of In[2288]:= 0.990385 0.000285708
During evaluation of In[2288]:= 0.992002 0.000291784
During evaluation of In[2288]:= 0.99375 0.000298485
During evaluation of In[2288]:= 0.995498 0.000305369
During evaluation of In[2288]:= 0.997115 0.00031197
During evaluation of In[2288]:= 0.998464 0.000317724
During evaluation of In[2288]:= 0.999414 0.000322009
During evaluation of In[2288]:= 0.999901 0.000324346
Out[2294]= 0.0000312806

During evaluation of In[2288]:= Integrate using compiled modules [[I print
integration variable, myint, cmyint, and the difference between the latter
two]]
During evaluation of In[2288]:= CompiledFunction::cfsa: Argument r at
position 1 should be a machine-size real number. >>
During evaluation of In[2288]:= 0.950796  0.000159148  0.000159148
-1.29508*10^-16
During evaluation of In[2288]:= 0.954691  0.000169447  0.000169447
2.46385*10^-17
During evaluation of In[2288]:= 0.962292  0.000191239  0.000191239
-3.69713*10^-17
During evaluation of In[2288]:= 0.973077  0.000225252  0.000225252
-2.46927*10^-17
During evaluation of In[2288]:= 0.986018  0.000269762  0.000269762
6.17995*10^-18
During evaluation of In[2288]:= 1.  0.00032485  0.00032485  6.17995*10^-18
During evaluation of In[2288]:= 1.01398  0.000381631  0.000381631
6.17995*10^-18
During evaluation of In[2288]:= 1.02692  0.000407476  0.000407476
3.08456*10^-17
During evaluation of In[2288]:= 1.03771  0.000412319  0.000412319
-1.11022*10^-16
During evaluation of In[2288]:= 1.04531  0.00040786  0.00040786
-7.39968*10^-17
During evaluation of In[2288]:= 1.0492  0.000403319  0.000403319
4.93312*10^-17
During evaluation of In[2288]:= 0.950398  0.000158134  0.000158134
-4.31784*10^-17
During evaluation of In[2288]:= 0.952346  0.000163168  0.000163168
-6.1664*10^-17
During evaluation of In[2288]:= 0.956146  0.000173454  0.000173454
6.17995*10^-18
During evaluation of In[2288]:= 0.961538  0.000188989  0.000188989
2.46927*10^-17
During evaluation of In[2288]:= 0.968009  0.000208878  0.000208878
-1.23599*10^-17
During evaluation of In[2288]:= 0.975  0.000231631  0.000231631
6.20706*10^-18
During evaluation of In[2288]:= 0.981991  0.000255507  0.000255507
-8.63567*10^-17
During evaluation of In[2288]:= 0.988462  0.000278613  0.000278613
-6.16911*10^-17
During evaluation of In[2288]:= 0.993854  0.000298889  0.000298889
5.55654*10^-17
During evaluation of In[2288]:= 0.997654  0.000314236  0.000314236
-7.4051*10^-17
During evaluation of In[2288]:= 0.999602  0.000322896  0.000322896
-8.64109*10^-17
During evaluation of In[2288]:= 1.0004  0.000326926  0.000326926
-8.63025*10^-17
During evaluation of In[2288]:= 1.00235  0.000336681  0.000336681
-7.39968*10^-17
During evaluation of In[2288]:= 1.00615  0.000353817  0.000353817
-1.84856*10^-17
During evaluation of In[2288]:= 1.01154  0.000374006  0.000374006
-8.01225*10^-17
During evaluation of In[2288]:= 1.01801  0.000392226  0.000392226
6.17995*10^-18
During evaluation of In[2288]:= 1.025  0.000405096  0.000405096
1.04951*10^-16
During evaluation of In[2288]:= 1.03199  0.00041149  0.00041149
-3.70255*10^-17
During evaluation of In[2288]:= 1.03846  0.000412151  0.000412151
-1.84856*10^-17
During evaluation of In[2288]:= 1.04385  0.000409176  0.000409176
-1.23382*10^-16
During evaluation of In[2288]:= 1.04765  0.000405299  0.000405299
6.16369*10^-17
During evaluation of In[2288]:= 1.0496  0.000402774  0.000402774
7.39968*10^-17

During evaluation of In[2288]:= NIntegrate::izero: Integral and error
estimates are 0 on all integration subregions. Try increasing the value of
the MinRecursion option. If value of integral may be 0, specify a finite
value for the AccuracyGoal option. >>

Out[2296]= 0.

The routines for myint and cmyint are below -- they both call other
noncompiled and compiled routines. myintdet and cmyintdet.  

ClearAll[myintdet,cmyintdet,myint,cmyint];

myintdet[t1_,t2_,a_,dt_,to_,j_?IntegerQ]:=Module[
{at2,at1,lat},
at2=Sqrt[-1. a^2+t2^2] ; at1= Sqrt[-1. a^2+t1^2];lat=
Log[(at2+t2)/(at1+t1)];
Switch[j,
1,0.05555555555555555` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+9.` (a^2 lat-at1
t1+at2 t2) (dt+0.5` to)+3.` lat (dt+to) (2.` dt+to) (3.` dt+to)+3.`
(at1-at2) (11.` dt^2+12.` dt to+3.` to^2)),
2,0.16666666666666666` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
(dt+to) (2.` dt+to)+1.5` (-1.` a^2 lat+at1 t1-at2 t2) (2.` dt+3.` to)+3.`
(at1-at2) (dt^2-4.` dt to-3.` to^2)),
3,0.16666666666666666` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+3.` lat (dt-to)
(2.` dt-to) (dt+to)+3.` (a^2 lat-at1 t1+at2 t2) (-1.` dt+1.5` to)+3.`
(-at1+at2) (dt^2+4.` dt to-3.` to^2)),
4,0.05555555555555555` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
(2.` dt-to) (3.` dt-to)+9.` (a^2 lat-at1 t1+at2 t2) (dt-0.5` to)-3.`
(at1-at2) (11.` dt^2-12.` dt to+3.` to^2)) ,
_,
0.
]];

cmyintdet=Compile[
{{t1, _Real},{t2, _Real}, {a, _Real},{dt, _Real},{to,_Real},{j, _Integer}},
Module[
{at2,at1,lat},
at2=Sqrt[-1. a^2+t2^2] ; at1= Sqrt[-1. a^2+t1^2];lat=
Log[(at2+t2)/(at1+t1)];
Switch[j,
1,0.05555555555555555` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+9.` (a^2 lat-at1
t1+at2 t2) (dt+0.5` to)+3.` lat (dt+to) (2.` dt+to) (3.` dt+to)+3.`
(at1-at2) (11.` dt^2+12.` dt to+3.` to^2)),
2,0.16666666666666666` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
(dt+to) (2.` dt+to)+1.5` (-1.` a^2 lat+at1 t1-at2 t2) (2.` dt+3.` to)+3.`
(at1-at2) (dt^2-4.` dt to-3.` to^2)),
3,0.16666666666666666` (2.` a^2 (at1-at2)+at1 t1^2-at2 t2^2+3.` lat (dt-to)
(2.` dt-to) (dt+to)+3.` (a^2 lat-at1 t1+at2 t2) (-1.` dt+1.5` to)+3.`
(-at1+at2) (dt^2+4.` dt to-3.` to^2)),
4,0.05555555555555555` (2.` a^2 (-at1+at2)-at1 t1^2+at2 t2^2+3.` lat (dt-to)
(2.` dt-to) (3.` dt-to)+9.` (a^2 lat-at1 t1+at2 t2) (dt-0.5` to)-3.`
(at1-at2) (11.` dt^2-12.` dt to+3.` to^2)),
_,
0.
]]
];

myint[r_?NumericQ,c_?NumericQ,dt_?NumericQ,i_?NumericQ]:=Module[
{a,to,limits,tab,t1,t2,j},
a=r/c;
to=i dt;
limits={{to,dt+to},{-dt+to,to},{-2 dt+to,-dt+to},{-3 dt+to,-2 dt+to}};
tab=Table[
If[a>limits[[j,2]],
0.,
{t1,t2}={Max[limits[[j,1]],a],limits[[j,2]]};
myintdet[t1,t2,a,dt,to,j]
],
{j,1,4}];
Total[tab]
];

cmyint=Compile[
{{r,_Real},{c, _Real},{dt, _Real},{i, _Integer}},
Module[
{a,to,limits,tab,t1,t2,j,u=0.},
a=r/c;
to=i dt;
limits={{to,dt+to},{-dt+to,to},{-2 dt+to,-dt+to},{-3 dt+to,-2 dt+to}};
tab=Table[
If[a>limits[[j,2]],
u=0.,
{t1,t2}={Max[limits[[j,1]],a],limits[[j,2]]};
u=cmyintdet[t1,t2,a,dt,to,j]   (***  FIX to cmyintdet ***)
];
u,
{j,1,4}];
Return[Total[tab]]
]];

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
Follow-Ups:
- Re: Problem with NIntegrating a compiled function
  - From: Daniel Lichtblau <danl@wolfram.com>
Prev by Date: Re: CUDA, Nvidia and Linux issue
Next by Date: Re: decoding inbuilt function
Previous by thread: Re: TransformedDistribution -- odd problem
Next by thread: Re: Problem with NIntegrating a compiled function