src/numerical/slatec/dqelg.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 dqelg (n epstab result abserr res3la nres)
  21   (declare (type (array double-float (*)) res3la)
  22            (type (double-float) abserr result)
  23            (type (array double-float (*)) epstab)
  24            (type (f2cl-lib:integer4) nres n))
  25   (f2cl-lib:with-multi-array-data
  26       ((epstab double-float epstab-%data% epstab-%offset%)
  27        (res3la double-float res3la-%data% res3la-%offset%))
  28     (prog ((i 0) (ib 0) (ib2 0) (ie 0) (indx 0) (k1 0) (k2 0) (k3 0) (limexp 0)
  29            (newelm 0) (num 0) (delta1 0.0) (delta2 0.0) (delta3 0.0)
  30            (epmach 0.0) (epsinf 0.0) (error$ 0.0) (err1 0.0) (err2 0.0)
  31            (err3 0.0) (e0 0.0) (e1 0.0) (e1abs 0.0) (e2 0.0) (e3 0.0)
  32            (oflow 0.0) (res 0.0) (ss 0.0) (tol1 0.0) (tol2 0.0) (tol3 0.0))
  33       (declare (type (double-float) tol3 tol2 tol1 ss res oflow e3 e2 e1abs e1
  34                                     e0 err3 err2 err1 error$ epsinf epmach
  35                                     delta3 delta2 delta1)
  36                (type (f2cl-lib:integer4) num newelm limexp k3 k2 k1 indx ie ib2
  37                                          ib i))
  38       (setf epmach (f2cl-lib:d1mach 4))
  39       (setf oflow (f2cl-lib:d1mach 2))
  40       (setf nres (f2cl-lib:int-add nres 1))
  41       (setf abserr oflow)
  42       (setf result (f2cl-lib:fref epstab-%data% (n) ((1 52)) epstab-%offset%))
  43       (if (< n 3) (go label100))
  44       (setf limexp 50)
  45       (setf (f2cl-lib:fref epstab-%data%
  46                            ((f2cl-lib:int-add n 2))
  47                            ((1 52))
  48                            epstab-%offset%)
  49               (f2cl-lib:fref epstab-%data% (n) ((1 52)) epstab-%offset%))
  50       (setf newelm (the f2cl-lib:integer4 (truncate (- n 1) 2)))
  51       (setf (f2cl-lib:fref epstab-%data% (n) ((1 52)) epstab-%offset%) oflow)
  52       (setf num n)
  53       (setf k1 n)
  54       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  55                     ((> i newelm) nil)
  56         (tagbody
  57           (setf k2 (f2cl-lib:int-sub k1 1))
  58           (setf k3 (f2cl-lib:int-sub k1 2))
  59           (setf res
  60                   (f2cl-lib:fref epstab-%data%
  61                                  ((f2cl-lib:int-add k1 2))
  62                                  ((1 52))
  63                                  epstab-%offset%))
  64           (setf e0 (f2cl-lib:fref epstab-%data% (k3) ((1 52)) epstab-%offset%))
  65           (setf e1 (f2cl-lib:fref epstab-%data% (k2) ((1 52)) epstab-%offset%))
  66           (setf e2 res)
  67           (setf e1abs (abs e1))
  68           (setf delta2 (- e2 e1))
  69           (setf err2 (abs delta2))
  70           (setf tol2 (* (max (abs e2) e1abs) epmach))
  71           (setf delta3 (- e1 e0))
  72           (setf err3 (abs delta3))
  73           (setf tol3 (* (max e1abs (abs e0)) epmach))
  74           (if (or (> err2 tol2) (> err3 tol3)) (go label10))
  75           (setf result res)
  76           (setf abserr (+ err2 err3))
  77           (go label100)
  78          label10
  79           (setf e3 (f2cl-lib:fref epstab-%data% (k1) ((1 52)) epstab-%offset%))
  80           (setf (f2cl-lib:fref epstab-%data% (k1) ((1 52)) epstab-%offset%) e1)
  81           (setf delta1 (- e1 e3))
  82           (setf err1 (abs delta1))
  83           (setf tol1 (* (max e1abs (abs e3)) epmach))
  84           (if (or (<= err1 tol1) (<= err2 tol2) (<= err3 tol3)) (go label20))
  85           (setf ss (+ (/ 1.0 delta1) (/ 1.0 delta2) (/ -1.0 delta3)))
  86           (setf epsinf (abs (* ss e1)))
  87           (if (> epsinf 1.0e-4) (go label30))
  88          label20
  89           (setf n (f2cl-lib:int-sub (f2cl-lib:int-add i i) 1))
  90           (go label50)
  91          label30
  92           (setf res (+ e1 (/ 1.0 ss)))
  93           (setf (f2cl-lib:fref epstab-%data% (k1) ((1 52)) epstab-%offset%)
  94                   res)
  95           (setf k1 (f2cl-lib:int-sub k1 2))
  96           (setf error$ (+ err2 (abs (- res e2)) err3))
  97           (if (> error$ abserr) (go label40))
  98           (setf abserr error$)
  99           (setf result res)
 100          label40))
 101      label50
 102       (if (= n limexp)
 103           (setf n (- (* 2 (the f2cl-lib:integer4 (truncate limexp 2))) 1)))
 104       (setf ib 1)
 105       (if (= (* (the f2cl-lib:integer4 (truncate num 2)) 2) num) (setf ib 2))
 106       (setf ie (f2cl-lib:int-add newelm 1))
 107       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 108                     ((> i ie) nil)
 109         (tagbody
 110           (setf ib2 (f2cl-lib:int-add ib 2))
 111           (setf (f2cl-lib:fref epstab-%data% (ib) ((1 52)) epstab-%offset%)
 112                   (f2cl-lib:fref epstab-%data% (ib2) ((1 52)) epstab-%offset%))
 113           (setf ib ib2)
 114          label60))
 115       (if (= num n) (go label80))
 116       (setf indx (f2cl-lib:int-add (f2cl-lib:int-sub num n) 1))
 117       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 118                     ((> i n) nil)
 119         (tagbody
 120           (setf (f2cl-lib:fref epstab-%data% (i) ((1 52)) epstab-%offset%)
 121                   (f2cl-lib:fref epstab-%data%
 122                                  (indx)
 123                                  ((1 52))
 124                                  epstab-%offset%))
 125           (setf indx (f2cl-lib:int-add indx 1))
 126          label70))
 127      label80
 128       (if (>= nres 4) (go label90))
 129       (setf (f2cl-lib:fref res3la-%data% (nres) ((1 3)) res3la-%offset%)
 130               result)
 131       (setf abserr oflow)
 132       (go label100)
 133      label90
 134       (setf abserr
 135               (+
 136                (abs
 137                 (- result
 138                    (f2cl-lib:fref res3la-%data% (3) ((1 3)) res3la-%offset%)))
 139                (abs
 140                 (- result
 141                    (f2cl-lib:fref res3la-%data% (2) ((1 3)) res3la-%offset%)))
 142                (abs
 143                 (- result
 144                    (f2cl-lib:fref res3la-%data%
 145                                   (1)
 146                                   ((1 3))
 147                                   res3la-%offset%)))))
 148       (setf (f2cl-lib:fref res3la-%data% (1) ((1 3)) res3la-%offset%)
 149               (f2cl-lib:fref res3la-%data% (2) ((1 3)) res3la-%offset%))
 150       (setf (f2cl-lib:fref res3la-%data% (2) ((1 3)) res3la-%offset%)
 151               (f2cl-lib:fref res3la-%data% (3) ((1 3)) res3la-%offset%))
 152       (setf (f2cl-lib:fref res3la-%data% (3) ((1 3)) res3la-%offset%) result)
 153      label100
 154       (setf abserr (max abserr (* 5.0 epmach (abs result))))
 155       (go end_label)
 156      end_label
 157       (return (values n nil result abserr nil nres)))))
 158
 159 (in-package #:cl-user)
 160 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 161 (eval-when (:load-toplevel :compile-toplevel :execute)
 162   (setf (gethash 'fortran-to-lisp::dqelg fortran-to-lisp::*f2cl-function-info*)
 163           (fortran-to-lisp::make-f2cl-finfo
 164            :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
 165                         (double-float) (double-float) (array double-float (*))
 166                         (fortran-to-lisp::integer4))
 167            :return-values '(fortran-to-lisp::n nil fortran-to-lisp::result
 168                             fortran-to-lisp::abserr nil fortran-to-lisp::nres)
 169            :calls '(fortran-to-lisp::d1mach))))
 170