share/hompack/lisp/stepds.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-2020-04 (21D 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 "HOMPACK")
  18
  19
  20 (let ((two
  21        (make-array 13
  22                    :element-type 'double-float
  23                    :initial-contents '(2.0 4.0 8.0 16.0 32.0 64.0 128.0 256.0
  24                                        512.0 1024.0 2048.0 4096.0 8192.0)))
  25       (gstr
  26        (make-array 13
  27                    :element-type 'double-float
  28                    :initial-contents '(0.5 0.0833 0.0417 0.0264 0.0188 0.0143
  29                                        0.0114 0.00936 0.00789 0.00679 0.00592
  30                                        0.00524 0.00468)))
  31       (psi (make-array 12 :element-type 'double-float))
  32       (beta (make-array 12 :element-type 'double-float))
  33       (sig (make-array 13 :element-type 'double-float))
  34       (v (make-array 12 :element-type 'double-float))
  35       (phase1 nil)
  36       (nornd nil)
  37       (i 0)
  38       (ifail 0)
  39       (im1 0)
  40       (ip1 0)
  41       (iq 0)
  42       (j 0)
  43       (jv 0)
  44       (km1 0)
  45       (km2 0)
  46       (knew 0)
  47       (kp1 0)
  48       (kp2 0)
  49       (kprev 0)
  50       (l 0)
  51       (limit1 0)
  52       (limit2 0)
  53       (ns 0)
  54       (nsm2 0)
  55       (nsp1 0)
  56       (nsp2 0)
  57       (absh 0.0)
  58       (erk 0.0)
  59       (erkm1 0.0)
  60       (erkm2 0.0)
  61       (erkp1 0.0)
  62       (err 0.0)
  63       (fouru 0.0)
  64       (hnew 0.0)
  65       (p5eps 0.0)
  66       (r 0.0)
  67       (reali 0.0)
  68       (realns 0.0)
  69       (rho 0.0)
  70       (round$ 0.0)
  71       (sum 0.0)
  72       (tau 0.0)
  73       (temp1 0.0)
  74       (temp2 0.0)
  75       (temp3 0.0)
  76       (temp4 0.0)
  77       (temp5 0.0)
  78       (temp6 0.0)
  79       (twou 0.0))
  80   (declare (type f2cl-lib:logical phase1 nornd)
  81            (type (f2cl-lib:integer4) i ifail im1 ip1 iq j jv km1 km2 knew kp1
  82                                      kp2 kprev l limit1 limit2 ns nsm2 nsp1
  83                                      nsp2)
  84            (type (double-float) absh erk erkm1 erkm2 erkp1 err fouru hnew p5eps
  85                                 r reali realns rho round$ sum tau temp1 temp2
  86                                 temp3 temp4 temp5 temp6 twou)
  87            (type (array double-float (12)) psi beta v)
  88            (type (array double-float (13)) two gstr sig))
  89   (defun stepds
  90          (f neqn y x h eps wt start hold k kold crash phi p yp alpha w g ksteps
  91           xold ivc iv kgi gi fpwa1 fpwa2 fpwa3 ifpc1 ifpwa1 fpwa4 fpwa5 ifpc2
  92           ifpc3 par ipar)
  93     (declare (type (array f2cl-lib:integer4 (*)) ipar)
  94              (type (array double-float (*)) par)
  95              (type (array f2cl-lib:integer4 (*)) ifpwa1)
  96              (type (array double-float (*)) gi)
  97              (type (array f2cl-lib:integer4 (*)) iv)
  98              (type (array double-float (*)) g)
  99              (type (array double-float (*)) w alpha)
 100              (type f2cl-lib:logical crash start)
 101              (type (double-float) xold hold eps h x)
 102              (type (array double-float (*)) fpwa5 fpwa4 fpwa3 fpwa2 fpwa1 yp p
 103                                             phi wt y)
 104              (type (f2cl-lib:integer4) ifpc3 ifpc2 ifpc1 kgi ivc ksteps kold k
 105                                        neqn))
 106     (f2cl-lib:with-multi-array-data
 107         ((y double-float y-%data% y-%offset%)
 108          (wt double-float wt-%data% wt-%offset%)
 109          (phi double-float phi-%data% phi-%offset%)
 110          (p double-float p-%data% p-%offset%)
 111          (yp double-float yp-%data% yp-%offset%)
 112          (fpwa1 double-float fpwa1-%data% fpwa1-%offset%)
 113          (fpwa2 double-float fpwa2-%data% fpwa2-%offset%)
 114          (fpwa3 double-float fpwa3-%data% fpwa3-%offset%)
 115          (fpwa4 double-float fpwa4-%data% fpwa4-%offset%)
 116          (fpwa5 double-float fpwa5-%data% fpwa5-%offset%)
 117          (alpha double-float alpha-%data% alpha-%offset%)
 118          (w double-float w-%data% w-%offset%)
 119          (g double-float g-%data% g-%offset%)
 120          (iv f2cl-lib:integer4 iv-%data% iv-%offset%)
 121          (gi double-float gi-%data% gi-%offset%)
 122          (ifpwa1 f2cl-lib:integer4 ifpwa1-%data% ifpwa1-%offset%)
 123          (par double-float par-%data% par-%offset%)
 124          (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
 125       (prog ()
 126         (declare)
 127         (setf twou (* 2.0f0 (f2cl-lib:d1mach 4)))
 128         (setf fouru (+ twou twou))
 129         (setf crash f2cl-lib:%true%)
 130         (if (>= (abs h) (* fouru (abs x))) (go label5))
 131         (setf h (f2cl-lib:sign (* fouru (abs x)) h))
 132         (go end_label)
 133        label5
 134         (setf p5eps (* 0.5f0 eps))
 135         (setf round$ (coerce 0.0f0 'double-float))
 136         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 137                       ((> l neqn) nil)
 138           (tagbody
 139            label10
 140             (setf round$
 141                     (+ round$
 142                        (expt
 143                         (/ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 144                            (f2cl-lib:fref wt-%data%
 145                                           (l)
 146                                           ((1 neqn))
 147                                           wt-%offset%))
 148                         2)))))
 149         (setf round$ (* twou (f2cl-lib:fsqrt round$)))
 150         (if (>= p5eps round$) (go label15))
 151         (setf eps (* 2.0f0 round$ (+ 1.0f0 fouru)))
 152         (go end_label)
 153        label15
 154         (setf crash f2cl-lib:%false%)
 155         (setf (f2cl-lib:fref g-%data% (1) ((1 13)) g-%offset%)
 156                 (coerce 1.0f0 'double-float))
 157         (setf (f2cl-lib:fref g-%data% (2) ((1 13)) g-%offset%)
 158                 (coerce 0.5f0 'double-float))
 159         (setf (f2cl-lib:fref sig (1) ((1 13))) (coerce 1.0f0 'double-float))
 160         (if (not start) (go label99))
 161         (multiple-value-bind
 162               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 163                var-10 var-11 var-12 var-13 var-14)
 164             (funcall f
 165                      x
 166                      y
 167                      yp
 168                      fpwa1
 169                      fpwa2
 170                      fpwa3
 171                      ifpc1
 172                      ifpwa1
 173                      fpwa4
 174                      fpwa5
 175                      ifpc2
 176                      (f2cl-lib:int-sub neqn 1)
 177                      ifpc3
 178                      par
 179                      ipar)
 180           (declare (ignore var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9
 181                            var-11 var-13 var-14))
 182           (when var-0
 183             (setf x var-0))
 184           (when var-6
 185             (setf ifpc1 var-6))
 186           (when var-10
 187             (setf ifpc2 var-10))
 188           (when var-12
 189             (setf ifpc3 var-12)))
 190         (if (> ifpc3 0) (go end_label))
 191         (setf sum (coerce 0.0f0 'double-float))
 192         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 193                       ((> l neqn) nil)
 194           (tagbody
 195             (setf (f2cl-lib:fref phi-%data%
 196                                  (l 1)
 197                                  ((1 neqn) (1 16))
 198                                  phi-%offset%)
 199                     (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%))
 200             (setf (f2cl-lib:fref phi-%data%
 201                                  (l 2)
 202                                  ((1 neqn) (1 16))
 203                                  phi-%offset%)
 204                     (coerce 0.0f0 'double-float))
 205            label20
 206             (setf sum
 207                     (+ sum
 208                        (expt
 209                         (/ (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
 210                            (f2cl-lib:fref wt-%data%
 211                                           (l)
 212                                           ((1 neqn))
 213                                           wt-%offset%))
 214                         2)))))
 215         (setf sum (f2cl-lib:fsqrt sum))
 216         (setf absh (abs h))
 217         (if (< eps (* 16.0f0 sum h h))
 218             (setf absh (* 0.25f0 (f2cl-lib:fsqrt (/ eps sum)))))
 219         (setf h (f2cl-lib:sign (max absh (* fouru (abs x))) h))
 220         (setf hold (coerce 0.0f0 'double-float))
 221         (setf k 1)
 222         (setf kold 0)
 223         (setf kprev 0)
 224         (setf start f2cl-lib:%false%)
 225         (setf phase1 f2cl-lib:%true%)
 226         (setf nornd f2cl-lib:%true%)
 227         (if (> p5eps (* 100.0f0 round$)) (go label99))
 228         (setf nornd f2cl-lib:%false%)
 229         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 230                       ((> l neqn) nil)
 231           (tagbody
 232            label25
 233             (setf (f2cl-lib:fref phi-%data%
 234                                  (l 15)
 235                                  ((1 neqn) (1 16))
 236                                  phi-%offset%)
 237                     (coerce 0.0f0 'double-float))))
 238        label99
 239         (setf ifail 0)
 240        label100
 241         (setf kp1 (f2cl-lib:int-add k 1))
 242         (setf kp2 (f2cl-lib:int-add k 2))
 243         (setf km1 (f2cl-lib:int-sub k 1))
 244         (setf km2 (f2cl-lib:int-sub k 2))
 245         (if (/= h hold) (setf ns 0))
 246         (if (<= ns kold) (setf ns (f2cl-lib:int-add ns 1)))
 247         (setf nsp1 (f2cl-lib:int-add ns 1))
 248         (if (< k ns) (go label199))
 249         (setf (f2cl-lib:fref beta (ns) ((1 12))) (coerce 1.0f0 'double-float))
 250         (setf realns (coerce (the f2cl-lib:integer4 ns) 'double-float))
 251         (setf (f2cl-lib:fref alpha-%data% (ns) ((1 12)) alpha-%offset%)
 252                 (/ 1.0f0 realns))
 253         (setf temp1 (* h realns))
 254         (setf (f2cl-lib:fref sig (nsp1) ((1 13))) (coerce 1.0f0 'double-float))
 255         (if (< k nsp1) (go label110))
 256         (f2cl-lib:fdo (i nsp1 (f2cl-lib:int-add i 1))
 257                       ((> i k) nil)
 258           (tagbody
 259             (setf im1 (f2cl-lib:int-sub i 1))
 260             (setf temp2 (f2cl-lib:fref psi (im1) ((1 12))))
 261             (setf (f2cl-lib:fref psi (im1) ((1 12))) temp1)
 262             (setf (f2cl-lib:fref beta (i) ((1 12)))
 263                     (/
 264                      (* (f2cl-lib:fref beta (im1) ((1 12)))
 265                         (f2cl-lib:fref psi (im1) ((1 12))))
 266                      temp2))
 267             (setf temp1 (+ temp2 h))
 268             (setf (f2cl-lib:fref alpha-%data% (i) ((1 12)) alpha-%offset%)
 269                     (/ h temp1))
 270             (setf reali (coerce (the f2cl-lib:integer4 i) 'double-float))
 271            label105
 272             (setf (f2cl-lib:fref sig ((f2cl-lib:int-add i 1)) ((1 13)))
 273                     (* reali
 274                        (f2cl-lib:fref alpha-%data% (i) ((1 12)) alpha-%offset%)
 275                        (f2cl-lib:fref sig (i) ((1 13)))))))
 276        label110
 277         (setf (f2cl-lib:fref psi (k) ((1 12))) temp1)
 278         (if (> ns 1) (go label120))
 279         (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
 280                       ((> iq k) nil)
 281           (tagbody
 282             (setf temp3
 283                     (coerce
 284                      (the f2cl-lib:integer4
 285                           (f2cl-lib:int-mul iq (f2cl-lib:int-add iq 1)))
 286                      'double-float))
 287             (setf (f2cl-lib:fref v (iq) ((1 12))) (/ 1.0f0 temp3))
 288            label115
 289             (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
 290                     (f2cl-lib:fref v (iq) ((1 12))))))
 291         (setf ivc 0)
 292         (setf kgi 0)
 293         (if (= k 1) (go label140))
 294         (setf kgi 1)
 295         (setf (f2cl-lib:fref gi-%data% (1) ((1 11)) gi-%offset%)
 296                 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
 297         (go label140)
 298        label120
 299         (if (<= k kprev) (go label130))
 300         (if (= ivc 0) (go label122))
 301         (setf jv
 302                 (f2cl-lib:int-sub kp1
 303                                   (f2cl-lib:fref iv-%data%
 304                                                  (ivc)
 305                                                  ((1 10))
 306                                                  iv-%offset%)))
 307         (setf ivc (f2cl-lib:int-sub ivc 1))
 308         (go label123)
 309        label122
 310         (setf jv 1)
 311         (setf temp4
 312                 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-mul k kp1))
 313                         'double-float))
 314         (setf (f2cl-lib:fref v (k) ((1 12))) (/ 1.0f0 temp4))
 315         (setf (f2cl-lib:fref w-%data% (k) ((1 12)) w-%offset%)
 316                 (f2cl-lib:fref v (k) ((1 12))))
 317         (if (/= k 2) (go label123))
 318         (setf kgi 1)
 319         (setf (f2cl-lib:fref gi-%data% (1) ((1 11)) gi-%offset%)
 320                 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
 321        label123
 322         (setf nsm2 (f2cl-lib:int-sub ns 2))
 323         (if (< nsm2 jv) (go label130))
 324         (f2cl-lib:fdo (j jv (f2cl-lib:int-add j 1))
 325                       ((> j nsm2) nil)
 326           (tagbody
 327             (setf i (f2cl-lib:int-sub k j))
 328             (setf (f2cl-lib:fref v (i) ((1 12)))
 329                     (- (f2cl-lib:fref v (i) ((1 12)))
 330                        (*
 331                         (f2cl-lib:fref alpha-%data%
 332                                        ((f2cl-lib:int-add j 1))
 333                                        ((1 12))
 334                                        alpha-%offset%)
 335                         (f2cl-lib:fref v ((f2cl-lib:int-add i 1)) ((1 12))))))
 336            label125
 337             (setf (f2cl-lib:fref w-%data% (i) ((1 12)) w-%offset%)
 338                     (f2cl-lib:fref v (i) ((1 12))))))
 339         (if (/= i 2) (go label130))
 340         (setf kgi (f2cl-lib:int-sub ns 1))
 341         (setf (f2cl-lib:fref gi-%data% (kgi) ((1 11)) gi-%offset%)
 342                 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
 343        label130
 344         (setf limit1 (f2cl-lib:int-sub kp1 ns))
 345         (setf temp5 (f2cl-lib:fref alpha-%data% (ns) ((1 12)) alpha-%offset%))
 346         (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
 347                       ((> iq limit1) nil)
 348           (tagbody
 349             (setf (f2cl-lib:fref v (iq) ((1 12)))
 350                     (- (f2cl-lib:fref v (iq) ((1 12)))
 351                        (* temp5
 352                           (f2cl-lib:fref v
 353                                          ((f2cl-lib:int-add iq 1))
 354                                          ((1 12))))))
 355            label135
 356             (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
 357                     (f2cl-lib:fref v (iq) ((1 12))))))
 358         (setf (f2cl-lib:fref g-%data% (nsp1) ((1 13)) g-%offset%)
 359                 (f2cl-lib:fref w-%data% (1) ((1 12)) w-%offset%))
 360         (if (= limit1 1) (go label137))
 361         (setf kgi ns)
 362         (setf (f2cl-lib:fref gi-%data% (kgi) ((1 11)) gi-%offset%)
 363                 (f2cl-lib:fref w-%data% (2) ((1 12)) w-%offset%))
 364        label137
 365         (setf (f2cl-lib:fref w-%data%
 366                              ((f2cl-lib:int-add limit1 1))
 367                              ((1 12))
 368                              w-%offset%)
 369                 (f2cl-lib:fref v ((f2cl-lib:int-add limit1 1)) ((1 12))))
 370         (if (>= k kold) (go label140))
 371         (setf ivc (f2cl-lib:int-add ivc 1))
 372         (setf (f2cl-lib:fref iv-%data% (ivc) ((1 10)) iv-%offset%)
 373                 (f2cl-lib:int-add limit1 2))
 374        label140
 375         (setf nsp2 (f2cl-lib:int-add ns 2))
 376         (setf kprev k)
 377         (if (< kp1 nsp2) (go label199))
 378         (f2cl-lib:fdo (i nsp2 (f2cl-lib:int-add i 1))
 379                       ((> i kp1) nil)
 380           (tagbody
 381             (setf limit2 (f2cl-lib:int-sub kp2 i))
 382             (setf temp6
 383                     (f2cl-lib:fref alpha-%data%
 384                                    ((f2cl-lib:int-sub i 1))
 385                                    ((1 12))
 386                                    alpha-%offset%))
 387             (f2cl-lib:fdo (iq 1 (f2cl-lib:int-add iq 1))
 388                           ((> iq limit2) nil)
 389               (tagbody
 390                label145
 391                 (setf (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
 392                         (- (f2cl-lib:fref w-%data% (iq) ((1 12)) w-%offset%)
 393                            (* temp6
 394                               (f2cl-lib:fref w-%data%
 395                                              ((f2cl-lib:int-add iq 1))
 396                                              ((1 12))
 397                                              w-%offset%))))))
 398            label150
 399             (setf (f2cl-lib:fref g-%data% (i) ((1 13)) g-%offset%)
 400                     (f2cl-lib:fref w-%data% (1) ((1 12)) w-%offset%))))
 401        label199
 402         (setf ksteps (f2cl-lib:int-add ksteps 1))
 403         (if (< k nsp1) (go label215))
 404         (f2cl-lib:fdo (i nsp1 (f2cl-lib:int-add i 1))
 405                       ((> i k) nil)
 406           (tagbody
 407             (setf temp1 (f2cl-lib:fref beta (i) ((1 12))))
 408             (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 409                           ((> l neqn) nil)
 410               (tagbody
 411                label205
 412                 (setf (f2cl-lib:fref phi-%data%
 413                                      (l i)
 414                                      ((1 neqn) (1 16))
 415                                      phi-%offset%)
 416                         (* temp1
 417                            (f2cl-lib:fref phi-%data%
 418                                           (l i)
 419                                           ((1 neqn) (1 16))
 420                                           phi-%offset%)))))
 421            label210))
 422        label215
 423         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 424                       ((> l neqn) nil)
 425           (tagbody
 426             (setf (f2cl-lib:fref phi-%data%
 427                                  (l kp2)
 428                                  ((1 neqn) (1 16))
 429                                  phi-%offset%)
 430                     (f2cl-lib:fref phi-%data%
 431                                    (l kp1)
 432                                    ((1 neqn) (1 16))
 433                                    phi-%offset%))
 434             (setf (f2cl-lib:fref phi-%data%
 435                                  (l kp1)
 436                                  ((1 neqn) (1 16))
 437                                  phi-%offset%)
 438                     (coerce 0.0f0 'double-float))
 439            label220
 440             (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 441                     (coerce 0.0f0 'double-float))))
 442         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 443                       ((> j k) nil)
 444           (tagbody
 445             (setf i (f2cl-lib:int-sub kp1 j))
 446             (setf ip1 (f2cl-lib:int-add i 1))
 447             (setf temp2 (f2cl-lib:fref g-%data% (i) ((1 13)) g-%offset%))
 448             (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 449                           ((> l neqn) nil)
 450               (tagbody
 451                 (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 452                         (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 453                            (* temp2
 454                               (f2cl-lib:fref phi-%data%
 455                                              (l i)
 456                                              ((1 neqn) (1 16))
 457                                              phi-%offset%))))
 458                label225
 459                 (setf (f2cl-lib:fref phi-%data%
 460                                      (l i)
 461                                      ((1 neqn) (1 16))
 462                                      phi-%offset%)
 463                         (+
 464                          (f2cl-lib:fref phi-%data%
 465                                         (l i)
 466                                         ((1 neqn) (1 16))
 467                                         phi-%offset%)
 468                          (f2cl-lib:fref phi-%data%
 469                                         (l ip1)
 470                                         ((1 neqn) (1 16))
 471                                         phi-%offset%)))))
 472            label230))
 473         (if nornd (go label240))
 474         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 475                       ((> l neqn) nil)
 476           (tagbody
 477             (setf tau
 478                     (- (* h (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%))
 479                        (f2cl-lib:fref phi-%data%
 480                                       (l 15)
 481                                       ((1 neqn) (1 16))
 482                                       phi-%offset%)))
 483             (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 484                     (+ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%) tau))
 485            label235
 486             (setf (f2cl-lib:fref phi-%data%
 487                                  (l 16)
 488                                  ((1 neqn) (1 16))
 489                                  phi-%offset%)
 490                     (- (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 491                        (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 492                        tau))))
 493         (go label250)
 494        label240
 495         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 496                       ((> l neqn) nil)
 497           (tagbody
 498            label245
 499             (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 500                     (+ (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 501                        (* h
 502                           (f2cl-lib:fref p-%data%
 503                                          (l)
 504                                          ((1 neqn))
 505                                          p-%offset%))))))
 506        label250
 507         (setf xold x)
 508         (setf x (+ x h))
 509         (setf absh (abs h))
 510         (multiple-value-bind
 511               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 512                var-10 var-11 var-12 var-13 var-14)
 513             (funcall f
 514                      x
 515                      p
 516                      yp
 517                      fpwa1
 518                      fpwa2
 519                      fpwa3
 520                      ifpc1
 521                      ifpwa1
 522                      fpwa4
 523                      fpwa5
 524                      ifpc2
 525                      (f2cl-lib:int-sub neqn 1)
 526                      ifpc3
 527                      par
 528                      ipar)
 529           (declare (ignore var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9
 530                            var-11 var-13 var-14))
 531           (when var-0
 532             (setf x var-0))
 533           (when var-6
 534             (setf ifpc1 var-6))
 535           (when var-10
 536             (setf ifpc2 var-10))
 537           (when var-12
 538             (setf ifpc3 var-12)))
 539         (if (> ifpc3 0) (go end_label))
 540         (setf erkm2 (coerce 0.0f0 'double-float))
 541         (setf erkm1 (coerce 0.0f0 'double-float))
 542         (setf erk (coerce 0.0f0 'double-float))
 543         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 544                       ((> l neqn) nil)
 545           (tagbody
 546             (setf temp3
 547                     (/ 1.0f0
 548                        (f2cl-lib:fref wt-%data% (l) ((1 neqn)) wt-%offset%)))
 549             (setf temp4
 550                     (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
 551                        (f2cl-lib:fref phi-%data%
 552                                       (l 1)
 553                                       ((1 neqn) (1 16))
 554                                       phi-%offset%)))
 555             (f2cl-lib:arithmetic-if km2
 556                                     (go label265)
 557                                     (go label260)
 558                                     (go label255))
 559            label255
 560             (setf erkm2
 561                     (+ erkm2
 562                        (expt
 563                         (*
 564                          (+
 565                           (f2cl-lib:fref phi-%data%
 566                                          (l km1)
 567                                          ((1 neqn) (1 16))
 568                                          phi-%offset%)
 569                           temp4)
 570                          temp3)
 571                         2)))
 572            label260
 573             (setf erkm1
 574                     (+ erkm1
 575                        (expt
 576                         (*
 577                          (+
 578                           (f2cl-lib:fref phi-%data%
 579                                          (l k)
 580                                          ((1 neqn) (1 16))
 581                                          phi-%offset%)
 582                           temp4)
 583                          temp3)
 584                         2)))
 585            label265
 586             (setf erk (+ erk (expt (* temp4 temp3) 2)))))
 587         (f2cl-lib:arithmetic-if km2 (go label280) (go label275) (go label270))
 588        label270
 589         (setf erkm2
 590                 (* absh
 591                    (f2cl-lib:fref sig (km1) ((1 13)))
 592                    (f2cl-lib:fref gstr (km2) ((1 13)))
 593                    (f2cl-lib:fsqrt erkm2)))
 594        label275
 595         (setf erkm1
 596                 (* absh
 597                    (f2cl-lib:fref sig (k) ((1 13)))
 598                    (f2cl-lib:fref gstr (km1) ((1 13)))
 599                    (f2cl-lib:fsqrt erkm1)))
 600        label280
 601         (setf temp5 (* absh (f2cl-lib:fsqrt erk)))
 602         (setf err
 603                 (* temp5
 604                    (- (f2cl-lib:fref g-%data% (k) ((1 13)) g-%offset%)
 605                       (f2cl-lib:fref g-%data% (kp1) ((1 13)) g-%offset%))))
 606         (setf erk
 607                 (* temp5
 608                    (f2cl-lib:fref sig (kp1) ((1 13)))
 609                    (f2cl-lib:fref gstr (k) ((1 13)))))
 610         (setf knew k)
 611         (f2cl-lib:arithmetic-if km2 (go label299) (go label290) (go label285))
 612        label285
 613         (if (<= (max erkm1 erkm2) erk) (setf knew km1))
 614         (go label299)
 615        label290
 616         (if (<= erkm1 (* 0.5f0 erk)) (setf knew km1))
 617        label299
 618         (if (<= err eps) (go label400))
 619         (setf phase1 f2cl-lib:%false%)
 620         (setf x xold)
 621         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 622                       ((> i k) nil)
 623           (tagbody
 624             (setf temp1 (/ 1.0f0 (f2cl-lib:fref beta (i) ((1 12)))))
 625             (setf ip1 (f2cl-lib:int-add i 1))
 626             (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 627                           ((> l neqn) nil)
 628               (tagbody
 629                label305
 630                 (setf (f2cl-lib:fref phi-%data%
 631                                      (l i)
 632                                      ((1 neqn) (1 16))
 633                                      phi-%offset%)
 634                         (* temp1
 635                            (-
 636                             (f2cl-lib:fref phi-%data%
 637                                            (l i)
 638                                            ((1 neqn) (1 16))
 639                                            phi-%offset%)
 640                             (f2cl-lib:fref phi-%data%
 641                                            (l ip1)
 642                                            ((1 neqn) (1 16))
 643                                            phi-%offset%))))))
 644            label310))
 645         (if (< k 2) (go label320))
 646         (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
 647                       ((> i k) nil)
 648           (tagbody
 649            label315
 650             (setf (f2cl-lib:fref psi ((f2cl-lib:int-sub i 1)) ((1 12)))
 651                     (- (f2cl-lib:fref psi (i) ((1 12))) h))))
 652        label320
 653         (setf ifail (f2cl-lib:int-add ifail 1))
 654         (setf temp2 (coerce 0.5f0 'double-float))
 655         (f2cl-lib:arithmetic-if (f2cl-lib:int-sub ifail 3)
 656                                 (go label335)
 657                                 (go label330)
 658                                 (go label325))
 659        label325
 660         (if (< p5eps (* 0.25f0 erk))
 661             (setf temp2 (f2cl-lib:fsqrt (/ p5eps erk))))
 662        label330
 663         (setf knew 1)
 664        label335
 665         (setf h (* temp2 h))
 666         (setf k knew)
 667         (setf ns 0)
 668         (if (>= (abs h) (* fouru (abs x))) (go label340))
 669         (setf crash f2cl-lib:%true%)
 670         (setf h (f2cl-lib:sign (* fouru (abs x)) h))
 671         (setf eps (+ eps eps))
 672         (go end_label)
 673        label340
 674         (go label100)
 675        label400
 676         (setf kold k)
 677         (setf hold h)
 678         (setf temp1 (* h (f2cl-lib:fref g-%data% (kp1) ((1 13)) g-%offset%)))
 679         (if nornd (go label410))
 680         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 681                       ((> l neqn) nil)
 682           (tagbody
 683             (setf temp3 (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%))
 684             (setf rho
 685                     (-
 686                      (* temp1
 687                         (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
 688                            (f2cl-lib:fref phi-%data%
 689                                           (l 1)
 690                                           ((1 neqn) (1 16))
 691                                           phi-%offset%)))
 692                      (f2cl-lib:fref phi-%data%
 693                                     (l 16)
 694                                     ((1 neqn) (1 16))
 695                                     phi-%offset%)))
 696             (setf (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 697                     (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) rho))
 698             (setf (f2cl-lib:fref phi-%data%
 699                                  (l 15)
 700                                  ((1 neqn) (1 16))
 701                                  phi-%offset%)
 702                     (- (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 703                        (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 704                        rho))
 705            label405
 706             (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) temp3)))
 707         (go label420)
 708        label410
 709         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 710                       ((> l neqn) nil)
 711           (tagbody
 712             (setf temp3 (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%))
 713             (setf (f2cl-lib:fref y-%data% (l) ((1 neqn)) y-%offset%)
 714                     (+ (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%)
 715                        (* temp1
 716                           (-
 717                            (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
 718                            (f2cl-lib:fref phi-%data%
 719                                           (l 1)
 720                                           ((1 neqn) (1 16))
 721                                           phi-%offset%)))))
 722            label415
 723             (setf (f2cl-lib:fref p-%data% (l) ((1 neqn)) p-%offset%) temp3)))
 724        label420
 725         (multiple-value-bind
 726               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 727                var-10 var-11 var-12 var-13 var-14)
 728             (funcall f
 729                      x
 730                      y
 731                      yp
 732                      fpwa1
 733                      fpwa2
 734                      fpwa3
 735                      ifpc1
 736                      ifpwa1
 737                      fpwa4
 738                      fpwa5
 739                      ifpc2
 740                      (f2cl-lib:int-sub neqn 1)
 741                      ifpc3
 742                      par
 743                      ipar)
 744           (declare (ignore var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9
 745                            var-11 var-13 var-14))
 746           (when var-0
 747             (setf x var-0))
 748           (when var-6
 749             (setf ifpc1 var-6))
 750           (when var-10
 751             (setf ifpc2 var-10))
 752           (when var-12
 753             (setf ifpc3 var-12)))
 754         (if (> ifpc3 0) (go end_label))
 755         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 756                       ((> l neqn) nil)
 757           (tagbody
 758             (setf (f2cl-lib:fref phi-%data%
 759                                  (l kp1)
 760                                  ((1 neqn) (1 16))
 761                                  phi-%offset%)
 762                     (- (f2cl-lib:fref yp-%data% (l) ((1 neqn)) yp-%offset%)
 763                        (f2cl-lib:fref phi-%data%
 764                                       (l 1)
 765                                       ((1 neqn) (1 16))
 766                                       phi-%offset%)))
 767            label425
 768             (setf (f2cl-lib:fref phi-%data%
 769                                  (l kp2)
 770                                  ((1 neqn) (1 16))
 771                                  phi-%offset%)
 772                     (-
 773                      (f2cl-lib:fref phi-%data%
 774                                     (l kp1)
 775                                     ((1 neqn) (1 16))
 776                                     phi-%offset%)
 777                      (f2cl-lib:fref phi-%data%
 778                                     (l kp2)
 779                                     ((1 neqn) (1 16))
 780                                     phi-%offset%)))))
 781         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 782                       ((> i k) nil)
 783           (tagbody
 784             (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 785                           ((> l neqn) nil)
 786               (tagbody
 787                label430
 788                 (setf (f2cl-lib:fref phi-%data%
 789                                      (l i)
 790                                      ((1 neqn) (1 16))
 791                                      phi-%offset%)
 792                         (+
 793                          (f2cl-lib:fref phi-%data%
 794                                         (l i)
 795                                         ((1 neqn) (1 16))
 796                                         phi-%offset%)
 797                          (f2cl-lib:fref phi-%data%
 798                                         (l kp1)
 799                                         ((1 neqn) (1 16))
 800                                         phi-%offset%)))))
 801            label435))
 802         (setf erkp1 (coerce 0.0f0 'double-float))
 803         (if (or (= knew km1) (= k 12)) (setf phase1 f2cl-lib:%false%))
 804         (if phase1 (go label450))
 805         (if (= knew km1) (go label455))
 806         (if (> kp1 ns) (go label460))
 807         (f2cl-lib:fdo (l 1 (f2cl-lib:int-add l 1))
 808                       ((> l neqn) nil)
 809           (tagbody
 810            label440
 811             (setf erkp1
 812                     (+ erkp1
 813                        (expt
 814                         (/
 815                          (f2cl-lib:fref phi-%data%
 816                                         (l kp2)
 817                                         ((1 neqn) (1 16))
 818                                         phi-%offset%)
 819                          (f2cl-lib:fref wt-%data% (l) ((1 neqn)) wt-%offset%))
 820                         2)))))
 821         (setf erkp1
 822                 (* absh
 823                    (f2cl-lib:fref gstr (kp1) ((1 13)))
 824                    (f2cl-lib:fsqrt erkp1)))
 825         (if (> k 1) (go label445))
 826         (if (>= erkp1 (* 0.5f0 erk)) (go label460))
 827         (go label450)
 828        label445
 829         (if (<= erkm1 (min erk erkp1)) (go label455))
 830         (if (or (>= erkp1 erk) (= k 12)) (go label460))
 831        label450
 832         (setf k kp1)
 833         (setf erk erkp1)
 834         (go label460)
 835        label455
 836         (setf k km1)
 837         (setf erk erkm1)
 838        label460
 839         (setf hnew (+ h h))
 840         (if phase1 (go label465))
 841         (if
 842          (>= p5eps
 843              (* erk (f2cl-lib:fref two ((f2cl-lib:int-add k 1)) ((1 13)))))
 844          (go label465))
 845         (setf hnew h)
 846         (if (>= p5eps erk) (go label465))
 847         (setf temp2
 848                 (coerce (the f2cl-lib:integer4 (f2cl-lib:int-add k 1))
 849                         'double-float))
 850         (setf r (expt (/ p5eps erk) (/ 1.0f0 temp2)))
 851         (setf hnew (* absh (max 0.5 (min 0.9 r))))
 852         (setf hnew (f2cl-lib:sign (max hnew (* fouru (abs x))) h))
 853        label465
 854         (setf h hnew)
 855         (go end_label)
 856        end_label
 857         (return
 858          (values nil
 859                  nil
 860                  nil
 861                  x
 862                  h
 863                  eps
 864                  nil
 865                  start
 866                  hold
 867                  k
 868                  kold
 869                  crash
 870                  nil
 871                  nil
 872                  nil
 873                  nil
 874                  nil
 875                  nil
 876                  ksteps
 877                  xold
 878                  ivc
 879                  nil
 880                  kgi
 881                  nil
 882                  nil
 883                  nil
 884                  nil
 885                  ifpc1
 886                  nil
 887                  nil
 888                  nil
 889                  ifpc2
 890                  ifpc3
 891                  nil
 892                  nil))))))
 893
 894 (in-package #:cl-user)
 895 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 896 (eval-when (:load-toplevel :compile-toplevel :execute)
 897   (setf (gethash 'fortran-to-lisp::stepds
 898                  fortran-to-lisp::*f2cl-function-info*)
 899           (fortran-to-lisp::make-f2cl-finfo
 900            :arg-types '(t (fortran-to-lisp::integer4) (array double-float (*))
 901                         (double-float) (double-float) (double-float)
 902                         (array double-float (*)) fortran-to-lisp::logical
 903                         (double-float) (fortran-to-lisp::integer4)
 904                         (fortran-to-lisp::integer4) fortran-to-lisp::logical
 905                         (array double-float (*)) (array double-float (*))
 906                         (array double-float (*)) (array double-float (*))
 907                         (array double-float (*)) (array double-float (*))
 908                         (fortran-to-lisp::integer4) (double-float)
 909                         (fortran-to-lisp::integer4)
 910                         (array fortran-to-lisp::integer4 (*))
 911                         (fortran-to-lisp::integer4) (array double-float (*))
 912                         (array double-float (*)) (array double-float (*))
 913                         (array double-float (*)) (fortran-to-lisp::integer4)
 914                         (array fortran-to-lisp::integer4 (*))
 915                         (array double-float (*)) (array double-float (*))
 916                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 917                         (array double-float (*))
 918                         (array fortran-to-lisp::integer4 (*)))
 919            :return-values '(nil nil nil fortran-to-lisp::x fortran-to-lisp::h
 920                             fortran-to-lisp::eps nil fortran-to-lisp::start
 921                             fortran-to-lisp::hold fortran-to-lisp::k
 922                             fortran-to-lisp::kold fortran-to-lisp::crash nil
 923                             nil nil nil nil nil fortran-to-lisp::ksteps
 924                             fortran-to-lisp::xold fortran-to-lisp::ivc nil
 925                             fortran-to-lisp::kgi nil nil nil nil
 926                             fortran-to-lisp::ifpc1 nil nil nil
 927                             fortran-to-lisp::ifpc2 fortran-to-lisp::ifpc3 nil
 928                             nil)
 929            :calls '(fortran-to-lisp::d1mach))))
 930