00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include "Quad.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034 #include "quit.h"
00035 #include "sun-utils.h"
00036
00037 static integrand_fcn user_fcn;
00038
00039
00040
00041
00042
00043 int quad_integration_error = 0;
00044
00045 typedef int (*quad_fcn_ptr) (double*, int&, double*);
00046
00047 extern "C"
00048 {
00049 F77_RET_T
00050 F77_FUNC (dqagp, DQAGP) (quad_fcn_ptr, const double&, const double&,
00051 const int&, const double*, const double&,
00052 const double&, double&, double&, int&,
00053 int&, const int&, const int&, int&, int*,
00054 double*);
00055
00056 F77_RET_T
00057 F77_FUNC (dqagi, DQAGI) (quad_fcn_ptr, const double&, const int&,
00058 const double&, const double&, double&,
00059 double&, int&, int&, const int&,
00060 const int&, int&, int*, double*);
00061 }
00062
00063 static int
00064 user_function (double *x, int& ierr, double *result)
00065 {
00066 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00067
00068 #if defined (sun) && defined (__GNUC__)
00069 double xx = access_double (x);
00070 #else
00071 double xx = *x;
00072 #endif
00073
00074 quad_integration_error = 0;
00075
00076 double xresult = (*user_fcn) (xx);
00077
00078 #if defined (sun) && defined (__GNUC__)
00079 assign_double (result, xresult);
00080 #else
00081 *result = xresult;
00082 #endif
00083
00084 if (quad_integration_error)
00085 ierr = -1;
00086
00087 END_INTERRUPT_WITH_EXCEPTIONS;
00088
00089 return 0;
00090 }
00091
00092 double
00093 DefQuad::do_integrate (int& ier, int& neval, double& abserr)
00094 {
00095 int npts = singularities.capacity () + 2;
00096 double *points = singularities.fortran_vec ();
00097 double result = 0.0;
00098
00099 int leniw = 183*npts - 122;
00100 Array<int> iwork (leniw);
00101 int *piwork = iwork.fortran_vec ();
00102
00103 int lenw = 2*leniw - npts;
00104 Array<double> work (lenw);
00105 double *pwork = work.fortran_vec ();
00106
00107 user_fcn = f;
00108 int last;
00109
00110 double abs_tol = absolute_tolerance ();
00111 double rel_tol = relative_tolerance ();
00112
00113 F77_XFCN (dqagp, DQAGP, (user_function, lower_limit, upper_limit,
00114 npts, points, abs_tol, rel_tol, result,
00115 abserr, neval, ier, leniw, lenw, last,
00116 piwork, pwork));
00117
00118 if (f77_exception_encountered)
00119 (*current_liboctave_error_handler) ("unrecoverable error in dqagp");
00120
00121 return result;
00122 }
00123
00124 double
00125 IndefQuad::do_integrate (int& ier, int& neval, double& abserr)
00126 {
00127 double result = 0.0;
00128
00129 int leniw = 128;
00130 Array<int> iwork (leniw);
00131 int *piwork = iwork.fortran_vec ();
00132
00133 int lenw = 8*leniw;
00134 Array<double> work (lenw);
00135 double *pwork = work.fortran_vec ();
00136
00137 user_fcn = f;
00138 int last;
00139
00140 int inf;
00141 switch (type)
00142 {
00143 case bound_to_inf:
00144 inf = 1;
00145 break;
00146
00147 case neg_inf_to_bound:
00148 inf = -1;
00149 break;
00150
00151 case doubly_infinite:
00152 inf = 2;
00153 break;
00154
00155 default:
00156 assert (0);
00157 break;
00158 }
00159
00160 double abs_tol = absolute_tolerance ();
00161 double rel_tol = relative_tolerance ();
00162
00163 F77_XFCN (dqagi, DQAGI, (user_function, bound, inf, abs_tol, rel_tol,
00164 result, abserr, neval, ier, leniw, lenw,
00165 last, piwork, pwork));
00166
00167 if (f77_exception_encountered)
00168 (*current_liboctave_error_handler) ("unrecoverable error in dqagi");
00169
00170 return result;
00171 }
00172
00173
00174
00175
00176
00177