src/numerical/slatec/dqk15i.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 8
  22                    :element-type 'double-float
  23                    :initial-contents '(0.0 0.1294849661688697 0.0
  24                                        0.27970539148927664 0.0
  25                                        0.3818300505051189 0.0
  26                                        0.4179591836734694)))
  27       (xgk
  28        (make-array 8
  29                    :element-type 'double-float
  30                    :initial-contents '(0.9914553711208126 0.9491079123427585
  31                                        0.8648644233597691 0.7415311855993945
  32                                        0.5860872354676911 0.4058451513773972
  33                                        0.20778495500789848 0.0)))
  34       (wgk
  35        (make-array 8
  36                    :element-type 'double-float
  37                    :initial-contents '(0.022935322010529224 0.06309209262997856
  38                                        0.10479001032225019 0.14065325971552592
  39                                        0.1690047266392679 0.19035057806478542
  40                                        0.20443294007529889
  41                                        0.20948214108472782))))
  42   (declare (type (array double-float (8)) wg xgk wgk))
  43   (defun dqk15i (f boun inf a b result abserr resabs resasc)
  44     (declare (type (f2cl-lib:integer4) inf)
  45              (type (double-float) resasc resabs abserr result b a boun))
  46     (prog ((fv1 (make-array 7 :element-type 'double-float))
  47            (fv2 (make-array 7 :element-type 'double-float)) (j 0) (absc 0.0)
  48            (absc1 0.0) (absc2 0.0) (centr 0.0) (dinf 0.0) (epmach 0.0) (fc 0.0)
  49            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  50            (reskh 0.0) (tabsc1 0.0) (tabsc2 0.0) (uflow 0.0))
  51       (declare (type (array double-float (7)) fv2 fv1)
  52                (type (double-float) uflow tabsc2 tabsc1 reskh resk resg hlgth
  53                                     fval2 fval1 fsum fc epmach dinf centr absc2
  54                                     absc1 absc)
  55                (type (f2cl-lib:integer4) j))
  56       (setf epmach (f2cl-lib:d1mach 4))
  57       (setf uflow (f2cl-lib:d1mach 1))
  58       (setf dinf
  59               (coerce
  60                (the f2cl-lib:integer4
  61                     (min (the f2cl-lib:integer4 1)
  62                          (the f2cl-lib:integer4 inf)))
  63                'double-float))
  64       (setf centr (* 0.5 (+ a b)))
  65       (setf hlgth (* 0.5 (- b a)))
  66       (setf tabsc1 (+ boun (/ (* dinf (- 1.0 centr)) centr)))
  67       (setf fval1
  68               (multiple-value-bind (ret-val var-0)
  69                   (funcall f tabsc1)
  70                 (declare (ignore))
  71                 (when var-0
  72                   (setf tabsc1 var-0))
  73                 ret-val))
  74       (if (= inf 2) (setf fval1 (+ fval1 (funcall f (- tabsc1)))))
  75       (setf fc (/ (/ fval1 centr) centr))
  76       (setf resg (* (f2cl-lib:fref wg (8) ((1 8))) fc))
  77       (setf resk (* (f2cl-lib:fref wgk (8) ((1 8))) fc))
  78       (setf resabs (abs resk))
  79       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  80                     ((> j 7) nil)
  81         (tagbody
  82           (setf absc (* hlgth (f2cl-lib:fref xgk (j) ((1 8)))))
  83           (setf absc1 (- centr absc))
  84           (setf absc2 (+ centr absc))
  85           (setf tabsc1 (+ boun (/ (* dinf (- 1.0 absc1)) absc1)))
  86           (setf tabsc2 (+ boun (/ (* dinf (- 1.0 absc2)) absc2)))
  87           (setf fval1
  88                   (multiple-value-bind (ret-val var-0)
  89                       (funcall f tabsc1)
  90                     (declare (ignore))
  91                     (when var-0
  92                       (setf tabsc1 var-0))
  93                     ret-val))
  94           (setf fval2
  95                   (multiple-value-bind (ret-val var-0)
  96                       (funcall f tabsc2)
  97                     (declare (ignore))
  98                     (when var-0
  99                       (setf tabsc2 var-0))
 100                     ret-val))
 101           (if (= inf 2) (setf fval1 (+ fval1 (funcall f (- tabsc1)))))
 102           (if (= inf 2) (setf fval2 (+ fval2 (funcall f (- tabsc2)))))
 103           (setf fval1 (/ (/ fval1 absc1) absc1))
 104           (setf fval2 (/ (/ fval2 absc2) absc2))
 105           (setf (f2cl-lib:fref fv1 (j) ((1 7))) fval1)
 106           (setf (f2cl-lib:fref fv2 (j) ((1 7))) fval2)
 107           (setf fsum (+ fval1 fval2))
 108           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 8))) fsum)))
 109           (setf resk (+ resk (* (f2cl-lib:fref wgk (j) ((1 8))) fsum)))
 110           (setf resabs
 111                   (+ resabs
 112                      (* (f2cl-lib:fref wgk (j) ((1 8)))
 113                         (+ (abs fval1) (abs fval2)))))
 114          label10))
 115       (setf reskh (* resk 0.5))
 116       (setf resasc (* (f2cl-lib:fref wgk (8) ((1 8))) (abs (- fc reskh))))
 117       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 118                     ((> j 7) nil)
 119         (tagbody
 120           (setf resasc
 121                   (+ resasc
 122                      (* (f2cl-lib:fref wgk (j) ((1 8)))
 123                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 7))) reskh))
 124                            (abs (- (f2cl-lib:fref fv2 (j) ((1 7))) reskh))))))
 125          label20))
 126       (setf result (* resk hlgth))
 127       (setf resasc (* resasc hlgth))
 128       (setf resabs (* resabs hlgth))
 129       (setf abserr (abs (* (- resk resg) hlgth)))
 130       (if (and (/= resasc 0.0) (/= abserr 0.0))
 131           (setf abserr
 132                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 133       (if (> resabs (/ uflow (* 50.0 epmach)))
 134           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 135       (go end_label)
 136      end_label
 137       (return (values nil nil nil nil nil result abserr resabs resasc)))))
 138
 139 (in-package #:cl-user)
 140 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 141 (eval-when (:load-toplevel :compile-toplevel :execute)
 142   (setf (gethash 'fortran-to-lisp::dqk15i
 143                  fortran-to-lisp::*f2cl-function-info*)
 144           (fortran-to-lisp::make-f2cl-finfo
 145            :arg-types '(t (double-float) (fortran-to-lisp::integer4)
 146                         (double-float) (double-float) (double-float)
 147                         (double-float) (double-float) (double-float))
 148            :return-values '(nil nil nil nil nil fortran-to-lisp::result
 149                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 150                             fortran-to-lisp::resasc)
 151            :calls '(fortran-to-lisp::d1mach))))
 152