# HG changeset patch # User David Bateman # Date 1211185358 -7200 # Node ID 62affb34e6487a40de790346e0aced86da958cc4 # Parent a0c550b22e61d86a6dbf6253208780e390f372c9 Make quad work with single precision diff -r a0c550b22e61 -r 62affb34e648 liboctave/ChangeLog --- a/liboctave/ChangeLog Wed May 21 19:25:08 2008 +0200 +++ b/liboctave/ChangeLog Mon May 19 10:22:38 2008 +0200 @@ -1,3 +1,16 @@ +2008-05-21 David Bateman + + * Quad-opts.in: Handle single precision tolerances. + * Quad.cc (float_user_fcn): New static variable. + (quad_float_fcn_ptr): New typedef. + (qagp, quagi): New QUADPACK decls. + (float_user_function): New function. + (DefQuad::do_integrate, IndefQuad::do_integrate): Float versions. + (FloatDefQuad::do_integrate, FloatIndefQuad::do_integrate): + New functions. + * Quad.h (class Quad): Handle float type. + (class FloatDefQuad, class FloatIndefQuad): New classes. + 2008-05-21 Jaroslav Hajek * fCMatrix.h (xgemm): Provide decl. diff -r a0c550b22e61 -r 62affb34e648 liboctave/Quad-opts.in --- a/liboctave/Quad-opts.in Wed May 21 19:25:08 2008 +0200 +++ b/liboctave/Quad-opts.in Mon May 19 10:22:38 2008 +0200 @@ -39,3 +39,26 @@ INIT_VALUE = "::sqrt (DBL_EPSILON)" SET_EXPR = "val" END_OPTION + +OPTION + NAME = "single precision absolute tolerance" + DOC_ITEM +Absolute tolerance for single precision; may be zero for pure relative +error test. + END_DOC_ITEM + TYPE = "float" + INIT_VALUE = "::sqrt (FLT_EPSILON)" + SET_EXPR = "val" +END_OPTION + +OPTION + NAME = "single precision relative tolerance" + DOC_ITEM +Nonnegative relative tolerance for single precision. If the absolute +tolerance is zero, the relative tolerance must be greater than or equal to +@code{max (50*eps, 0.5e-28)}. + END_DOC_ITEM + TYPE = "float" + INIT_VALUE = "::sqrt (FLT_EPSILON)" + SET_EXPR = "val" +END_OPTION diff -r a0c550b22e61 -r 62affb34e648 liboctave/Quad.cc --- a/liboctave/Quad.cc Wed May 21 19:25:08 2008 +0200 +++ b/liboctave/Quad.cc Mon May 19 10:22:38 2008 +0200 @@ -32,6 +32,7 @@ #include "sun-utils.h" static integrand_fcn user_fcn; +static float_integrand_fcn float_user_fcn; // FIXME -- would be nice to not have to have this global // variable. @@ -40,6 +41,7 @@ int quad_integration_error = 0; typedef octave_idx_type (*quad_fcn_ptr) (double*, int&, double*); +typedef octave_idx_type (*quad_float_fcn_ptr) (float*, int&, float*); extern "C" { @@ -55,6 +57,19 @@ const double&, const double&, double&, double&, octave_idx_type&, octave_idx_type&, const octave_idx_type&, const octave_idx_type&, octave_idx_type&, octave_idx_type*, double*); + + F77_RET_T + F77_FUNC (qagp, QAGP) (quad_float_fcn_ptr, const float&, const float&, + const octave_idx_type&, const float*, const float&, + const float&, float&, float&, octave_idx_type&, + octave_idx_type&, const octave_idx_type&, const octave_idx_type&, octave_idx_type&, octave_idx_type*, + float*); + + F77_RET_T + F77_FUNC (qagi, QAGI) (quad_float_fcn_ptr, const float&, const octave_idx_type&, + const float&, const float&, float&, + float&, octave_idx_type&, octave_idx_type&, const octave_idx_type&, + const octave_idx_type&, octave_idx_type&, octave_idx_type*, float*); } static octave_idx_type @@ -86,6 +101,35 @@ return 0; } +static octave_idx_type +float_user_function (float *x, int& ierr, float *result) +{ + BEGIN_INTERRUPT_WITH_EXCEPTIONS; + +#if defined (sun) && defined (__GNUC__) + float xx = access_float (x); +#else + float xx = *x; +#endif + + quad_integration_error = 0; + + float xresult = (*float_user_fcn) (xx); + +#if defined (sun) && defined (__GNUC__) + assign_float (result, xresult); +#else + *result = xresult; +#endif + + if (quad_integration_error) + ierr = -1; + + END_INTERRUPT_WITH_EXCEPTIONS; + + return 0; +} + double DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) { @@ -115,6 +159,13 @@ return result; } +float +DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) +{ + (*current_liboctave_error_handler) ("incorrect integration function called"); + return 0.0; +} + double IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) { @@ -161,6 +212,102 @@ return result; } +float +IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) +{ + (*current_liboctave_error_handler) ("incorrect integration function called"); + return 0.0; +} + +double +FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) +{ + (*current_liboctave_error_handler) ("incorrect integration function called"); + return 0.0; +} + +float +FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) +{ + octave_idx_type npts = singularities.capacity () + 2; + float *points = singularities.fortran_vec (); + float result = 0.0; + + octave_idx_type leniw = 183*npts - 122; + Array iwork (leniw); + octave_idx_type *piwork = iwork.fortran_vec (); + + octave_idx_type lenw = 2*leniw - npts; + Array work (lenw); + float *pwork = work.fortran_vec (); + + float_user_fcn = ff; + octave_idx_type last; + + float abs_tol = single_precision_absolute_tolerance (); + float rel_tol = single_precision_relative_tolerance (); + + F77_XFCN (qagp, QAGP, (float_user_function, lower_limit, upper_limit, + npts, points, abs_tol, rel_tol, result, + abserr, neval, ier, leniw, lenw, last, + piwork, pwork)); + + return result; +} + +double +FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) +{ + (*current_liboctave_error_handler) ("incorrect integration function called"); + return 0.0; +} + +float +FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) +{ + float result = 0.0; + + octave_idx_type leniw = 128; + Array iwork (leniw); + octave_idx_type *piwork = iwork.fortran_vec (); + + octave_idx_type lenw = 8*leniw; + Array work (lenw); + float *pwork = work.fortran_vec (); + + float_user_fcn = ff; + octave_idx_type last; + + octave_idx_type inf; + switch (type) + { + case bound_to_inf: + inf = 1; + break; + + case neg_inf_to_bound: + inf = -1; + break; + + case doubly_infinite: + inf = 2; + break; + + default: + assert (0); + break; + } + + float abs_tol = single_precision_absolute_tolerance (); + float rel_tol = single_precision_relative_tolerance (); + + F77_XFCN (qagi, QAGI, (float_user_function, bound, inf, abs_tol, rel_tol, + result, abserr, neval, ier, leniw, lenw, + last, piwork, pwork)); + + return result; +} + /* ;;; Local Variables: *** ;;; mode: C++ *** diff -r a0c550b22e61 -r 62affb34e648 liboctave/Quad.h --- a/liboctave/Quad.h Wed May 21 19:25:08 2008 +0200 +++ b/liboctave/Quad.h Mon May 19 10:22:38 2008 +0200 @@ -27,12 +27,14 @@ #include #include "dColVector.h" +#include "fColVector.h" #include "lo-math.h" #if !defined (octave_Quad_typedefs) #define octave_Quad_typedefs 1 typedef double (*integrand_fcn) (double x); +typedef float (*float_integrand_fcn) (float x); #endif @@ -53,6 +55,9 @@ Quad (integrand_fcn fcn) : Quad_options (), f (fcn) { } + Quad (float_integrand_fcn fcn) + : Quad_options (), ff (fcn) { } + virtual ~Quad (void) { } virtual double integrate (void) @@ -62,6 +67,13 @@ return do_integrate (ier, neval, abserr); } + virtual float float_integrate (void) + { + octave_idx_type ier, neval; + float abserr; + return do_integrate (ier, neval, abserr); + } + virtual double integrate (octave_idx_type& ier) { octave_idx_type neval; @@ -69,22 +81,43 @@ return do_integrate (ier, neval, abserr); } + virtual float float_integrate (octave_idx_type& ier) + { + octave_idx_type neval; + float abserr; + return do_integrate (ier, neval, abserr); + } + virtual double integrate (octave_idx_type& ier, octave_idx_type& neval) { double abserr; return do_integrate (ier, neval, abserr); } + virtual float float_integrate (octave_idx_type& ier, octave_idx_type& neval) + { + float abserr; + return do_integrate (ier, neval, abserr); + } + virtual double integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) { return do_integrate (ier, neval, abserr); } + virtual float float_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) + { + return do_integrate (ier, neval, abserr); + } + virtual double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) = 0; + virtual float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) = 0; + protected: integrand_fcn f; + float_integrand_fcn ff; }; class @@ -112,6 +145,8 @@ double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr); + float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr); + private: double lower_limit; @@ -138,6 +173,8 @@ double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr); + float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr); + private: double bound; @@ -145,6 +182,68 @@ int integration_error; }; +class +OCTAVE_API +FloatDefQuad : public Quad +{ + public: + + FloatDefQuad (float_integrand_fcn fcn) + : Quad (fcn), lower_limit (0.0), upper_limit (1.0), singularities () { } + + FloatDefQuad (float_integrand_fcn fcn, float ll, float ul) + : Quad (fcn), lower_limit (ll), upper_limit (ul), singularities () { } + + FloatDefQuad (float_integrand_fcn fcn, float ll, float ul, + const FloatColumnVector& sing) + : Quad (fcn), lower_limit (ll), upper_limit (ul), + singularities (sing) { } + + FloatDefQuad (float_integrand_fcn fcn, const FloatColumnVector& sing) + : Quad (fcn), lower_limit (0.0), upper_limit (1.0), + singularities (sing) { } + + ~FloatDefQuad (void) { } + + double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr); + + float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr); + + private: + + float lower_limit; + float upper_limit; + + FloatColumnVector singularities; +}; + +class +OCTAVE_API +FloatIndefQuad : public Quad +{ + public: + + enum IntegralType { bound_to_inf, neg_inf_to_bound, doubly_infinite }; + + FloatIndefQuad (float_integrand_fcn fcn) + : Quad (fcn), bound (0.0), type (bound_to_inf) { } + + FloatIndefQuad (float_integrand_fcn fcn, double b, IntegralType t) + : Quad (fcn), bound (b), type (t) { } + + ~FloatIndefQuad (void) { } + + double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr); + + float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr); + + private: + + float bound; + IntegralType type; + int integration_error; +}; + #endif /* diff -r a0c550b22e61 -r 62affb34e648 src/ChangeLog --- a/src/ChangeLog Wed May 21 19:25:08 2008 +0200 +++ b/src/ChangeLog Mon May 19 10:22:38 2008 +0200 @@ -1,3 +1,9 @@ +2008-05-21 David Bateman + + * DLD-FUNCTIONS/quad.cc (quad_float_user_function): New function. + (Fquad): Handle float type. + New tests. + 2008-05-21 Jaroslav Hajek * OPERATORS/op-fcm-fcm.cc (trans_mul, mul_trans, herm_mul, mul_herm): diff -r a0c550b22e61 -r 62affb34e648 src/DLD-FUNCTIONS/quad.cc --- a/src/DLD-FUNCTIONS/quad.cc Wed May 21 19:25:08 2008 +0200 +++ b/src/DLD-FUNCTIONS/quad.cc Mon May 19 10:22:38 2008 +0200 @@ -103,6 +103,51 @@ return retval; } +float +quad_float_user_function (float x) +{ + float retval = 0.0; + + octave_value_list args; + args(0) = x; + + if (quad_fcn) + { + octave_value_list tmp = quad_fcn->do_multi_index_op (1, args); + + if (error_state) + { + quad_integration_error = 1; // FIXME + gripe_user_supplied_eval ("quad"); + return retval; + } + + if (tmp.length () && tmp(0).is_defined ()) + { + if (! warned_imaginary && tmp(0).is_complex_type ()) + { + warning ("quad: ignoring imaginary part returned from user-supplied function"); + warned_imaginary = true; + } + + retval = tmp(0).float_value (); + + if (error_state) + { + quad_integration_error = 1; // FIXME + gripe_user_supplied_eval ("quad"); + } + } + else + { + quad_integration_error = 1; // FIXME + gripe_user_supplied_eval ("quad"); + } + } + + return retval; +} + #define QUAD_ABORT() \ do \ { \ @@ -203,110 +248,222 @@ if (! quad_fcn) QUAD_ABORT (); - double a = args(1).double_value (); - - if (error_state) - QUAD_ABORT1 ("expecting second argument to be a scalar"); - - double b = args(2).double_value (); - - if (error_state) - QUAD_ABORT1 ("expecting third argument to be a scalar"); - - int indefinite = 0; - IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; - double bound = 0.0; - if (xisinf (a) && xisinf (b)) - { - indefinite = 1; - indef_type = IndefQuad::doubly_infinite; - } - else if (xisinf (a)) + if (args(1).is_single_type () || args(2).is_single_type ()) { - indefinite = 1; - bound = b; - indef_type = IndefQuad::neg_inf_to_bound; - } - else if (xisinf (b)) - { - indefinite = 1; - bound = a; - indef_type = IndefQuad::bound_to_inf; - } + float a = args(1).float_value (); - octave_idx_type ier = 0; - octave_idx_type nfun = 0; - double abserr = 0.0; - double val = 0.0; - bool have_sing = false; - ColumnVector sing; - ColumnVector tol; + if (error_state) + QUAD_ABORT1 ("expecting second argument to be a scalar"); - switch (nargin) - { - case 5: - if (indefinite) - QUAD_ABORT1 ("singularities not allowed on infinite intervals"); - - have_sing = true; - - sing = ColumnVector (args(4).vector_value ()); + float b = args(2).float_value (); if (error_state) - QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); + QUAD_ABORT1 ("expecting third argument to be a scalar"); - case 4: - tol = ColumnVector (args(3).vector_value ()); + int indefinite = 0; + FloatIndefQuad::IntegralType indef_type = FloatIndefQuad::doubly_infinite; + float bound = 0.0; + if (xisinf (a) && xisinf (b)) + { + indefinite = 1; + indef_type = FloatIndefQuad::doubly_infinite; + } + else if (xisinf (a)) + { + indefinite = 1; + bound = b; + indef_type = FloatIndefQuad::neg_inf_to_bound; + } + else if (xisinf (b)) + { + indefinite = 1; + bound = a; + indef_type = FloatIndefQuad::bound_to_inf; + } - if (error_state) - QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); + octave_idx_type ier = 0; + octave_idx_type nfun = 0; + float abserr = 0.0; + float val = 0.0; + bool have_sing = false; + FloatColumnVector sing; + FloatColumnVector tol; + + switch (nargin) + { + case 5: + if (indefinite) + QUAD_ABORT1 ("singularities not allowed on infinite intervals"); + + have_sing = true; + + sing = FloatColumnVector (args(4).float_vector_value ()); - switch (tol.capacity ()) - { - case 2: - quad_opts.set_relative_tolerance (tol (1)); + if (error_state) + QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); + + case 4: + tol = FloatColumnVector (args(3).float_vector_value ()); + + if (error_state) + QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); + + switch (tol.capacity ()) + { + case 2: + quad_opts.set_single_precision_relative_tolerance (tol (1)); + + case 1: + quad_opts.set_single_precision_absolute_tolerance (tol (0)); + break; + + default: + QUAD_ABORT1 ("expecting tol to contain no more than two values"); + } - case 1: - quad_opts.set_absolute_tolerance (tol (0)); + case 3: + if (indefinite) + { + FloatIndefQuad iq (quad_float_user_function, bound, + indef_type); + iq.set_options (quad_opts); + val = iq.float_integrate (ier, nfun, abserr); + } + else + { + if (have_sing) + { + FloatDefQuad dq (quad_float_user_function, a, b, sing); + dq.set_options (quad_opts); + val = dq.float_integrate (ier, nfun, abserr); + } + else + { + FloatDefQuad dq (quad_float_user_function, a, b); + dq.set_options (quad_opts); + val = dq.float_integrate (ier, nfun, abserr); + } + } break; default: - QUAD_ABORT1 ("expecting tol to contain no more than two values"); + panic_impossible (); + break; + } + + retval(3) = abserr; + retval(2) = static_cast (nfun); + retval(1) = static_cast (ier); + retval(0) = val; + + } + else + { + double a = args(1).double_value (); + + if (error_state) + QUAD_ABORT1 ("expecting second argument to be a scalar"); + + double b = args(2).double_value (); + + if (error_state) + QUAD_ABORT1 ("expecting third argument to be a scalar"); + + int indefinite = 0; + IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite; + double bound = 0.0; + if (xisinf (a) && xisinf (b)) + { + indefinite = 1; + indef_type = IndefQuad::doubly_infinite; + } + else if (xisinf (a)) + { + indefinite = 1; + bound = b; + indef_type = IndefQuad::neg_inf_to_bound; + } + else if (xisinf (b)) + { + indefinite = 1; + bound = a; + indef_type = IndefQuad::bound_to_inf; } - case 3: - if (indefinite) + octave_idx_type ier = 0; + octave_idx_type nfun = 0; + double abserr = 0.0; + double val = 0.0; + bool have_sing = false; + ColumnVector sing; + ColumnVector tol; + + switch (nargin) { - IndefQuad iq (quad_user_function, bound, indef_type); - iq.set_options (quad_opts); - val = iq.integrate (ier, nfun, abserr); - } - else - { - if (have_sing) + case 5: + if (indefinite) + QUAD_ABORT1 ("singularities not allowed on infinite intervals"); + + have_sing = true; + + sing = ColumnVector (args(4).vector_value ()); + + if (error_state) + QUAD_ABORT1 ("expecting vector of singularities as fourth argument"); + + case 4: + tol = ColumnVector (args(3).vector_value ()); + + if (error_state) + QUAD_ABORT1 ("expecting vector of tolerances as fifth argument"); + + switch (tol.capacity ()) { - DefQuad dq (quad_user_function, a, b, sing); - dq.set_options (quad_opts); - val = dq.integrate (ier, nfun, abserr); + case 2: + quad_opts.set_relative_tolerance (tol (1)); + + case 1: + quad_opts.set_absolute_tolerance (tol (0)); + break; + + default: + QUAD_ABORT1 ("expecting tol to contain no more than two values"); + } + + case 3: + if (indefinite) + { + IndefQuad iq (quad_user_function, bound, indef_type); + iq.set_options (quad_opts); + val = iq.integrate (ier, nfun, abserr); } else { - DefQuad dq (quad_user_function, a, b); - dq.set_options (quad_opts); - val = dq.integrate (ier, nfun, abserr); + if (have_sing) + { + DefQuad dq (quad_user_function, a, b, sing); + dq.set_options (quad_opts); + val = dq.integrate (ier, nfun, abserr); + } + else + { + DefQuad dq (quad_user_function, a, b); + dq.set_options (quad_opts); + val = dq.integrate (ier, nfun, abserr); + } } - } - break; + break; - default: - panic_impossible (); - break; - } + default: + panic_impossible (); + break; + } - retval(3) = abserr; - retval(2) = static_cast (nfun); - retval(1) = static_cast (ier); - retval(0) = val; + retval(3) = abserr; + retval(2) = static_cast (nfun); + retval(1) = static_cast (ier); + retval(0) = val; + } if (fcn_name.length()) clear_function (fcn_name); @@ -327,12 +484,19 @@ %! [v, ier, nfun, err] = quad ("f", 0, 5); %! assert(ier == 0 && abs (v - 17.5) < sqrt (eps) && nfun > 0 && %! err < sqrt (eps)) +%!test +%! [v, ier, nfun, err] = quad ("f", single(0), single(5)); +%! assert(ier == 0 && abs (v - 17.5) < sqrt (eps ("single")) && nfun > 0 && +%! err < sqrt (eps ("single"))) %!function y = f (x) %! y = x .* sin (1 ./ x) .* sqrt (abs (1 - x)); %!test %! [v, ier, nfun, err] = quad ("f", 0.001, 3); %! assert((ier == 0 || ier == 1) && abs (v - 1.98194120273598) < sqrt (eps) && nfun > 0); +%!test +%! [v, ier, nfun, err] = quad ("f", single(0.001), single(3)); +%! assert((ier == 0 || ier == 1) && abs (v - 1.98194120273598) < sqrt (eps ("single")) && nfun > 0); %!error quad ();