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