src/numerical/slatec/dqng.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 (let ((x1
  21        (make-array 5
  22                    :element-type 'double-float
  23                    :initial-contents '(0.9739065285171717 0.8650633666889845
  24                                        0.6794095682990244 0.4333953941292472
  25                                        0.14887433898163122)))
  26       (w10
  27        (make-array 5
  28                    :element-type 'double-float
  29                    :initial-contents '(0.06667134430868814 0.1494513491505806
  30                                        0.21908636251598204 0.26926671930999635
  31                                        0.29552422471475287)))
  32       (x2
  33        (make-array 5
  34                    :element-type 'double-float
  35                    :initial-contents '(0.9956571630258081 0.9301574913557082
  36                                        0.7808177265864169 0.5627571346686047
  37                                        0.2943928627014602)))
  38       (w21a
  39        (make-array 5
  40                    :element-type 'double-float
  41                    :initial-contents '(0.032558162307964725 0.07503967481091996
  42                                        0.10938715880229764 0.13470921731147334
  43                                        0.14773910490133849)))
  44       (w21b
  45        (make-array 6
  46                    :element-type 'double-float
  47                    :initial-contents '(0.011694638867371874
  48                                        0.054755896574351995 0.0931254545836976
  49                                        0.12349197626206584 0.14277593857706009
  50                                        0.1494455540029169)))
  51       (x3
  52        (make-array 11
  53                    :element-type 'double-float
  54                    :initial-contents '(0.999333360901932 0.9874334029080889
  55                                        0.9548079348142663 0.9001486957483283
  56                                        0.8251983149831141 0.732148388989305
  57                                        0.6228479705377252 0.4994795740710565
  58                                        0.36490166134658075 0.2222549197766013
  59                                        0.07465061746138332)))
  60       (w43a
  61        (make-array 10
  62                    :element-type 'double-float
  63                    :initial-contents '(0.016296734289666565 0.0375228761208695
  64                                        0.05469490205825544 0.06735541460947808
  65                                        0.07387019963239395 0.005768556059769796
  66                                        0.027371890593248842 0.04656082691042883
  67                                        0.06174499520144257
  68                                        0.07138726726869339)))
  69       (w43b
  70        (make-array 12
  71                    :element-type 'double-float
  72                    :initial-contents '(0.001844477640212414
  73                                        0.010798689585891651
  74                                        0.021895363867795427
  75                                        0.032597463975345686 0.04216313793519181
  76                                        0.050741939600184575 0.05837939554261925
  77                                        0.06474640495144589 0.06956619791235648
  78                                        0.07282444147183322 0.07450775101417512
  79                                        0.07472214751740301)))
  80       (x4
  81        (make-array 22
  82                    :element-type 'double-float
  83                    :initial-contents '(0.9999029772627293 0.9979898959866788
  84                                        0.9921754978606873 0.9813581635727128
  85                                        0.9650576238583847 0.9431676131336706
  86                                        0.9158064146855072 0.8832216577713164
  87                                        0.8457107484624157 0.8035576580352309
  88                                        0.7570057306854956 0.7062732097873218
  89                                        0.6515894665011779 0.5932233740579611
  90                                        0.531493605970832 0.46676362304202285
  91                                        0.3994248478592188 0.3298748771061883
  92                                        0.25850355920216156 0.18569539656834666
  93                                        0.11184221317990747
  94                                        0.03735212339461987)))
  95       (w87a
  96        (make-array 21
  97                    :element-type 'double-float
  98                    :initial-contents '(0.008148377384149173
  99                                        0.018761438201562824
 100                                        0.027347451050052287 0.03367770731163793
 101                                        0.036935099820427905
 102                                        0.0028848724302115306
 103                                        0.013685946022712702 0.02328041350288831
 104                                        0.03087249761171336 0.03569363363941877
 105                                        9.152833452022414e-4
 106                                        0.005399280219300471 0.01094767960111893
 107                                        0.016298731696787336
 108                                        0.021081568889203834
 109                                        0.025370969769253827
 110                                        0.029189697756475754 0.03237320246720279
 111                                        0.034783098950365146 0.03641222073135179
 112                                        0.037253875503047706)))
 113       (w87b
 114        (make-array 23
 115                    :element-type 'double-float
 116                    :initial-contents '(2.7414556376207234e-4
 117                                        0.0018071241550579428
 118                                        0.0040968692827591646
 119                                        0.006758290051847379
 120                                        0.009549957672201646
 121                                        0.012329447652244854
 122                                        0.015010447346388952 0.01754896798624319
 123                                        0.019938037786440887
 124                                        0.022194935961012286
 125                                        0.024339147126000805
 126                                        0.026374505414839208 0.0282869107887712
 127                                        0.030052581128092695 0.03164675137143993
 128                                        0.033050413419978504
 129                                        0.034255099704226064 0.03526241266015668
 130                                        0.0360769896228887 0.03669860449845609
 131                                        0.037120549269832576 0.03733422875193504
 132                                        0.037361073762679026))))
 133   (declare (type (array double-float (5)) x1 w10 x2 w21a)
 134            (type (array double-float (6)) w21b)
 135            (type (array double-float (11)) x3)
 136            (type (array double-float (10)) w43a)
 137            (type (array double-float (12)) w43b)
 138            (type (array double-float (22)) x4)
 139            (type (array double-float (21)) w87a)
 140            (type (array double-float (23)) w87b))
 141   (defun dqng (f a b epsabs epsrel result abserr neval ier)
 142     (declare (type (f2cl-lib:integer4) ier neval)
 143              (type (double-float) abserr result epsrel epsabs b a))
 144     (prog ((fv1 (make-array 5 :element-type 'double-float))
 145            (fv2 (make-array 5 :element-type 'double-float))
 146            (fv3 (make-array 5 :element-type 'double-float))
 147            (fv4 (make-array 5 :element-type 'double-float))
 148            (savfun (make-array 21 :element-type 'double-float)) (ipx 0) (k 0)
 149            (l 0) (absc 0.0) (centr 0.0) (dhlgth 0.0) (epmach 0.0) (fcentr 0.0)
 150            (fval 0.0) (fval1 0.0) (fval2 0.0) (hlgth 0.0) (res10 0.0)
 151            (res21 0.0) (res43 0.0) (res87 0.0) (resabs 0.0) (resasc 0.0)
 152            (reskh 0.0) (uflow 0.0))
 153       (declare (type (array double-float (21)) savfun)
 154                (type (array double-float (5)) fv4 fv3 fv2 fv1)
 155                (type (double-float) uflow reskh resasc resabs res87 res43 res21
 156                                     res10 hlgth fval2 fval1 fval fcentr epmach
 157                                     dhlgth centr absc)
 158                (type (f2cl-lib:integer4) l k ipx))
 159       (setf epmach (f2cl-lib:d1mach 4))
 160       (setf uflow (f2cl-lib:d1mach 1))
 161       (setf result 0.0)
 162       (setf abserr 0.0)
 163       (setf neval 0)
 164       (setf ier 6)
 165       (if (and (<= epsabs 0.0) (< epsrel (max (* 50.0 epmach) 5.0e-29)))
 166           (go label80))
 167       (setf hlgth (* 0.5 (- b a)))
 168       (setf dhlgth (abs hlgth))
 169       (setf centr (* 0.5 (+ b a)))
 170       (setf fcentr
 171               (multiple-value-bind (ret-val var-0)
 172                   (funcall f centr)
 173                 (declare (ignore))
 174                 (when var-0
 175                   (setf centr var-0))
 176                 ret-val))
 177       (setf neval 21)
 178       (setf ier 1)
 179       (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 180                     ((> l 3) nil)
 181         (tagbody
 182           (f2cl-lib:computed-goto (label5 label25 label45) l)
 183          label5
 184           (setf res10 0.0)
 185           (setf res21 (* (f2cl-lib:fref w21b (6) ((1 6))) fcentr))
 186           (setf resabs (* (f2cl-lib:fref w21b (6) ((1 6))) (abs fcentr)))
 187           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 188                         ((> k 5) nil)
 189             (tagbody
 190               (setf absc (* hlgth (f2cl-lib:fref x1 (k) ((1 5)))))
 191               (setf fval1 (funcall f (+ centr absc)))
 192               (setf fval2 (funcall f (- centr absc)))
 193               (setf fval (+ fval1 fval2))
 194               (setf res10 (+ res10 (* (f2cl-lib:fref w10 (k) ((1 5))) fval)))
 195               (setf res21 (+ res21 (* (f2cl-lib:fref w21a (k) ((1 5))) fval)))
 196               (setf resabs
 197                       (+ resabs
 198                          (* (f2cl-lib:fref w21a (k) ((1 5)))
 199                             (+ (abs fval1) (abs fval2)))))
 200               (setf (f2cl-lib:fref savfun (k) ((1 21))) fval)
 201               (setf (f2cl-lib:fref fv1 (k) ((1 5))) fval1)
 202               (setf (f2cl-lib:fref fv2 (k) ((1 5))) fval2)
 203              label10))
 204           (setf ipx 5)
 205           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 206                         ((> k 5) nil)
 207             (tagbody
 208               (setf ipx (f2cl-lib:int-add ipx 1))
 209               (setf absc (* hlgth (f2cl-lib:fref x2 (k) ((1 5)))))
 210               (setf fval1 (funcall f (+ centr absc)))
 211               (setf fval2 (funcall f (- centr absc)))
 212               (setf fval (+ fval1 fval2))
 213               (setf res21 (+ res21 (* (f2cl-lib:fref w21b (k) ((1 6))) fval)))
 214               (setf resabs
 215                       (+ resabs
 216                          (* (f2cl-lib:fref w21b (k) ((1 6)))
 217                             (+ (abs fval1) (abs fval2)))))
 218               (setf (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
 219               (setf (f2cl-lib:fref fv3 (k) ((1 5))) fval1)
 220               (setf (f2cl-lib:fref fv4 (k) ((1 5))) fval2)
 221              label15))
 222           (setf result (* res21 hlgth))
 223           (setf resabs (* resabs dhlgth))
 224           (setf reskh (* 0.5 res21))
 225           (setf resasc
 226                   (* (f2cl-lib:fref w21b (6) ((1 6))) (abs (- fcentr reskh))))
 227           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 228                         ((> k 5) nil)
 229             (tagbody
 230               (setf resasc
 231                       (+ resasc
 232                          (* (f2cl-lib:fref w21a (k) ((1 5)))
 233                             (+ (abs (- (f2cl-lib:fref fv1 (k) ((1 5))) reskh))
 234                                (abs
 235                                 (- (f2cl-lib:fref fv2 (k) ((1 5))) reskh))))
 236                          (* (f2cl-lib:fref w21b (k) ((1 6)))
 237                             (+ (abs (- (f2cl-lib:fref fv3 (k) ((1 5))) reskh))
 238                                (abs
 239                                 (- (f2cl-lib:fref fv4 (k) ((1 5))) reskh))))))
 240              label20))
 241           (setf abserr (abs (* (- res21 res10) hlgth)))
 242           (setf resasc (* resasc dhlgth))
 243           (go label65)
 244          label25
 245           (setf res43 (* (f2cl-lib:fref w43b (12) ((1 12))) fcentr))
 246           (setf neval 43)
 247           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 248                         ((> k 10) nil)
 249             (tagbody
 250               (setf res43
 251                       (+ res43
 252                          (* (f2cl-lib:fref savfun (k) ((1 21)))
 253                             (f2cl-lib:fref w43a (k) ((1 10))))))
 254              label30))
 255           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 256                         ((> k 11) nil)
 257             (tagbody
 258               (setf ipx (f2cl-lib:int-add ipx 1))
 259               (setf absc (* hlgth (f2cl-lib:fref x3 (k) ((1 11)))))
 260               (setf fval
 261                       (+ (funcall f (+ absc centr))
 262                          (funcall f (- centr absc))))
 263               (setf res43 (+ res43 (* fval (f2cl-lib:fref w43b (k) ((1 12))))))
 264               (setf (f2cl-lib:fref savfun (ipx) ((1 21))) fval)
 265              label40))
 266           (setf result (* res43 hlgth))
 267           (setf abserr (abs (* (- res43 res21) hlgth)))
 268           (go label65)
 269          label45
 270           (setf res87 (* (f2cl-lib:fref w87b (23) ((1 23))) fcentr))
 271           (setf neval 87)
 272           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 273                         ((> k 21) nil)
 274             (tagbody
 275               (setf res87
 276                       (+ res87
 277                          (* (f2cl-lib:fref savfun (k) ((1 21)))
 278                             (f2cl-lib:fref w87a (k) ((1 21))))))
 279              label50))
 280           (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
 281                         ((> k 22) nil)
 282             (tagbody
 283               (setf absc (* hlgth (f2cl-lib:fref x4 (k) ((1 22)))))
 284               (setf res87
 285                       (+ res87
 286                          (* (f2cl-lib:fref w87b (k) ((1 23)))
 287                             (+ (funcall f (+ absc centr))
 288                                (funcall f (- centr absc))))))
 289              label60))
 290           (setf result (* res87 hlgth))
 291           (setf abserr (abs (* (- res87 res43) hlgth)))
 292          label65
 293           (if (and (/= resasc 0.0) (/= abserr 0.0))
 294               (setf abserr
 295                       (* resasc
 296                          (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 297           (if (> resabs (/ uflow (* 50.0 epmach)))
 298               (setf abserr (max (* epmach 50.0 resabs) abserr)))
 299           (if (<= abserr (max epsabs (* epsrel (abs result)))) (setf ier 0))
 300           (if (= ier 0) (go label999))
 301          label70))
 302      label80
 303       (xermsg "SLATEC" "DQNG" "ABNORMAL RETURN" ier 0)
 304      label999
 305       (go end_label)
 306      end_label
 307       (return (values nil nil nil nil nil result abserr neval ier)))))
 308
 309 (in-package #:cl-user)
 310 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 311 (eval-when (:load-toplevel :compile-toplevel :execute)
 312   (setf (gethash 'fortran-to-lisp::dqng fortran-to-lisp::*f2cl-function-info*)
 313           (fortran-to-lisp::make-f2cl-finfo
 314            :arg-types '(t (double-float) (double-float) (double-float)
 315                         (double-float) (double-float) (double-float)
 316                         (fortran-to-lisp::integer4)
 317                         (fortran-to-lisp::integer4))
 318            :return-values '(nil nil nil nil nil fortran-to-lisp::result
 319                             fortran-to-lisp::abserr fortran-to-lisp::neval
 320                             fortran-to-lisp::ier)
 321            :calls '(fortran-to-lisp::xermsg fortran-to-lisp::d1mach))))
 322