share/odepack/src/dlhin.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
   3 ;;;  "f2cl2.l,v 95098eb54f13 2013/04/01 00:45:16 toy $"
   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 95098eb54f13 2013/04/01 00:45:16 toy $"
   7 ;;;  "f2cl6.l,v 1d5cbacbb977 2008/08/24 00:56:27 rtoy $"
   8 ;;;  "macros.l,v 1409c1352feb 2013/03/24 20:44:50 toy $")
   9
  10 ;;; Using Lisp CMU Common Lisp snapshot-2013-11 (20E 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 single-float))
  16
  17 (in-package "ODEPACK")
  18
  19
  20 (let ((half 0.5d0) (hun 100.0d0) (pt1 0.1d0) (two 2.0d0))
  21   (declare (type (double-float) half hun pt1 two))
  22   (defun dlhin
  23          (neq n t0 y0 ydot f tout uround ewt itol atol y temp h0 niter ier)
  24     (declare (type (array double-float (*)) temp y atol ewt ydot y0)
  25              (type (double-float) h0 uround tout t0)
  26              (type (f2cl-lib:integer4) ier niter itol n)
  27              (type (array f2cl-lib:integer4 (*)) neq))
  28     (f2cl-lib:with-multi-array-data
  29         ((neq f2cl-lib:integer4 neq-%data% neq-%offset%)
  30          (y0 double-float y0-%data% y0-%offset%)
  31          (ydot double-float ydot-%data% ydot-%offset%)
  32          (ewt double-float ewt-%data% ewt-%offset%)
  33          (atol double-float atol-%data% atol-%offset%)
  34          (y double-float y-%data% y-%offset%)
  35          (temp double-float temp-%data% temp-%offset%))
  36       (prog ((i 0) (iter 0) (afi 0.0d0) (atoli 0.0d0) (delyi 0.0d0) (hg 0.0d0)
  37              (hlb 0.0d0) (hnew 0.0d0) (hrat 0.0d0) (hub 0.0d0) (t1 0.0d0)
  38              (tdist 0.0d0) (tround 0.0d0) (yddnrm 0.0d0))
  39         (declare (type (double-float) yddnrm tround tdist t1 hub hrat hnew hlb
  40                                       hg delyi atoli afi)
  41                  (type (f2cl-lib:integer4) iter i))
  42         (setf niter 0)
  43         (setf tdist (abs (- tout t0)))
  44         (setf tround (* uround (max (abs t0) (abs tout))))
  45         (if (< tdist (* two tround)) (go label100))
  46         (setf hlb (* hun tround))
  47         (setf hub (* pt1 tdist))
  48         (setf atoli (f2cl-lib:fref atol-%data% (1) ((1 *)) atol-%offset%))
  49         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  50                       ((> i n) nil)
  51           (tagbody
  52             (if (or (= itol 2) (= itol 4))
  53                 (setf atoli
  54                         (f2cl-lib:fref atol-%data% (i) ((1 *)) atol-%offset%)))
  55             (setf delyi
  56                     (+
  57                      (* pt1
  58                         (abs
  59                          (f2cl-lib:fref y0-%data% (i) ((1 *)) y0-%offset%)))
  60                      atoli))
  61             (setf afi
  62                     (abs
  63                      (f2cl-lib:fref ydot-%data% (i) ((1 *)) ydot-%offset%)))
  64             (if (> (* afi hub) delyi) (setf hub (/ delyi afi)))
  65            label10))
  66         (setf iter 0)
  67         (setf hg (f2cl-lib:fsqrt (* hlb hub)))
  68         (cond
  69           ((< hub hlb)
  70            (setf h0 hg)
  71            (go label90)))
  72        label50
  73         (setf t1 (+ t0 hg))
  74         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  75                       ((> i n) nil)
  76           (tagbody
  77            label60
  78             (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
  79                     (+ (f2cl-lib:fref y0-%data% (i) ((1 *)) y0-%offset%)
  80                        (* hg
  81                           (f2cl-lib:fref ydot-%data%
  82                                          (i)
  83                                          ((1 *))
  84                                          ydot-%offset%))))))
  85         (multiple-value-bind (var-0 var-1 var-2 var-3)
  86             (funcall f neq t1 y temp)
  87           (declare (ignore var-0 var-2 var-3))
  88           (when var-1
  89             (setf t1 var-1)))
  90         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  91                       ((> i n) nil)
  92           (tagbody
  93            label70
  94             (setf (f2cl-lib:fref temp-%data% (i) ((1 *)) temp-%offset%)
  95                     (/
  96                      (- (f2cl-lib:fref temp-%data% (i) ((1 *)) temp-%offset%)
  97                         (f2cl-lib:fref ydot-%data% (i) ((1 *)) ydot-%offset%))
  98                      hg))))
  99         (setf yddnrm (dvnorm n temp ewt))
 100         (cond
 101           ((> (* yddnrm hub hub) two)
 102            (setf hnew (f2cl-lib:fsqrt (/ two yddnrm))))
 103           (t
 104            (setf hnew (f2cl-lib:fsqrt (* hg hub)))))
 105         (setf iter (f2cl-lib:int-add iter 1))
 106         (if (>= iter 4) (go label80))
 107         (setf hrat (/ hnew hg))
 108         (if (and (> hrat half) (< hrat two)) (go label80))
 109         (cond
 110           ((and (>= iter 2) (> hnew (* two hg)))
 111            (setf hnew hg)
 112            (go label80)))
 113         (setf hg hnew)
 114         (go label50)
 115        label80
 116         (setf h0 (* hnew half))
 117         (if (< h0 hlb) (setf h0 hlb))
 118         (if (> h0 hub) (setf h0 hub))
 119        label90
 120         (setf h0 (f2cl-lib:sign h0 (- tout t0)))
 121         (dcopy n y0 1 y 1)
 122         (setf niter iter)
 123         (setf ier 0)
 124         (go end_label)
 125        label100
 126         (setf ier -1)
 127         (go end_label)
 128        end_label
 129         (return
 130          (values nil
 131                  nil
 132                  nil
 133                  nil
 134                  nil
 135                  nil
 136                  nil
 137                  nil
 138                  nil
 139                  nil
 140                  nil
 141                  nil
 142                  nil
 143                  h0
 144                  niter
 145                  ier))))))
 146
 147 (in-package #:cl-user)
 148 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 149 (eval-when (:load-toplevel :compile-toplevel :execute)
 150   (setf (gethash 'fortran-to-lisp::dlhin fortran-to-lisp::*f2cl-function-info*)
 151           (fortran-to-lisp::make-f2cl-finfo
 152            :arg-types '((array fortran-to-lisp::integer4 (*))
 153                         (fortran-to-lisp::integer4) (double-float)
 154                         (array double-float (*)) (array double-float (*)) t
 155                         (double-float) (double-float) (array double-float (*))
 156                         (fortran-to-lisp::integer4) (array double-float (*))
 157                         (array double-float (*)) (array double-float (*))
 158                         (double-float) (fortran-to-lisp::integer4)
 159                         (fortran-to-lisp::integer4))
 160            :return-values '(nil nil nil nil nil nil nil nil nil nil nil nil nil
 161                             fortran-to-lisp::h0 fortran-to-lisp::niter
 162                             fortran-to-lisp::ier)
 163            :calls '(fortran-to-lisp::dcopy fortran-to-lisp::dvnorm))))
 164