Mercurial > hg > octave-lyh

annotate liboctave/Quad.cc @ 8920:eb63fbe60fab

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update copyright notices
author John W. Eaton <jwe@octave.org>
date Sat, 07 Mar 2009 10:41:27 -0500
parents 4976f66d469b
children b1e82cc8a9f3
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
1 /*
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
2
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2003,
8920
eb63fbe60fab update copyright notices
John W. Eaton <jwe@octave.org>
parents: 7924
diff changeset
4 2004, 2005, 2006, 2007, 2008 John W. Eaton
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
5
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
6 This file is part of Octave.
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
7
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
8 Octave is free software; you can redistribute it and/or modify it
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
9 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5775
diff changeset
10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5775
diff changeset
11 option) any later version.
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
12
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
13 Octave is distributed in the hope that it will be useful, but WITHOUT
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
16 for more details.
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
17
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
18 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5775
diff changeset
19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5775
diff changeset
20 <http://www.gnu.org/licenses/>.
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
21
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
22 */
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
23
238
780cbbc57b7c [project @ 1993-11-30 20:23:04 by jwe]
jwe
parents: 3
diff changeset
24 #ifdef HAVE_CONFIG_H
1192
b6360f2d4fa6 [project @ 1995-03-30 21:38:35 by jwe]
jwe
parents: 1011
diff changeset
25 #include <config.h>
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
26 #endif
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
27
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
28 #include "Quad.h"
1847
2ffe49eb95a5 [project @ 1996-02-03 12:47:55 by jwe]
jwe
parents: 1546
diff changeset
29 #include "f77-fcn.h"
2292
a265f5e924e3 [project @ 1996-06-06 14:59:05 by jwe]
jwe
parents: 1993
diff changeset
30 #include "lo-error.h"
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
31 #include "quit.h"
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
32 #include "sun-utils.h"
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
33
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
34 static integrand_fcn user_fcn;
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
35 static float_integrand_fcn float_user_fcn;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
36
5775
ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe]
jwe
parents: 5307
diff changeset
37 // FIXME -- would be nice to not have to have this global
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
38 // variable.
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
39 // Nonzero means an error occurred in the calculation of the integrand
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
40 // function, and the user wants us to quit.
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
41 int quad_integration_error = 0;
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
42
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
43 typedef octave_idx_type (*quad_fcn_ptr) (double*, int&, double*);
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
44 typedef octave_idx_type (*quad_float_fcn_ptr) (float*, int&, float*);
3507
00fdd363c098 [project @ 2000-02-01 21:31:44 by jwe]
jwe
parents: 3136
diff changeset
45
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
46 extern "C"
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
47 {
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
48 F77_RET_T
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
49 F77_FUNC (dqagp, DQAGP) (quad_fcn_ptr, const double&, const double&,
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
50 const octave_idx_type&, const double*, const double&,
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
51 const double&, double&, double&, octave_idx_type&,
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
52 octave_idx_type&, const octave_idx_type&, const octave_idx_type&, octave_idx_type&, octave_idx_type*,
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
53 double*);
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
54
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
55 F77_RET_T
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
56 F77_FUNC (dqagi, DQAGI) (quad_fcn_ptr, const double&, const octave_idx_type&,
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
57 const double&, const double&, double&,
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
58 double&, octave_idx_type&, octave_idx_type&, const octave_idx_type&,
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
59 const octave_idx_type&, octave_idx_type&, octave_idx_type*, double*);
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
60
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
61 F77_RET_T
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
62 F77_FUNC (qagp, QAGP) (quad_float_fcn_ptr, const float&, const float&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
63 const octave_idx_type&, const float*, const float&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
64 const float&, float&, float&, octave_idx_type&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
65 octave_idx_type&, const octave_idx_type&, const octave_idx_type&, octave_idx_type&, octave_idx_type*,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
66 float*);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
67
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
68 F77_RET_T
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
69 F77_FUNC (qagi, QAGI) (quad_float_fcn_ptr, const float&, const octave_idx_type&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
70 const float&, const float&, float&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
71 float&, octave_idx_type&, octave_idx_type&, const octave_idx_type&,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
72 const octave_idx_type&, octave_idx_type&, octave_idx_type*, float*);
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
73 }
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
74
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
75 static octave_idx_type
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
76 user_function (double *x, int& ierr, double *result)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
77 {
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
78 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
79
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
80 #if defined (sun) && defined (__GNUC__)
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
81 double xx = access_double (x);
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
82 #else
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
83 double xx = *x;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
84 #endif
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
85
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
86 quad_integration_error = 0;
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
87
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
88 double xresult = (*user_fcn) (xx);
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
89
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
90 #if defined (sun) && defined (__GNUC__)
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
91 assign_double (result, xresult);
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
92 #else
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
93 *result = xresult;
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
94 #endif
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
95
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
96 if (quad_integration_error)
1251
97eac19837dc [project @ 1995-04-11 15:58:32 by jwe]
jwe
parents: 1192
diff changeset
97 ierr = -1;
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
98
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
99 END_INTERRUPT_WITH_EXCEPTIONS;
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
100
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
101 return 0;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
102 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
103
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
104 static octave_idx_type
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
105 float_user_function (float *x, int& ierr, float *result)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
106 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
107 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
108
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
109 #if defined (sun) && defined (__GNUC__)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
110 float xx = access_float (x);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
111 #else
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
112 float xx = *x;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
113 #endif
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
114
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
115 quad_integration_error = 0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
116
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
117 float xresult = (*float_user_fcn) (xx);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
118
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
119 #if defined (sun) && defined (__GNUC__)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
120 assign_float (result, xresult);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
121 #else
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
122 *result = xresult;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
123 #endif
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
124
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
125 if (quad_integration_error)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
126 ierr = -1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
127
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
128 END_INTERRUPT_WITH_EXCEPTIONS;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
129
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
130 return 0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
131 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
132
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
133 double
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
134 DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
135 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
136 octave_idx_type npts = singularities.capacity () + 2;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
137 double *points = singularities.fortran_vec ();
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
138 double result = 0.0;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
139
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
140 octave_idx_type leniw = 183*npts - 122;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
141 Array<octave_idx_type> iwork (leniw);
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
142 octave_idx_type *piwork = iwork.fortran_vec ();
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
143
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
144 octave_idx_type lenw = 2*leniw - npts;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
145 Array<double> work (lenw);
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
146 double *pwork = work.fortran_vec ();
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
147
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
148 user_fcn = f;
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
149 octave_idx_type last;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
150
289
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
151 double abs_tol = absolute_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
152 double rel_tol = relative_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
153
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
154 F77_XFCN (dqagp, DQAGP, (user_function, lower_limit, upper_limit,
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
155 npts, points, abs_tol, rel_tol, result,
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
156 abserr, neval, ier, leniw, lenw, last,
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
157 piwork, pwork));
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
158
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
159 return result;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
160 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
161
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
162 float
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
163 DefQuad::do_integrate (octave_idx_type&, octave_idx_type&, float&)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
164 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
165 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
166 return 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
167 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
168
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
169 double
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
170 IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
171 {
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
172 double result = 0.0;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
173
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
174 octave_idx_type leniw = 128;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
175 Array<octave_idx_type> iwork (leniw);
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
176 octave_idx_type *piwork = iwork.fortran_vec ();
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
177
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
178 octave_idx_type lenw = 8*leniw;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
179 Array<double> work (lenw);
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
180 double *pwork = work.fortran_vec ();
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
181
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
182 user_fcn = f;
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
183 octave_idx_type last;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
184
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
185 octave_idx_type inf;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
186 switch (type)
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
187 {
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
188 case bound_to_inf:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
189 inf = 1;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
190 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
191
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
192 case neg_inf_to_bound:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
193 inf = -1;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
194 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
195
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
196 case doubly_infinite:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
197 inf = 2;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
198 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
199
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
200 default:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
201 assert (0);
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
202 break;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
203 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
204
289
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
205 double abs_tol = absolute_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
206 double rel_tol = relative_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
207
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
208 F77_XFCN (dqagi, DQAGI, (user_function, bound, inf, abs_tol, rel_tol,
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
209 result, abserr, neval, ier, leniw, lenw,
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
210 last, piwork, pwork));
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
211
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
212 return result;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
213 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
214
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
215 float
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
216 IndefQuad::do_integrate (octave_idx_type&, octave_idx_type&, float&)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
217 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
218 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
219 return 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
220 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
221
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
222 double
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
223 FloatDefQuad::do_integrate (octave_idx_type&, octave_idx_type&, double&)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
224 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
225 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
226 return 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
227 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
228
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
229 float
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
230 FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
231 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
232 octave_idx_type npts = singularities.capacity () + 2;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
233 float *points = singularities.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
234 float result = 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
235
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
236 octave_idx_type leniw = 183*npts - 122;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
237 Array<octave_idx_type> iwork (leniw);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
238 octave_idx_type *piwork = iwork.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
239
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
240 octave_idx_type lenw = 2*leniw - npts;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
241 Array<float> work (lenw);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
242 float *pwork = work.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
243
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
244 float_user_fcn = ff;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
245 octave_idx_type last;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
246
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
247 float abs_tol = single_precision_absolute_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
248 float rel_tol = single_precision_relative_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
249
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
250 F77_XFCN (qagp, QAGP, (float_user_function, lower_limit, upper_limit,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
251 npts, points, abs_tol, rel_tol, result,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
252 abserr, neval, ier, leniw, lenw, last,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
253 piwork, pwork));
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
254
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
255 return result;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
256 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
257
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
258 double
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
259 FloatIndefQuad::do_integrate (octave_idx_type&, octave_idx_type&, double&)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
260 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
261 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
262 return 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
263 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
264
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
265 float
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
266 FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
267 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
268 float result = 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
269
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
270 octave_idx_type leniw = 128;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
271 Array<octave_idx_type> iwork (leniw);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
272 octave_idx_type *piwork = iwork.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
273
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
274 octave_idx_type lenw = 8*leniw;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
275 Array<float> work (lenw);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
276 float *pwork = work.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
277
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
278 float_user_fcn = ff;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
279 octave_idx_type last;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
280
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
281 octave_idx_type inf;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
282 switch (type)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
283 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
284 case bound_to_inf:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
285 inf = 1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
286 break;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
287
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
288 case neg_inf_to_bound:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
289 inf = -1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
290 break;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
291
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
292 case doubly_infinite:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
293 inf = 2;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
294 break;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
295
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
296 default:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
297 assert (0);
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
298 break;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
299 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
300
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
301 float abs_tol = single_precision_absolute_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
302 float rel_tol = single_precision_relative_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
303
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
304 F77_XFCN (qagi, QAGI, (float_user_function, bound, inf, abs_tol, rel_tol,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
305 result, abserr, neval, ier, leniw, lenw,
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
306 last, piwork, pwork));
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
307
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
308 return result;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
309 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
310
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
311 /*
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
312 ;;; Local Variables: ***
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
313 ;;; mode: C++ ***
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
314 ;;; End: ***
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
315 */