Menu
▾
▴
[Maxima-commits] Maxima, A Computer Algebra System branch master updated. cd72ad4d4b9d1602014f1bc48f9f35e807b4e808
|
From: Raymond T. <rt...@us...> - 2011-10-16 05:38:57
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "Maxima, A Computer Algebra System".
The branch, master has been updated
via cd72ad4d4b9d1602014f1bc48f9f35e807b4e808 (commit)
via 82a6070c0b73c772111ebd4829a3b4e14985f7d0 (commit)
via ffbfbd964574a6e725534f2518fb77f08c489ec4 (commit)
via 4e6576139f3f769f7249f6d75ad44a4e488faafb (commit)
via 9067349233cc217786983b07c23a2b78cf7d27bc (commit)
from ce45cc56bdcf6a9fcfeb24db46ddd8bd2553144a (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit cd72ad4d4b9d1602014f1bc48f9f35e807b4e808
Merge: 82a6070 ce45cc5
Author: Raymond Toy <toy...@gm...>
Date: Sat Oct 15 22:37:59 2011 -0700
Merge branch 'master' of ssh://maxima.git.sourceforge.net/gitroot/maxima/maxima
commit 82a6070c0b73c772111ebd4829a3b4e14985f7d0
Author: Raymond Toy <toy...@gm...>
Date: Sat Oct 15 22:31:09 2011 -0700
Document quad_qagp.
diff --git a/doc/info/Integration.texi b/doc/info/Integration.texi
index 2618efc..48c0ed8 100644
--- a/doc/info/Integration.texi
+++ b/doc/info/Integration.texi
@@ -943,8 +943,14 @@ interval @math{(a, b)} and specified @math{c}.
The strategy is globally adaptive, and modified
Clenshaw-Curtis integration is used on the subranges
which contain the point @math{x = c}.
+
+@item quad_qagp
+Basically the same as @code{quad_qags} but points of singularity or
+discontinuity of the integrand must be supplied. This makes it easier
+for the integrator to produce a good solution.
@end table
+
@opencatbox
@category{Integral calculus} @category{Numerical methods} @category{Share packages} @category{Package quadpack}
@closecatbox
@@ -1748,3 +1754,103 @@ pos;
@closecatbox
@end deffn
+@c THERE ARE OPTIONAL ARGUMENTS WHICH MAKES LISTING THE VARIANTS A LITTLE TEDIOUS
+@c NEED A MORE CONVENIENT (AND NONAMBIGUOUS) NOTATION FOR OPTIONAL ARGUMENTS
+
+@c -----------------------------------------------------------------------------
+@anchor{quad_qagp}
+@deffn {Function} quad_qagp (@var{f(x)}, @var{x}, @var{a}, @var{b}, @var{points}, [@var{epsrel}, @var{epsabs}, @var{limit}])
+@deffnx {Function} quad_qagp (@var{f}, @var{x}, @var{a}, @var{b}, @var{points}, [@var{epsrel}, @var{epsabs}, @var{limit}])
+
+Integration of a general function over a finite interval.
+@code{quad_qagp} implements globally adaptive interval subdivision with
+extrapolation (de Doncker, 1978) by the Epsilon algorithm (Wynn, 1956).
+
+@code{quad_qagp} computes the integral
+
+@ifnottex
+@math{integrate (f(x), x, a, b)}
+@end ifnottex
+@tex
+$$\int_a^b {f(x) \, dx}$$
+@end tex
+
+The function to be integrated is @var{f(x)}, with
+dependent variable @var{x}, and the function is to be integrated
+between the limits @var{a} and @var{b}.
+
+The integrand may be specified as the name of a Maxima or Lisp function or
+operator, a Maxima lambda expression, or a general Maxima expression.
+
+To help the integrator, the user must supply a list of points where
+the integrand is singular or discontinous.
+
+The keyword arguments are optional and may be specified in any order.
+They all take the form @code{key=val}. The keyword arguments are:
+
+@table @var
+@item epsrel
+Desired relative error of approximation. Default is 1d-8.
+@item epsabs
+Desired absolute error of approximation. Default is 0.
+@item limit
+Size of internal work array. @var{limit} is the
+maximum number of subintervals to use. Default is 200.
+@end table
+
+@code{quad_qagp} returns a list of four elements:
+
+@itemize
+@item
+an approximation to the integral,
+@item
+the estimated absolute error of the approximation,
+@item
+the number integrand evaluations,
+@item
+an error code.
+@end itemize
+
+The error code (fourth element of the return value) can have the values:
+
+@table @code
+@item 0
+no problems were encountered;
+@item 1
+too many sub-intervals were done;
+@item 2
+excessive roundoff error is detected;
+@item 3
+extremely bad integrand behavior occurs;
+@item 4
+failed to converge
+@item 5
+integral is probably divergent or slowly convergent
+@item 6
+if the input is invalid.
+@end table
+
+@c NEED CROSS REFS HERE -- EITHER CROSS REF A QUADPACK OVERVIEW, OR CROSS REF EACH OF THE quad_* FUNCTIONS
+
+Examples:
+
+@example
+@group
+(%i1) quad_qagp(x^3*log(abs((x^2-1)*(x^2-2))), x, 0, 3, [1,sqrt(2)]);
+(%o1) [52.74074838347143, 2.6247632689546663e-7, 1029, 0]
+@end group
+@group
+(%i2) quad_qags(x^3*log(abs((x^2-1)*(x^2-2))), x, 0, 3);
+(%o2) [52.74074847951494, 4.088443219529836e-7, 1869, 0]
+@end group
+@end example
+
+The integrand has singularities at 1 and sqrt(2) so we supply these
+points to @code{quad_qagp}. We also note that @code{quad_qagp} is
+more accurate and more efficient that @code{quad_qags}.
+
+@opencatbox
+@category{Numerical methods} @category{Package quadpack}
+@closecatbox
+@end deffn
+
commit ffbfbd964574a6e725534f2518fb77f08c489ec4
Author: Raymond Toy <toy...@gm...>
Date: Sat Oct 15 22:13:42 2011 -0700
Update dependencies to include dqagp.lisp.
diff --git a/src/clisp-depends.mk b/src/clisp-depends.mk
index 79bfd00..67001c8 100644
--- a/src/clisp-depends.mk
+++ b/src/clisp-depends.mk
@@ -148,6 +148,7 @@ binary-clisp/maxima.mem : numerical/slatec/dqc25f.lisp
binary-clisp/maxima.mem : numerical/slatec/dqage.lisp
binary-clisp/maxima.mem : numerical/slatec/dqagie.lisp
binary-clisp/maxima.mem : numerical/slatec/dqagpe.lisp
+binary-clisp/maxima.mem : numerical/slatec/dqagp.lisp
binary-clisp/maxima.mem : numerical/slatec/dqagse.lisp
binary-clisp/maxima.mem : numerical/slatec/dqawoe.lisp
binary-clisp/maxima.mem : numerical/slatec/dqawfe.lisp
diff --git a/src/cmucl-depends.mk b/src/cmucl-depends.mk
index 5d9f7f5..d1f16dc 100644
--- a/src/cmucl-depends.mk
+++ b/src/cmucl-depends.mk
@@ -147,6 +147,7 @@ binary-cmucl/maxima.core : numerical/slatec/dqc25f.lisp
binary-cmucl/maxima.core : numerical/slatec/dqage.lisp
binary-cmucl/maxima.core : numerical/slatec/dqagie.lisp
binary-cmucl/maxima.core : numerical/slatec/dqagpe.lisp
+binary-cmucl/maxima.core : numerical/slatec/dqagp.lisp
binary-cmucl/maxima.core : numerical/slatec/dqagse.lisp
binary-cmucl/maxima.core : numerical/slatec/dqawoe.lisp
binary-cmucl/maxima.core : numerical/slatec/dqawfe.lisp
diff --git a/src/ecl-depends.mk b/src/ecl-depends.mk
index 04f5968..17a1c95 100644
--- a/src/ecl-depends.mk
+++ b/src/ecl-depends.mk
@@ -149,6 +149,7 @@ binary-ecl/maxima : numerical/slatec/dqc25f.lisp
binary-ecl/maxima : numerical/slatec/dqage.lisp
binary-ecl/maxima : numerical/slatec/dqagie.lisp
binary-ecl/maxima : numerical/slatec/dqagpe.lisp
+binary-ecl/maxima : numerical/slatec/dqagp.lisp
binary-ecl/maxima : numerical/slatec/dqagse.lisp
binary-ecl/maxima : numerical/slatec/dqawoe.lisp
binary-ecl/maxima : numerical/slatec/dqawfe.lisp
diff --git a/src/openmcl-depends.mk b/src/openmcl-depends.mk
index f6d00d0..5abfed1 100644
--- a/src/openmcl-depends.mk
+++ b/src/openmcl-depends.mk
@@ -148,6 +148,7 @@ binary-openmcl/maxima.image : numerical/slatec/dqc25f.lisp
binary-openmcl/maxima.image : numerical/slatec/dqage.lisp
binary-openmcl/maxima.image : numerical/slatec/dqagie.lisp
binary-openmcl/maxima.image : numerical/slatec/dqagpe.lisp
+binary-openmcl/maxima.image : numerical/slatec/dqagp.lisp
binary-openmcl/maxima.image : numerical/slatec/dqagse.lisp
binary-openmcl/maxima.image : numerical/slatec/dqawoe.lisp
binary-openmcl/maxima.image : numerical/slatec/dqawfe.lisp
diff --git a/src/sbcl-depends.mk b/src/sbcl-depends.mk
index acd2f2a..2151242 100644
--- a/src/sbcl-depends.mk
+++ b/src/sbcl-depends.mk
@@ -148,6 +148,7 @@ binary-sbcl/maxima.core : numerical/slatec/dqc25f.lisp
binary-sbcl/maxima.core : numerical/slatec/dqage.lisp
binary-sbcl/maxima.core : numerical/slatec/dqagie.lisp
binary-sbcl/maxima.core : numerical/slatec/dqagpe.lisp
+binary-sbcl/maxima.core : numerical/slatec/dqagp.lisp
binary-sbcl/maxima.core : numerical/slatec/dqagse.lisp
binary-sbcl/maxima.core : numerical/slatec/dqawoe.lisp
binary-sbcl/maxima.core : numerical/slatec/dqawfe.lisp
@@ -292,5 +293,6 @@ binary-sbcl/maxima.core : gnuplot_def.lisp
binary-sbcl/maxima.core : xmaxima_def.lisp
binary-sbcl/maxima.core : autol.lisp
binary-sbcl/maxima.core : max_ext.lisp
+binary-sbcl/maxima.core : share-subdirs.lisp
binary-sbcl/maxima.core : init-cl.lisp
commit 4e6576139f3f769f7249f6d75ad44a4e488faafb
Author: Raymond Toy <toy...@gm...>
Date: Sat Oct 15 22:11:12 2011 -0700
Fix a few compiler warnings.
diff --git a/src/numerical/slatec/quadpack.lisp b/src/numerical/slatec/quadpack.lisp
index fc6ade3..c6fd0a1 100644
--- a/src/numerical/slatec/quadpack.lisp
+++ b/src/numerical/slatec/quadpack.lisp
@@ -107,7 +107,7 @@
(work (make-array lenw :element-type 'flonum))
(iwork (make-array limit :element-type 'f2cl-lib:integer4))
(f (get-integrand fun var))
- (infinity (case inf-type
+ (infinity (ecase inf-type
((1 $inf)
;; Interval is [bound, infinity]
1)
@@ -299,8 +299,8 @@
(float-or-lose epsrel)
0.0 0.0 0 0
leniw lenw 0 iwork work)
- (declare (ignore junk z-a z-b z-int z-epsabs z-epsrel
- z-limit z-lenw last))
+ (declare (ignore junk z-a z-b z-npts z-points z-epsabs z-epsrel
+ z-leniw z-lenw last))
(list '(mlist) result abserr neval ier))
(error ()
`(($quad_qagp) ,fun ,var ,a ,b ,points
commit 9067349233cc217786983b07c23a2b78cf7d27bc
Author: Raymond Toy <toy...@gm...>
Date: Sat Oct 15 21:59:52 2011 -0700
Implement interface to quadpack QAGP
src/maxima.system:
o Add dqapg
src/numerical/f2cl-lib.lisp:
o Update to newer version. Should be backwared compatible with older
version.
src/numerical/slatec.lisp:
o Export dqagp
src/numerical/slatec/dqagp.lisp:
o New file containing translation of dqagp.f
src/numerical/slatec/quadpack.lisp:
o Implement interface to qagp.
src/numerical/slatec/quadpack.system:
o Add dqagp.
diff --git a/src/maxima.system b/src/maxima.system
index e89ec87..4e3c69b 100644
--- a/src/maxima.system
+++ b/src/maxima.system
@@ -430,6 +430,8 @@
:depends-on ("dqelg"
"dqk15i"
"dqpsrt"))
+ (:file "dqagp"
+ :depends-on ("dqagpe"))
(:file "dqagpe"
:depends-on ("dqelg"
"dqpsrt"
diff --git a/src/numerical/f2cl-lib.lisp b/src/numerical/f2cl-lib.lisp
index d6987b0..4705118 100644
--- a/src/numerical/f2cl-lib.lisp
+++ b/src/numerical/f2cl-lib.lisp
@@ -175,7 +175,7 @@ is not included")
;;
;; This is done by making a displaced array to VNAME with the
;; appropriate offset.
-(defmacro array-slice (vname type indices bounds)
+(defmacro array-slice (vname type indices bounds &optional offset)
;; To figure the size of the sliced array, use ARRAY-TOTAL-SIZE
;; instead of the f2cl derived/declared BOUNDS, just in case we
;; screwed up or in case we changed the size of the array in some
@@ -197,10 +197,16 @@ is not included")
;;
;; This seems somewhat reasonable, so let's do that for array
;; slices.
- `(make-array (max 0 (- (array-total-size ,vname) ,(col-major-index indices bounds)))
+ `(make-array (max 0 (- (array-total-size ,vname)
+ (the fixnum
+ (+ ,(col-major-index indices bounds)
+ (or ,offset 0)))))
:element-type ',type
:displaced-to ,vname
- :displaced-index-offset (min (array-total-size ,vname) ,(col-major-index indices bounds))))
+ :displaced-index-offset (min (array-total-size ,vname)
+ (the fixnum
+ (+ ,(col-major-index indices bounds)
+ (or ,offset 0))))))
;; Compute an initializer for make-array given the data in the list
;; DATA. The array has en element type of TYPE and has dimensions of
@@ -993,7 +999,7 @@ causing all pending operations to be flushed"
(open file :direction :io :if-exists :supersede
:if-does-not-exist :create))
((string-equal s "old")
- (open file :direction :io :if-does-not-exist nil))
+ (open file :direction :io :if-does-not-exist nil :if-exists :overwrite))
((string-equal s "new")
(open file :direction :io :if-exists nil))
(t
@@ -1461,9 +1467,64 @@ causing all pending operations to be flushed"
(defun stop (&optional arg)
(when arg
(format cl::*error-output* "~A~%" arg))
- (unless *stop-signals-error-p*
+ (when *stop-signals-error-p*
(cerror "Continue anyway" "STOP reached")))
+(defmacro f2cl-copy-seq (dst src dst-type src-type)
+ (flet ((copy-error ()
+ (error "F2CL cannot copy arrays of element type ~A to ~A~%"
+ src-type dst-type)))
+ (cond ((subtypep dst-type 'float)
+ ;; Copy to float array
+ (cond ((subtypep src-type 'float)
+ `(replace ,dst ,src))
+ ((subtypep src-type 'complex)
+ ;; Copy complex to float by putting each real and
+ ;; imaginary part into the float array, in order.
+ (let ((idx (gensym "IDX-"))
+ (el (gensym "EL-")))
+ `(loop for ,idx of-type fixnum from 0 by 2 below (length ,dst)
+ for ,el of-type ,src-type across ,src
+ do
+ (progn
+ (setf (aref ,dst ,idx) (realpart ,el))
+ (setf (aref ,dst (1+ ,idx)) (imagpart ,el))))))
+ (t
+ (copy-error))))
+ ((subtypep dst-type 'complex)
+ ;; Copy to complex array
+ (cond ((subtypep src-type 'float)
+ (let ((idx (gensym "IDX-"))
+ (dst-idx (gensym "DST-IDX-")))
+ `(loop for ,idx of-type fixnum from 0 by 2 below (length ,src)
+ for ,dst-idx of-type fixnum from 0 below (length ,dst)
+ do
+ (setf (aref ,dst ,dst-idx) (complex (aref ,src ,idx)
+ (aref ,src (1+ ,idx)))))))
+ ((subtypep src-type 'complex)
+ `(replace ,dst ,src))
+ (t
+ (copy-error))))
+ (t
+ (copy-error)))))
+
+(defmacro make-compatible-seq (type array array-type)
+ (let ((element-type (second type))
+ (array-type (second array-type)))
+ (cond ((subtypep element-type 'float)
+ (cond ((subtypep array-type 'complex)
+ `(make-array (* 2 (length ,array)) :element-type ',element-type))
+ (t
+ `(make-array (length ,array) :element-type ',element-type))))
+ ((subtypep element-type 'complex)
+ (cond ((subtypep array-type 'complex)
+ `(make-array (length ,array) :element-type ',element-type))
+ (t
+ `(make-array (ceiling (length ,array) 2) :element-type ',element-type))))
+ (t
+ (error "Don't know how to make an array with element-type ~A~%" element-type)))))
+
+
;;;-------------------------------------------------------------------------
;;; end of macros.l
;;;
diff --git a/src/numerical/slatec.lisp b/src/numerical/slatec.lisp
index 451cb83..1f63190 100644
--- a/src/numerical/slatec.lisp
+++ b/src/numerical/slatec.lisp
@@ -30,6 +30,6 @@
#:dspenc
;; Quadpack routines
- #:dqag #:dqags #:dqagi #:dqawc #:dqawf #:dqawo #:dqaws
+ #:dqag #:dqags #:dqagi #:dqawc #:dqawf #:dqawo #:dqaws #:dqagp
)
(:documentation "Package for the Fortran routines we need from SLATEC"))
diff --git a/src/numerical/slatec/dqagp.lisp b/src/numerical/slatec/dqagp.lisp
new file mode 100644
index 0000000..76690e5
--- /dev/null
+++ b/src/numerical/slatec/dqagp.lisp
@@ -0,0 +1,151 @@
+;;; Compiled by f2cl version:
+;;; ("f2cl1.l,v f0f149e72999 2010/10/08 03:05:30 rtoy $"
+;;; "f2cl2.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
+;;; "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
+;;; "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
+;;; "f2cl5.l,v 11bea7dae5a0 2011/06/11 17:53:39 toy $"
+;;; "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
+;;; "macros.l,v 11bea7dae5a0 2011/06/11 17:53:39 toy $")
+
+;;; Using Lisp CMU Common Lisp Snapshot 2011-10 ce1db1c (20B Unicode)
+;;;
+;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
+;;; (:coerce-assigns :as-needed) (:array-type ':array)
+;;; (:array-slicing t) (:declare-common nil)
+;;; (:float-format double-float))
+
+(in-package :slatec)
+
+
+(defun dqagp
+ (f a b npts2 points epsabs epsrel result abserr neval ier leniw lenw
+ last$ iwork work)
+ (declare (type (array f2cl-lib:integer4 (*)) iwork)
+ (type (array double-float (*)) work points)
+ (type (f2cl-lib:integer4) last$ lenw leniw ier neval npts2)
+ (type (double-float) abserr result epsrel epsabs b a))
+ (f2cl-lib:with-multi-array-data
+ ((points double-float points-%data% points-%offset%)
+ (work double-float work-%data% work-%offset%)
+ (iwork f2cl-lib:integer4 iwork-%data% iwork-%offset%))
+ (prog ((limit 0) (lvl 0) (l1 0) (l2 0) (l3 0) (l4 0))
+ (declare (type (f2cl-lib:integer4) l4 l3 l2 l1 lvl limit))
+ (setf ier 6)
+ (setf neval 0)
+ (setf last$ 0)
+ (setf result 0.0)
+ (setf abserr 0.0)
+ (if
+ (or (< leniw (f2cl-lib:int-sub (f2cl-lib:int-mul 3 npts2) 2))
+ (< lenw (f2cl-lib:int-sub (f2cl-lib:int-mul leniw 2) npts2))
+ (< npts2 2))
+ (go label10))
+ (setf limit (the f2cl-lib:integer4 (truncate (- leniw npts2) 2)))
+ (setf l1 (f2cl-lib:int-add limit 1))
+ (setf l2 (f2cl-lib:int-add limit l1))
+ (setf l3 (f2cl-lib:int-add limit l2))
+ (setf l4 (f2cl-lib:int-add limit l3))
+ (multiple-value-bind
+ (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
+ var-11 var-12 var-13 var-14 var-15 var-16 var-17 var-18 var-19
+ var-20)
+ (dqagpe f a b npts2 points epsabs epsrel limit result abserr neval
+ ier
+ (f2cl-lib:array-slice work-%data%
+ double-float
+ (1)
+ ((1 *))
+ work-%offset%)
+ (f2cl-lib:array-slice work-%data%
+ double-float
+ (l1)
+ ((1 *))
+ work-%offset%)
+ (f2cl-lib:array-slice work-%data%
+ double-float
+ (l2)
+ ((1 *))
+ work-%offset%)
+ (f2cl-lib:array-slice work-%data%
+ double-float
+ (l3)
+ ((1 *))
+ work-%offset%)
+ (f2cl-lib:array-slice work-%data%
+ double-float
+ (l4)
+ ((1 *))
+ work-%offset%)
+ (f2cl-lib:array-slice iwork-%data%
+ f2cl-lib:integer4
+ (1)
+ ((1 *))
+ iwork-%offset%)
+ (f2cl-lib:array-slice iwork-%data%
+ f2cl-lib:integer4
+ (l1)
+ ((1 *))
+ iwork-%offset%)
+ (f2cl-lib:array-slice iwork-%data%
+ f2cl-lib:integer4
+ (l2)
+ ((1 *))
+ iwork-%offset%)
+ last$)
+ (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-12
+ var-13 var-14 var-15 var-16 var-17 var-18 var-19))
+ (setf result var-8)
+ (setf abserr var-9)
+ (setf neval var-10)
+ (setf ier var-11)
+ (setf last$ var-20))
+ (setf lvl 0)
+ label10
+ (if (= ier 6) (setf lvl 1))
+ (if (/= ier 0)
+ (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
+ (xermsg "SLATEC" "DQAGP" "ABNORMAL RETURN" ier lvl)
+ (declare (ignore var-0 var-1 var-2))
+ (when var-3
+ (setf ier var-3))
+ (when var-4
+ (setf lvl var-4))))
+ (go end_label)
+ end_label
+ (return
+ (values nil
+ nil
+ nil
+ nil
+ nil
+ nil
+ nil
+ result
+ abserr
+ neval
+ ier
+ nil
+ nil
+ last$
+ nil
+ nil)))))
+
+(in-package #-gcl #:cl-user #+gcl "CL-USER")
+#+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
+(eval-when (:load-toplevel :compile-toplevel :execute)
+ (setf (gethash 'fortran-to-lisp::dqagp fortran-to-lisp::*f2cl-function-info*)
+ (fortran-to-lisp::make-f2cl-finfo
+ :arg-types '(t (double-float) (double-float)
+ (fortran-to-lisp::integer4) (array double-float (*))
+ (double-float) (double-float) (double-float)
+ (double-float) (fortran-to-lisp::integer4)
+ (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
+ (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
+ (array fortran-to-lisp::integer4 (*))
+ (array double-float (*)))
+ :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::result
+ fortran-to-lisp::abserr fortran-to-lisp::neval
+ fortran-to-lisp::ier nil nil fortran-to-lisp::last$
+ nil nil)
+ :calls '(fortran-to-lisp::dqagpe))))
+
diff --git a/src/numerical/slatec/quadpack.lisp b/src/numerical/slatec/quadpack.lisp
index 0c8f6f0..fc6ade3 100644
--- a/src/numerical/slatec/quadpack.lisp
+++ b/src/numerical/slatec/quadpack.lisp
@@ -274,6 +274,40 @@
((mequal) $epsabs ,epsabs)
((mequal) $limit ,limit))))))
+(defun quad-qagp (fun var a b points
+ &key (epsrel 1e-8) (epsabs 0.0) (limit 200))
+ (quad_argument_check fun var a b)
+ (let* ((npts2 (+ 2 (length (cdr points))))
+ (p (make-array npts2 :element-type 'flonum))
+ (leniw (max limit (- (* 3 npts2) 2)))
+ (lenw (- (* 2 leniw) npts2))
+ (work (make-array lenw :element-type 'flonum))
+ (iwork (make-array limit :element-type 'f2cl-lib:integer4))
+ (f (get-integrand fun var)))
+ (map-into p #'float-or-lose (cdr points))
+ (handler-case
+ (multiple-value-bind (junk z-a z-b z-npts z-points z-epsabs z-epsrel
+ result abserr neval ier
+ z-leniw z-lenw last)
+ (slatec:dqagp #'(lambda (x)
+ (float (funcall f x)))
+ (float-or-lose a)
+ (float-or-lose b)
+ npts2
+ p
+ (float-or-lose epsabs)
+ (float-or-lose epsrel)
+ 0.0 0.0 0 0
+ leniw lenw 0 iwork work)
+ (declare (ignore junk z-a z-b z-int z-epsabs z-epsrel
+ z-limit z-lenw last))
+ (list '(mlist) result abserr neval ier))
+ (error ()
+ `(($quad_qagp) ,fun ,var ,a ,b ,points
+ ((mequal) $epsrel ,epsrel)
+ ((mequal) $epsabs ,epsabs)
+ ((mequal) $limit ,limit))))))
+
;; error checking similar to that done by $defint
(defun quad_argument_check (exp var ll ul)
(setq exp (ratdisrep exp))
@@ -306,7 +340,8 @@
(frob $quad_qawc quad-qawc (fun var c a b) ($epsrel $limit $epsabs))
(frob $quad_qawf quad-qawf (fun var a omega trig) ($limit $epsabs $maxp1 $limlst))
(frob $quad_qawo quad-qawo (fun var a b omega trig) ($epsrel $limit $epsabs $maxp1))
- (frob $quad_qaws quad-qaws (fun var a b alfa beta wfun) ($epsrel $limit $epsabs)))
+ (frob $quad_qaws quad-qaws (fun var a b alfa beta wfun) ($epsrel $limit $epsabs))
+ (frob $quad_qagp quad-qagp (fun var a b points) ($epsrel $limit $epsabs)))
;; Tests
;;
diff --git a/src/numerical/slatec/quadpack.system b/src/numerical/slatec/quadpack.system
index 0a144cd..7fab907 100644
--- a/src/numerical/slatec/quadpack.system
+++ b/src/numerical/slatec/quadpack.system
@@ -31,6 +31,7 @@
:compiler-options (:float-format double-float :package "SLATEC")
:compile-only t
:binary-pathname (logical-pathname "slatec:")
+ :binary-extension "lisp"
:components
(
;; Support
@@ -78,6 +79,8 @@
:depends-on ("dqelg"
"dqk15i"
"dqpsrt"))
+ (:file "dqagp"
+ :depends-on ("dqagpe"))
(:file "dqagpe"
:depends-on ("dqelg"
"dqpsrt"
-----------------------------------------------------------------------
Summary of changes:
doc/info/Integration.texi | 106 ++++++++++++++++++++++++
src/clisp-depends.mk | 1 +
src/cmucl-depends.mk | 1 +
src/ecl-depends.mk | 1 +
src/maxima.system | 2 +
src/numerical/f2cl-lib.lisp | 71 +++++++++++++++-
src/numerical/slatec.lisp | 2 +-
src/numerical/slatec/dqagp.lisp | 151 ++++++++++++++++++++++++++++++++++
src/numerical/slatec/quadpack.lisp | 39 ++++++++-
src/numerical/slatec/quadpack.system | 3 +
src/openmcl-depends.mk | 1 +
src/sbcl-depends.mk | 2 +
12 files changed, 372 insertions(+), 8 deletions(-)
create mode 100644 src/numerical/slatec/dqagp.lisp
hooks/post-receive
--
Maxima, A Computer Algebra System
|