src/numerical/slatec/dqk61.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 15
  22                    :element-type 'double-float
  23                    :initial-contents '(0.007968192496166605 0.01846646831109096
  24                                        0.02878470788332337 0.03879919256962705
  25                                        0.04840267283059405 0.057493156217619065
  26                                        0.06597422988218049 0.0737559747377052
  27                                        0.08075589522942021 0.08689978720108298
  28                                        0.09212252223778612 0.09636873717464425
  29                                        0.09959342058679527 0.1017623897484055
  30                                        0.10285265289355884)))
  31       (xgk
  32        (make-array 31
  33                    :element-type 'double-float
  34                    :initial-contents '(0.9994844100504906 0.9968934840746495
  35                                        0.9916309968704046 0.9836681232797472
  36                                        0.9731163225011262 0.9600218649683075
  37                                        0.94437444474856 0.9262000474292743
  38                                        0.9055733076999078 0.8825605357920527
  39                                        0.8572052335460612 0.8295657623827684
  40                                        0.799727835821839 0.7677774321048262
  41                                        0.7337900624532268 0.6978504947933158
  42                                        0.6600610641266269 0.6205261829892429
  43                                        0.5793452358263617 0.5366241481420199
  44                                        0.49248046786177857 0.44703376953808915
  45                                        0.4004012548303944 0.3527047255308781
  46                                        0.30407320227362505 0.25463692616788985
  47                                        0.20452511668230988 0.15386991360858354
  48                                        0.10280693796673702 0.0514718425553177
  49                                        0.0)))
  50       (wgk
  51        (make-array 31
  52                    :element-type 'double-float
  53                    :initial-contents '(0.0013890136986770077
  54                                        0.003890461127099884
  55                                        0.0066307039159312926
  56                                        0.009273279659517764
  57                                        0.011823015253496341
  58                                        0.014369729507045804 0.01692088918905327
  59                                        0.019414141193942382
  60                                        0.021828035821609193 0.0241911620780806
  61                                        0.0265099548823331 0.02875404876504129
  62                                        0.030907257562387762 0.03298144705748372
  63                                        0.034979338028060025 0.03688236465182123
  64                                        0.038678945624727595
  65                                        0.040374538951535956
  66                                        0.041969810215164244 0.04345253970135607
  67                                        0.04481480013316266 0.04605923827100699
  68                                        0.04718554656929915 0.04818586175708713
  69                                        0.04905543455502978 0.04979568342707421
  70                                        0.05040592140278235 0.05088179589874961
  71                                        0.051221547849258774 0.05142612853745902
  72                                        0.05149472942945157))))
  73   (declare (type (array double-float (15)) wg)
  74            (type (array double-float (31)) xgk wgk))
  75   (defun dqk61 (f a b result abserr resabs resasc)
  76     (declare (type (double-float) resasc resabs abserr result b a))
  77     (prog ((fv1 (make-array 30 :element-type 'double-float))
  78            (fv2 (make-array 30 :element-type 'double-float)) (j 0) (jtw 0)
  79            (jtwm1 0) (dabsc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0)
  80            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  81            (reskh 0.0) (uflow 0.0))
  82       (declare (type (array double-float (30)) fv2 fv1)
  83                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  84                                     fsum fc epmach dhlgth centr dabsc)
  85                (type (f2cl-lib:integer4) jtwm1 jtw j))
  86       (setf epmach (f2cl-lib:d1mach 4))
  87       (setf uflow (f2cl-lib:d1mach 1))
  88       (setf centr (* 0.5 (+ b a)))
  89       (setf hlgth (* 0.5 (- b a)))
  90       (setf dhlgth (abs hlgth))
  91       (setf resg 0.0)
  92       (setf fc
  93               (multiple-value-bind (ret-val var-0)
  94                   (funcall f centr)
  95                 (declare (ignore))
  96                 (when var-0
  97                   (setf centr var-0))
  98                 ret-val))
  99       (setf resk (* (f2cl-lib:fref wgk (31) ((1 31))) fc))
 100       (setf resabs (abs resk))
 101       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 102                     ((> j 15) nil)
 103         (tagbody
 104           (setf jtw (f2cl-lib:int-mul j 2))
 105           (setf dabsc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 31)))))
 106           (setf fval1 (funcall f (- centr dabsc)))
 107           (setf fval2 (funcall f (+ centr dabsc)))
 108           (setf (f2cl-lib:fref fv1 (jtw) ((1 30))) fval1)
 109           (setf (f2cl-lib:fref fv2 (jtw) ((1 30))) fval2)
 110           (setf fsum (+ fval1 fval2))
 111           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 15))) fsum)))
 112           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 31))) fsum)))
 113           (setf resabs
 114                   (+ resabs
 115                      (* (f2cl-lib:fref wgk (jtw) ((1 31)))
 116                         (+ (abs fval1) (abs fval2)))))
 117          label10))
 118       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 119                     ((> j 15) nil)
 120         (tagbody
 121           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1))
 122           (setf dabsc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 31)))))
 123           (setf fval1 (funcall f (- centr dabsc)))
 124           (setf fval2 (funcall f (+ centr dabsc)))
 125           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 30))) fval1)
 126           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 30))) fval2)
 127           (setf fsum (+ fval1 fval2))
 128           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 31))) fsum)))
 129           (setf resabs
 130                   (+ resabs
 131                      (* (f2cl-lib:fref wgk (jtwm1) ((1 31)))
 132                         (+ (abs fval1) (abs fval2)))))
 133          label15))
 134       (setf reskh (* resk 0.5))
 135       (setf resasc (* (f2cl-lib:fref wgk (31) ((1 31))) (abs (- fc reskh))))
 136       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 137                     ((> j 30) nil)
 138         (tagbody
 139           (setf resasc
 140                   (+ resasc
 141                      (* (f2cl-lib:fref wgk (j) ((1 31)))
 142                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 30))) reskh))
 143                            (abs (- (f2cl-lib:fref fv2 (j) ((1 30))) reskh))))))
 144          label20))
 145       (setf result (* resk hlgth))
 146       (setf resabs (* resabs dhlgth))
 147       (setf resasc (* resasc dhlgth))
 148       (setf abserr (abs (* (- resk resg) hlgth)))
 149       (if (and (/= resasc 0.0) (/= abserr 0.0))
 150           (setf abserr
 151                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 152       (if (> resabs (/ uflow (* 50.0 epmach)))
 153           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 154       (go end_label)
 155      end_label
 156       (return (values nil nil nil result abserr resabs resasc)))))
 157
 158 (in-package #:cl-user)
 159 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 160 (eval-when (:load-toplevel :compile-toplevel :execute)
 161   (setf (gethash 'fortran-to-lisp::dqk61 fortran-to-lisp::*f2cl-function-info*)
 162           (fortran-to-lisp::make-f2cl-finfo
 163            :arg-types '(t (double-float) (double-float) (double-float)
 164                         (double-float) (double-float) (double-float))
 165            :return-values '(nil nil nil fortran-to-lisp::result
 166                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 167                             fortran-to-lisp::resasc)
 168            :calls '(fortran-to-lisp::d1mach))))
 169