Source code: riso/numerical/qagse.java
1 package riso.numerical;
2
3 /** qagse.java, qags.java, qelg.java, qk21.java, and qpsrt.java are
4 * derivative works (translations) of Fortran code by Robert Piessens
5 * and Elise de Doncker. These five files are released under GPL
6 * by permission of Robert Piessens. In response to my question,
7 * <pre>
8 * >I would like to have your permission to distribute my
9 * >Java translation of your QUADPACK routines under the
10 * >terms of the GPL. Do I have your permission to do so?
11 * </pre>
12 * Robert Piessens writes:
13 * <pre>
14 * Date: Mon, 28 Jan 2002 14:41:58 +0100
15 * To: "Robert Dodier" <robert_dodier@yahoo.com>
16 * From: "Robert Piessens" <Robert.Piessens@cs.kuleuven.ac.be>
17 * Subject: Re: Permission to redistribute QUADPACK translation?
18 *
19 * OK, You have my permission.
20 *
21 * Robert Piessens
22 * </pre>
23 */
24 public class qagse implements java.io.Serializable
25 {
26 public static void do_qagse ( Callback_1d f , double a , double b , double epsabs , double epsrel , int limit , double[] result , double[] abserr , int[] neval , int[] ier , double[] alist , double[] blist , double[] rlist , double[] elist , int[] iord , int[] last ) throws Exception // SHOULD USE ier EXCLUSIVELY OR EXCEPTIONS EXCLUSIVELY, NOT BOTH !!!
27 {
28 double area, area12, a1, a2, b1, b2, correc = -999, dres, epmach, erlarg = -999, erlast, errbnd, erro12, errsum, ertest = -999, oflow, small = -999, uflow;
29 double[] res3la = new double [ 3 ], rlist2 = new double [ 52 ];
30 double[] defabs = new double[1], resabs = new double[1];
31 double[] area1 = new double[1], area2 = new double[1];
32 double[] error1 = new double[1], error2 = new double[1];
33 double[] defab1 = new double[1], defab2 = new double[1];
34 double[] reseps = new double[1], abseps = new double[1];
35 double[] errmax = new double[1];
36 int id, ierro, iroff1, iroff2, iroff3, jupbnd, k, ksgn, ktmin;
37 int[] maxerr = new int[1], nrmax = new int[1], numrl2 = new int[1], nres = new int[1];
38 boolean extrap,noext;
39 epmach = qk21.D1MACH [ 4-1 ];
40 ier[0] = 0;
41 neval[0] = 0;
42 last[0] = 0;
43 result[0] = 0;
44 abserr[0] = 0;
45 alist [ 1 -1 ] = a;
46 blist [ 1 -1 ] = b;
47 rlist [ 1 -1 ] = 0;
48 elist [ 1 -1 ] = 0;
49 if ( epsabs <= 0 && epsrel < Math.max ( 50 * epmach , 0.5e-28 ) ) ier[0] = 6;
50 if ( ier[0] == 6 ) return;
51 uflow = qk21.D1MACH [ 1-1 ];
52 oflow = qk21.D1MACH [ 2-1 ];
53 ierro = 0;
54 qk21.do_qk21 ( f , a , b , result , abserr , defabs , resabs );
55 dres = Math.abs ( result[0] );
56 errbnd = Math.max ( epsabs , epsrel * dres );
57 last[0] = 1;
58 rlist [ 1 -1 ] = result[0];
59 elist [ 1 -1 ] = abserr[0];
60 iord [ 1 -1 ] = 1;
61 if ( abserr[0] <= 100 * epmach * defabs[0] && abserr[0] > errbnd ) ier[0] = 2;
62 if ( limit == 1 ) ier[0] = 1;
63 if ( ier[0] != 0 || ( abserr[0] <= errbnd && abserr[0] != resabs[0] ) || abserr[0] == 0 )
64 {
65 neval[0] = 42*last[0]-21;
66 // System.err.println( "do_qagse: return after 1st qk21; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
67 return;
68 }
69 rlist2 [ 1 -1 ] = result[0];
70 errmax[0] = abserr[0];
71 maxerr[0] = 1;
72 area = result[0];
73 errsum = abserr[0];
74 abserr[0] = oflow;
75 nrmax[0] = 1;
76 nres[0] = 0;
77 numrl2[0] = 2;
78 ktmin = 0;
79 extrap = false;
80 noext = false;
81 iroff1 = 0;
82 iroff2 = 0;
83 iroff3 = 0;
84 ksgn = -1;
85 if ( dres >= ( 1 - 50 * epmach ) * defabs[0] ) ksgn = 1;
86 for ( last[0] = 2 ; last[0] <= limit ; last[0] += 1 )
87 {
88 a1 = alist [ maxerr[0] -1 ];
89 b1 = 0.5 * ( alist [ maxerr[0] -1 ] + blist [ maxerr[0] -1 ] );
90 a2 = b1;
91 b2 = blist [ maxerr[0] -1 ];
92 erlast = errmax[0];
93 qk21.do_qk21 ( f , a1 , b1 , area1 , error1 , resabs , defab1 );
94 qk21.do_qk21 ( f , a2 , b2 , area2 , error2 , resabs , defab2 );
95 area12 = area1[0]+area2[0];
96 erro12 = error1[0]+error2[0];
97 errsum = errsum+erro12-errmax[0];
98 area = area + area12 - rlist [ maxerr[0] -1 ];
99 if ( ! ( defab1[0] == error1[0] || defab2[0] == error2[0] ) )
100 {
101 if ( ! ( Math.abs ( rlist [ maxerr[0] -1 ] - area12 ) > 1e-5 * Math.abs ( area12 ) || erro12 < 0.99 * errmax[0] ) )
102 {
103 if ( extrap ) iroff2 = iroff2 + 1;
104 if ( ! extrap ) iroff1 = iroff1 + 1;
105 }
106 if ( last[0] > 10 && erro12 > errmax[0] ) iroff3 = iroff3 + 1;
107 }
108 rlist [ maxerr[0] -1 ] = area1[0];
109 rlist [ last[0] -1 ] = area2[0];
110 errbnd = Math.max ( epsabs , epsrel * Math.abs ( area ) );
111 if ( iroff1 + iroff2 >= 10 || iroff3 >= 20 ) ier[0] = 2;
112 if ( iroff2 >= 5 ) ierro = 3;
113 if ( last[0] == limit ) ier[0] = 1;
114 if ( Math.max ( Math.abs ( a1 ) , Math.abs ( b2 ) ) <= ( 1 + 100 * epmach ) * ( Math.abs ( a2 ) + 1000 * uflow ) ) ier[0] = 4;
115 if ( ! ( error2[0] > error1[0] ) )
116 {
117 alist [ last[0] -1 ] = a2;
118 blist [ maxerr[0] -1 ] = b1;
119 blist [ last[0] -1 ] = b2;
120 elist [ maxerr[0] -1 ] = error1[0];
121 elist [ last[0] -1 ] = error2[0];
122 }
123 else
124 {
125 alist [ maxerr[0] -1 ] = a2;
126 alist [ last[0] -1 ] = a1;
127 blist [ last[0] -1 ] = b1;
128 rlist [ maxerr[0] -1 ] = area2[0];
129 rlist [ last[0] -1 ] = area1[0];
130 elist [ maxerr[0] -1 ] = error2[0];
131 elist [ last[0] -1 ] = error1[0];
132 }
133 qpsrt.do_qpsrt ( limit , last[0] , maxerr , errmax , elist , iord , nrmax );
134 if ( errsum <= errbnd )
135 {
136 result[0] = 0;
137 for ( k = 1 ; k <= last[0] ; k += 1 )
138 {
139 result[0] = result[0] + rlist [ k -1 ];
140 }
141 abserr[0] = errsum;
142 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
143 neval[0] = 42*last[0]-21;
144 // System.err.println( "do_qagse: return thru errsum <= errbnd; errsum: "+errsum+" errbnd: "+errbnd+" result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
145 return;
146 }
147 // System.err.println( "reached ier != 0 test; ier: "+ier[0] );
148 if ( ier[0] != 0 ) break;
149 // System.err.println( "reached last != 2 test; last: "+last[0] );
150 if ( ! ( last[0] == 2 ) )
151 {
152 // System.err.println( "reach noext test; noext: "+noext );
153 if ( noext ) continue;
154 if ( erlarg == -999 ) throw new Exception( "do_qagse: erlarg NOT DEFINED!!!" );
155 erlarg = erlarg-erlast;
156 if ( small == -999 ) throw new Exception( "do_qagse: small NOT DEFINED!!!" );
157 if ( Math.abs ( b1 - a1 ) > small ) erlarg = erlarg + erro12;
158 if ( ! ( extrap ) )
159 {
160 // System.err.println( "reach small test 1; blist[maxerr]: "+blist [ maxerr[0] -1 ]+" alist[maxerr]: "+alist [ maxerr[0] -1 ]+" small: "+small );
161 if ( Math.abs ( blist [ maxerr[0] -1 ] - alist [ maxerr[0] -1 ] ) > small ) continue;
162 // System.err.println( "passed small test 1" );
163 extrap = true;
164 nrmax[0] = 2;
165 }
166 if ( ertest == -999 ) throw new Exception( "do_qagse: ertest NOT DEFINED!!!" );
167 // System.err.println( "reached ierro test; ierro: "+ierro+" erlarg: "+erlarg+" ertest: "+ertest );
168 if ( ! ( ierro == 3 || erlarg <= ertest ) )
169 {
170 id = nrmax[0];
171 jupbnd = last[0];
172 if ( last[0] > ( 2 + limit / 2 ) ) jupbnd = limit + 3 - last[0];
173 boolean goto90 = false;
174 for ( k = id ; k <= jupbnd ; k += 1 )
175 {
176 maxerr[0] = iord [ nrmax[0] -1 ];
177 errmax[0] = elist [ maxerr[0] -1 ];
178 // System.err.println( "reach small test 2; small: "+small );
179 if ( Math.abs ( blist [ maxerr[0] -1 ] - alist [ maxerr[0] -1 ] ) > small )
180 {
181 // System.err.println( "set goto90 = true" );
182 goto90 = true;
183 break;
184 }
185 nrmax[0] = nrmax[0]+1;
186 }
187 if ( goto90 ) continue;
188 }
189 numrl2[0] = numrl2[0]+1;
190 // System.err.println( "do_qagse: assign area == "+area+" to rlist2["+(numrl2[0]-1)+"]" );
191 rlist2 [ numrl2[0] -1 ] = area;
192 qelg.do_qelg( numrl2 , rlist2 , reseps , abseps , res3la , nres );
193 ktmin = ktmin+1;
194 if ( ktmin > 5 && abserr[0] < 1e-3 * errsum ) ier[0] = 5;
195 // System.err.println( "reached abseps < abserr test; abseps: "+abseps[0]+", abserr: "+abserr[0] );
196 if ( ! ( abseps[0] >= abserr[0] ) )
197 {
198 ktmin = 0;
199 abserr[0] = abseps[0];
200 // System.err.println( "assign reseps ("+reseps[0]+") to result" );
201 result[0] = reseps[0];
202 correc = erlarg;
203 ertest = Math.max ( epsabs , epsrel * Math.abs ( reseps[0] ) );
204 if ( abserr[0] <= ertest ) break;
205 }
206 if ( numrl2[0] == 1 ) noext = true;
207 if ( ier[0] == 5 ) break;
208 maxerr[0] = iord [ 1 -1 ];
209 errmax[0] = elist [ maxerr[0] -1 ];
210 nrmax[0] = 1;
211 extrap = false;
212 small = small*0.5;
213 erlarg = errsum;
214 continue;
215 }
216 small = Math.abs ( b - a ) * 0.375;
217 erlarg = errsum;
218 ertest = errbnd;
219 // System.err.println( "do_qagse: assign (#2) area == "+area+" to rlist2[1]" );
220 rlist2 [ 2 -1 ] = area;
221 }
222 if ( abserr[0] == oflow )
223 {
224 // System.err.println( "do_qagse: abserr: "+abserr[0]+" oflow: "+oflow );
225 result[0] = 0;
226 for ( k = 1 ; k <= last[0] ; k += 1 )
227 {
228 result[0] = result[0] + rlist [ k -1 ];
229 }
230 abserr[0] = errsum;
231 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
232 neval[0] = 42*last[0]-21;
233 // System.err.println( "do_qagse: return thru abserr == oflow; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
234 return;
235 }
236 if ( ier[0] + ierro == 0 )
237 {
238 if ( correc == -999 ) throw new Exception( "do_qagse: correc NOT DEFINED!!!" );
239 if ( ierro == 3 ) abserr[0] = abserr[0] + correc;
240 if ( ier[0] == 0 ) ier[0] = 3;
241 if ( result[0] != 0 && area != 0 )
242 {
243 if ( abserr[0] / Math.abs ( result[0] ) > errsum / Math.abs ( area ) )
244 {
245 result[0] = 0;
246 for ( k = 1 ; k <= last[0] ; k += 1 )
247 {
248 result[0] = result[0] + rlist [ k -1 ];
249 }
250 abserr[0] = errsum;
251 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
252 neval[0] = 42*last[0]-21;
253 // System.err.println( "do_qagse: return thru rel. err. ineq.; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
254 return;
255 }
256 if ( ksgn == ( - 1 ) && Math.max ( Math.abs ( result[0] ) , Math.abs ( area ) ) <= defabs[0] * 0.01 )
257 {
258 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
259 neval[0] = 42*last[0]-21;
260 // System.err.println( "do_qagse: return thru 2nd rel. err. ineq.; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
261 return;
262 }
263 if ( 0.01 > ( result[0] / area ) || ( result[0] / area ) > 100 || errsum > Math.abs ( area ) ) ier[0] = 6;
264 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
265 neval[0] = 42*last[0]-21;
266 // System.err.println( "do_qagse: didn't return thru other ways; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
267 return;
268 }
269 if ( abserr[0] > errsum )
270 {
271 result[0] = 0;
272 for ( k = 1 ; k <= last[0] ; k += 1 )
273 {
274 result[0] = result[0] + rlist [ k -1 ];
275 }
276 abserr[0] = errsum;
277 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
278 neval[0] = 42*last[0]-21;
279 // System.err.println( "do_qagse: return thru abserr > errsum; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
280 return;
281 }
282 if ( area == 0 )
283 {
284 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
285 neval[0] = 42*last[0]-21;
286 // System.err.println( "do_qagse: return thru area == 0; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
287 return;
288 }
289 }
290 if ( ksgn == ( - 1 ) && Math.max ( Math.abs ( result[0] ) , Math.abs ( area ) ) <= defabs[0] * 0.01 )
291 {
292 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
293 neval[0] = 42*last[0]-21;
294 // System.err.println( "do_qagse: return thru ineq 3; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
295 return;
296 }
297 if ( 0.01 > ( result[0] / area ) || ( result[0] / area ) > 100 || errsum > Math.abs ( area ) ) ier[0] = 6;
298 if ( ier[0] > 2 ) ier[0] = ier[0] - 1;
299 neval[0] = 42*last[0]-21;
300 // System.err.println( "do_qagse: didn't return other ways #2; result: "+result[0]+" abserr: "+abserr[0]+" ier: "+ier[0] );
301 return;
302 }
303 }