src/numerical/slatec/dqage.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 dqage
  21        (f a b epsabs epsrel key limit result abserr neval ier alist blist rlist
  22         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 key)
  26            (type (double-float) abserr result epsrel epsabs 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 ((iroff1 0) (iroff2 0) (k 0) (keyf 0) (maxerr 0) (nrmax 0) (area 0.0)
  34            (area1 0.0) (area12 0.0) (area2 0.0) (a1 0.0) (a2 0.0) (b1 0.0)
  35            (b2 0.0) (defabs 0.0) (defab1 0.0) (defab2 0.0) (epmach 0.0)
  36            (errbnd 0.0) (errmax 0.0) (error1 0.0) (error2 0.0) (erro12 0.0)
  37            (errsum 0.0) (resabs 0.0) (uflow 0.0))
  38       (declare (type (double-float) uflow resabs errsum erro12 error2 error1
  39                                     errmax errbnd epmach defab2 defab1 defabs
  40                                     b2 b1 a2 a1 area2 area12 area1 area)
  41                (type (f2cl-lib:integer4) nrmax maxerr keyf k iroff2 iroff1))
  42       (setf epmach (f2cl-lib:d1mach 4))
  43       (setf uflow (f2cl-lib:d1mach 1))
  44       (setf ier 0)
  45       (setf neval 0)
  46       (setf last$ 0)
  47       (setf result 0.0)
  48       (setf abserr 0.0)
  49       (setf (f2cl-lib:fref alist-%data% (1) ((1 *)) alist-%offset%) a)
  50       (setf (f2cl-lib:fref blist-%data% (1) ((1 *)) blist-%offset%) b)
  51       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) 0.0)
  52       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) 0.0)
  53       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 0)
  54       (if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
  55           (setf ier 6))
  56       (if (= ier 6) (go label999))
  57       (setf keyf key)
  58       (if (<= key 0) (setf keyf 1))
  59       (if (>= key 7) (setf keyf 6))
  60       (setf neval 0)
  61       (if (= keyf 1)
  62           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
  63               (dqk15 f a b result abserr defabs resabs)
  64             (declare (ignore var-0 var-1 var-2))
  65             (setf result var-3)
  66             (setf abserr var-4)
  67             (setf defabs var-5)
  68             (setf resabs var-6)))
  69       (if (= keyf 2)
  70           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
  71               (dqk21 f a b result abserr defabs resabs)
  72             (declare (ignore var-0 var-1 var-2))
  73             (setf result var-3)
  74             (setf abserr var-4)
  75             (setf defabs var-5)
  76             (setf resabs var-6)))
  77       (if (= keyf 3)
  78           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
  79               (dqk31 f a b result abserr defabs resabs)
  80             (declare (ignore var-0 var-1 var-2))
  81             (setf result var-3)
  82             (setf abserr var-4)
  83             (setf defabs var-5)
  84             (setf resabs var-6)))
  85       (if (= keyf 4)
  86           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
  87               (dqk41 f a b result abserr defabs resabs)
  88             (declare (ignore var-0 var-1 var-2))
  89             (setf result var-3)
  90             (setf abserr var-4)
  91             (setf defabs var-5)
  92             (setf resabs var-6)))
  93       (if (= keyf 5)
  94           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
  95               (dqk51 f a b result abserr defabs resabs)
  96             (declare (ignore var-0 var-1 var-2))
  97             (setf result var-3)
  98             (setf abserr var-4)
  99             (setf defabs var-5)
 100             (setf resabs var-6)))
 101       (if (= keyf 6)
 102           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 103               (dqk61 f a b result abserr defabs resabs)
 104             (declare (ignore var-0 var-1 var-2))
 105             (setf result var-3)
 106             (setf abserr var-4)
 107             (setf defabs var-5)
 108             (setf resabs var-6)))
 109       (setf last$ 1)
 110       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) result)
 111       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) abserr)
 112       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 1)
 113       (setf errbnd (max epsabs (* epsrel (abs result))))
 114       (if (and (<= abserr (* 50.0 epmach defabs)) (> abserr errbnd))
 115           (setf ier 2))
 116       (if (= limit 1) (setf ier 1))
 117       (if
 118        (or (/= ier 0)
 119            (and (<= abserr errbnd) (/= abserr resabs))
 120            (= abserr 0.0))
 121        (go label60))
 122       (setf errmax abserr)
 123       (setf maxerr 1)
 124       (setf area result)
 125       (setf errsum abserr)
 126       (setf nrmax 1)
 127       (setf iroff1 0)
 128       (setf iroff2 0)
 129       (f2cl-lib:fdo (last$ 2 (f2cl-lib:int-add last$ 1))
 130                     ((> last$ limit) nil)
 131         (tagbody
 132           (setf a1
 133                   (f2cl-lib:fref alist-%data% (maxerr) ((1 *)) alist-%offset%))
 134           (setf b1
 135                   (* 0.5
 136                      (+
 137                       (f2cl-lib:fref alist-%data%
 138                                      (maxerr)
 139                                      ((1 *))
 140                                      alist-%offset%)
 141                       (f2cl-lib:fref blist-%data%
 142                                      (maxerr)
 143                                      ((1 *))
 144                                      blist-%offset%))))
 145           (setf a2 b1)
 146           (setf b2
 147                   (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%))
 148           (if (= keyf 1)
 149               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 150                   (dqk15 f a1 b1 area1 error1 resabs defab1)
 151                 (declare (ignore var-0 var-1 var-2))
 152                 (setf area1 var-3)
 153                 (setf error1 var-4)
 154                 (setf resabs var-5)
 155                 (setf defab1 var-6)))
 156           (if (= keyf 2)
 157               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 158                   (dqk21 f a1 b1 area1 error1 resabs defab1)
 159                 (declare (ignore var-0 var-1 var-2))
 160                 (setf area1 var-3)
 161                 (setf error1 var-4)
 162                 (setf resabs var-5)
 163                 (setf defab1 var-6)))
 164           (if (= keyf 3)
 165               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 166                   (dqk31 f a1 b1 area1 error1 resabs defab1)
 167                 (declare (ignore var-0 var-1 var-2))
 168                 (setf area1 var-3)
 169                 (setf error1 var-4)
 170                 (setf resabs var-5)
 171                 (setf defab1 var-6)))
 172           (if (= keyf 4)
 173               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 174                   (dqk41 f a1 b1 area1 error1 resabs defab1)
 175                 (declare (ignore var-0 var-1 var-2))
 176                 (setf area1 var-3)
 177                 (setf error1 var-4)
 178                 (setf resabs var-5)
 179                 (setf defab1 var-6)))
 180           (if (= keyf 5)
 181               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 182                   (dqk51 f a1 b1 area1 error1 resabs defab1)
 183                 (declare (ignore var-0 var-1 var-2))
 184                 (setf area1 var-3)
 185                 (setf error1 var-4)
 186                 (setf resabs var-5)
 187                 (setf defab1 var-6)))
 188           (if (= keyf 6)
 189               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 190                   (dqk61 f a1 b1 area1 error1 resabs defab1)
 191                 (declare (ignore var-0 var-1 var-2))
 192                 (setf area1 var-3)
 193                 (setf error1 var-4)
 194                 (setf resabs var-5)
 195                 (setf defab1 var-6)))
 196           (if (= keyf 1)
 197               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 198                   (dqk15 f a2 b2 area2 error2 resabs defab2)
 199                 (declare (ignore var-0 var-1 var-2))
 200                 (setf area2 var-3)
 201                 (setf error2 var-4)
 202                 (setf resabs var-5)
 203                 (setf defab2 var-6)))
 204           (if (= keyf 2)
 205               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 206                   (dqk21 f a2 b2 area2 error2 resabs defab2)
 207                 (declare (ignore var-0 var-1 var-2))
 208                 (setf area2 var-3)
 209                 (setf error2 var-4)
 210                 (setf resabs var-5)
 211                 (setf defab2 var-6)))
 212           (if (= keyf 3)
 213               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 214                   (dqk31 f a2 b2 area2 error2 resabs defab2)
 215                 (declare (ignore var-0 var-1 var-2))
 216                 (setf area2 var-3)
 217                 (setf error2 var-4)
 218                 (setf resabs var-5)
 219                 (setf defab2 var-6)))
 220           (if (= keyf 4)
 221               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 222                   (dqk41 f a2 b2 area2 error2 resabs defab2)
 223                 (declare (ignore var-0 var-1 var-2))
 224                 (setf area2 var-3)
 225                 (setf error2 var-4)
 226                 (setf resabs var-5)
 227                 (setf defab2 var-6)))
 228           (if (= keyf 5)
 229               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 230                   (dqk51 f a2 b2 area2 error2 resabs defab2)
 231                 (declare (ignore var-0 var-1 var-2))
 232                 (setf area2 var-3)
 233                 (setf error2 var-4)
 234                 (setf resabs var-5)
 235                 (setf defab2 var-6)))
 236           (if (= keyf 6)
 237               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 238                   (dqk61 f a2 b2 area2 error2 resabs defab2)
 239                 (declare (ignore var-0 var-1 var-2))
 240                 (setf area2 var-3)
 241                 (setf error2 var-4)
 242                 (setf resabs var-5)
 243                 (setf defab2 var-6)))
 244           (setf neval (f2cl-lib:int-add neval 1))
 245           (setf area12 (+ area1 area2))
 246           (setf erro12 (+ error1 error2))
 247           (setf errsum (- (+ errsum erro12) errmax))
 248           (setf area
 249                   (- (+ area area12)
 250                      (f2cl-lib:fref rlist-%data%
 251                                     (maxerr)
 252                                     ((1 *))
 253                                     rlist-%offset%)))
 254           (if (or (= defab1 error1) (= defab2 error2)) (go label5))
 255           (if
 256            (and
 257             (<=
 258              (abs
 259               (- (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 260                  area12))
 261              (* 1.0e-5 (abs area12)))
 262             (>= erro12 (* 0.99 errmax)))
 263            (setf iroff1 (f2cl-lib:int-add iroff1 1)))
 264           (if (and (> last$ 10) (> erro12 errmax))
 265               (setf iroff2 (f2cl-lib:int-add iroff2 1)))
 266          label5
 267           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 268                   area1)
 269           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 270                   area2)
 271           (setf errbnd (max epsabs (* epsrel (abs area))))
 272           (if (<= errsum errbnd) (go label8))
 273           (if (or (>= iroff1 6) (>= iroff2 20)) (setf ier 2))
 274           (if (= last$ limit) (setf ier 1))
 275           (if
 276            (<= (max (abs a1) (abs b2))
 277                (* (+ 1.0 (* 100.0 epmach)) (+ (abs a2) (* 1000.0 uflow))))
 278            (setf ier 3))
 279          label8
 280           (if (> error2 error1) (go label10))
 281           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a2)
 282           (setf (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%)
 283                   b1)
 284           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b2)
 285           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 286                   error1)
 287           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 288                   error2)
 289           (go label20)
 290          label10
 291           (setf (f2cl-lib:fref alist-%data% (maxerr) ((1 *)) alist-%offset%)
 292                   a2)
 293           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a1)
 294           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b1)
 295           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 296                   area2)
 297           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 298                   area1)
 299           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 300                   error2)
 301           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 302                   error1)
 303          label20
 304           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 305               (dqpsrt limit last$ maxerr errmax elist iord nrmax)
 306             (declare (ignore var-0 var-1 var-4 var-5))
 307             (setf maxerr var-2)
 308             (setf errmax var-3)
 309             (setf nrmax var-6))
 310           (if (or (/= ier 0) (<= errsum errbnd)) (go label40))
 311          label30))
 312      label40
 313       (setf result 0.0)
 314       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 315                     ((> k last$) nil)
 316         (tagbody
 317           (setf result
 318                   (+ result
 319                      (f2cl-lib:fref rlist-%data% (k) ((1 *)) rlist-%offset%)))
 320          label50))
 321       (setf abserr errsum)
 322      label60
 323       (if (/= keyf 1)
 324           (setf neval
 325                   (f2cl-lib:int-mul
 326                    (f2cl-lib:int-add (f2cl-lib:int-mul 10 keyf) 1)
 327                    (f2cl-lib:int-add (f2cl-lib:int-mul 2 neval) 1))))
 328       (if (= keyf 1)
 329           (setf neval (f2cl-lib:int-add (f2cl-lib:int-mul 30 neval) 15)))
 330      label999
 331       (go end_label)
 332      end_label
 333       (return
 334        (values nil
 335                nil
 336                nil
 337                nil
 338                nil
 339                nil
 340                nil
 341                result
 342                abserr
 343                neval
 344                ier
 345                nil
 346                nil
 347                nil
 348                nil
 349                nil
 350                last$)))))
 351
 352 (in-package #:cl-user)
 353 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 354 (eval-when (:load-toplevel :compile-toplevel :execute)
 355   (setf (gethash 'fortran-to-lisp::dqage fortran-to-lisp::*f2cl-function-info*)
 356           (fortran-to-lisp::make-f2cl-finfo
 357            :arg-types '(t (double-float) (double-float) (double-float)
 358                         (double-float) (fortran-to-lisp::integer4)
 359                         (fortran-to-lisp::integer4) (double-float)
 360                         (double-float) (fortran-to-lisp::integer4)
 361                         (fortran-to-lisp::integer4) (array double-float (*))
 362                         (array double-float (*)) (array double-float (*))
 363                         (array double-float (*))
 364                         (array fortran-to-lisp::integer4 (*))
 365                         (fortran-to-lisp::integer4))
 366            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::result
 367                             fortran-to-lisp::abserr fortran-to-lisp::neval
 368                             fortran-to-lisp::ier nil nil nil nil nil
 369                             fortran-to-lisp::last$)
 370            :calls '(fortran-to-lisp::dqpsrt fortran-to-lisp::dqk61
 371                     fortran-to-lisp::dqk51 fortran-to-lisp::dqk41
 372                     fortran-to-lisp::dqk31 fortran-to-lisp::dqk21
 373                     fortran-to-lisp::dqk15 fortran-to-lisp::d1mach))))
 374