src/numerical/slatec/dqk21.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 5
  22                    :element-type 'double-float
  23                    :initial-contents '(0.06667134430868814 0.1494513491505806
  24                                        0.21908636251598204 0.26926671930999635
  25                                        0.29552422471475287)))
  26       (xgk
  27        (make-array 11
  28                    :element-type 'double-float
  29                    :initial-contents '(0.9956571630258081 0.9739065285171717
  30                                        0.9301574913557082 0.8650633666889845
  31                                        0.7808177265864169 0.6794095682990244
  32                                        0.5627571346686047 0.4333953941292472
  33                                        0.2943928627014602 0.14887433898163122
  34                                        0.0)))
  35       (wgk
  36        (make-array 11
  37                    :element-type 'double-float
  38                    :initial-contents '(0.011694638867371874
  39                                        0.032558162307964725
  40                                        0.054755896574351995 0.07503967481091996
  41                                        0.0931254545836976 0.10938715880229764
  42                                        0.12349197626206584 0.13470921731147334
  43                                        0.14277593857706009 0.14773910490133849
  44                                        0.1494455540029169))))
  45   (declare (type (array double-float (5)) wg)
  46            (type (array double-float (11)) xgk wgk))
  47   (defun dqk21 (f a b result abserr resabs resasc)
  48     (declare (type (double-float) resasc resabs abserr result b a))
  49     (prog ((fv1 (make-array 10 :element-type 'double-float))
  50            (fv2 (make-array 10 :element-type 'double-float)) (j 0) (jtw 0)
  51            (jtwm1 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fc 0.0)
  52            (fsum 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0)
  53            (reskh 0.0) (uflow 0.0))
  54       (declare (type (array double-float (10)) fv2 fv1)
  55                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  56                                     fsum fc epmach dhlgth centr absc)
  57                (type (f2cl-lib:integer4) jtwm1 jtw j))
  58       (setf epmach (f2cl-lib:d1mach 4))
  59       (setf uflow (f2cl-lib:d1mach 1))
  60       (setf centr (* 0.5 (+ a b)))
  61       (setf hlgth (* 0.5 (- b a)))
  62       (setf dhlgth (abs hlgth))
  63       (setf resg 0.0)
  64       (setf fc
  65               (multiple-value-bind (ret-val var-0)
  66                   (funcall f centr)
  67                 (declare (ignore))
  68                 (when var-0
  69                   (setf centr var-0))
  70                 ret-val))
  71       (setf resk (* (f2cl-lib:fref wgk (11) ((1 11))) fc))
  72       (setf resabs (abs resk))
  73       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  74                     ((> j 5) nil)
  75         (tagbody
  76           (setf jtw (f2cl-lib:int-mul 2 j))
  77           (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 11)))))
  78           (setf fval1 (funcall f (- centr absc)))
  79           (setf fval2 (funcall f (+ centr absc)))
  80           (setf (f2cl-lib:fref fv1 (jtw) ((1 10))) fval1)
  81           (setf (f2cl-lib:fref fv2 (jtw) ((1 10))) fval2)
  82           (setf fsum (+ fval1 fval2))
  83           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 5))) fsum)))
  84           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 11))) fsum)))
  85           (setf resabs
  86                   (+ resabs
  87                      (* (f2cl-lib:fref wgk (jtw) ((1 11)))
  88                         (+ (abs fval1) (abs fval2)))))
  89          label10))
  90       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  91                     ((> j 5) nil)
  92         (tagbody
  93           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul 2 j) 1))
  94           (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 11)))))
  95           (setf fval1 (funcall f (- centr absc)))
  96           (setf fval2 (funcall f (+ centr absc)))
  97           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 10))) fval1)
  98           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 10))) fval2)
  99           (setf fsum (+ fval1 fval2))
 100           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 11))) fsum)))
 101           (setf resabs
 102                   (+ resabs
 103                      (* (f2cl-lib:fref wgk (jtwm1) ((1 11)))
 104                         (+ (abs fval1) (abs fval2)))))
 105          label15))
 106       (setf reskh (* resk 0.5))
 107       (setf resasc (* (f2cl-lib:fref wgk (11) ((1 11))) (abs (- fc reskh))))
 108       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 109                     ((> j 10) nil)
 110         (tagbody
 111           (setf resasc
 112                   (+ resasc
 113                      (* (f2cl-lib:fref wgk (j) ((1 11)))
 114                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 10))) reskh))
 115                            (abs (- (f2cl-lib:fref fv2 (j) ((1 10))) reskh))))))
 116          label20))
 117       (setf result (* resk hlgth))
 118       (setf resabs (* resabs dhlgth))
 119       (setf resasc (* resasc dhlgth))
 120       (setf abserr (abs (* (- resk resg) hlgth)))
 121       (if (and (/= resasc 0.0) (/= abserr 0.0))
 122           (setf abserr
 123                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 124       (if (> resabs (/ uflow (* 50.0 epmach)))
 125           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 126       (go end_label)
 127      end_label
 128       (return (values nil nil nil result abserr resabs resasc)))))
 129
 130 (in-package #:cl-user)
 131 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 132 (eval-when (:load-toplevel :compile-toplevel :execute)
 133   (setf (gethash 'fortran-to-lisp::dqk21 fortran-to-lisp::*f2cl-function-info*)
 134           (fortran-to-lisp::make-f2cl-finfo
 135            :arg-types '(t (double-float) (double-float) (double-float)
 136                         (double-float) (double-float) (double-float))
 137            :return-values '(nil nil nil fortran-to-lisp::result
 138                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 139                             fortran-to-lisp::resasc)
 140            :calls '(fortran-to-lisp::d1mach))))
 141