share/hompack/lisp/stepnf.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* ((litfh 4))
  21   (declare (type (f2cl-lib:integer4 4 4) litfh) (ignorable litfh))
  22   (defun stepnf
  23          (n nfe iflag start crash hold h relerr abserr s y yp yold ypold a qr
  24           alpha tz pivot w wp z0 z1 sspar par ipar)
  25     (declare (type (array f2cl-lib:integer4 (*)) ipar)
  26              (type (array double-float (*)) par)
  27              (type (array double-float (*)) sspar)
  28              (type (array f2cl-lib:integer4 (*)) pivot)
  29              (type (array double-float (*)) z1 z0 wp w tz alpha qr a ypold yold
  30                                             yp y)
  31              (type (double-float) s abserr relerr h hold)
  32              (type f2cl-lib:logical crash start)
  33              (type (f2cl-lib:integer4) iflag nfe n))
  34     (f2cl-lib:with-multi-array-data
  35         ((y double-float y-%data% y-%offset%)
  36          (yp double-float yp-%data% yp-%offset%)
  37          (yold double-float yold-%data% yold-%offset%)
  38          (ypold double-float ypold-%data% ypold-%offset%)
  39          (a double-float a-%data% a-%offset%)
  40          (qr double-float qr-%data% qr-%offset%)
  41          (alpha double-float alpha-%data% alpha-%offset%)
  42          (tz double-float tz-%data% tz-%offset%)
  43          (w double-float w-%data% w-%offset%)
  44          (wp double-float wp-%data% wp-%offset%)
  45          (z0 double-float z0-%data% z0-%offset%)
  46          (z1 double-float z1-%data% z1-%offset%)
  47          (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  48          (sspar double-float sspar-%data% sspar-%offset%)
  49          (par double-float par-%data% par-%offset%)
  50          (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
  51       (labels ((dd01 (f0 f1 dels)
  52                  (f2cl-lib:f2cl/ (+ f1 (- f0)) dels))
  53                (dd001 (f0 fp0 f1 dels)
  54                  (/ (- (dd01 f0 f1 dels) fp0) dels))
  55                (dd011 (f0 f1 fp1 dels)
  56                  (/ (- fp1 (dd01 f0 f1 dels)) dels))
  57                (dd0011 (f0 fp0 f1 fp1 dels)
  58                  (/ (- (dd011 f0 f1 fp1 dels) (dd001 f0 fp0 f1 dels)) dels))
  59                (qofs (f0 fp0 f1 fp1 dels s)
  60                  (+
  61                   (*
  62                    (+
  63                     (*
  64                      (+ (* (dd0011 f0 fp0 f1 fp1 dels) (- s dels))
  65                         (dd001 f0 fp0 f1 dels))
  66                      s)
  67                     fp0)
  68                    s)
  69                   f0)))
  70         (declare (ftype (function (double-float double-float double-float)
  71                          (values double-float &rest t))
  72                         dd01))
  73         (declare (ftype (function
  74                          (double-float double-float double-float double-float)
  75                          (values double-float &rest t))
  76                         dd001))
  77         (declare (ftype (function
  78                          (double-float double-float double-float double-float)
  79                          (values double-float &rest t))
  80                         dd011))
  81         (declare (ftype (function
  82                          (double-float double-float double-float double-float
  83                           double-float)
  84                          (values double-float &rest t))
  85                         dd0011))
  86         (declare (ftype (function
  87                          (double-float double-float double-float double-float
  88                           double-float double-float)
  89                          (values double-float &rest t))
  90                         qofs))
  91         (prog ((fail nil) (itnum 0) (j 0) (judy 0) (np1 0) (dcalc 0.0)
  92                (dels 0.0) (f0 0.0) (f1 0.0) (fouru 0.0) (fp0 0.0) (fp1 0.0)
  93                (hfail 0.0) (ht 0.0) (lcalc 0.0) (rcalc 0.0) (rholen 0.0)
  94                (temp 0.0) (twou 0.0))
  95           (declare (type f2cl-lib:logical fail)
  96                    (type (f2cl-lib:integer4) itnum j judy np1)
  97                    (type (double-float) dcalc dels f0 f1 fouru fp0 fp1 hfail ht
  98                                         lcalc rcalc rholen temp twou))
  99           (setf twou (* 2.0f0 (f2cl-lib:d1mach 4)))
 100           (setf fouru (+ twou twou))
 101           (setf np1 (f2cl-lib:int-add n 1))
 102           (setf crash f2cl-lib:%true%)
 103           (if (< s 0.0f0) (go end_label))
 104           (cond
 105             ((< h (* fouru (+ 1.0f0 s)))
 106              (setf h (* fouru (+ 1.0f0 s)))
 107              (go end_label)))
 108           (setf temp (dnrm2 np1 y 1))
 109           (if (>= (* 0.5f0 (+ (* relerr temp) abserr)) (* twou temp))
 110               (go label40))
 111           (cond
 112             ((/= relerr 0.0f0)
 113              (setf relerr (* fouru (+ 1.0f0 fouru)))
 114              (setf abserr (max abserr 0.0)))
 115             (t
 116              (setf abserr (* fouru temp))))
 117           (go end_label)
 118          label40
 119           (setf crash f2cl-lib:%false%)
 120           (if (not start) (go label300))
 121           (setf fail f2cl-lib:%false%)
 122           (setf start f2cl-lib:%false%)
 123           (setf h
 124                   (min h
 125                        0.1
 126                        (f2cl-lib:fsqrt
 127                         (f2cl-lib:fsqrt (+ (* relerr temp) abserr)))))
 128           (setf (f2cl-lib:fref ypold-%data%
 129                                (1)
 130                                ((1 (f2cl-lib:int-add n 1)))
 131                                ypold-%offset%)
 132                   (coerce 1.0f0 'double-float))
 133           (f2cl-lib:fdo (j 2 (f2cl-lib:int-add j 1))
 134                         ((> j np1) nil)
 135             (tagbody
 136               (setf (f2cl-lib:fref ypold-%data%
 137                                    (j)
 138                                    ((1 (f2cl-lib:int-add n 1)))
 139                                    ypold-%offset%)
 140                       (coerce 0.0f0 'double-float))
 141              label50))
 142           (multiple-value-bind
 143                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 144                  var-10 var-11 var-12 var-13)
 145               (tangnf s y yp ypold a qr alpha tz pivot nfe n iflag par ipar)
 146             (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 147                              var-10 var-12 var-13))
 148             (setf s var-0)
 149             (setf nfe var-9)
 150             (setf iflag var-11))
 151           (if (> iflag 0) (go end_label))
 152          label70
 153           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 154                         ((> j np1) nil)
 155             (tagbody
 156               (setf temp
 157                       (+
 158                        (f2cl-lib:fref y-%data%
 159                                       (j)
 160                                       ((1 (f2cl-lib:int-add n 1)))
 161                                       y-%offset%)
 162                        (* h
 163                           (f2cl-lib:fref yp-%data%
 164                                          (j)
 165                                          ((1 (f2cl-lib:int-add n 1)))
 166                                          yp-%offset%))))
 167               (setf (f2cl-lib:fref w-%data%
 168                                    (j)
 169                                    ((1 (f2cl-lib:int-add n 1)))
 170                                    w-%offset%)
 171                       temp)
 172               (setf (f2cl-lib:fref z0-%data%
 173                                    (j)
 174                                    ((1 (f2cl-lib:int-add n 1)))
 175                                    z0-%offset%)
 176                       temp)
 177              label80))
 178           (f2cl-lib:fdo (judy 1 (f2cl-lib:int-add judy 1))
 179                         ((> judy litfh) nil)
 180             (tagbody
 181               (setf rholen (coerce -1.0f0 'double-float))
 182               (multiple-value-bind
 183                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 184                      var-9 var-10 var-11 var-12 var-13)
 185                   (tangnf rholen w wp ypold a qr alpha tz pivot nfe n iflag par
 186                    ipar)
 187                 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7
 188                                  var-8 var-10 var-12 var-13))
 189                 (setf rholen var-0)
 190                 (setf nfe var-9)
 191                 (setf iflag var-11))
 192               (if (> iflag 0) (go end_label))
 193               (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 194                             ((> j np1) nil)
 195                 (tagbody
 196                   (setf (f2cl-lib:fref w-%data%
 197                                        (j)
 198                                        ((1 (f2cl-lib:int-add n 1)))
 199                                        w-%offset%)
 200                           (+
 201                            (f2cl-lib:fref w-%data%
 202                                           (j)
 203                                           ((1 (f2cl-lib:int-add n 1)))
 204                                           w-%offset%)
 205                            (f2cl-lib:fref tz-%data%
 206                                           (j)
 207                                           ((1 (f2cl-lib:int-add n 1)))
 208                                           tz-%offset%)))
 209                  label90))
 210               (setf itnum judy)
 211               (cond
 212                 ((= judy 1)
 213                  (setf lcalc (dnrm2 np1 tz 1))
 214                  (setf rcalc rholen)
 215                  (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 216                                ((> j np1) nil)
 217                    (tagbody
 218                      (setf (f2cl-lib:fref z1-%data%
 219                                           (j)
 220                                           ((1 (f2cl-lib:int-add n 1)))
 221                                           z1-%offset%)
 222                              (f2cl-lib:fref w-%data%
 223                                             (j)
 224                                             ((1 (f2cl-lib:int-add n 1)))
 225                                             w-%offset%))
 226                     label110)))
 227                 ((= judy 2)
 228                  (setf lcalc (/ (dnrm2 np1 tz 1) lcalc))
 229                  (setf rcalc (/ rholen rcalc))))
 230               (if (<= (dnrm2 np1 tz 1) (+ (* relerr (dnrm2 np1 w 1)) abserr))
 231                   (go label600))
 232              label200))
 233           (cond
 234             ((<= h (* fouru (+ 1.0f0 s)))
 235              (setf iflag 6)
 236              (go end_label)))
 237           (setf h (* 0.5f0 h))
 238           (go label70)
 239          label300
 240           (setf fail f2cl-lib:%false%)
 241          label320
 242           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 243                         ((> j np1) nil)
 244             (tagbody
 245               (setf temp
 246                       (qofs
 247                        (f2cl-lib:fref yold-%data%
 248                                       (j)
 249                                       ((1 (f2cl-lib:int-add n 1)))
 250                                       yold-%offset%)
 251                        (f2cl-lib:fref ypold-%data%
 252                                       (j)
 253                                       ((1 (f2cl-lib:int-add n 1)))
 254                                       ypold-%offset%)
 255                        (f2cl-lib:fref y-%data%
 256                                       (j)
 257                                       ((1 (f2cl-lib:int-add n 1)))
 258                                       y-%offset%)
 259                        (f2cl-lib:fref yp-%data%
 260                                       (j)
 261                                       ((1 (f2cl-lib:int-add n 1)))
 262                                       yp-%offset%)
 263                        hold (+ hold h)))
 264               (setf (f2cl-lib:fref w-%data%
 265                                    (j)
 266                                    ((1 (f2cl-lib:int-add n 1)))
 267                                    w-%offset%)
 268                       temp)
 269               (setf (f2cl-lib:fref z0-%data%
 270                                    (j)
 271                                    ((1 (f2cl-lib:int-add n 1)))
 272                                    z0-%offset%)
 273                       temp)
 274              label330))
 275           (f2cl-lib:fdo (judy 1 (f2cl-lib:int-add judy 1))
 276                         ((> judy litfh) nil)
 277             (tagbody
 278               (setf rholen (coerce -1.0f0 'double-float))
 279               (multiple-value-bind
 280                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 281                      var-9 var-10 var-11 var-12 var-13)
 282                   (tangnf rholen w wp yp a qr alpha tz pivot nfe n iflag par
 283                    ipar)
 284                 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7
 285                                  var-8 var-10 var-12 var-13))
 286                 (setf rholen var-0)
 287                 (setf nfe var-9)
 288                 (setf iflag var-11))
 289               (if (> iflag 0) (go end_label))
 290               (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 291                             ((> j np1) nil)
 292                 (tagbody
 293                   (setf (f2cl-lib:fref w-%data%
 294                                        (j)
 295                                        ((1 (f2cl-lib:int-add n 1)))
 296                                        w-%offset%)
 297                           (+
 298                            (f2cl-lib:fref w-%data%
 299                                           (j)
 300                                           ((1 (f2cl-lib:int-add n 1)))
 301                                           w-%offset%)
 302                            (f2cl-lib:fref tz-%data%
 303                                           (j)
 304                                           ((1 (f2cl-lib:int-add n 1)))
 305                                           tz-%offset%)))
 306                  label420))
 307               (setf itnum judy)
 308               (cond
 309                 ((= judy 1)
 310                  (setf lcalc (dnrm2 np1 tz 1))
 311                  (setf rcalc rholen)
 312                  (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 313                                ((> j np1) nil)
 314                    (tagbody
 315                      (setf (f2cl-lib:fref z1-%data%
 316                                           (j)
 317                                           ((1 (f2cl-lib:int-add n 1)))
 318                                           z1-%offset%)
 319                              (f2cl-lib:fref w-%data%
 320                                             (j)
 321                                             ((1 (f2cl-lib:int-add n 1)))
 322                                             w-%offset%))
 323                     label440)))
 324                 ((= judy 2)
 325                  (setf lcalc (/ (dnrm2 np1 tz 1) lcalc))
 326                  (setf rcalc (/ rholen rcalc))))
 327               (if (<= (dnrm2 np1 tz 1) (+ (* relerr (dnrm2 np1 w 1)) abserr))
 328                   (go label600))
 329              label500))
 330           (setf fail f2cl-lib:%true%)
 331           (setf hfail h)
 332           (cond
 333             ((<= h (* fouru (+ 1.0f0 s)))
 334              (setf iflag 6)
 335              (go end_label)))
 336           (setf h (* 0.5f0 h))
 337           (go label320)
 338          label600
 339           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 340                         ((> j np1) nil)
 341             (tagbody
 342               (setf (f2cl-lib:fref yold-%data%
 343                                    (j)
 344                                    ((1 (f2cl-lib:int-add n 1)))
 345                                    yold-%offset%)
 346                       (f2cl-lib:fref y-%data%
 347                                      (j)
 348                                      ((1 (f2cl-lib:int-add n 1)))
 349                                      y-%offset%))
 350               (setf (f2cl-lib:fref ypold-%data%
 351                                    (j)
 352                                    ((1 (f2cl-lib:int-add n 1)))
 353                                    ypold-%offset%)
 354                       (f2cl-lib:fref yp-%data%
 355                                      (j)
 356                                      ((1 (f2cl-lib:int-add n 1)))
 357                                      yp-%offset%))
 358               (setf (f2cl-lib:fref y-%data%
 359                                    (j)
 360                                    ((1 (f2cl-lib:int-add n 1)))
 361                                    y-%offset%)
 362                       (f2cl-lib:fref w-%data%
 363                                      (j)
 364                                      ((1 (f2cl-lib:int-add n 1)))
 365                                      w-%offset%))
 366               (setf (f2cl-lib:fref yp-%data%
 367                                    (j)
 368                                    ((1 (f2cl-lib:int-add n 1)))
 369                                    yp-%offset%)
 370                       (f2cl-lib:fref wp-%data%
 371                                      (j)
 372                                      ((1 (f2cl-lib:int-add n 1)))
 373                                      wp-%offset%))
 374               (setf (f2cl-lib:fref w-%data%
 375                                    (j)
 376                                    ((1 (f2cl-lib:int-add n 1)))
 377                                    w-%offset%)
 378                       (-
 379                        (f2cl-lib:fref y-%data%
 380                                       (j)
 381                                       ((1 (f2cl-lib:int-add n 1)))
 382                                       y-%offset%)
 383                        (f2cl-lib:fref yold-%data%
 384                                       (j)
 385                                       ((1 (f2cl-lib:int-add n 1)))
 386                                       yold-%offset%)))
 387              label620))
 388           (setf hold (dnrm2 np1 w 1))
 389           (setf s (+ s hold))
 390          label700
 391           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 392                         ((> j np1) nil)
 393             (tagbody
 394               (setf (f2cl-lib:fref tz-%data%
 395                                    (j)
 396                                    ((1 (f2cl-lib:int-add n 1)))
 397                                    tz-%offset%)
 398                       (-
 399                        (f2cl-lib:fref z0-%data%
 400                                       (j)
 401                                       ((1 (f2cl-lib:int-add n 1)))
 402                                       z0-%offset%)
 403                        (f2cl-lib:fref y-%data%
 404                                       (j)
 405                                       ((1 (f2cl-lib:int-add n 1)))
 406                                       y-%offset%)))
 407               (setf (f2cl-lib:fref w-%data%
 408                                    (j)
 409                                    ((1 (f2cl-lib:int-add n 1)))
 410                                    w-%offset%)
 411                       (-
 412                        (f2cl-lib:fref z1-%data%
 413                                       (j)
 414                                       ((1 (f2cl-lib:int-add n 1)))
 415                                       z1-%offset%)
 416                        (f2cl-lib:fref y-%data%
 417                                       (j)
 418                                       ((1 (f2cl-lib:int-add n 1)))
 419                                       y-%offset%)))
 420              label710))
 421           (setf dcalc (dnrm2 np1 tz 1))
 422           (if (/= dcalc 0.0f0) (setf dcalc (/ (dnrm2 np1 w 1) dcalc)))
 423           (if (= itnum 1) (setf lcalc (coerce 0.0f0 'double-float)))
 424           (cond
 425             ((= (+ lcalc rcalc dcalc) 0.0f0)
 426              (setf ht
 427                      (* (f2cl-lib:fref sspar-%data% (7) ((1 8)) sspar-%offset%)
 428                         hold)))
 429             (t
 430              (setf ht
 431                      (*
 432                       (expt
 433                        (/ 1.0f0
 434                           (max
 435                            (/ lcalc
 436                               (f2cl-lib:fref sspar-%data%
 437                                              (1)
 438                                              ((1 8))
 439                                              sspar-%offset%))
 440                            (/ rcalc
 441                               (f2cl-lib:fref sspar-%data%
 442                                              (2)
 443                                              ((1 8))
 444                                              sspar-%offset%))
 445                            (/ dcalc
 446                               (f2cl-lib:fref sspar-%data%
 447                                              (3)
 448                                              ((1 8))
 449                                              sspar-%offset%))))
 450                        (/ 1.0f0
 451                           (f2cl-lib:fref sspar-%data%
 452                                          (8)
 453                                          ((1 8))
 454                                          sspar-%offset%)))
 455                       hold))))
 456           (setf h
 457                   (min
 458                    (max ht
 459                         (*
 460                          (f2cl-lib:fref sspar-%data%
 461                                         (6)
 462                                         ((1 8))
 463                                         sspar-%offset%)
 464                          hold)
 465                         (f2cl-lib:fref sspar-%data%
 466                                        (4)
 467                                        ((1 8))
 468                                        sspar-%offset%))
 469                    (* (f2cl-lib:fref sspar-%data% (7) ((1 8)) sspar-%offset%)
 470                       hold)
 471                    (f2cl-lib:fref sspar-%data% (5) ((1 8)) sspar-%offset%)))
 472           (cond
 473             ((= itnum 1)
 474              (setf h (max h hold)))
 475             ((= itnum litfh)
 476              (setf h (min h hold))))
 477           (if fail (setf h (min h hfail)))
 478           (go end_label)
 479          end_label
 480           (return
 481            (values nil
 482                    nfe
 483                    iflag
 484                    start
 485                    crash
 486                    hold
 487                    h
 488                    relerr
 489                    abserr
 490                    s
 491                    nil
 492                    nil
 493                    nil
 494                    nil
 495                    nil
 496                    nil
 497                    nil
 498                    nil
 499                    nil
 500                    nil
 501                    nil
 502                    nil
 503                    nil
 504                    nil
 505                    nil
 506                    nil)))))))
 507
 508 (in-package #:cl-user)
 509 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 510 (eval-when (:load-toplevel :compile-toplevel :execute)
 511   (setf (gethash 'fortran-to-lisp::stepnf
 512                  fortran-to-lisp::*f2cl-function-info*)
 513           (fortran-to-lisp::make-f2cl-finfo
 514            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 515                         (fortran-to-lisp::integer4) fortran-to-lisp::logical
 516                         fortran-to-lisp::logical (double-float) (double-float)
 517                         (double-float) (double-float) (double-float)
 518                         (array double-float (*)) (array double-float (*))
 519                         (array double-float (*)) (array double-float (*))
 520                         (array double-float (*)) (array double-float (*))
 521                         (array double-float (*)) (array double-float (*))
 522                         (array fortran-to-lisp::integer4 (*))
 523                         (array double-float (*)) (array double-float (*))
 524                         (array double-float (*)) (array double-float (*))
 525                         (array double-float (*)) (array double-float (*))
 526                         (array fortran-to-lisp::integer4 (*)))
 527            :return-values '(nil fortran-to-lisp::nfe fortran-to-lisp::iflag
 528                             fortran-to-lisp::start fortran-to-lisp::crash
 529                             fortran-to-lisp::hold fortran-to-lisp::h
 530                             fortran-to-lisp::relerr fortran-to-lisp::abserr
 531                             fortran-to-lisp::s nil nil nil nil nil nil nil nil
 532                             nil nil nil nil nil nil nil nil)
 533            :calls '(fortran-to-lisp::dnrm2 fortran-to-lisp::tangnf
 534                     fortran-to-lisp::d1mach))))
 535