src/numerical/slatec/dqk31.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.03075324199611727 0.07036604748810812
  24                                        0.10715922046717194 0.13957067792615432
  25                                        0.16626920581699392 0.1861610000155622
  26                                        0.19843148532711158
  27                                        0.2025782419255613)))
  28       (xgk
  29        (make-array 16
  30                    :element-type 'double-float
  31                    :initial-contents '(0.9980022986933971 0.9879925180204854
  32                                        0.9677390756791391 0.937273392400706
  33                                        0.8972645323440819 0.8482065834104272
  34                                        0.790418501442466 0.7244177313601701
  35                                        0.650996741297417 0.5709721726085388
  36                                        0.4850818636402397 0.3941513470775634
  37                                        0.29918000715316884 0.20119409399743451
  38                                        0.1011420669187175 0.0)))
  39       (wgk
  40        (make-array 16
  41                    :element-type 'double-float
  42                    :initial-contents '(0.005377479872923349
  43                                        0.015007947329316122 0.02546084732671532
  44                                        0.03534636079137585 0.04458975132476488
  45                                        0.05348152469092809 0.06200956780067064
  46                                        0.06985412131872826 0.07684968075772038
  47                                        0.08308050282313302 0.08856444305621176
  48                                        0.09312659817082532 0.09664272698362368
  49                                        0.09917359872179196 0.10076984552387559
  50                                        0.10133000701479154))))
  51   (declare (type (array double-float (8)) wg)
  52            (type (array double-float (16)) xgk wgk))
  53   (defun dqk31 (f a b result abserr resabs resasc)
  54     (declare (type (double-float) resasc resabs abserr result b a))
  55     (prog ((fv1 (make-array 15 :element-type 'double-float))
  56            (fv2 (make-array 15 :element-type 'double-float)) (j 0) (jtw 0)
  57            (jtwm1 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0)
  58            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  59            (reskh 0.0) (uflow 0.0))
  60       (declare (type (array double-float (15)) fv2 fv1)
  61                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  62                                     fsum fc epmach dhlgth centr absc)
  63                (type (f2cl-lib:integer4) jtwm1 jtw j))
  64       (setf epmach (f2cl-lib:d1mach 4))
  65       (setf uflow (f2cl-lib:d1mach 1))
  66       (setf centr (* 0.5 (+ a b)))
  67       (setf hlgth (* 0.5 (- b a)))
  68       (setf dhlgth (abs hlgth))
  69       (setf fc
  70               (multiple-value-bind (ret-val var-0)
  71                   (funcall f centr)
  72                 (declare (ignore))
  73                 (when var-0
  74                   (setf centr var-0))
  75                 ret-val))
  76       (setf resg (* (f2cl-lib:fref wg (8) ((1 8))) fc))
  77       (setf resk (* (f2cl-lib:fref wgk (16) ((1 16))) 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 jtw (f2cl-lib:int-mul j 2))
  83           (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 16)))))
  84           (setf fval1 (funcall f (- centr absc)))
  85           (setf fval2 (funcall f (+ centr absc)))
  86           (setf (f2cl-lib:fref fv1 (jtw) ((1 15))) fval1)
  87           (setf (f2cl-lib:fref fv2 (jtw) ((1 15))) fval2)
  88           (setf fsum (+ fval1 fval2))
  89           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 8))) fsum)))
  90           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 16))) fsum)))
  91           (setf resabs
  92                   (+ resabs
  93                      (* (f2cl-lib:fref wgk (jtw) ((1 16)))
  94                         (+ (abs fval1) (abs fval2)))))
  95          label10))
  96       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  97                     ((> j 8) nil)
  98         (tagbody
  99           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1))
 100           (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 16)))))
 101           (setf fval1 (funcall f (- centr absc)))
 102           (setf fval2 (funcall f (+ centr absc)))
 103           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 15))) fval1)
 104           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 15))) fval2)
 105           (setf fsum (+ fval1 fval2))
 106           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 16))) fsum)))
 107           (setf resabs
 108                   (+ resabs
 109                      (* (f2cl-lib:fref wgk (jtwm1) ((1 16)))
 110                         (+ (abs fval1) (abs fval2)))))
 111          label15))
 112       (setf reskh (* resk 0.5))
 113       (setf resasc (* (f2cl-lib:fref wgk (16) ((1 16))) (abs (- fc reskh))))
 114       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 115                     ((> j 15) nil)
 116         (tagbody
 117           (setf resasc
 118                   (+ resasc
 119                      (* (f2cl-lib:fref wgk (j) ((1 16)))
 120                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 15))) reskh))
 121                            (abs (- (f2cl-lib:fref fv2 (j) ((1 15))) reskh))))))
 122          label20))
 123       (setf result (* resk hlgth))
 124       (setf resabs (* resabs dhlgth))
 125       (setf resasc (* resasc dhlgth))
 126       (setf abserr (abs (* (- resk resg) hlgth)))
 127       (if (and (/= resasc 0.0) (/= abserr 0.0))
 128           (setf abserr
 129                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 130       (if (> resabs (/ uflow (* 50.0 epmach)))
 131           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 132       (go end_label)
 133      end_label
 134       (return (values nil nil nil result abserr resabs resasc)))))
 135
 136 (in-package #:cl-user)
 137 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 138 (eval-when (:load-toplevel :compile-toplevel :execute)
 139   (setf (gethash 'fortran-to-lisp::dqk31 fortran-to-lisp::*f2cl-function-info*)
 140           (fortran-to-lisp::make-f2cl-finfo
 141            :arg-types '(t (double-float) (double-float) (double-float)
 142                         (double-float) (double-float) (double-float))
 143            :return-values '(nil nil nil fortran-to-lisp::result
 144                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 145                             fortran-to-lisp::resasc)
 146            :calls '(fortran-to-lisp::d1mach))))
 147