src/numerical/slatec/dqk51.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 13
  22                    :element-type 'double-float
  23                    :initial-contents '(0.011393798501026288
  24                                        0.026354986615032137
  25                                        0.040939156701306316
  26                                        0.054904695975835194 0.06803833381235691
  27                                        0.08014070033500102 0.09102826198296365
  28                                        0.10053594906705064 0.10851962447426365
  29                                        0.11485825914571164 0.11945576353578477
  30                                        0.12224244299031004
  31                                        0.12317605372671545)))
  32       (xgk
  33        (make-array 26
  34                    :element-type 'double-float
  35                    :initial-contents '(0.9992621049926098 0.9955569697904981
  36                                        0.9880357945340772 0.9766639214595175
  37                                        0.9616149864258425 0.9429745712289743
  38                                        0.9207471152817016 0.8949919978782753
  39                                        0.8658470652932756 0.833442628760834
  40                                        0.7978737979985001 0.7592592630373576
  41                                        0.7177664068130843 0.6735663684734684
  42                                        0.6268100990103174 0.577662930241223
  43                                        0.5263252843347191 0.473002731445715
  44                                        0.4178853821930377 0.36117230580938786
  45                                        0.30308953893110785 0.24386688372098844
  46                                        0.1837189394210489 0.1228646926107104
  47                                        0.06154448300568508 0.0)))
  48       (wgk
  49        (make-array 26
  50                    :element-type 'double-float
  51                    :initial-contents '(0.001987383892330316
  52                                        0.005561932135356714
  53                                        0.009473973386174152
  54                                        0.013236229195571676 0.0168478177091283
  55                                        0.020435371145882834
  56                                        0.024009945606953215 0.02747531758785174
  57                                        0.030792300167387487
  58                                        0.034002130274329335 0.03711627148341554
  59                                        0.04008382550403238 0.04287284502017005
  60                                        0.04550291304992179 0.04798253713883671
  61                                        0.05027767908071567 0.05236288580640747
  62                                        0.05425112988854549 0.055950811220412316
  63                                        0.057437116361567835
  64                                        0.058689680022394206 0.05972034032417406
  65                                        0.06053945537604586 0.061128509717053046
  66                                        0.061471189871425316
  67                                        0.061580818067832936))))
  68   (declare (type (array double-float (13)) wg)
  69            (type (array double-float (26)) xgk wgk))
  70   (defun dqk51 (f a b result abserr resabs resasc)
  71     (declare (type (double-float) resasc resabs abserr result b a))
  72     (prog ((fv1 (make-array 25 :element-type 'double-float))
  73            (fv2 (make-array 25 :element-type 'double-float)) (j 0) (jtw 0)
  74            (jtwm1 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0)
  75            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  76            (reskh 0.0) (uflow 0.0))
  77       (declare (type (array double-float (25)) fv2 fv1)
  78                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  79                                     fsum fc epmach dhlgth centr absc)
  80                (type (f2cl-lib:integer4) jtwm1 jtw j))
  81       (setf epmach (f2cl-lib:d1mach 4))
  82       (setf uflow (f2cl-lib:d1mach 1))
  83       (setf centr (* 0.5 (+ a b)))
  84       (setf hlgth (* 0.5 (- b a)))
  85       (setf dhlgth (abs hlgth))
  86       (setf fc
  87               (multiple-value-bind (ret-val var-0)
  88                   (funcall f centr)
  89                 (declare (ignore))
  90                 (when var-0
  91                   (setf centr var-0))
  92                 ret-val))
  93       (setf resg (* (f2cl-lib:fref wg (13) ((1 13))) fc))
  94       (setf resk (* (f2cl-lib:fref wgk (26) ((1 26))) fc))
  95       (setf resabs (abs resk))
  96       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  97                     ((> j 12) nil)
  98         (tagbody
  99           (setf jtw (f2cl-lib:int-mul j 2))
 100           (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 26)))))
 101           (setf fval1 (funcall f (- centr absc)))
 102           (setf fval2 (funcall f (+ centr absc)))
 103           (setf (f2cl-lib:fref fv1 (jtw) ((1 25))) fval1)
 104           (setf (f2cl-lib:fref fv2 (jtw) ((1 25))) fval2)
 105           (setf fsum (+ fval1 fval2))
 106           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 13))) fsum)))
 107           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 26))) fsum)))
 108           (setf resabs
 109                   (+ resabs
 110                      (* (f2cl-lib:fref wgk (jtw) ((1 26)))
 111                         (+ (abs fval1) (abs fval2)))))
 112          label10))
 113       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 114                     ((> j 13) nil)
 115         (tagbody
 116           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1))
 117           (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 26)))))
 118           (setf fval1 (funcall f (- centr absc)))
 119           (setf fval2 (funcall f (+ centr absc)))
 120           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 25))) fval1)
 121           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 25))) fval2)
 122           (setf fsum (+ fval1 fval2))
 123           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 26))) fsum)))
 124           (setf resabs
 125                   (+ resabs
 126                      (* (f2cl-lib:fref wgk (jtwm1) ((1 26)))
 127                         (+ (abs fval1) (abs fval2)))))
 128          label15))
 129       (setf reskh (* resk 0.5))
 130       (setf resasc (* (f2cl-lib:fref wgk (26) ((1 26))) (abs (- fc reskh))))
 131       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 132                     ((> j 25) nil)
 133         (tagbody
 134           (setf resasc
 135                   (+ resasc
 136                      (* (f2cl-lib:fref wgk (j) ((1 26)))
 137                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 25))) reskh))
 138                            (abs (- (f2cl-lib:fref fv2 (j) ((1 25))) reskh))))))
 139          label20))
 140       (setf result (* resk hlgth))
 141       (setf resabs (* resabs dhlgth))
 142       (setf resasc (* resasc dhlgth))
 143       (setf abserr (abs (* (- resk resg) hlgth)))
 144       (if (and (/= resasc 0.0) (/= abserr 0.0))
 145           (setf abserr
 146                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 147       (if (> resabs (/ uflow (* 50.0 epmach)))
 148           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 149       (go end_label)
 150      end_label
 151       (return (values nil nil nil result abserr resabs resasc)))))
 152
 153 (in-package #:cl-user)
 154 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 155 (eval-when (:load-toplevel :compile-toplevel :execute)
 156   (setf (gethash 'fortran-to-lisp::dqk51 fortran-to-lisp::*f2cl-function-info*)
 157           (fortran-to-lisp::make-f2cl-finfo
 158            :arg-types '(t (double-float) (double-float) (double-float)
 159                         (double-float) (double-float) (double-float))
 160            :return-values '(nil nil nil fortran-to-lisp::result
 161                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 162                             fortran-to-lisp::resasc)
 163            :calls '(fortran-to-lisp::d1mach))))
 164