src/numerical/slatec/dqawfe.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 nil) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((p 0.9) (pi$ 3.141592653589793))
  21   (declare (type (double-float) p pi$))
  22   (defun dqawfe
  23          (f a omega integr epsabs limlst limit maxp1 result abserr neval ier
  24           rslst erlst ierlst lst alist blist rlist elist iord nnlog chebmo)
  25     (declare (type (array f2cl-lib:integer4 (*)) nnlog iord ierlst)
  26              (type (array double-float (*)) chebmo elist rlist blist alist
  27                                             erlst rslst)
  28              (type (f2cl-lib:integer4) lst ier neval maxp1 limit limlst integr)
  29              (type (double-float) abserr result epsabs omega a))
  30     (f2cl-lib:with-multi-array-data
  31         ((rslst double-float rslst-%data% rslst-%offset%)
  32          (erlst double-float erlst-%data% erlst-%offset%)
  33          (alist double-float alist-%data% alist-%offset%)
  34          (blist double-float blist-%data% blist-%offset%)
  35          (rlist double-float rlist-%data% rlist-%offset%)
  36          (elist double-float elist-%data% elist-%offset%)
  37          (chebmo double-float chebmo-%data% chebmo-%offset%)
  38          (ierlst f2cl-lib:integer4 ierlst-%data% ierlst-%offset%)
  39          (iord f2cl-lib:integer4 iord-%data% iord-%offset%)
  40          (nnlog f2cl-lib:integer4 nnlog-%data% nnlog-%offset%))
  41       (prog ((psum (make-array 52 :element-type 'double-float))
  42              (res3la (make-array 3 :element-type 'double-float)) (ktmin 0)
  43              (l 0) (ll 0) (momcom 0) (nev 0) (nres 0) (numrl2 0) (abseps 0.0)
  44              (correc 0.0) (cycle 0.0) (c1 0.0) (c2 0.0) (dl 0.0) (drl 0.0)
  45              (ep 0.0) (eps 0.0) (epsa 0.0) (errsum 0.0) (fact 0.0) (p1 0.0)
  46              (reseps 0.0) (uflow 0.0) (last$ 0))
  47         (declare (type (array double-float (3)) res3la)
  48                  (type (array double-float (52)) psum)
  49                  (type (double-float) uflow reseps p1 fact errsum epsa eps ep
  50                                       drl dl c2 c1 cycle correc abseps)
  51                  (type (f2cl-lib:integer4) last$ numrl2 nres nev momcom ll l
  52                                            ktmin))
  53         (setf result 0.0)
  54         (setf abserr 0.0)
  55         (setf neval 0)
  56         (setf lst 0)
  57         (setf ier 0)
  58         (if (or (and (/= integr 1) (/= integr 2)) (<= epsabs 0.0) (< limlst 3))
  59             (setf ier 6))
  60         (if (= ier 6) (go label999))
  61         (if (/= omega 0.0) (go label10))
  62         (if (= integr 1)
  63             (multiple-value-bind
  64                   (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  65                    var-10 var-11 var-12 var-13 var-14 var-15)
  66                 (dqagie f a 1 epsabs 0.0 limit result abserr neval ier alist
  67                  blist rlist elist iord last$)
  68               (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-10
  69                                var-11 var-12 var-13 var-14))
  70               (setf result var-6)
  71               (setf abserr var-7)
  72               (setf neval var-8)
  73               (setf ier var-9)
  74               (setf last$ var-15)))
  75         (setf (f2cl-lib:fref rslst-%data% (1) ((1 *)) rslst-%offset%) result)
  76         (setf (f2cl-lib:fref erlst-%data% (1) ((1 *)) erlst-%offset%) abserr)
  77         (setf (f2cl-lib:fref ierlst-%data% (1) ((1 *)) ierlst-%offset%) ier)
  78         (setf lst 1)
  79         (go label999)
  80        label10
  81         (setf l (f2cl-lib:int (abs omega)))
  82         (setf dl
  83                 (coerce
  84                  (the f2cl-lib:integer4
  85                       (f2cl-lib:int-add (f2cl-lib:int-mul 2 l) 1))
  86                  'double-float))
  87         (setf cycle (/ (* dl pi$) (abs omega)))
  88         (setf ier 0)
  89         (setf ktmin 0)
  90         (setf neval 0)
  91         (setf numrl2 0)
  92         (setf nres 0)
  93         (setf c1 a)
  94         (setf c2 (+ cycle a))
  95         (setf p1 (- 1.0 p))
  96         (setf uflow (f2cl-lib:d1mach 1))
  97         (setf eps epsabs)
  98         (if (> epsabs (/ uflow p1)) (setf eps (* epsabs p1)))
  99         (setf ep eps)
 100         (setf fact 1.0)
 101         (setf correc 0.0)
 102         (setf abserr 0.0)
 103         (setf errsum 0.0)
 104         (f2cl-lib:fdo (lst 1 (f2cl-lib:int-add lst 1))
 105                       ((> lst limlst) nil)
 106           (tagbody
 107             (setf epsa (* eps fact))
 108             (multiple-value-bind
 109                   (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 110                    var-10 var-11 var-12 var-13 var-14 var-15 var-16 var-17
 111                    var-18 var-19 var-20 var-21 var-22)
 112                 (dqawoe f c1 c2 omega integr epsa 0.0 limit lst maxp1
 113                  (f2cl-lib:fref rslst-%data% (lst) ((1 *)) rslst-%offset%)
 114                  (f2cl-lib:fref erlst-%data% (lst) ((1 *)) erlst-%offset%) nev
 115                  (f2cl-lib:fref ierlst-%data% (lst) ((1 *)) ierlst-%offset%)
 116                  last$ alist blist rlist elist iord nnlog momcom chebmo)
 117               (declare (ignore var-0 var-1 var-2 var-3 var-5 var-6 var-7 var-8
 118                                var-9 var-15 var-16 var-17 var-18 var-19 var-20
 119                                var-22))
 120               (setf integr var-4)
 121               (setf (f2cl-lib:fref rslst-%data% (lst) ((1 *)) rslst-%offset%)
 122                       var-10)
 123               (setf (f2cl-lib:fref erlst-%data% (lst) ((1 *)) erlst-%offset%)
 124                       var-11)
 125               (setf nev var-12)
 126               (setf (f2cl-lib:fref ierlst-%data% (lst) ((1 *)) ierlst-%offset%)
 127                       var-13)
 128               (setf last$ var-14)
 129               (setf momcom var-21))
 130             (setf neval (f2cl-lib:int-add neval nev))
 131             (setf fact (* fact p))
 132             (setf errsum
 133                     (+ errsum
 134                        (f2cl-lib:fref erlst-%data%
 135                                       (lst)
 136                                       ((1 *))
 137                                       erlst-%offset%)))
 138             (setf drl
 139                     (* 50.0
 140                        (abs
 141                         (f2cl-lib:fref rslst-%data%
 142                                        (lst)
 143                                        ((1 *))
 144                                        rslst-%offset%))))
 145             (if (and (<= (+ errsum drl) epsabs) (>= lst 6)) (go label80))
 146             (setf correc
 147                     (max correc
 148                          (f2cl-lib:fref erlst-%data%
 149                                         (lst)
 150                                         ((1 *))
 151                                         erlst-%offset%)))
 152             (if
 153              (/= (f2cl-lib:fref ierlst-%data% (lst) ((1 *)) ierlst-%offset%) 0)
 154              (setf eps (max ep (* correc p1))))
 155             (if
 156              (/= (f2cl-lib:fref ierlst-%data% (lst) ((1 *)) ierlst-%offset%) 0)
 157              (setf ier 7))
 158             (if (and (= ier 7) (<= (+ errsum drl) (* correc 10.0)) (> lst 5))
 159                 (go label80))
 160             (setf numrl2 (f2cl-lib:int-add numrl2 1))
 161             (if (> lst 1) (go label20))
 162             (setf (f2cl-lib:fref psum (1) ((1 52)))
 163                     (f2cl-lib:fref rslst-%data% (1) ((1 *)) rslst-%offset%))
 164             (go label40)
 165            label20
 166             (setf (f2cl-lib:fref psum (numrl2) ((1 52)))
 167                     (+ (f2cl-lib:fref psum (ll) ((1 52)))
 168                        (f2cl-lib:fref rslst-%data%
 169                                       (lst)
 170                                       ((1 *))
 171                                       rslst-%offset%)))
 172             (if (= lst 2) (go label40))
 173             (if (= lst limlst) (setf ier 1))
 174             (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 175                 (dqelg numrl2 psum reseps abseps res3la nres)
 176               (declare (ignore var-1 var-4))
 177               (setf numrl2 var-0)
 178               (setf reseps var-2)
 179               (setf abseps var-3)
 180               (setf nres var-5))
 181             (setf ktmin (f2cl-lib:int-add ktmin 1))
 182             (if (and (>= ktmin 15) (<= abserr (* 0.001 (+ errsum drl))))
 183                 (setf ier 4))
 184             (if (and (> abseps abserr) (/= lst 3)) (go label30))
 185             (setf abserr abseps)
 186             (setf result reseps)
 187             (setf ktmin 0)
 188             (if
 189              (or (<= (+ abserr (* 10.0 correc)) epsabs)
 190                  (and (<= abserr epsabs) (>= (* 10.0 correc) epsabs)))
 191              (go label60))
 192            label30
 193             (if (and (/= ier 0) (/= ier 7)) (go label60))
 194            label40
 195             (setf ll numrl2)
 196             (setf c1 c2)
 197             (setf c2 (+ c2 cycle))
 198            label50))
 199        label60
 200         (setf abserr (+ abserr (* 10.0 correc)))
 201         (if (= ier 0) (go label999))
 202         (if
 203          (and (/= result 0.0) (/= (f2cl-lib:fref psum (numrl2) ((1 52))) 0.0))
 204          (go label70))
 205         (if (> abserr errsum) (go label80))
 206         (if (= (f2cl-lib:fref psum (numrl2) ((1 52))) 0.0) (go label999))
 207        label70
 208         (if
 209          (> (/ abserr (abs result))
 210             (/ (+ errsum drl) (abs (f2cl-lib:fref psum (numrl2) ((1 52))))))
 211          (go label80))
 212         (if (and (>= ier 1) (/= ier 7)) (setf abserr (+ abserr drl)))
 213         (go label999)
 214        label80
 215         (setf result (f2cl-lib:fref psum (numrl2) ((1 52))))
 216         (setf abserr (+ errsum drl))
 217        label999
 218         (go end_label)
 219        end_label
 220         (return
 221          (values nil
 222                  nil
 223                  nil
 224                  integr
 225                  nil
 226                  nil
 227                  nil
 228                  nil
 229                  result
 230                  abserr
 231                  neval
 232                  ier
 233                  nil
 234                  nil
 235                  nil
 236                  lst
 237                  nil
 238                  nil
 239                  nil
 240                  nil
 241                  nil
 242                  nil
 243                  nil))))))
 244
 245 (in-package #:cl-user)
 246 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 247 (eval-when (:load-toplevel :compile-toplevel :execute)
 248   (setf (gethash 'fortran-to-lisp::dqawfe
 249                  fortran-to-lisp::*f2cl-function-info*)
 250           (fortran-to-lisp::make-f2cl-finfo
 251            :arg-types '(t (double-float) (double-float)
 252                         (fortran-to-lisp::integer4) (double-float)
 253                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 254                         (fortran-to-lisp::integer4) (double-float)
 255                         (double-float) (fortran-to-lisp::integer4)
 256                         (fortran-to-lisp::integer4) (array double-float (*))
 257                         (array double-float (*))
 258                         (array fortran-to-lisp::integer4 (*))
 259                         (fortran-to-lisp::integer4) (array double-float (*))
 260                         (array double-float (*)) (array double-float (*))
 261                         (array double-float (*))
 262                         (array fortran-to-lisp::integer4 (*))
 263                         (array fortran-to-lisp::integer4 (*))
 264                         (array double-float (*)))
 265            :return-values '(nil nil nil fortran-to-lisp::integr nil nil nil nil
 266                             fortran-to-lisp::result fortran-to-lisp::abserr
 267                             fortran-to-lisp::neval fortran-to-lisp::ier nil nil
 268                             nil fortran-to-lisp::lst nil nil nil nil nil nil
 269                             nil)
 270            :calls '(fortran-to-lisp::dqelg fortran-to-lisp::dqawoe
 271                     fortran-to-lisp::d1mach fortran-to-lisp::dqagie))))
 272