share/odepack/src/droots.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 ((zero 0.0d0))
  21   (declare (type (double-float) zero))
  22   (defun droots (ng hmin jflag x0 x1 g0 g1 gx x jroot)
  23     (declare (type (array f2cl-lib:integer4 (*)) jroot)
  24              (type (array double-float (*)) gx g1 g0)
  25              (type (double-float) x x1 x0 hmin)
  26              (type (f2cl-lib:integer4) jflag ng))
  27     (let ()
  28       (symbol-macrolet ((alpha (aref (dlsr01-part-0 *dlsr01-common-block*) 0))
  29                         (x2 (aref (dlsr01-part-0 *dlsr01-common-block*) 1))
  30                         (imax (aref (dlsr01-part-1 *dlsr01-common-block*) 3))
  31                         (last$ (aref (dlsr01-part-1 *dlsr01-common-block*) 4)))
  32         (f2cl-lib:with-multi-array-data
  33             ((g0 double-float g0-%data% g0-%offset%)
  34              (g1 double-float g1-%data% g1-%offset%)
  35              (gx double-float gx-%data% gx-%offset%)
  36              (jroot f2cl-lib:integer4 jroot-%data% jroot-%offset%))
  37           (prog ((nxlast 0) (imxold 0) (i 0) (tmax 0.0d0) (t2 0.0d0)
  38                  (xroot nil) (sgnchg nil) (zroot nil))
  39             (declare (type f2cl-lib:logical zroot sgnchg xroot)
  40                      (type (double-float) t2 tmax)
  41                      (type (f2cl-lib:integer4) i imxold nxlast))
  42             (if (= jflag 1) (go label200))
  43             (setf imax 0)
  44             (setf tmax zero)
  45             (setf zroot f2cl-lib:%false%)
  46             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  47                           ((> i ng) nil)
  48               (tagbody
  49                 (if
  50                  (> (abs (f2cl-lib:fref g1-%data% (i) ((1 ng)) g1-%offset%))
  51                     zero)
  52                  (go label110))
  53                 (setf zroot f2cl-lib:%true%)
  54                 (go label120)
  55                label110
  56                 (if
  57                  (=
  58                   (f2cl-lib:sign 1.0d0
  59                                  (f2cl-lib:fref g0-%data%
  60                                                 (i)
  61                                                 ((1 ng))
  62                                                 g0-%offset%))
  63                   (f2cl-lib:sign 1.0d0
  64                                  (f2cl-lib:fref g1-%data%
  65                                                 (i)
  66                                                 ((1 ng))
  67                                                 g1-%offset%)))
  68                  (go label120))
  69                 (setf t2
  70                         (abs
  71                          (/ (f2cl-lib:fref g1-%data% (i) ((1 ng)) g1-%offset%)
  72                             (-
  73                              (f2cl-lib:fref g1-%data% (i) ((1 ng)) g1-%offset%)
  74                              (f2cl-lib:fref g0-%data%
  75                                             (i)
  76                                             ((1 ng))
  77                                             g0-%offset%)))))
  78                 (if (<= t2 tmax) (go label120))
  79                 (setf tmax t2)
  80                 (setf imax i)
  81                label120))
  82             (if (> imax 0) (go label130))
  83             (setf sgnchg f2cl-lib:%false%)
  84             (go label140)
  85            label130
  86             (setf sgnchg f2cl-lib:%true%)
  87            label140
  88             (if (not sgnchg) (go label400))
  89             (setf xroot f2cl-lib:%false%)
  90             (setf nxlast 0)
  91             (setf last$ 1)
  92            label150
  93             (if xroot (go label300))
  94             (if (= nxlast last$) (go label160))
  95             (setf alpha 1.0d0)
  96             (go label180)
  97            label160
  98             (if (= last$ 0) (go label170))
  99             (setf alpha (* 0.5d0 alpha))
 100             (go label180)
 101            label170
 102             (setf alpha (* 2.0d0 alpha))
 103            label180
 104             (setf x2
 105                     (+ x1
 106                        (/
 107                         (* (- (- x1 x0))
 108                            (f2cl-lib:fref g1-%data%
 109                                           (imax)
 110                                           ((1 ng))
 111                                           g1-%offset%))
 112                         (-
 113                          (f2cl-lib:fref g1-%data% (imax) ((1 ng)) g1-%offset%)
 114                          (* alpha
 115                             (f2cl-lib:fref g0-%data%
 116                                            (imax)
 117                                            ((1 ng))
 118                                            g0-%offset%))))))
 119             (if
 120              (and (< (abs (- x2 x0)) hmin) (> (abs (- x1 x0)) (* 10.0d0 hmin)))
 121              (setf x2 (+ x0 (* 0.1d0 (- x1 x0)))))
 122             (setf jflag 1)
 123             (setf x x2)
 124             (go end_label)
 125            label200
 126             (setf imxold imax)
 127             (setf imax 0)
 128             (setf tmax zero)
 129             (setf zroot f2cl-lib:%false%)
 130             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 131                           ((> i ng) nil)
 132               (tagbody
 133                 (if
 134                  (> (abs (f2cl-lib:fref gx-%data% (i) ((1 ng)) gx-%offset%))
 135                     zero)
 136                  (go label210))
 137                 (setf zroot f2cl-lib:%true%)
 138                 (go label220)
 139                label210
 140                 (if
 141                  (=
 142                   (f2cl-lib:sign 1.0d0
 143                                  (f2cl-lib:fref g0-%data%
 144                                                 (i)
 145                                                 ((1 ng))
 146                                                 g0-%offset%))
 147                   (f2cl-lib:sign 1.0d0
 148                                  (f2cl-lib:fref gx-%data%
 149                                                 (i)
 150                                                 ((1 ng))
 151                                                 gx-%offset%)))
 152                  (go label220))
 153                 (setf t2
 154                         (abs
 155                          (/ (f2cl-lib:fref gx-%data% (i) ((1 ng)) gx-%offset%)
 156                             (-
 157                              (f2cl-lib:fref gx-%data% (i) ((1 ng)) gx-%offset%)
 158                              (f2cl-lib:fref g0-%data%
 159                                             (i)
 160                                             ((1 ng))
 161                                             g0-%offset%)))))
 162                 (if (<= t2 tmax) (go label220))
 163                 (setf tmax t2)
 164                 (setf imax i)
 165                label220))
 166             (if (> imax 0) (go label230))
 167             (setf sgnchg f2cl-lib:%false%)
 168             (setf imax imxold)
 169             (go label240)
 170            label230
 171             (setf sgnchg f2cl-lib:%true%)
 172            label240
 173             (setf nxlast last$)
 174             (if (not sgnchg) (go label250))
 175             (setf x1 x2)
 176             (dcopy ng gx 1 g1 1)
 177             (setf last$ 1)
 178             (setf xroot f2cl-lib:%false%)
 179             (go label270)
 180            label250
 181             (if (not zroot) (go label260))
 182             (setf x1 x2)
 183             (dcopy ng gx 1 g1 1)
 184             (setf xroot f2cl-lib:%true%)
 185             (go label270)
 186            label260
 187             (dcopy ng gx 1 g0 1)
 188             (setf x0 x2)
 189             (setf last$ 0)
 190             (setf xroot f2cl-lib:%false%)
 191            label270
 192             (if (<= (abs (- x1 x0)) hmin) (setf xroot f2cl-lib:%true%))
 193             (go label150)
 194            label300
 195             (setf jflag 2)
 196             (setf x x1)
 197             (dcopy ng g1 1 gx 1)
 198             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 199                           ((> i ng) nil)
 200               (tagbody
 201                 (setf (f2cl-lib:fref jroot-%data% (i) ((1 ng)) jroot-%offset%)
 202                         0)
 203                 (if
 204                  (> (abs (f2cl-lib:fref g1-%data% (i) ((1 ng)) g1-%offset%))
 205                     zero)
 206                  (go label310))
 207                 (setf (f2cl-lib:fref jroot-%data% (i) ((1 ng)) jroot-%offset%)
 208                         1)
 209                 (go label320)
 210                label310
 211                 (if
 212                  (/=
 213                   (f2cl-lib:sign 1.0d0
 214                                  (f2cl-lib:fref g0-%data%
 215                                                 (i)
 216                                                 ((1 ng))
 217                                                 g0-%offset%))
 218                   (f2cl-lib:sign 1.0d0
 219                                  (f2cl-lib:fref g1-%data%
 220                                                 (i)
 221                                                 ((1 ng))
 222                                                 g1-%offset%)))
 223                  (setf (f2cl-lib:fref jroot-%data% (i) ((1 ng)) jroot-%offset%)
 224                          1))
 225                label320))
 226             (go end_label)
 227            label400
 228             (if (not zroot) (go label420))
 229             (setf x x1)
 230             (dcopy ng g1 1 gx 1)
 231             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 232                           ((> i ng) nil)
 233               (tagbody
 234                 (setf (f2cl-lib:fref jroot-%data% (i) ((1 ng)) jroot-%offset%)
 235                         0)
 236                 (if
 237                  (<= (abs (f2cl-lib:fref g1-%data% (i) ((1 ng)) g1-%offset%))
 238                      zero)
 239                  (setf (f2cl-lib:fref jroot-%data% (i) ((1 ng)) jroot-%offset%)
 240                          1))
 241                label410))
 242             (setf jflag 3)
 243             (go end_label)
 244            label420
 245             (dcopy ng g1 1 gx 1)
 246             (setf x x1)
 247             (setf jflag 4)
 248             (go end_label)
 249            end_label
 250             (return (values nil nil jflag x0 x1 nil nil nil x nil))))))))
 251
 252 (in-package #:cl-user)
 253 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 254 (eval-when (:load-toplevel :compile-toplevel :execute)
 255   (setf (gethash 'fortran-to-lisp::droots
 256                  fortran-to-lisp::*f2cl-function-info*)
 257           (fortran-to-lisp::make-f2cl-finfo
 258            :arg-types '((fortran-to-lisp::integer4) (double-float)
 259                         (fortran-to-lisp::integer4) (double-float)
 260                         (double-float) (array double-float (*))
 261                         (array double-float (*)) (array double-float (*))
 262                         (double-float) (array fortran-to-lisp::integer4 (*)))
 263            :return-values '(nil nil fortran-to-lisp::jflag fortran-to-lisp::x0
 264                             fortran-to-lisp::x1 nil nil nil fortran-to-lisp::x
 265                             nil)
 266            :calls '(fortran-to-lisp::dcopy))))
 267