src/numerical/slatec/dqawse.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 (defun dqawse
  21        (f a b alfa beta integr epsabs epsrel limit result abserr neval ier
  22         alist blist rlist elist iord last$)
  23   (declare (type (array f2cl-lib:integer4 (*)) iord)
  24            (type (array double-float (*)) elist rlist blist alist)
  25            (type (f2cl-lib:integer4) last$ ier neval limit integr)
  26            (type (double-float) abserr result epsrel epsabs beta alfa b a))
  27   (f2cl-lib:with-multi-array-data
  28       ((alist double-float alist-%data% alist-%offset%)
  29        (blist double-float blist-%data% blist-%offset%)
  30        (rlist double-float rlist-%data% rlist-%offset%)
  31        (elist double-float elist-%data% elist-%offset%)
  32        (iord f2cl-lib:integer4 iord-%data% iord-%offset%))
  33     (prog ((ri (make-array 25 :element-type 'double-float))
  34            (rj (make-array 25 :element-type 'double-float))
  35            (rh (make-array 25 :element-type 'double-float))
  36            (rg (make-array 25 :element-type 'double-float)) (iroff1 0)
  37            (iroff2 0) (k 0) (maxerr 0) (nev 0) (nrmax 0) (area 0.0) (area1 0.0)
  38            (area12 0.0) (area2 0.0) (a1 0.0) (a2 0.0) (b1 0.0) (b2 0.0)
  39            (centre 0.0) (epmach 0.0) (errbnd 0.0) (errmax 0.0) (error1 0.0)
  40            (erro12 0.0) (error2 0.0) (errsum 0.0) (resas1 0.0) (resas2 0.0)
  41            (uflow 0.0))
  42       (declare (type (array double-float (25)) rj ri rh rg)
  43                (type (double-float) uflow resas2 resas1 errsum error2 erro12
  44                                     error1 errmax errbnd epmach centre b2 b1 a2
  45                                     a1 area2 area12 area1 area)
  46                (type (f2cl-lib:integer4) nrmax nev maxerr k iroff2 iroff1))
  47       (setf epmach (f2cl-lib:d1mach 4))
  48       (setf uflow (f2cl-lib:d1mach 1))
  49       (setf ier 6)
  50       (setf neval 0)
  51       (setf last$ 0)
  52       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) 0.0)
  53       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) 0.0)
  54       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 0)
  55       (setf result 0.0)
  56       (setf abserr 0.0)
  57       (if
  58        (or (<= b a)
  59            (and (= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
  60            (<= alfa -1.0)
  61            (<= beta -1.0)
  62            (< integr 1)
  63            (> integr 4)
  64            (< limit 2))
  65        (go label999))
  66       (setf ier 0)
  67       (dqmomo alfa beta ri rj rg rh integr)
  68       (setf centre (* 0.5 (+ b a)))
  69       (multiple-value-bind
  70             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  71              var-11 var-12 var-13 var-14 var-15)
  72           (dqc25s f a b a centre alfa beta ri rj rg rh area1 error1 resas1
  73            integr nev)
  74         (declare (ignore var-0 var-3 var-4 var-7 var-8 var-9 var-10))
  75         (setf a var-1)
  76         (setf b var-2)
  77         (setf alfa var-5)
  78         (setf beta var-6)
  79         (setf area1 var-11)
  80         (setf error1 var-12)
  81         (setf resas1 var-13)
  82         (setf integr var-14)
  83         (setf nev var-15))
  84       (setf neval nev)
  85       (multiple-value-bind
  86             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  87              var-11 var-12 var-13 var-14 var-15)
  88           (dqc25s f a b centre b alfa beta ri rj rg rh area2 error2 resas2
  89            integr nev)
  90         (declare (ignore var-0 var-3 var-4 var-7 var-8 var-9 var-10))
  91         (setf a var-1)
  92         (setf b var-2)
  93         (setf alfa var-5)
  94         (setf beta var-6)
  95         (setf area2 var-11)
  96         (setf error2 var-12)
  97         (setf resas2 var-13)
  98         (setf integr var-14)
  99         (setf nev var-15))
 100       (setf last$ 2)
 101       (setf neval (f2cl-lib:int-add neval nev))
 102       (setf result (+ area1 area2))
 103       (setf abserr (+ error1 error2))
 104       (setf errbnd (max epsabs (* epsrel (abs result))))
 105       (if (> error2 error1) (go label10))
 106       (setf (f2cl-lib:fref alist-%data% (1) ((1 *)) alist-%offset%) a)
 107       (setf (f2cl-lib:fref alist-%data% (2) ((1 *)) alist-%offset%) centre)
 108       (setf (f2cl-lib:fref blist-%data% (1) ((1 *)) blist-%offset%) centre)
 109       (setf (f2cl-lib:fref blist-%data% (2) ((1 *)) blist-%offset%) b)
 110       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) area1)
 111       (setf (f2cl-lib:fref rlist-%data% (2) ((1 *)) rlist-%offset%) area2)
 112       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) error1)
 113       (setf (f2cl-lib:fref elist-%data% (2) ((1 *)) elist-%offset%) error2)
 114       (go label20)
 115      label10
 116       (setf (f2cl-lib:fref alist-%data% (1) ((1 *)) alist-%offset%) centre)
 117       (setf (f2cl-lib:fref alist-%data% (2) ((1 *)) alist-%offset%) a)
 118       (setf (f2cl-lib:fref blist-%data% (1) ((1 *)) blist-%offset%) b)
 119       (setf (f2cl-lib:fref blist-%data% (2) ((1 *)) blist-%offset%) centre)
 120       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) area2)
 121       (setf (f2cl-lib:fref rlist-%data% (2) ((1 *)) rlist-%offset%) area1)
 122       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) error2)
 123       (setf (f2cl-lib:fref elist-%data% (2) ((1 *)) elist-%offset%) error1)
 124      label20
 125       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 1)
 126       (setf (f2cl-lib:fref iord-%data% (2) ((1 *)) iord-%offset%) 2)
 127       (if (= limit 2) (setf ier 1))
 128       (if (or (<= abserr errbnd) (= ier 1)) (go label999))
 129       (setf errmax (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%))
 130       (setf maxerr 1)
 131       (setf nrmax 1)
 132       (setf area result)
 133       (setf errsum abserr)
 134       (setf iroff1 0)
 135       (setf iroff2 0)
 136       (f2cl-lib:fdo (last$ 3 (f2cl-lib:int-add last$ 1))
 137                     ((> last$ limit) nil)
 138         (tagbody
 139           (setf a1
 140                   (f2cl-lib:fref alist-%data% (maxerr) ((1 *)) alist-%offset%))
 141           (setf b1
 142                   (* 0.5
 143                      (+
 144                       (f2cl-lib:fref alist-%data%
 145                                      (maxerr)
 146                                      ((1 *))
 147                                      alist-%offset%)
 148                       (f2cl-lib:fref blist-%data%
 149                                      (maxerr)
 150                                      ((1 *))
 151                                      blist-%offset%))))
 152           (setf a2 b1)
 153           (setf b2
 154                   (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%))
 155           (multiple-value-bind
 156                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 157                  var-10 var-11 var-12 var-13 var-14 var-15)
 158               (dqc25s f a b a1 b1 alfa beta ri rj rg rh area1 error1 resas1
 159                integr nev)
 160             (declare (ignore var-0 var-3 var-4 var-7 var-8 var-9 var-10))
 161             (setf a var-1)
 162             (setf b var-2)
 163             (setf alfa var-5)
 164             (setf beta var-6)
 165             (setf area1 var-11)
 166             (setf error1 var-12)
 167             (setf resas1 var-13)
 168             (setf integr var-14)
 169             (setf nev var-15))
 170           (setf neval (f2cl-lib:int-add neval nev))
 171           (multiple-value-bind
 172                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 173                  var-10 var-11 var-12 var-13 var-14 var-15)
 174               (dqc25s f a b a2 b2 alfa beta ri rj rg rh area2 error2 resas2
 175                integr nev)
 176             (declare (ignore var-0 var-3 var-4 var-7 var-8 var-9 var-10))
 177             (setf a var-1)
 178             (setf b var-2)
 179             (setf alfa var-5)
 180             (setf beta var-6)
 181             (setf area2 var-11)
 182             (setf error2 var-12)
 183             (setf resas2 var-13)
 184             (setf integr var-14)
 185             (setf nev var-15))
 186           (setf neval (f2cl-lib:int-add neval nev))
 187           (setf area12 (+ area1 area2))
 188           (setf erro12 (+ error1 error2))
 189           (setf errsum (- (+ errsum erro12) errmax))
 190           (setf area
 191                   (- (+ area area12)
 192                      (f2cl-lib:fref rlist-%data%
 193                                     (maxerr)
 194                                     ((1 *))
 195                                     rlist-%offset%)))
 196           (if (or (= a a1) (= b b2)) (go label30))
 197           (if (or (= resas1 error1) (= resas2 error2)) (go label30))
 198           (if
 199            (and
 200             (<
 201              (abs
 202               (- (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 203                  area12))
 204              (* 1.0e-5 (abs area12)))
 205             (>= erro12 (* 0.99 errmax)))
 206            (setf iroff1 (f2cl-lib:int-add iroff1 1)))
 207           (if (and (> last$ 10) (> erro12 errmax))
 208               (setf iroff2 (f2cl-lib:int-add iroff2 1)))
 209          label30
 210           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 211                   area1)
 212           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 213                   area2)
 214           (setf errbnd (max epsabs (* epsrel (abs area))))
 215           (if (<= errsum errbnd) (go label35))
 216           (if (= last$ limit) (setf ier 1))
 217           (if (or (>= iroff1 6) (>= iroff2 20)) (setf ier 2))
 218           (if
 219            (<= (max (abs a1) (abs b2))
 220                (* (+ 1.0 (* 100.0 epmach)) (+ (abs a2) (* 1000.0 uflow))))
 221            (setf ier 3))
 222          label35
 223           (if (> error2 error1) (go label40))
 224           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a2)
 225           (setf (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%)
 226                   b1)
 227           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b2)
 228           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 229                   error1)
 230           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 231                   error2)
 232           (go label50)
 233          label40
 234           (setf (f2cl-lib:fref alist-%data% (maxerr) ((1 *)) alist-%offset%)
 235                   a2)
 236           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a1)
 237           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b1)
 238           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 239                   area2)
 240           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 241                   area1)
 242           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 243                   error2)
 244           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 245                   error1)
 246          label50
 247           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 248               (dqpsrt limit last$ maxerr errmax elist iord nrmax)
 249             (declare (ignore var-0 var-1 var-4 var-5))
 250             (setf maxerr var-2)
 251             (setf errmax var-3)
 252             (setf nrmax var-6))
 253           (if (or (/= ier 0) (<= errsum errbnd)) (go label70))
 254          label60))
 255      label70
 256       (setf result 0.0)
 257       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 258                     ((> k last$) nil)
 259         (tagbody
 260           (setf result
 261                   (+ result
 262                      (f2cl-lib:fref rlist-%data% (k) ((1 *)) rlist-%offset%)))
 263          label80))
 264       (setf abserr errsum)
 265      label999
 266       (go end_label)
 267      end_label
 268       (return
 269        (values nil
 270                a
 271                b
 272                alfa
 273                beta
 274                integr
 275                nil
 276                nil
 277                nil
 278                result
 279                abserr
 280                neval
 281                ier
 282                nil
 283                nil
 284                nil
 285                nil
 286                nil
 287                last$)))))
 288
 289 (in-package #:cl-user)
 290 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 291 (eval-when (:load-toplevel :compile-toplevel :execute)
 292   (setf (gethash 'fortran-to-lisp::dqawse
 293                  fortran-to-lisp::*f2cl-function-info*)
 294           (fortran-to-lisp::make-f2cl-finfo
 295            :arg-types '(t (double-float) (double-float) (double-float)
 296                         (double-float) (fortran-to-lisp::integer4)
 297                         (double-float) (double-float)
 298                         (fortran-to-lisp::integer4) (double-float)
 299                         (double-float) (fortran-to-lisp::integer4)
 300                         (fortran-to-lisp::integer4) (array double-float (*))
 301                         (array double-float (*)) (array double-float (*))
 302                         (array double-float (*))
 303                         (array fortran-to-lisp::integer4 (*))
 304                         (fortran-to-lisp::integer4))
 305            :return-values '(nil fortran-to-lisp::a fortran-to-lisp::b
 306                             fortran-to-lisp::alfa fortran-to-lisp::beta
 307                             fortran-to-lisp::integr nil nil nil
 308                             fortran-to-lisp::result fortran-to-lisp::abserr
 309                             fortran-to-lisp::neval fortran-to-lisp::ier nil nil
 310                             nil nil nil fortran-to-lisp::last$)
 311            :calls '(fortran-to-lisp::dqpsrt fortran-to-lisp::dqc25s
 312                     fortran-to-lisp::dqmomo fortran-to-lisp::d1mach))))
 313