src/numerical/slatec/dqk15w.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 ((xgk
  21        (make-array 8
  22                    :element-type 'double-float
  23                    :initial-contents '(0.9914553711208126 0.9491079123427585
  24                                        0.8648644233597691 0.7415311855993943
  25                                        0.5860872354676911 0.4058451513773972
  26                                        0.2077849550078985 0.0)))
  27       (wgk
  28        (make-array 8
  29                    :element-type 'double-float
  30                    :initial-contents '(0.02293532201052922 0.06309209262997854
  31                                        0.1047900103222502 0.1406532597155259
  32                                        0.1690047266392679 0.1903505780647854
  33                                        0.2044329400752989 0.2094821410847278)))
  34       (wg
  35        (make-array 4
  36                    :element-type 'double-float
  37                    :initial-contents '(0.1294849661688697 0.2797053914892767
  38                                        0.3818300505051889 0.4179591836734694))))
  39   (declare (type (array double-float (8)) xgk wgk)
  40            (type (array double-float (4)) wg))
  41   (defun dqk15w (f w p1 p2 p3 p4 kp a b result abserr resabs resasc)
  42     (declare (type (f2cl-lib:integer4) kp)
  43              (type (double-float) resasc resabs abserr result b a p4 p3 p2 p1))
  44     (prog ((fv1 (make-array 7 :element-type 'double-float))
  45            (fv2 (make-array 7 :element-type 'double-float)) (j 0) (jtw 0)
  46            (jtwm1 0) (absc 0.0) (absc1 0.0) (absc2 0.0) (centr 0.0)
  47            (dhlgth 0.0) (epmach 0.0) (fc 0.0) (fsum 0.0) (fval1 0.0)
  48            (fval2 0.0) (hlgth 0.0) (resg 0.0) (resk 0.0) (reskh 0.0)
  49            (uflow 0.0))
  50       (declare (type (array double-float (7)) fv2 fv1)
  51                (type (double-float) uflow reskh resk resg hlgth fval2 fval1
  52                                     fsum fc epmach dhlgth centr absc2 absc1
  53                                     absc)
  54                (type (f2cl-lib:integer4) jtwm1 jtw j))
  55       (setf epmach (f2cl-lib:d1mach 4))
  56       (setf uflow (f2cl-lib:d1mach 1))
  57       (setf centr (* 0.5 (+ a b)))
  58       (setf hlgth (* 0.5 (- b a)))
  59       (setf dhlgth (abs hlgth))
  60       (setf fc
  61               (*
  62                (multiple-value-bind (ret-val var-0)
  63                    (funcall f centr)
  64                  (declare (ignore))
  65                  (when var-0
  66                    (setf centr var-0))
  67                  ret-val)
  68                (multiple-value-bind
  69                      (ret-val var-0 var-1 var-2 var-3 var-4 var-5)
  70                    (funcall w centr p1 p2 p3 p4 kp)
  71                  (declare (ignore))
  72                  (when var-0
  73                    (setf centr var-0))
  74                  (when var-1
  75                    (setf p1 var-1))
  76                  (when var-2
  77                    (setf p2 var-2))
  78                  (when var-3
  79                    (setf p3 var-3))
  80                  (when var-4
  81                    (setf p4 var-4))
  82                  (when var-5
  83                    (setf kp var-5))
  84                  ret-val)))
  85       (setf resg (* (f2cl-lib:fref wg (4) ((1 4))) fc))
  86       (setf resk (* (f2cl-lib:fref wgk (8) ((1 8))) fc))
  87       (setf resabs (abs resk))
  88       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  89                     ((> j 3) nil)
  90         (tagbody
  91           (setf jtw (f2cl-lib:int-mul j 2))
  92           (setf absc (* hlgth (f2cl-lib:fref xgk (jtw) ((1 8)))))
  93           (setf absc1 (- centr absc))
  94           (setf absc2 (+ centr absc))
  95           (setf fval1
  96                   (*
  97                    (multiple-value-bind (ret-val var-0)
  98                        (funcall f absc1)
  99                      (declare (ignore))
 100                      (when var-0
 101                        (setf absc1 var-0))
 102                      ret-val)
 103                    (multiple-value-bind
 104                          (ret-val var-0 var-1 var-2 var-3 var-4 var-5)
 105                        (funcall w absc1 p1 p2 p3 p4 kp)
 106                      (declare (ignore))
 107                      (when var-0
 108                        (setf absc1 var-0))
 109                      (when var-1
 110                        (setf p1 var-1))
 111                      (when var-2
 112                        (setf p2 var-2))
 113                      (when var-3
 114                        (setf p3 var-3))
 115                      (when var-4
 116                        (setf p4 var-4))
 117                      (when var-5
 118                        (setf kp var-5))
 119                      ret-val)))
 120           (setf fval2
 121                   (*
 122                    (multiple-value-bind (ret-val var-0)
 123                        (funcall f absc2)
 124                      (declare (ignore))
 125                      (when var-0
 126                        (setf absc2 var-0))
 127                      ret-val)
 128                    (multiple-value-bind
 129                          (ret-val var-0 var-1 var-2 var-3 var-4 var-5)
 130                        (funcall w absc2 p1 p2 p3 p4 kp)
 131                      (declare (ignore))
 132                      (when var-0
 133                        (setf absc2 var-0))
 134                      (when var-1
 135                        (setf p1 var-1))
 136                      (when var-2
 137                        (setf p2 var-2))
 138                      (when var-3
 139                        (setf p3 var-3))
 140                      (when var-4
 141                        (setf p4 var-4))
 142                      (when var-5
 143                        (setf kp var-5))
 144                      ret-val)))
 145           (setf (f2cl-lib:fref fv1 (jtw) ((1 7))) fval1)
 146           (setf (f2cl-lib:fref fv2 (jtw) ((1 7))) fval2)
 147           (setf fsum (+ fval1 fval2))
 148           (setf resg (+ resg (* (f2cl-lib:fref wg (j) ((1 4))) fsum)))
 149           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtw) ((1 8))) fsum)))
 150           (setf resabs
 151                   (+ resabs
 152                      (* (f2cl-lib:fref wgk (jtw) ((1 8)))
 153                         (+ (abs fval1) (abs fval2)))))
 154          label10))
 155       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 156                     ((> j 4) nil)
 157         (tagbody
 158           (setf jtwm1 (f2cl-lib:int-sub (f2cl-lib:int-mul j 2) 1))
 159           (setf absc (* hlgth (f2cl-lib:fref xgk (jtwm1) ((1 8)))))
 160           (setf absc1 (- centr absc))
 161           (setf absc2 (+ centr absc))
 162           (setf fval1
 163                   (*
 164                    (multiple-value-bind (ret-val var-0)
 165                        (funcall f absc1)
 166                      (declare (ignore))
 167                      (when var-0
 168                        (setf absc1 var-0))
 169                      ret-val)
 170                    (multiple-value-bind
 171                          (ret-val var-0 var-1 var-2 var-3 var-4 var-5)
 172                        (funcall w absc1 p1 p2 p3 p4 kp)
 173                      (declare (ignore))
 174                      (when var-0
 175                        (setf absc1 var-0))
 176                      (when var-1
 177                        (setf p1 var-1))
 178                      (when var-2
 179                        (setf p2 var-2))
 180                      (when var-3
 181                        (setf p3 var-3))
 182                      (when var-4
 183                        (setf p4 var-4))
 184                      (when var-5
 185                        (setf kp var-5))
 186                      ret-val)))
 187           (setf fval2
 188                   (*
 189                    (multiple-value-bind (ret-val var-0)
 190                        (funcall f absc2)
 191                      (declare (ignore))
 192                      (when var-0
 193                        (setf absc2 var-0))
 194                      ret-val)
 195                    (multiple-value-bind
 196                          (ret-val var-0 var-1 var-2 var-3 var-4 var-5)
 197                        (funcall w absc2 p1 p2 p3 p4 kp)
 198                      (declare (ignore))
 199                      (when var-0
 200                        (setf absc2 var-0))
 201                      (when var-1
 202                        (setf p1 var-1))
 203                      (when var-2
 204                        (setf p2 var-2))
 205                      (when var-3
 206                        (setf p3 var-3))
 207                      (when var-4
 208                        (setf p4 var-4))
 209                      (when var-5
 210                        (setf kp var-5))
 211                      ret-val)))
 212           (setf (f2cl-lib:fref fv1 (jtwm1) ((1 7))) fval1)
 213           (setf (f2cl-lib:fref fv2 (jtwm1) ((1 7))) fval2)
 214           (setf fsum (+ fval1 fval2))
 215           (setf resk (+ resk (* (f2cl-lib:fref wgk (jtwm1) ((1 8))) fsum)))
 216           (setf resabs
 217                   (+ resabs
 218                      (* (f2cl-lib:fref wgk (jtwm1) ((1 8)))
 219                         (+ (abs fval1) (abs fval2)))))
 220          label15))
 221       (setf reskh (* resk 0.5))
 222       (setf resasc (* (f2cl-lib:fref wgk (8) ((1 8))) (abs (- fc reskh))))
 223       (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 224                     ((> j 7) nil)
 225         (tagbody
 226           (setf resasc
 227                   (+ resasc
 228                      (* (f2cl-lib:fref wgk (j) ((1 8)))
 229                         (+ (abs (- (f2cl-lib:fref fv1 (j) ((1 7))) reskh))
 230                            (abs (- (f2cl-lib:fref fv2 (j) ((1 7))) reskh))))))
 231          label20))
 232       (setf result (* resk hlgth))
 233       (setf resabs (* resabs dhlgth))
 234       (setf resasc (* resasc dhlgth))
 235       (setf abserr (abs (* (- resk resg) hlgth)))
 236       (if (and (/= resasc 0.0) (/= abserr 0.0))
 237           (setf abserr
 238                   (* resasc (min 1.0 (expt (/ (* 200.0 abserr) resasc) 1.5)))))
 239       (if (> resabs (/ uflow (* 50.0 epmach)))
 240           (setf abserr (max (* epmach 50.0 resabs) abserr)))
 241       (go end_label)
 242      end_label
 243       (return
 244        (values nil nil p1 p2 p3 p4 kp nil nil result abserr resabs resasc)))))
 245
 246 (in-package #:cl-user)
 247 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 248 (eval-when (:load-toplevel :compile-toplevel :execute)
 249   (setf (gethash 'fortran-to-lisp::dqk15w
 250                  fortran-to-lisp::*f2cl-function-info*)
 251           (fortran-to-lisp::make-f2cl-finfo
 252            :arg-types '(t t (double-float) (double-float) (double-float)
 253                         (double-float) (fortran-to-lisp::integer4)
 254                         (double-float) (double-float) (double-float)
 255                         (double-float) (double-float) (double-float))
 256            :return-values '(nil nil fortran-to-lisp::p1 fortran-to-lisp::p2
 257                             fortran-to-lisp::p3 fortran-to-lisp::p4
 258                             fortran-to-lisp::kp nil nil fortran-to-lisp::result
 259                             fortran-to-lisp::abserr fortran-to-lisp::resabs
 260                             fortran-to-lisp::resasc)
 261            :calls '(fortran-to-lisp::d1mach))))
 262