src/numerical/slatec/dqk15.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $"
   3 ;;;  "f2cl2.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   4 ;;;  "f2cl3.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   5 ;;;  "f2cl4.l,v 96616d88fb7e 2008/02/22 22:19:34 rtoy $"
   6 ;;;  "f2cl5.l,v 46c1f6a93b0d 2012/05/03 04:40:28 toy $"
   7 ;;;  "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
   8 ;;;  "macros.l,v fceac530ef0c 2011/11/26 04:02:26 toy $")
   9
  10 ;;; Using Lisp CMU Common Lisp snapshot-2012-04 (20C Unicode)
  11 ;;;
  12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  13 ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((wg
  21        (make-array 4
  22                    :element-type 'double-float
  23                    :initial-contents '(0.1294849661688697 0.27970539148927664
  24                                        0.3818300505051189 0.4179591836734694)))
  25       (xgk
  26        (make-array 8
  27                    :element-type 'double-float
  28                    :initial-contents '(0.9914553711208126 0.9491079123427585
  29                                        0.8648644233597691 0.7415311855993945
  30                                        0.5860872354676911 0.4058451513773972
  31                                        0.20778495500789848 0.0)))
  32       (wgk
  33        (make-array 8
  34                    :element-type 'double-float
  35                    :initial-contents '(0.022935322010529224 0.06309209262997856
  36                                        0.10479001032225019 0.14065325971552592
  37                                        0.1690047266392679 0.19035057806478542
  38                                        0.20443294007529889
  39                                        0.20948214108472782))))
  40   (declare (type (array double-float (4)) wg)
  41            (type (array double-float (8)) xgk wgk))
  42   (defun dqk15 (f a b result abserr resabs resasc)
  43     (declare (type (double-float) resasc resabs abserr result b a))
  44     (prog ((fv1 (make-array 7 :element-type 'double-float))
  45            (fv2 (make-array 7 :element-type 'double-float)) (j 0) (jtw 0)
  46            (jtwm1 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0)
  47            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  48            (reskh 0.0) (uflow 0.0))
  49       (declare (type (array double-float (7)) fv2 fv1)
  50                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  51                                     fsum fc epmach dhlgth centr absc)
  52                (type (f2cl-lib:integer4) jtwm1 jtw j))
  53       (setf epmach (f2cl-lib:d1mach 4))
  54       (setf uflow (f2cl-lib:d1mach 1))
  55       (setf centr (* 0.5 (+ a b)))
  56       (setf hlgth (* 0.5 (- b a)))
  57       (setf dhlgth (abs hlgth))
  58       (setf fc
  59               (multiple-value-bind (ret-val var-0)
  60                   (funcall f centr)
  61                 (declare (ignore))
  62                 (when var-0
  63                   (setf centr var-0))
  64                 ret-val))
  65       (setf resg (* fc (f2cl-lib:fref wg (4) ((1 4)))))
  66       (setf resk (* fc (f2cl-lib:fref wgk (8) ((1 8)))))
  67       (setf resabs (abs resk))
  68       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  69                     ((> j 3) nil)
  70         (tagbody
  71           (setf jtw (f2cl-lib:int-mul j 2))
  72           (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 8)))))
  73           (setf fval1 (funcall f (- centr absc)))
  74           (setf fval2 (funcall f (+ centr absc)))
  75           (setf (f2cl-lib:fref fv1 (jtw) ((1 7))) fval1)
  76           (setf (f2cl-lib:fref fv2 (jtw) ((1 7))) fval2)
  77           (setf fsum (+ fval1 fval2))
  78           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 4))) fsum)))
  79           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 8))) fsum)))
  80           (setf resabs
  81                   (+ resabs
  82                      (* (f2cl-lib:fref wgk (jtw) ((1 8)))
  83                         (+ (abs fval1) (abs fval2)))))
  84          label10))
  85       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  86                     ((> j 4) nil)
  87         (tagbody
  88           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1))
  89           (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 8)))))
  90           (setf fval1 (funcall f (- centr absc)))
  91           (setf fval2 (funcall f (+ centr absc)))
  92           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 7))) fval1)
  93           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 7))) fval2)
  94           (setf fsum (+ fval1 fval2))
  95           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 8))) fsum)))
  96           (setf resabs
  97                   (+ resabs
  98                      (* (f2cl-lib:fref wgk (jtwm1) ((1 8)))
  99                         (+ (abs fval1) (abs fval2)))))
 100          label15))
 101       (setf reskh (* resk 0.5))
 102       (setf resasc (* (f2cl-lib:fref wgk (8) ((1 8))) (abs (- fc reskh))))
 103       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 104                     ((> j 7) nil)
 105         (tagbody
 106           (setf resasc
 107                   (+ resasc
 108                      (* (f2cl-lib:fref wgk (j) ((1 8)))
 109                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 7))) reskh))
 110                            (abs (- (f2cl-lib:fref fv2 (j) ((1 7))) reskh))))))
 111          label20))
 112       (setf result (* resk hlgth))
 113       (setf resabs (* resabs dhlgth))
 114       (setf resasc (* resasc dhlgth))
 115       (setf abserr (abs (* (- resk resg) hlgth)))
 116       (if (and (/= resasc 0.0) (/= abserr 0.0))
 117           (setf abserr
 118                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 119       (if (> resabs (/ uflow (* 50.0 epmach)))
 120           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 121       (go end_label)
 122      end_label
 123       (return (values nil nil nil result abserr resabs resasc)))))
 124
 125 (in-package #:cl-user)
 126 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 127 (eval-when (:load-toplevel :compile-toplevel :execute)
 128   (setf (gethash 'fortran-to-lisp::dqk15 fortran-to-lisp::*f2cl-function-info*)
 129           (fortran-to-lisp::make-f2cl-finfo
 130            :arg-types '(t (double-float) (double-float) (double-float)
 131                         (double-float) (double-float) (double-float))
 132            :return-values '(nil nil nil fortran-to-lisp::result
 133                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 134                             fortran-to-lisp::resasc)
 135            :calls '(fortran-to-lisp::d1mach))))
 136