GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Quad.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1993-2013 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "Quad.h"
28 #include "f77-fcn.h"
29 #include "lo-error.h"
30 #include "quit.h"
31 #include "sun-utils.h"
32 
33 static integrand_fcn user_fcn;
34 static float_integrand_fcn float_user_fcn;
35 
36 // FIXME: would be nice to not have to have this global variable.
37 // Nonzero means an error occurred in the calculation of the integrand
38 // function, and the user wants us to quit.
39 int quad_integration_error = 0;
40 
41 typedef octave_idx_type (*quad_fcn_ptr) (double*, int&, double*);
42 typedef octave_idx_type (*quad_float_fcn_ptr) (float*, int&, float*);
43 
44 extern "C"
45 {
46  F77_RET_T
47  F77_FUNC (dqagp, DQAGP) (quad_fcn_ptr, const double&, const double&,
48  const octave_idx_type&, const double*,
49  const double&, const double&, double&,
50  double&, octave_idx_type&, octave_idx_type&,
51  const octave_idx_type&, const octave_idx_type&,
52  octave_idx_type&, octave_idx_type*, double*);
53 
54  F77_RET_T
55  F77_FUNC (dqagi, DQAGI) (quad_fcn_ptr, const double&,
56  const octave_idx_type&, const double&,
57  const double&, double&, double&,
58  octave_idx_type&, octave_idx_type&,
59  const octave_idx_type&, const octave_idx_type&,
60  octave_idx_type&, octave_idx_type*, double*);
61 
62  F77_RET_T
63  F77_FUNC (qagp, QAGP) (quad_float_fcn_ptr, const float&, const float&,
64  const octave_idx_type&, const float*, const float&,
65  const float&, float&, float&, octave_idx_type&,
66  octave_idx_type&, const octave_idx_type&,
67  const octave_idx_type&, octave_idx_type&,
68  octave_idx_type*, float*);
69 
70  F77_RET_T
71  F77_FUNC (qagi, QAGI) (quad_float_fcn_ptr, const float&,
72  const octave_idx_type&, const float&,
73  const float&, float&, float&, octave_idx_type&,
74  octave_idx_type&, const octave_idx_type&,
75  const octave_idx_type&, octave_idx_type&,
76  octave_idx_type*, float*);
77 }
78 
79 static octave_idx_type
80 user_function (double *x, int& ierr, double *result)
81 {
82  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
83 
84 #if defined (__sparc) && defined (__GNUC__)
85  double xx = access_double (x);
86 #else
87  double xx = *x;
88 #endif
89 
90  quad_integration_error = 0;
91 
92  double xresult = (*user_fcn) (xx);
93 
94 #if defined (__sparc) && defined (__GNUC__)
95  assign_double (result, xresult);
96 #else
97  *result = xresult;
98 #endif
99 
100  if (quad_integration_error)
101  ierr = -1;
102 
103  END_INTERRUPT_WITH_EXCEPTIONS;
104 
105  return 0;
106 }
107 
108 static octave_idx_type
109 float_user_function (float *x, int& ierr, float *result)
110 {
111  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
112 
113  quad_integration_error = 0;
114 
115  *result = (*float_user_fcn) (*x);
116 
117  if (quad_integration_error)
118  ierr = -1;
119 
120  END_INTERRUPT_WITH_EXCEPTIONS;
121 
122  return 0;
123 }
124 
125 double
126 DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
127  double& abserr)
128 {
129  octave_idx_type npts = singularities.capacity () + 2;
130  double *points = singularities.fortran_vec ();
131  double result = 0.0;
132 
133  octave_idx_type leniw = 183*npts - 122;
134  Array<octave_idx_type> iwork (dim_vector (leniw, 1));
135  octave_idx_type *piwork = iwork.fortran_vec ();
136 
137  octave_idx_type lenw = 2*leniw - npts;
138  Array<double> work (dim_vector (lenw, 1));
139  double *pwork = work.fortran_vec ();
140 
141  user_fcn = f;
142  octave_idx_type last;
143 
144  double abs_tol = absolute_tolerance ();
145  double rel_tol = relative_tolerance ();
146 
147  F77_XFCN (dqagp, DQAGP, (user_function, lower_limit, upper_limit,
148  npts, points, abs_tol, rel_tol, result,
149  abserr, neval, ier, leniw, lenw, last,
150  piwork, pwork));
151 
152  return result;
153 }
154 
155 float
156 DefQuad::do_integrate (octave_idx_type&, octave_idx_type&, float&)
157 {
158  (*current_liboctave_error_handler) ("incorrect integration function called");
159  return 0.0;
160 }
161 
162 double
163 IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
164  double& abserr)
165 {
166  double result = 0.0;
167 
168  octave_idx_type leniw = 128;
169  Array<octave_idx_type> iwork (dim_vector (leniw, 1));
170  octave_idx_type *piwork = iwork.fortran_vec ();
171 
172  octave_idx_type lenw = 8*leniw;
173  Array<double> work (dim_vector (lenw, 1));
174  double *pwork = work.fortran_vec ();
175 
176  user_fcn = f;
177  octave_idx_type last;
178 
179  octave_idx_type inf;
180  switch (type)
181  {
182  case bound_to_inf:
183  inf = 1;
184  break;
185 
186  case neg_inf_to_bound:
187  inf = -1;
188  break;
189 
190  case doubly_infinite:
191  inf = 2;
192  break;
193 
194  default:
195  assert (0);
196  break;
197  }
198 
199  double abs_tol = absolute_tolerance ();
200  double rel_tol = relative_tolerance ();
201 
202  F77_XFCN (dqagi, DQAGI, (user_function, bound, inf, abs_tol, rel_tol,
203  result, abserr, neval, ier, leniw, lenw,
204  last, piwork, pwork));
205 
206  return result;
207 }
208 
209 float
210 IndefQuad::do_integrate (octave_idx_type&, octave_idx_type&, float&)
211 {
212  (*current_liboctave_error_handler) ("incorrect integration function called");
213  return 0.0;
214 }
215 
216 double
217 FloatDefQuad::do_integrate (octave_idx_type&, octave_idx_type&, double&)
218 {
219  (*current_liboctave_error_handler) ("incorrect integration function called");
220  return 0.0;
221 }
222 
223 float
224 FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
225  float& abserr)
226 {
227  octave_idx_type npts = singularities.capacity () + 2;
228  float *points = singularities.fortran_vec ();
229  float result = 0.0;
230 
231  octave_idx_type leniw = 183*npts - 122;
232  Array<octave_idx_type> iwork (dim_vector (leniw, 1));
233  octave_idx_type *piwork = iwork.fortran_vec ();
234 
235  octave_idx_type lenw = 2*leniw - npts;
236  Array<float> work (dim_vector (lenw, 1));
237  float *pwork = work.fortran_vec ();
238 
239  float_user_fcn = ff;
240  octave_idx_type last;
241 
242  float abs_tol = single_precision_absolute_tolerance ();
243  float rel_tol = single_precision_relative_tolerance ();
244 
245  F77_XFCN (qagp, QAGP, (float_user_function, lower_limit, upper_limit,
246  npts, points, abs_tol, rel_tol, result,
247  abserr, neval, ier, leniw, lenw, last,
248  piwork, pwork));
249 
250  return result;
251 }
252 
253 double
254 FloatIndefQuad::do_integrate (octave_idx_type&, octave_idx_type&, double&)
255 {
256  (*current_liboctave_error_handler) ("incorrect integration function called");
257  return 0.0;
258 }
259 
260 float
261 FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
262  float& abserr)
263 {
264  float result = 0.0;
265 
266  octave_idx_type leniw = 128;
267  Array<octave_idx_type> iwork (dim_vector (leniw, 1));
268  octave_idx_type *piwork = iwork.fortran_vec ();
269 
270  octave_idx_type lenw = 8*leniw;
271  Array<float> work (dim_vector (lenw, 1));
272  float *pwork = work.fortran_vec ();
273 
274  float_user_fcn = ff;
275  octave_idx_type last;
276 
277  octave_idx_type inf;
278  switch (type)
279  {
280  case bound_to_inf:
281  inf = 1;
282  break;
283 
284  case neg_inf_to_bound:
285  inf = -1;
286  break;
287 
288  case doubly_infinite:
289  inf = 2;
290  break;
291 
292  default:
293  assert (0);
294  break;
295  }
296 
297  float abs_tol = single_precision_absolute_tolerance ();
298  float rel_tol = single_precision_relative_tolerance ();
299 
300  F77_XFCN (qagi, QAGI, (float_user_function, bound, inf, abs_tol, rel_tol,
301  result, abserr, neval, ier, leniw, lenw,
302  last, piwork, pwork));
303 
304  return result;
305 }
Generated on Mon Dec 30 2013 03:04:50 for GNU Octave by   doxygen 1.8.1.2