src/numerical/slatec/zmlri.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 ':simple-array)
  14 ;;;           (:array-slicing nil) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :slatec)
  18
  19
  20 (let ((zeror 0.0) (zeroi 0.0) (coner 1.0) (conei 0.0))
  21   (declare (type (double-float) zeror zeroi coner conei))
  22   (defun zmlri (zr zi fnu kode n yr yi nz tol)
  23     (declare (type (simple-array double-float (*)) yi yr)
  24              (type (f2cl-lib:integer4) nz n kode)
  25              (type (double-float) tol fnu zi zr))
  26     (prog ((i 0) (iaz 0) (idum 0) (ifnu 0) (inu 0) (itime 0) (k 0) (kk 0)
  27            (km 0) (m 0) (ack 0.0) (ak 0.0) (ap 0.0) (at 0.0) (az 0.0) (bk 0.0)
  28            (cki 0.0) (ckr 0.0) (cnormi 0.0) (cnormr 0.0) (fkap 0.0) (fkk 0.0)
  29            (flam 0.0) (fnf 0.0) (pti 0.0) (ptr 0.0) (p1i 0.0) (p1r 0.0)
  30            (p2i 0.0) (p2r 0.0) (raz 0.0) (rho 0.0) (rho2 0.0) (rzi 0.0)
  31            (rzr 0.0) (scle 0.0) (sti 0.0) (str 0.0) (sumi 0.0) (sumr 0.0)
  32            (tfnf 0.0) (tst 0.0))
  33       (declare (type (double-float) tst tfnf sumr sumi str sti scle rzr rzi
  34                                     rho2 rho raz p2r p2i p1r p1i ptr pti fnf
  35                                     flam fkk fkap cnormr cnormi ckr cki bk az
  36                                     at ap ak ack)
  37                (type (f2cl-lib:integer4) m km kk k itime inu ifnu idum iaz i))
  38       (setf scle (/ (f2cl-lib:d1mach 1) tol))
  39       (setf nz 0)
  40       (setf az (coerce (realpart (zabs zr zi)) 'double-float))
  41       (setf iaz (f2cl-lib:int az))
  42       (setf ifnu (f2cl-lib:int fnu))
  43       (setf inu (f2cl-lib:int-sub (f2cl-lib:int-add ifnu n) 1))
  44       (setf at (+ iaz 1.0))
  45       (setf raz (/ 1.0 az))
  46       (setf str (* zr raz))
  47       (setf sti (* (- zi) raz))
  48       (setf ckr (* str at raz))
  49       (setf cki (* sti at raz))
  50       (setf rzr (* (+ str str) raz))
  51       (setf rzi (* (+ sti sti) raz))
  52       (setf p1r zeror)
  53       (setf p1i zeroi)
  54       (setf p2r coner)
  55       (setf p2i conei)
  56       (setf ack (* (+ at 1.0) raz))
  57       (setf rho (+ ack (f2cl-lib:fsqrt (- (* ack ack) 1.0))))
  58       (setf rho2 (* rho rho))
  59       (setf tst (/ (+ rho2 rho2) (* (- rho2 1.0) (- rho 1.0))))
  60       (setf tst (/ tst tol))
  61       (setf ak at)
  62       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  63                     ((> i 80) nil)
  64         (tagbody
  65           (setf ptr p2r)
  66           (setf pti p2i)
  67           (setf p2r (- p1r (- (* ckr ptr) (* cki pti))))
  68           (setf p2i (- p1i (+ (* cki ptr) (* ckr pti))))
  69           (setf p1r ptr)
  70           (setf p1i pti)
  71           (setf ckr (+ ckr rzr))
  72           (setf cki (+ cki rzi))
  73           (setf ap (coerce (realpart (zabs p2r p2i)) 'double-float))
  74           (if (> ap (* tst ak ak)) (go label20))
  75           (setf ak (+ ak 1.0))
  76          label10))
  77       (go label110)
  78      label20
  79       (setf i (f2cl-lib:int-add i 1))
  80       (setf k 0)
  81       (if (< inu iaz) (go label40))
  82       (setf p1r zeror)
  83       (setf p1i zeroi)
  84       (setf p2r coner)
  85       (setf p2i conei)
  86       (setf at (+ inu 1.0))
  87       (setf str (* zr raz))
  88       (setf sti (* (- zi) raz))
  89       (setf ckr (* str at raz))
  90       (setf cki (* sti at raz))
  91       (setf ack (* at raz))
  92       (setf tst (f2cl-lib:fsqrt (/ ack tol)))
  93       (setf itime 1)
  94       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  95                     ((> k 80) nil)
  96         (tagbody
  97           (setf ptr p2r)
  98           (setf pti p2i)
  99           (setf p2r (- p1r (- (* ckr ptr) (* cki pti))))
 100           (setf p2i (- p1i (+ (* ckr pti) (* cki ptr))))
 101           (setf p1r ptr)
 102           (setf p1i pti)
 103           (setf ckr (+ ckr rzr))
 104           (setf cki (+ cki rzi))
 105           (setf ap (coerce (realpart (zabs p2r p2i)) 'double-float))
 106           (if (< ap tst) (go label30))
 107           (if (= itime 2) (go label40))
 108           (setf ack (coerce (realpart (zabs ckr cki)) 'double-float))
 109           (setf flam (+ ack (f2cl-lib:fsqrt (- (* ack ack) 1.0))))
 110           (setf fkap (coerce (realpart (/ ap (zabs p1r p1i))) 'double-float))
 111           (setf rho (min flam fkap))
 112           (setf tst (* tst (f2cl-lib:fsqrt (/ rho (- (* rho rho) 1.0)))))
 113           (setf itime 2)
 114          label30))
 115       (go label110)
 116      label40
 117       (setf k (f2cl-lib:int-add k 1))
 118       (setf kk
 119               (max (the f2cl-lib:integer4 (f2cl-lib:int-add i iaz))
 120                    (the f2cl-lib:integer4 (f2cl-lib:int-add k inu))))
 121       (setf fkk (coerce (the f2cl-lib:integer4 kk) 'double-float))
 122       (setf p1r zeror)
 123       (setf p1i zeroi)
 124       (setf p2r scle)
 125       (setf p2i zeroi)
 126       (setf fnf (- fnu ifnu))
 127       (setf tfnf (+ fnf fnf))
 128       (setf bk
 129               (-
 130                (multiple-value-bind (ret-val var-0 var-1)
 131                    (dgamln (+ fkk tfnf 1.0) idum)
 132                  (declare (ignore var-0))
 133                  (setf idum var-1)
 134                  ret-val)
 135                (multiple-value-bind (ret-val var-0 var-1)
 136                    (dgamln (+ fkk 1.0) idum)
 137                  (declare (ignore var-0))
 138                  (setf idum var-1)
 139                  ret-val)
 140                (multiple-value-bind (ret-val var-0 var-1)
 141                    (dgamln (+ tfnf 1.0) idum)
 142                  (declare (ignore var-0))
 143                  (setf idum var-1)
 144                  ret-val)))
 145       (setf bk (exp bk))
 146       (setf sumr zeror)
 147       (setf sumi zeroi)
 148       (setf km (f2cl-lib:int-sub kk inu))
 149       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 150                     ((> i km) nil)
 151         (tagbody
 152           (setf ptr p2r)
 153           (setf pti p2i)
 154           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
 155           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzi ptr) (* rzr pti)))))
 156           (setf p1r ptr)
 157           (setf p1i pti)
 158           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
 159           (setf ack (* bk ak))
 160           (setf sumr (+ sumr (* (+ ack bk) p1r)))
 161           (setf sumi (+ sumi (* (+ ack bk) p1i)))
 162           (setf bk ack)
 163           (setf fkk (- fkk 1.0))
 164          label50))
 165       (setf (f2cl-lib:fref yr (n) ((1 n))) p2r)
 166       (setf (f2cl-lib:fref yi (n) ((1 n))) p2i)
 167       (if (= n 1) (go label70))
 168       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 169                     ((> i n) nil)
 170         (tagbody
 171           (setf ptr p2r)
 172           (setf pti p2i)
 173           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
 174           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzi ptr) (* rzr pti)))))
 175           (setf p1r ptr)
 176           (setf p1i pti)
 177           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
 178           (setf ack (* bk ak))
 179           (setf sumr (+ sumr (* (+ ack bk) p1r)))
 180           (setf sumi (+ sumi (* (+ ack bk) p1i)))
 181           (setf bk ack)
 182           (setf fkk (- fkk 1.0))
 183           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub n i) 1))
 184           (setf (f2cl-lib:fref yr (m) ((1 n))) p2r)
 185           (setf (f2cl-lib:fref yi (m) ((1 n))) p2i)
 186          label60))
 187      label70
 188       (if (<= ifnu 0) (go label90))
 189       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 190                     ((> i ifnu) nil)
 191         (tagbody
 192           (setf ptr p2r)
 193           (setf pti p2i)
 194           (setf p2r (+ p1r (* (+ fkk fnf) (- (* rzr ptr) (* rzi pti)))))
 195           (setf p2i (+ p1i (* (+ fkk fnf) (+ (* rzr pti) (* rzi ptr)))))
 196           (setf p1r ptr)
 197           (setf p1i pti)
 198           (setf ak (+ 1.0 (/ (- tfnf) (+ fkk tfnf))))
 199           (setf ack (* bk ak))
 200           (setf sumr (+ sumr (* (+ ack bk) p1r)))
 201           (setf sumi (+ sumi (* (+ ack bk) p1i)))
 202           (setf bk ack)
 203           (setf fkk (- fkk 1.0))
 204          label80))
 205      label90
 206       (setf ptr zr)
 207       (setf pti zi)
 208       (if (= kode 2) (setf ptr zeror))
 209       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4)
 210           (zlog rzr rzi str sti idum)
 211         (declare (ignore var-0 var-1))
 212         (setf str var-2)
 213         (setf sti var-3)
 214         (setf idum var-4))
 215       (setf p1r (+ (* (- fnf) str) ptr))
 216       (setf p1i (+ (* (- fnf) sti) pti))
 217       (setf ap
 218               (multiple-value-bind (ret-val var-0 var-1)
 219                   (dgamln (+ 1.0 fnf) idum)
 220                 (declare (ignore var-0))
 221                 (setf idum var-1)
 222                 ret-val))
 223       (setf ptr (- p1r ap))
 224       (setf pti p1i)
 225       (setf p2r (+ p2r sumr))
 226       (setf p2i (+ p2i sumi))
 227       (setf ap (coerce (realpart (zabs p2r p2i)) 'double-float))
 228       (setf p1r (/ 1.0 ap))
 229       (multiple-value-bind (var-0 var-1 var-2 var-3)
 230           (zexp ptr pti str sti)
 231         (declare (ignore var-0 var-1))
 232         (setf str var-2)
 233         (setf sti var-3))
 234       (setf ckr (* str p1r))
 235       (setf cki (* sti p1r))
 236       (setf ptr (* p2r p1r))
 237       (setf pti (* (- p2i) p1r))
 238       (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
 239           (zmlt ckr cki ptr pti cnormr cnormi)
 240         (declare (ignore var-0 var-1 var-2 var-3))
 241         (setf cnormr var-4)
 242         (setf cnormi var-5))
 243       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 244                     ((> i n) nil)
 245         (tagbody
 246           (setf str
 247                   (- (* (f2cl-lib:fref yr (i) ((1 n))) cnormr)
 248                      (* (f2cl-lib:fref yi (i) ((1 n))) cnormi)))
 249           (setf (f2cl-lib:fref yi (i) ((1 n)))
 250                   (+ (* (f2cl-lib:fref yr (i) ((1 n))) cnormi)
 251                      (* (f2cl-lib:fref yi (i) ((1 n))) cnormr)))
 252           (setf (f2cl-lib:fref yr (i) ((1 n))) str)
 253          label100))
 254       (go end_label)
 255      label110
 256       (setf nz -2)
 257       (go end_label)
 258      end_label
 259       (return (values nil nil nil nil nil nil nil nz nil)))))
 260
 261 (in-package #:cl-user)
 262 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 263 (eval-when (:load-toplevel :compile-toplevel :execute)
 264   (setf (gethash 'fortran-to-lisp::zmlri fortran-to-lisp::*f2cl-function-info*)
 265           (fortran-to-lisp::make-f2cl-finfo
 266            :arg-types '((double-float) (double-float) (double-float)
 267                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 268                         (simple-array double-float (*))
 269                         (simple-array double-float (*))
 270                         (fortran-to-lisp::integer4) (double-float))
 271            :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::nz
 272                             nil)
 273            :calls '(fortran-to-lisp::zmlt fortran-to-lisp::zexp
 274                     fortran-to-lisp::zlog fortran-to-lisp::dgamln
 275                     fortran-to-lisp::zabs fortran-to-lisp::d1mach))))
 276