home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 9 / FreshFishVol9-CD2.bin / bbs / gnu / octave-1.1.1-src.lha / octave-1.1.1 / liboctave / Quad.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1995-01-04  |  6.4 KB  |  330 lines
  1. // Quad.cc                                                -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include <math.h>
  29. #include <float.h>
  30.  
  31. #include "Quad.h"
  32. #include "f77-uscore.h"
  33. #include "sun-utils.h"
  34.  
  35. static integrand_fcn user_fcn;
  36.  
  37. // XXX FIXME XXX -- would be nice to not have to have this global
  38. // variable.
  39. // Nonzero means an error occurred in the calculation of the integrand
  40. // function, and the user wants us to quit.
  41. int quad_integration_error = 0;
  42.  
  43. extern "C"
  44. {
  45.   int F77_FCN (dqagp) (const double (*)(double*, int*), const double*,
  46.                const double*, const int*, const double*,
  47.                const double*, const double*, double*, double*,
  48.                int*, int*, const int*, const int*, int*, int*,
  49.                double*);
  50.  
  51.   int F77_FCN (dqagi) (const double (*)(double*, int*), const double*,
  52.                const int*, const double*, const double*,
  53.                double*, double*, int*, int*, const int*,
  54.                const int*, int*, int*, double*);
  55. }
  56.  
  57. Quad::Quad (integrand_fcn fcn)
  58. {
  59.   f = fcn;
  60. }
  61.  
  62. Quad::Quad (integrand_fcn fcn, double abs, double rel)
  63. {
  64.   f = fcn;
  65. }
  66.  
  67. double
  68. Quad::integrate (void)
  69. {
  70.   int ier, neval;
  71.   double abserr;
  72.   return integrate (ier, neval, abserr);
  73. }
  74.  
  75. double
  76. Quad::integrate (int& ier)
  77. {
  78.   int neval;
  79.   double abserr;
  80.   return integrate (ier, neval, abserr);
  81. }
  82.  
  83. double
  84. Quad::integrate (int& ier, int& neval)
  85. {
  86.   double abserr;
  87.   return integrate (ier, neval, abserr);
  88. }
  89.  
  90. static double
  91. user_function (double *x, int *ierr)
  92. {
  93. #if defined (sun) && defined (__GNUC__)
  94.   double xx = access_double (x);
  95. #else
  96.   double xx = *x;
  97. #endif
  98.  
  99.   quad_integration_error = 0;
  100.  
  101.   double retval = (*user_fcn) (xx);
  102.  
  103.   if (quad_integration_error)
  104.     *ierr = -1;
  105.  
  106.   return retval;
  107. }
  108.  
  109. DefQuad::DefQuad (integrand_fcn fcn) : Quad (fcn)
  110. {
  111.   lower_limit = 0.0;
  112.   upper_limit = 1.0;
  113. }
  114.  
  115. DefQuad::DefQuad (integrand_fcn fcn, double ll, double ul)
  116.   : Quad (fcn) 
  117. {
  118.   lower_limit = ll;
  119.   upper_limit = ul;
  120. }
  121.  
  122. DefQuad::DefQuad (integrand_fcn fcn, double ll, double ul,
  123.           double abs, double rel) : Quad (fcn, abs, rel)
  124. {
  125.   lower_limit = ll;
  126.   upper_limit = ul;
  127. }
  128.  
  129. DefQuad::DefQuad (integrand_fcn fcn, double ll, double ul,
  130.           const Vector& sing) : Quad (fcn)
  131. {
  132.   lower_limit = ll;
  133.   upper_limit = ul;
  134.   singularities = sing;
  135. }
  136.  
  137. DefQuad::DefQuad (integrand_fcn fcn, const Vector& sing,
  138.           double abs, double rel) : Quad (fcn, abs, rel)
  139. {
  140.   lower_limit = 0.0;
  141.   upper_limit = 1.0;
  142.   singularities = sing;
  143. }
  144.  
  145. DefQuad::DefQuad (integrand_fcn fcn, const Vector& sing)
  146.   : Quad (fcn)
  147. {
  148.   lower_limit = 0.0;
  149.   upper_limit = 1.0;
  150.   singularities = sing;
  151. }
  152.  
  153. DefQuad::DefQuad (integrand_fcn fcn, double ll, double ul,
  154.           const Vector& sing, double abs, double rel)
  155.   : Quad (fcn, abs, rel)
  156. {
  157.   lower_limit = ll;
  158.   upper_limit = ul;
  159.   singularities = sing;
  160. }
  161.  
  162. double
  163. DefQuad::integrate (int& ier, int& neval, double& abserr)
  164. {
  165.   int npts = singularities.capacity () + 2;
  166.   double *points = singularities.fortran_vec ();
  167.   double result = 0.0;
  168.   int leniw = 183*npts - 122;
  169.   int lenw = 2*leniw - npts;
  170.   int *iwork = new int [leniw];
  171.   double *work = new double [lenw];
  172.   user_fcn = f;
  173.   int last;
  174.  
  175.   double abs_tol = absolute_tolerance ();
  176.   double rel_tol = relative_tolerance ();
  177.  
  178.   F77_FCN (dqagp) (user_function, &lower_limit, &upper_limit, &npts,
  179.            points, &abs_tol, &rel_tol, &result, &abserr,
  180.            &neval, &ier, &leniw, &lenw, &last, iwork, work);
  181.  
  182.   delete [] iwork;
  183.   delete [] work;
  184.  
  185.   return result;
  186. }
  187.  
  188. IndefQuad::IndefQuad (integrand_fcn fcn) : Quad (fcn)
  189. {
  190.   bound = 0.0;
  191.   type = bound_to_inf;
  192. }
  193.  
  194. IndefQuad::IndefQuad (integrand_fcn fcn, double b, IntegralType t)
  195.   : Quad (fcn)
  196. {
  197.   bound = b;
  198.   type = t;
  199. }
  200.  
  201. IndefQuad::IndefQuad (integrand_fcn fcn, double b, IntegralType t,
  202.               double abs, double rel) : Quad (fcn, abs, rel)
  203. {
  204.   bound = b;
  205.   type = t;
  206. }
  207.  
  208. IndefQuad::IndefQuad (integrand_fcn fcn, double abs, double rel)
  209.   : Quad (fcn, abs, rel)
  210. {
  211.   bound = 0.0;
  212.   type = bound_to_inf;
  213. }
  214.  
  215. double
  216. IndefQuad::integrate (int& ier, int& neval, double& abserr)
  217. {
  218.   double result = 0.0;
  219.   int leniw = 128;
  220.   int lenw = 8*leniw;
  221.   int *iwork = new int [leniw];
  222.   double *work = new double [lenw];
  223.   user_fcn = f;
  224.   int last;
  225.  
  226.   int inf;
  227.   switch (type)
  228.     {
  229.     case bound_to_inf:
  230.       inf = 1;
  231.       break;
  232.     case neg_inf_to_bound:
  233.       inf = -1;
  234.       break;
  235.     case doubly_infinite:
  236.       inf = 2;
  237.       break;
  238.     default:
  239.       assert (0);
  240.       break;
  241.     }
  242.  
  243.   double abs_tol = absolute_tolerance ();
  244.   double rel_tol = relative_tolerance ();
  245.  
  246.   F77_FCN (dqagi) (user_function, &bound, &inf, &abs_tol, &rel_tol,
  247.            &result, &abserr, &neval, &ier, &leniw, &lenw,
  248.            &last, iwork, work);
  249.  
  250.   delete [] iwork;
  251.   delete [] work;
  252.  
  253.   return result;
  254. }
  255.  
  256. Quad_options::Quad_options (void)
  257. {
  258.   init ();
  259. }
  260.  
  261. Quad_options::Quad_options (const Quad_options& opt)
  262. {
  263.   copy (opt);
  264. }
  265.  
  266. Quad_options&
  267. Quad_options::operator = (const Quad_options& opt)
  268. {
  269.   if (this != &opt)
  270.     copy (opt);
  271.  
  272.   return *this;
  273. }
  274.  
  275. Quad_options::~Quad_options (void)
  276. {
  277. }
  278.  
  279. void
  280. Quad_options::init (void)
  281. {
  282.   double sqrt_eps = sqrt (DBL_EPSILON);
  283.   x_absolute_tolerance = sqrt_eps;
  284.   x_relative_tolerance = sqrt_eps;
  285. }
  286.  
  287. void
  288. Quad_options::copy (const Quad_options& opt)
  289. {
  290.   x_absolute_tolerance = opt.x_absolute_tolerance;
  291.   x_relative_tolerance = opt.x_relative_tolerance;
  292. }
  293.  
  294. void
  295. Quad_options::set_default_options (void)
  296. {
  297.   init ();
  298. }
  299.  
  300. void
  301. Quad_options::set_absolute_tolerance (double val)
  302. {
  303.   x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  304. }
  305.  
  306. void
  307. Quad_options::set_relative_tolerance (double val)
  308. {
  309.   x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON);
  310. }
  311.  
  312. double
  313. Quad_options::absolute_tolerance (void)
  314. {
  315.   return x_absolute_tolerance;
  316. }
  317.  
  318. double
  319. Quad_options::relative_tolerance (void)
  320. {
  321.   return x_relative_tolerance;
  322. }
  323.  
  324. /*
  325. ;;; Local Variables: ***
  326. ;;; mode: C++ ***
  327. ;;; page-delimiter: "^/\\*" ***
  328. ;;; End: ***
  329. */
  330.