src/numerical/slatec/dqawoe.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 dqawoe
  21        (f a b omega integr epsabs epsrel limit icall maxp1 result abserr neval
  22         ier last$ alist blist rlist elist iord nnlog momcom chebmo)
  23   (declare (type (array f2cl-lib:integer4 (*)) nnlog iord)
  24            (type (array double-float (*)) chebmo elist rlist blist alist)
  25            (type (f2cl-lib:integer4) momcom last$ ier neval maxp1 icall limit
  26                                      integr)
  27            (type (double-float) abserr result epsrel epsabs omega b a))
  28   (f2cl-lib:with-multi-array-data
  29       ((alist double-float alist-%data% alist-%offset%)
  30        (blist double-float blist-%data% blist-%offset%)
  31        (rlist double-float rlist-%data% rlist-%offset%)
  32        (elist double-float elist-%data% elist-%offset%)
  33        (chebmo double-float chebmo-%data% chebmo-%offset%)
  34        (iord f2cl-lib:integer4 iord-%data% iord-%offset%)
  35        (nnlog f2cl-lib:integer4 nnlog-%data% nnlog-%offset%))
  36     (prog ((rlist2 (make-array 52 :element-type 'double-float))
  37            (res3la (make-array 3 :element-type 'double-float)) (extrap nil)
  38            (noext nil) (extall nil) (id 0) (ierro 0) (iroff1 0) (iroff2 0)
  39            (iroff3 0) (jupbnd 0) (k 0) (ksgn 0) (ktmin 0) (maxerr 0) (nev 0)
  40            (nres 0) (nrmax 0) (nrmom 0) (numrl2 0) (abseps 0.0) (area 0.0)
  41            (area1 0.0) (area12 0.0) (area2 0.0) (a1 0.0) (a2 0.0) (b1 0.0)
  42            (b2 0.0) (correc 0.0) (defab1 0.0) (defab2 0.0) (defabs 0.0)
  43            (domega 0.0) (dres 0.0) (epmach 0.0) (erlarg 0.0) (erlast 0.0)
  44            (errbnd 0.0) (errmax 0.0) (error1 0.0) (erro12 0.0) (error2 0.0)
  45            (errsum 0.0) (ertest 0.0) (oflow 0.0) (resabs 0.0) (reseps 0.0)
  46            (small 0.0) (uflow 0.0) (width 0.0))
  47       (declare (type (array double-float (52)) rlist2)
  48                (type (array double-float (3)) res3la)
  49                (type (double-float) width uflow small reseps resabs oflow
  50                                     ertest errsum error2 erro12 error1 errmax
  51                                     errbnd erlast erlarg epmach dres domega
  52                                     defabs defab2 defab1 correc b2 b1 a2 a1
  53                                     area2 area12 area1 area abseps)
  54                (type (f2cl-lib:integer4) numrl2 nrmom nrmax nres nev maxerr
  55                                          ktmin ksgn k jupbnd iroff3 iroff2
  56                                          iroff1 ierro id)
  57                (type f2cl-lib:logical extall noext extrap))
  58       (setf epmach (f2cl-lib:d1mach 4))
  59       (setf ier 0)
  60       (setf neval 0)
  61       (setf last$ 0)
  62       (setf result 0.0)
  63       (setf abserr 0.0)
  64       (setf (f2cl-lib:fref alist-%data% (1) ((1 *)) alist-%offset%) a)
  65       (setf (f2cl-lib:fref blist-%data% (1) ((1 *)) blist-%offset%) b)
  66       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) 0.0)
  67       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) 0.0)
  68       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 0)
  69       (setf (f2cl-lib:fref nnlog-%data% (1) ((1 *)) nnlog-%offset%) 0)
  70       (if
  71        (or (and (/= integr 1) (/= integr 2))
  72            (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
  73            (< icall 1)
  74            (< maxp1 1))
  75        (setf ier 6))
  76       (if (= ier 6) (go label999))
  77       (setf domega (abs omega))
  78       (setf nrmom 0)
  79       (if (> icall 1) (go label5))
  80       (setf momcom 0)
  81      label5
  82       (multiple-value-bind
  83             (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  84              var-11 var-12 var-13 var-14)
  85           (dqc25f f a b domega integr nrmom maxp1 0 result abserr neval defabs
  86            resabs momcom chebmo)
  87         (declare (ignore var-0 var-1 var-2 var-5 var-6 var-7 var-14))
  88         (setf domega var-3)
  89         (setf integr var-4)
  90         (setf result var-8)
  91         (setf abserr var-9)
  92         (setf neval var-10)
  93         (setf defabs var-11)
  94         (setf resabs var-12)
  95         (setf momcom var-13))
  96       (setf dres (abs result))
  97       (setf errbnd (max epsabs (* epsrel dres)))
  98       (setf (f2cl-lib:fref rlist-%data% (1) ((1 *)) rlist-%offset%) result)
  99       (setf (f2cl-lib:fref elist-%data% (1) ((1 *)) elist-%offset%) abserr)
 100       (setf (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%) 1)
 101       (if (and (<= abserr (* 100.0 epmach defabs)) (> abserr errbnd))
 102           (setf ier 2))
 103       (if (= limit 1) (setf ier 1))
 104       (if (or (/= ier 0) (<= abserr errbnd)) (go label200))
 105       (setf uflow (f2cl-lib:d1mach 1))
 106       (setf oflow (f2cl-lib:d1mach 2))
 107       (setf errmax abserr)
 108       (setf maxerr 1)
 109       (setf area result)
 110       (setf errsum abserr)
 111       (setf abserr oflow)
 112       (setf nrmax 1)
 113       (setf extrap f2cl-lib:%false%)
 114       (setf noext f2cl-lib:%false%)
 115       (setf ierro 0)
 116       (setf iroff1 0)
 117       (setf iroff2 0)
 118       (setf iroff3 0)
 119       (setf ktmin 0)
 120       (setf small (* (abs (- b a)) 0.75))
 121       (setf nres 0)
 122       (setf numrl2 0)
 123       (setf extall f2cl-lib:%false%)
 124       (if (> (* 0.5 (abs (- b a)) domega) 2.0) (go label10))
 125       (setf numrl2 1)
 126       (setf extall f2cl-lib:%true%)
 127       (setf (f2cl-lib:fref rlist2 (1) ((1 52))) result)
 128      label10
 129       (if (<= (* 0.25 (abs (- b a)) domega) 2.0) (setf extall f2cl-lib:%true%))
 130       (setf ksgn -1)
 131       (if (>= dres (* (- 1.0 (* 50.0 epmach)) defabs)) (setf ksgn 1))
 132       (f2cl-lib:fdo (last$ 2 (f2cl-lib:int-add last$ 1))
 133                     ((> last$ limit) nil)
 134         (tagbody
 135           (setf nrmom
 136                   (f2cl-lib:int-add
 137                    (f2cl-lib:fref nnlog-%data% (maxerr) ((1 *)) nnlog-%offset%)
 138                    1))
 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           (setf erlast errmax)
 156           (multiple-value-bind
 157                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 158                  var-10 var-11 var-12 var-13 var-14)
 159               (dqc25f f a1 b1 domega integr nrmom maxp1 0 area1 error1 nev
 160                resabs defab1 momcom chebmo)
 161             (declare (ignore var-0 var-1 var-2 var-5 var-6 var-7 var-14))
 162             (setf domega var-3)
 163             (setf integr var-4)
 164             (setf area1 var-8)
 165             (setf error1 var-9)
 166             (setf nev var-10)
 167             (setf resabs var-11)
 168             (setf defab1 var-12)
 169             (setf momcom var-13))
 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)
 174               (dqc25f f a2 b2 domega integr nrmom maxp1 1 area2 error2 nev
 175                resabs defab2 momcom chebmo)
 176             (declare (ignore var-0 var-1 var-2 var-5 var-6 var-7 var-14))
 177             (setf domega var-3)
 178             (setf integr var-4)
 179             (setf area2 var-8)
 180             (setf error2 var-9)
 181             (setf nev var-10)
 182             (setf resabs var-11)
 183             (setf defab2 var-12)
 184             (setf momcom var-13))
 185           (setf neval (f2cl-lib:int-add neval nev))
 186           (setf area12 (+ area1 area2))
 187           (setf erro12 (+ error1 error2))
 188           (setf errsum (- (+ errsum erro12) errmax))
 189           (setf area
 190                   (- (+ area area12)
 191                      (f2cl-lib:fref rlist-%data%
 192                                     (maxerr)
 193                                     ((1 *))
 194                                     rlist-%offset%)))
 195           (if (or (= defab1 error1) (= defab2 error2)) (go label25))
 196           (if
 197            (or
 198             (>
 199              (abs
 200               (- (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 201                  area12))
 202              (* 1.0e-5 (abs area12)))
 203             (< erro12 (* 0.99 errmax)))
 204            (go label20))
 205           (if extrap (setf iroff2 (f2cl-lib:int-add iroff2 1)))
 206           (if (not extrap) (setf iroff1 (f2cl-lib:int-add iroff1 1)))
 207          label20
 208           (if (and (> last$ 10) (> erro12 errmax))
 209               (setf iroff3 (f2cl-lib:int-add iroff3 1)))
 210          label25
 211           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 212                   area1)
 213           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 214                   area2)
 215           (setf (f2cl-lib:fref nnlog-%data% (maxerr) ((1 *)) nnlog-%offset%)
 216                   nrmom)
 217           (setf (f2cl-lib:fref nnlog-%data% (last$) ((1 *)) nnlog-%offset%)
 218                   nrmom)
 219           (setf errbnd (max epsabs (* epsrel (abs area))))
 220           (if (or (>= (f2cl-lib:int-add iroff1 iroff2) 10) (>= iroff3 20))
 221               (setf ier 2))
 222           (if (>= iroff2 5) (setf ierro 3))
 223           (if (= last$ limit) (setf ier 1))
 224           (if
 225            (<= (max (abs a1) (abs b2))
 226                (* (+ 1.0 (* 100.0 epmach)) (+ (abs a2) (* 1000.0 uflow))))
 227            (setf ier 4))
 228           (if (> error2 error1) (go label30))
 229           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a2)
 230           (setf (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%)
 231                   b1)
 232           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b2)
 233           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 234                   error1)
 235           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 236                   error2)
 237           (go label40)
 238          label30
 239           (setf (f2cl-lib:fref alist-%data% (maxerr) ((1 *)) alist-%offset%)
 240                   a2)
 241           (setf (f2cl-lib:fref alist-%data% (last$) ((1 *)) alist-%offset%) a1)
 242           (setf (f2cl-lib:fref blist-%data% (last$) ((1 *)) blist-%offset%) b1)
 243           (setf (f2cl-lib:fref rlist-%data% (maxerr) ((1 *)) rlist-%offset%)
 244                   area2)
 245           (setf (f2cl-lib:fref rlist-%data% (last$) ((1 *)) rlist-%offset%)
 246                   area1)
 247           (setf (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%)
 248                   error2)
 249           (setf (f2cl-lib:fref elist-%data% (last$) ((1 *)) elist-%offset%)
 250                   error1)
 251          label40
 252           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 253               (dqpsrt limit last$ maxerr errmax elist iord nrmax)
 254             (declare (ignore var-0 var-1 var-4 var-5))
 255             (setf maxerr var-2)
 256             (setf errmax var-3)
 257             (setf nrmax var-6))
 258           (if (<= errsum errbnd) (go label170))
 259           (if (/= ier 0) (go label150))
 260           (if (and (= last$ 2) extall) (go label120))
 261           (if noext (go label140))
 262           (if (not extall) (go label50))
 263           (setf erlarg (- erlarg erlast))
 264           (if (> (abs (- b1 a1)) small) (setf erlarg (+ erlarg erro12)))
 265           (if extrap (go label70))
 266          label50
 267           (setf width
 268                   (abs
 269                    (-
 270                     (f2cl-lib:fref blist-%data%
 271                                    (maxerr)
 272                                    ((1 *))
 273                                    blist-%offset%)
 274                     (f2cl-lib:fref alist-%data%
 275                                    (maxerr)
 276                                    ((1 *))
 277                                    alist-%offset%))))
 278           (if (> width small) (go label140))
 279           (if extall (go label60))
 280           (setf small (* small 0.5))
 281           (if (> (* 0.25 width domega) 2.0) (go label140))
 282           (setf extall f2cl-lib:%true%)
 283           (go label130)
 284          label60
 285           (setf extrap f2cl-lib:%true%)
 286           (setf nrmax 2)
 287          label70
 288           (if (or (= ierro 3) (<= erlarg ertest)) (go label90))
 289           (setf jupbnd last$)
 290           (if (> last$ (+ (the f2cl-lib:integer4 (truncate limit 2)) 2))
 291               (setf jupbnd
 292                       (f2cl-lib:int-sub (f2cl-lib:int-add limit 3) last$)))
 293           (setf id nrmax)
 294           (f2cl-lib:fdo (k id (f2cl-lib:int-add k 1))
 295                         ((> k jupbnd) nil)
 296             (tagbody
 297               (setf maxerr
 298                       (f2cl-lib:fref iord-%data%
 299                                      (nrmax)
 300                                      ((1 *))
 301                                      iord-%offset%))
 302               (setf errmax
 303                       (f2cl-lib:fref elist-%data%
 304                                      (maxerr)
 305                                      ((1 *))
 306                                      elist-%offset%))
 307               (if
 308                (>
 309                 (abs
 310                  (-
 311                   (f2cl-lib:fref blist-%data% (maxerr) ((1 *)) blist-%offset%)
 312                   (f2cl-lib:fref alist-%data%
 313                                  (maxerr)
 314                                  ((1 *))
 315                                  alist-%offset%)))
 316                 small)
 317                (go label140))
 318               (setf nrmax (f2cl-lib:int-add nrmax 1))
 319              label80))
 320          label90
 321           (setf numrl2 (f2cl-lib:int-add numrl2 1))
 322           (setf (f2cl-lib:fref rlist2 (numrl2) ((1 52))) area)
 323           (if (< numrl2 3) (go label110))
 324           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 325               (dqelg numrl2 rlist2 reseps abseps res3la nres)
 326             (declare (ignore var-1 var-4))
 327             (setf numrl2 var-0)
 328             (setf reseps var-2)
 329             (setf abseps var-3)
 330             (setf nres var-5))
 331           (setf ktmin (f2cl-lib:int-add ktmin 1))
 332           (if (and (> ktmin 5) (< abserr (* 0.001 errsum))) (setf ier 5))
 333           (if (>= abseps abserr) (go label100))
 334           (setf ktmin 0)
 335           (setf abserr abseps)
 336           (setf result reseps)
 337           (setf correc erlarg)
 338           (setf ertest (max epsabs (* epsrel (abs reseps))))
 339           (if (<= abserr ertest) (go label150))
 340          label100
 341           (if (= numrl2 1) (setf noext f2cl-lib:%true%))
 342           (if (= ier 5) (go label150))
 343          label110
 344           (setf maxerr (f2cl-lib:fref iord-%data% (1) ((1 *)) iord-%offset%))
 345           (setf errmax
 346                   (f2cl-lib:fref elist-%data% (maxerr) ((1 *)) elist-%offset%))
 347           (setf nrmax 1)
 348           (setf extrap f2cl-lib:%false%)
 349           (setf small (* small 0.5))
 350           (setf erlarg errsum)
 351           (go label140)
 352          label120
 353           (setf small (* small 0.5))
 354           (setf numrl2 (f2cl-lib:int-add numrl2 1))
 355           (setf (f2cl-lib:fref rlist2 (numrl2) ((1 52))) area)
 356          label130
 357           (setf ertest errbnd)
 358           (setf erlarg errsum)
 359          label140))
 360      label150
 361       (if (or (= abserr oflow) (= nres 0)) (go label170))
 362       (if (= (f2cl-lib:int-add ier ierro) 0) (go label165))
 363       (if (= ierro 3) (setf abserr (+ abserr correc)))
 364       (if (= ier 0) (setf ier 3))
 365       (if (and (/= result 0.0) (/= area 0.0)) (go label160))
 366       (if (> abserr errsum) (go label170))
 367       (if (= area 0.0) (go label190))
 368       (go label165)
 369      label160
 370       (if (> (/ abserr (abs result)) (/ errsum (abs area))) (go label170))
 371      label165
 372       (if (and (= ksgn -1) (<= (max (abs result) (abs area)) (* defabs 0.01)))
 373           (go label190))
 374       (if
 375        (or (> 0.01 (/ result area))
 376            (> (/ result area) 100.0)
 377            (>= errsum (abs area)))
 378        (setf ier 6))
 379       (go label190)
 380      label170
 381       (setf result 0.0)
 382       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 383                     ((> k last$) nil)
 384         (tagbody
 385           (setf result
 386                   (+ result
 387                      (f2cl-lib:fref rlist-%data% (k) ((1 *)) rlist-%offset%)))
 388          label180))
 389       (setf abserr errsum)
 390      label190
 391       (if (> ier 2) (setf ier (f2cl-lib:int-sub ier 1)))
 392      label200
 393       (if (and (= integr 2) (< omega 0.0)) (setf result (- result)))
 394      label999
 395       (go end_label)
 396      end_label
 397       (return
 398        (values nil
 399                nil
 400                nil
 401                nil
 402                integr
 403                nil
 404                nil
 405                nil
 406                nil
 407                nil
 408                result
 409                abserr
 410                neval
 411                ier
 412                last$
 413                nil
 414                nil
 415                nil
 416                nil
 417                nil
 418                nil
 419                momcom
 420                nil)))))
 421
 422 (in-package #:cl-user)
 423 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 424 (eval-when (:load-toplevel :compile-toplevel :execute)
 425   (setf (gethash 'fortran-to-lisp::dqawoe
 426                  fortran-to-lisp::*f2cl-function-info*)
 427           (fortran-to-lisp::make-f2cl-finfo
 428            :arg-types '(t (double-float) (double-float) (double-float)
 429                         (fortran-to-lisp::integer4) (double-float)
 430                         (double-float) (fortran-to-lisp::integer4)
 431                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 432                         (double-float) (double-float)
 433                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 434                         (fortran-to-lisp::integer4) (array double-float (*))
 435                         (array double-float (*)) (array double-float (*))
 436                         (array double-float (*))
 437                         (array fortran-to-lisp::integer4 (*))
 438                         (array fortran-to-lisp::integer4 (*))
 439                         (fortran-to-lisp::integer4) (array double-float (*)))
 440            :return-values '(nil nil nil nil fortran-to-lisp::integr nil nil nil
 441                             nil nil fortran-to-lisp::result
 442                             fortran-to-lisp::abserr fortran-to-lisp::neval
 443                             fortran-to-lisp::ier fortran-to-lisp::last$ nil nil
 444                             nil nil nil nil fortran-to-lisp::momcom nil)
 445            :calls '(fortran-to-lisp::dqelg fortran-to-lisp::dqpsrt
 446                     fortran-to-lisp::dqc25f fortran-to-lisp::d1mach))))
 447