Mercurial > octave

annotate liboctave/numeric/Quad.cc @ 22323:bac0d6f07a3e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Update copyright notices for 2016.
author John W. Eaton <jwe@octave.org>
date Wed, 17 Aug 2016 01:05:19 -0400
parents 59cadee1c74b
children 4caa7b28d183
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
22323
bac0d6f07a3e maint: Update copyright notices for 2016.
John W. Eaton <jwe@octave.org>
parents: 22133
diff changeset
3 Copyright (C) 1993-2016 John W. Eaton
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
4
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
6
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
7 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
8 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
9 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
10 option) any later version.
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
11
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
12 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
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
14 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
15 for more details.
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
16
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
17 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
18 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
19 <http://www.gnu.org/licenses/>.
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
20
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
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
23 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21202
diff changeset
24 # include "config.h"
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
25 #endif
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
26
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
27 #include "Quad.h"
1847
2ffe49eb95a5 [project @ 1996-02-03 12:47:55 by jwe]
jwe
parents: 1546
diff changeset
28 #include "f77-fcn.h"
2292
a265f5e924e3 [project @ 1996-06-06 14:59:05 by jwe]
jwe
parents: 1993
diff changeset
29 #include "lo-error.h"
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
30 #include "quit.h"
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
31 #include "sun-utils.h"
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
32
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
33 static integrand_fcn user_fcn;
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
34 static float_integrand_fcn float_user_fcn;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
35
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
36 // FIXME: would be nice to not have to have this global variable.
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
37 // 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
38 // function, and the user wants us to quit.
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
39 int quad_integration_error = 0;
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
40
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
41 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
42 typedef octave_idx_type (*quad_float_fcn_ptr) (float*, int&, float*);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
43
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
44 extern "C"
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
45 {
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
46 F77_RET_T
22133
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
47 F77_FUNC (dqagp, DQAGP) (quad_fcn_ptr, const F77_DBLE&, const F77_DBLE&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
48 const F77_INT&, const F77_DBLE*,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
49 const F77_DBLE&, const F77_DBLE&, F77_DBLE&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
50 F77_DBLE&, F77_INT&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
51 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
52 F77_INT&, F77_INT*, F77_DBLE*);
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
53
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
54 F77_RET_T
22133
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
55 F77_FUNC (dqagi, DQAGI) (quad_fcn_ptr, const F77_DBLE&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
56 const F77_INT&, const F77_DBLE&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
57 const F77_DBLE&, F77_DBLE&, F77_DBLE&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
58 F77_INT&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
59 const F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
60 F77_INT&, F77_INT*, F77_DBLE*);
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
61
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
62 F77_RET_T
22133
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
63 F77_FUNC (qagp, QAGP) (quad_float_fcn_ptr, const F77_REAL&, const F77_REAL&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
64 const F77_INT&, const F77_REAL*, const F77_REAL&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
65 const F77_REAL&, F77_REAL&, F77_REAL&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
66 F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
67 const F77_INT&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
68 F77_INT*, F77_REAL*);
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
69
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
70 F77_RET_T
22133
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
71 F77_FUNC (qagi, QAGI) (quad_float_fcn_ptr, const F77_REAL&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
72 const F77_INT&, const F77_REAL&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
73 const F77_REAL&, F77_REAL&, F77_REAL&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
74 F77_INT&, const F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
75 const F77_INT&, F77_INT&,
59cadee1c74b new macros for F77 data types
John W. Eaton <jwe@octave.org>
parents: 21724
diff changeset
76 F77_INT*, F77_REAL*);
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
77 }
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
78
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
79 static octave_idx_type
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
80 user_function (double *x, int& ierr, double *result)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
81 {
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
82 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
83
9179
5be2e6696772 use access_double and assign_double on sparc only
Carsten Clark <tantumquantum+gnuoctave@gmail.com>
parents: 9169
diff changeset
84 #if defined (__sparc) && defined (__GNUC__)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
85 double xx = access_double (x);
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
86 #else
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
87 double xx = *x;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
88 #endif
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
89
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
90 quad_integration_error = 0;
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
91
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
92 double xresult = (*user_fcn) (xx);
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
93
9179
5be2e6696772 use access_double and assign_double on sparc only
Carsten Clark <tantumquantum+gnuoctave@gmail.com>
parents: 9169
diff changeset
94 #if defined (__sparc) && defined (__GNUC__)
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
95 assign_double (result, xresult);
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
96 #else
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
97 *result = xresult;
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
98 #endif
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
99
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
100 if (quad_integration_error)
1251
97eac19837dc [project @ 1995-04-11 15:58:32 by jwe]
jwe
parents: 1192
diff changeset
101 ierr = -1;
260
74d73a4b3fc7 [project @ 1993-12-14 23:57:42 by jwe]
jwe
parents: 238
diff changeset
102
4180
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
103 END_INTERRUPT_WITH_EXCEPTIONS;
84fe3ca3a246 [project @ 2002-11-15 04:47:01 by jwe]
jwe
parents: 4066
diff changeset
104
3136
af7ec9d3a5e6 [project @ 1998-02-01 20:11:06 by jwe]
jwe
parents: 2847
diff changeset
105 return 0;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
106 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
107
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
108 static octave_idx_type
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
109 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
110 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
111 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
112
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
113 quad_integration_error = 0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
114
9169
b1e82cc8a9f3 eliminate broken special case for copying floats on Sun systems
Carsten Clark <tantumquantum+gnuoctave@gmail.com>
parents: 8920
diff changeset
115 *result = (*float_user_fcn) (*x);
7805
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 if (quad_integration_error)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
118 ierr = -1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
119
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
120 END_INTERRUPT_WITH_EXCEPTIONS;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
121
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
122 return 0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
123 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
124
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
125 double
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
126 DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
127 double& abserr)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
128 {
20229
5dfaaaae784f Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents: 19697
diff changeset
129 octave_idx_type npts = singularities.numel () + 2;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
130 double *points = singularities.fortran_vec ();
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
131 double result = 0.0;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
132
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
133 octave_idx_type leniw = 183*npts - 122;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
134 Array<octave_idx_type> iwork (dim_vector (leniw, 1));
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
135 octave_idx_type *piwork = iwork.fortran_vec ();
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
136
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
137 octave_idx_type lenw = 2*leniw - npts;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
138 Array<double> work (dim_vector (lenw, 1));
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
139 double *pwork = work.fortran_vec ();
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
140
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
141 user_fcn = f;
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
142 octave_idx_type last;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
143
289
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
144 double abs_tol = absolute_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
145 double rel_tol = relative_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
146
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
147 F77_XFCN (dqagp, DQAGP, (user_function, lower_limit, upper_limit,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
148 npts, points, abs_tol, rel_tol, result,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
149 abserr, neval, ier, leniw, lenw, last,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
150 piwork, pwork));
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
151
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
152 return result;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
153 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
154
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
155 float
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
156 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
157 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
158 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
159 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
160
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
161 double
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
162 IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
163 double& abserr)
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
164 {
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
165 double result = 0.0;
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
166
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
167 octave_idx_type leniw = 128;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
168 Array<octave_idx_type> iwork (dim_vector (leniw, 1));
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
169 octave_idx_type *piwork = iwork.fortran_vec ();
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
170
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
171 octave_idx_type lenw = 8*leniw;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
172 Array<double> work (dim_vector (lenw, 1));
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
173 double *pwork = work.fortran_vec ();
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
174
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
175 user_fcn = f;
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
176 octave_idx_type last;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
177
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
178 octave_idx_type inf;
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
179 switch (type)
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
180 {
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
181 case bound_to_inf:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
182 inf = 1;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
183 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
184
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
185 case neg_inf_to_bound:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
186 inf = -1;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
187 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
188
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
189 case doubly_infinite:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
190 inf = 2;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
191 break;
1360
7eb93d12654c [project @ 1995-09-05 21:51:54 by jwe]
jwe
parents: 1315
diff changeset
192
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
193 default:
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
194 assert (0);
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
195 break;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
196 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
197
289
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
198 double abs_tol = absolute_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
199 double rel_tol = relative_tolerance ();
c23f50e61c58 [project @ 1994-01-13 06:25:58 by jwe]
jwe
parents: 260
diff changeset
200
1935
23c350d0cf9d [project @ 1996-02-12 04:03:05 by jwe]
jwe
parents: 1882
diff changeset
201 F77_XFCN (dqagi, DQAGI, (user_function, bound, inf, abs_tol, rel_tol,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
202 result, abserr, neval, ier, leniw, lenw,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
203 last, piwork, pwork));
3
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
204
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
205 return result;
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
206 }
9a4c07481e61 [project @ 1993-08-08 01:20:23 by jwe]
jwe
parents:
diff changeset
207
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
208 float
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
209 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
210 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
211 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
212 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
213
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
214 double
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
215 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
216 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
217 (*current_liboctave_error_handler) ("incorrect integration function called");
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
218 }
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
219
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
220 float
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
221 FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
222 float& abserr)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
223 {
20229
5dfaaaae784f Deprecate Array::capacity() and Sparse::capacity() for numel() and nzmax().
Carnë Draug <carandraug@octave.org>
parents: 19697
diff changeset
224 octave_idx_type npts = singularities.numel () + 2;
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
225 float *points = singularities.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
226 float result = 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 octave_idx_type leniw = 183*npts - 122;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
229 Array<octave_idx_type> iwork (dim_vector (leniw, 1));
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
230 octave_idx_type *piwork = iwork.fortran_vec ();
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 lenw = 2*leniw - npts;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
233 Array<float> work (dim_vector (lenw, 1));
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
234 float *pwork = work.fortran_vec ();
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 float_user_fcn = ff;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
237 octave_idx_type last;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
238
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
239 float abs_tol = single_precision_absolute_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
240 float rel_tol = single_precision_relative_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
241
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
242 F77_XFCN (qagp, QAGP, (float_user_function, lower_limit, upper_limit,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
243 npts, points, abs_tol, rel_tol, result,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
244 abserr, neval, ier, leniw, lenw, last,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
245 piwork, pwork));
7805
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 return result;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
248 }
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 double
7924
4976f66d469b miscellaneous cleanup
John W. Eaton <jwe@octave.org>
parents: 7805
diff changeset
251 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
252 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
253 (*current_liboctave_error_handler) ("incorrect integration function called");
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
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
256 float
17769
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
257 FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval,
49a5a4be04a1 maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents: 17744
diff changeset
258 float& abserr)
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
259 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
260 float result = 0.0;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
261
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
262 octave_idx_type leniw = 128;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
263 Array<octave_idx_type> iwork (dim_vector (leniw, 1));
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
264 octave_idx_type *piwork = iwork.fortran_vec ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
265
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
266 octave_idx_type lenw = 8*leniw;
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
267 Array<float> work (dim_vector (lenw, 1));
7805
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
268 float *pwork = work.fortran_vec ();
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 float_user_fcn = ff;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
271 octave_idx_type last;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
272
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
273 octave_idx_type inf;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
274 switch (type)
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
275 {
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
276 case bound_to_inf:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
277 inf = 1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
278 break;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
279
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
280 case neg_inf_to_bound:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
281 inf = -1;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
282 break;
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 doubly_infinite:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
285 inf = 2;
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 default:
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
289 assert (0);
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
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
293 float abs_tol = single_precision_absolute_tolerance ();
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
294 float rel_tol = single_precision_relative_tolerance ();
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 F77_XFCN (qagi, QAGI, (float_user_function, bound, inf, abs_tol, rel_tol,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
297 result, abserr, neval, ier, leniw, lenw,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
298 last, piwork, pwork));
7805
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 return result;
62affb34e648 Make quad work with single precision
David Bateman <dbateman@free.fr>
parents: 7482
diff changeset
301 }