share/hompack/lisp/fixpds.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* ((limitd 1000))
  21   (declare (type (f2cl-lib:integer4 1000 1000) limitd) (ignorable limitd))
  22   (let ((y1sout 0.0)
  23         (xold 0.0)
  24         (sqnp1 0.0)
  25         (sout 0.0)
  26         (sb 0.0)
  27         (sa 0.0)
  28         (s99 0.0)
  29         (s 0.0)
  30         (hold 0.0)
  31         (h 0.0)
  32         (epst 0.0)
  33         (epsstp 0.0)
  34         (curtol 0.0)
  35         (cursw 0.0)
  36         (aold 0.0)
  37         (np1 0)
  38         (nfec 0)
  39         (limit 0)
  40         (lcode 0)
  41         (ksteps 0)
  42         (kold 0)
  43         (kgi 0)
  44         (k 0)
  45         (jw 0)
  46         (j 0)
  47         (ivc 0)
  48         (iter 0)
  49         (iflagc 0)
  50         (st99 nil)
  51         (crash nil)
  52         (start nil)
  53         (gi (make-array 11 :element-type 'double-float))
  54         (g (make-array 13 :element-type 'double-float))
  55         (w (make-array 12 :element-type 'double-float))
  56         (alphas (make-array 12 :element-type 'double-float))
  57         (iv (make-array 10 :element-type 'f2cl-lib:integer4)))
  58     (declare (type (double-float) y1sout xold sqnp1 sout sb sa s99 s hold h
  59                                   epst epsstp curtol cursw aold)
  60              (type (f2cl-lib:integer4) np1 nfec limit lcode ksteps kold kgi k
  61                                        jw j ivc iter iflagc)
  62              (type f2cl-lib:logical st99 crash start)
  63              (type (array double-float (11)) gi)
  64              (type (array double-float (13)) g)
  65              (type (array double-float (12)) w alphas)
  66              (type (array f2cl-lib:integer4 (10)) iv))
  67     (defun fixpds
  68            (n y iflag arctol eps trace$ a ndima nfe arclen yp ypold qr lenqr
  69             pivot pp work wt phi p par ipar)
  70       (declare (type (array f2cl-lib:integer4 (*)) ipar)
  71                (type (array double-float (*)) par)
  72                (type (array f2cl-lib:integer4 (*)) pivot)
  73                (type (double-float) arclen eps arctol)
  74                (type (array double-float (*)) p phi wt work pp qr ypold yp a y)
  75                (type (f2cl-lib:integer4) lenqr nfe ndima trace$ iflag n))
  76       (f2cl-lib:with-multi-array-data
  77           ((y double-float y-%data% y-%offset%)
  78            (a double-float a-%data% a-%offset%)
  79            (yp double-float yp-%data% yp-%offset%)
  80            (ypold double-float ypold-%data% ypold-%offset%)
  81            (qr double-float qr-%data% qr-%offset%)
  82            (pp double-float pp-%data% pp-%offset%)
  83            (work double-float work-%data% work-%offset%)
  84            (wt double-float wt-%data% wt-%offset%)
  85            (phi double-float phi-%data% phi-%offset%)
  86            (p double-float p-%data% p-%offset%)
  87            (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  88            (par double-float par-%data% par-%offset%)
  89            (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
  90         (prog ()
  91           (declare)
  92           (if (or (<= n 0) (<= eps 0.0f0)) (setf iflag 7))
  93           (if (and (>= iflag -2) (<= iflag 0)) (go label10))
  94           (if (= iflag 2) (go label35))
  95           (if (= iflag 3) (go label30))
  96           (setf iflag 7)
  97           (go end_label)
  98          label10
  99           (setf arclen (coerce 0.0f0 'double-float))
 100           (setf s (coerce 0.0f0 'double-float))
 101           (if (<= arctol 0.0f0) (setf arctol (* 0.5f0 (f2cl-lib:fsqrt eps))))
 102           (setf nfec 0)
 103           (setf iflagc iflag)
 104           (setf np1 (f2cl-lib:int-add n 1))
 105           (setf sqnp1 (f2cl-lib:fsqrt (f2cl-lib:dble np1)))
 106           (setf cursw (coerce 10.0f0 'double-float))
 107           (setf st99 f2cl-lib:%false%)
 108           (setf start f2cl-lib:%true%)
 109           (setf crash f2cl-lib:%false%)
 110           (setf hold (coerce 1.0f0 'double-float))
 111           (setf h (coerce 0.1f0 'double-float))
 112           (setf epsstp arctol)
 113           (setf ksteps 0)
 114           (setf (f2cl-lib:fref ypold-%data%
 115                                (np1)
 116                                ((1 (f2cl-lib:int-add n 1)))
 117                                ypold-%offset%)
 118                   (coerce 1.0f0 'double-float))
 119           (setf (f2cl-lib:fref yp-%data%
 120                                (np1)
 121                                ((1 (f2cl-lib:int-add n 1)))
 122                                yp-%offset%)
 123                   (coerce 1.0f0 'double-float))
 124           (setf (f2cl-lib:fref y-%data%
 125                                (np1)
 126                                ((1 (f2cl-lib:int-add n 1)))
 127                                y-%offset%)
 128                   (coerce 0.0f0 'double-float))
 129           (setf (f2cl-lib:fref work-%data%
 130                                ((f2cl-lib:int-mul 2 np1))
 131                                ((1
 132                                  (f2cl-lib:int-add
 133                                   (f2cl-lib:int-mul 6 (f2cl-lib:int-add n 1))
 134                                   lenqr)))
 135                                work-%offset%)
 136                   (coerce 0.0f0 'double-float))
 137           (setf (f2cl-lib:fref work-%data%
 138                                ((f2cl-lib:int-mul 3 np1))
 139                                ((1
 140                                  (f2cl-lib:int-add
 141                                   (f2cl-lib:int-mul 6 (f2cl-lib:int-add n 1))
 142                                   lenqr)))
 143                                work-%offset%)
 144                   (coerce 0.0f0 'double-float))
 145           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 146                         ((> j n) nil)
 147             (tagbody
 148               (setf (f2cl-lib:fref ypold-%data%
 149                                    (j)
 150                                    ((1 (f2cl-lib:int-add n 1)))
 151                                    ypold-%offset%)
 152                       (coerce 0.0f0 'double-float))
 153               (setf (f2cl-lib:fref yp-%data%
 154                                    (j)
 155                                    ((1 (f2cl-lib:int-add n 1)))
 156                                    yp-%offset%)
 157                       (coerce 0.0f0 'double-float))
 158               (setf (f2cl-lib:fref work-%data%
 159                                    ((f2cl-lib:int-add np1 j))
 160                                    ((1
 161                                      (f2cl-lib:int-add
 162                                       (f2cl-lib:int-mul 6
 163                                                         (f2cl-lib:int-add n 1))
 164                                       lenqr)))
 165                                    work-%offset%)
 166                       (coerce 0.0f0 'double-float))
 167               (setf (f2cl-lib:fref work-%data%
 168                                    ((f2cl-lib:int-add (f2cl-lib:int-mul 2 np1)
 169                                                       j))
 170                                    ((1
 171                                      (f2cl-lib:int-add
 172                                       (f2cl-lib:int-mul 6
 173                                                         (f2cl-lib:int-add n 1))
 174                                       lenqr)))
 175                                    work-%offset%)
 176                       (coerce 0.0f0 'double-float))
 177              label20))
 178           (cond
 179             ((>= iflagc (f2cl-lib:int-sub 1))
 180              (dcopy n y 1 a 1)))
 181          label30
 182           (setf limit limitd)
 183          label35
 184           (f2cl-lib:fdo (iter 1 (f2cl-lib:int-add iter 1))
 185                         ((> iter limit) nil)
 186             (tagbody
 187               (cond
 188                 ((< (f2cl-lib:fref y (np1) ((1 (f2cl-lib:int-add n 1)))) 0.0f0)
 189                  (tagbody
 190                   label40
 191                    (setf arclen (+ arclen s))
 192                    (setf iflag 5)
 193                    (go end_label))))
 194              label50
 195               (if (<= s (* 7.0f0 sqnp1)) (go label80))
 196               (setf arclen (+ arclen s))
 197               (setf s (coerce 0.0f0 'double-float))
 198              label60
 199               (setf start f2cl-lib:%true%)
 200               (setf crash f2cl-lib:%false%)
 201               (cond
 202                 ((= iflagc (f2cl-lib:int-sub 2))
 203                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 204                                ((> jw ndima) nil)
 205                    (tagbody
 206                      (setf (f2cl-lib:fref qr-%data%
 207                                           (jw)
 208                                           ((1 lenqr))
 209                                           qr-%offset%)
 210                              (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%))
 211                     label63))
 212                  (rhoa a
 213                   (f2cl-lib:fref y-%data%
 214                                  (np1)
 215                                  ((1 (f2cl-lib:int-add n 1)))
 216                                  y-%offset%)
 217                   y par ipar)
 218                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 219                                ((> jw ndima) nil)
 220                    (tagbody
 221                      (setf aold
 222                              (f2cl-lib:fref qr-%data%
 223                                             (jw)
 224                                             ((1 lenqr))
 225                                             qr-%offset%))
 226                      (cond
 227                        ((> (abs (+ (f2cl-lib:fref a (jw) ((1 n))) (- aold)))
 228                            (+ 1.0f0 (abs aold)))
 229                         (setf arclen (+ arclen s))
 230                         (setf iflag 5)
 231                         (go end_label)))
 232                     label65)))
 233                 (t
 234                  (f y yp)
 235                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 236                                ((> jw n) nil)
 237                    (tagbody
 238                      (setf aold
 239                              (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%))
 240                      (cond
 241                        ((= iflagc (f2cl-lib:int-sub 1))
 242                         (setf (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%)
 243                                 (+
 244                                  (/
 245                                   (*
 246                                    (f2cl-lib:fref y-%data%
 247                                                   (np1)
 248                                                   ((1 (f2cl-lib:int-add n 1)))
 249                                                   y-%offset%)
 250                                    (f2cl-lib:fref yp-%data%
 251                                                   (jw)
 252                                                   ((1 (f2cl-lib:int-add n 1)))
 253                                                   yp-%offset%))
 254                                   (- 1.0f0
 255                                      (f2cl-lib:fref y-%data%
 256                                                     (np1)
 257                                                     ((1
 258                                                       (f2cl-lib:int-add n 1)))
 259                                                     y-%offset%)))
 260                                  (f2cl-lib:fref y-%data%
 261                                                 (jw)
 262                                                 ((1 (f2cl-lib:int-add n 1)))
 263                                                 y-%offset%))))
 264                        (t
 265                         (setf (f2cl-lib:fref a-%data% (jw) ((1 n)) a-%offset%)
 266                                 (/
 267                                  (-
 268                                   (f2cl-lib:fref y-%data%
 269                                                  (jw)
 270                                                  ((1 (f2cl-lib:int-add n 1)))
 271                                                  y-%offset%)
 272                                   (*
 273                                    (f2cl-lib:fref y-%data%
 274                                                   (np1)
 275                                                   ((1 (f2cl-lib:int-add n 1)))
 276                                                   y-%offset%)
 277                                    (f2cl-lib:fref yp-%data%
 278                                                   (jw)
 279                                                   ((1 (f2cl-lib:int-add n 1)))
 280                                                   yp-%offset%)))
 281                                  (- 1.0f0
 282                                     (f2cl-lib:fref y-%data%
 283                                                    (np1)
 284                                                    ((1 (f2cl-lib:int-add n 1)))
 285                                                    y-%offset%))))))
 286                      (cond
 287                        ((> (abs (+ (f2cl-lib:fref a (jw) ((1 n))) (- aold)))
 288                            (+ 1.0f0 (abs aold)))
 289                         (setf arclen (+ arclen s))
 290                         (setf iflag 5)
 291                         (go end_label)))
 292                     label70))))
 293               (go label100)
 294              label80
 295               (if
 296                (or
 297                 (<=
 298                  (f2cl-lib:fref y-%data%
 299                                 (np1)
 300                                 ((1 (f2cl-lib:int-add n 1)))
 301                                 y-%offset%)
 302                  0.99f0)
 303                 st99)
 304                (go label100))
 305              label90
 306               (setf st99 f2cl-lib:%true%)
 307               (setf epsstp eps)
 308               (setf arctol eps)
 309               (go label60)
 310              label100
 311               (setf curtol (* cursw hold))
 312               (setf epst (/ eps epsstp))
 313               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 314                             ((> jw np1) nil)
 315                 (tagbody
 316                   (cond
 317                     ((<=
 318                       (abs
 319                        (+ (f2cl-lib:fref yp (jw) ((1 (f2cl-lib:int-add n 1))))
 320                           (-
 321                            (f2cl-lib:fref ypold
 322                                           (jw)
 323                                           ((1 (f2cl-lib:int-add n 1)))))))
 324                       curtol)
 325                      (setf (f2cl-lib:fref wt-%data%
 326                                           (jw)
 327                                           ((1 (f2cl-lib:int-add n 1)))
 328                                           wt-%offset%)
 329                              (+
 330                               (abs
 331                                (f2cl-lib:fref y-%data%
 332                                               (jw)
 333                                               ((1 (f2cl-lib:int-add n 1)))
 334                                               y-%offset%))
 335                               1.0f0)))
 336                     (t
 337                      (setf (f2cl-lib:fref wt-%data%
 338                                           (jw)
 339                                           ((1 (f2cl-lib:int-add n 1)))
 340                                           wt-%offset%)
 341                              (*
 342                               (+
 343                                (abs
 344                                 (f2cl-lib:fref y-%data%
 345                                                (jw)
 346                                                ((1 (f2cl-lib:int-add n 1)))
 347                                                y-%offset%))
 348                                1.0f0)
 349                               epst))))
 350                  label110))
 351               (multiple-value-bind
 352                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 353                      var-9 var-10 var-11 var-12 var-13 var-14 var-15 var-16
 354                      var-17 var-18 var-19 var-20 var-21 var-22 var-23 var-24
 355                      var-25 var-26 var-27 var-28 var-29 var-30 var-31 var-32
 356                      var-33 var-34)
 357                   (stepds #'fodeds np1 y s h epsstp wt start hold k kold crash
 358                    phi p yp alphas w g ksteps xold ivc iv kgi gi ypold a qr
 359                    lenqr pivot pp work nfec iflagc par ipar)
 360                 (declare (ignore var-0 var-1 var-2 var-6 var-12 var-13 var-14
 361                                  var-15 var-16 var-17 var-21 var-23 var-24
 362                                  var-25 var-26 var-28 var-29 var-30 var-33
 363                                  var-34))
 364                 (setf s var-3)
 365                 (setf h var-4)
 366                 (setf epsstp var-5)
 367                 (setf start var-7)
 368                 (setf hold var-8)
 369                 (setf k var-9)
 370                 (setf kold var-10)
 371                 (setf crash var-11)
 372                 (setf ksteps var-18)
 373                 (setf xold var-19)
 374                 (setf ivc var-20)
 375                 (setf kgi var-22)
 376                 (setf lenqr var-27)
 377                 (setf nfec var-31)
 378                 (setf iflagc var-32))
 379               (cond
 380                 ((> trace$ 0)
 381                  (f2cl-lib:fformat trace
 382                                    ("~%" " STEP" 1 (("~5D")) "~3@T" "NFE =" 1
 383                                     (("~5D")) "~3@T" "ARC LENGTH =" 1
 384                                     (("~9,4,0,'*,F")) "~3@T" "LAMBDA =" 1
 385                                     (("~7,4,0,'*,F")) "~5@T" "X vector:" "~%" t
 386                                     ("~1@T" 6 (("~12,4,2,0,'*,,'EE"))) "~%")
 387                                    iter
 388                                    nfec
 389                                    s
 390                                    (f2cl-lib:fref y-%data%
 391                                                   (np1)
 392                                                   ((1 (f2cl-lib:int-add n 1)))
 393                                                   y-%offset%)
 394                                    (do ((jw 1 (f2cl-lib:int-add jw 1))
 395                                         (%ret nil))
 396                                        ((> jw n) (nreverse %ret))
 397                                      (declare (type f2cl-lib:integer4 jw))
 398                                      (push
 399                                       (f2cl-lib:fref y-%data%
 400                                                      (jw)
 401                                                      ((1
 402                                                        (f2cl-lib:int-add n 1)))
 403                                                      y-%offset%)
 404                                       %ret)))))
 405               (setf nfe nfec)
 406               (cond
 407                 ((= iflagc 4)
 408                  (setf arclen (+ arclen s))
 409                  (setf iflag 4)
 410                  (go end_label)))
 411              label120
 412               (cond
 413                 (crash
 414                  (setf iflag 2)
 415                  (setf eps epsstp)
 416                  (if (< arctol eps) (setf arctol eps))
 417                  (setf limit (f2cl-lib:int-sub limit iter))
 418                  (go end_label)))
 419              label130
 420               (cond
 421                 ((>= (f2cl-lib:fref y (np1) ((1 (f2cl-lib:int-add n 1)))) 1.0f0)
 422                  (tagbody
 423                    (if st99 (go label160))
 424                    (setf s99 (- s (* 0.5f0 hold)))
 425                   label135
 426                    (sintrp s y s99 wt yp np1 kold phi ivc iv kgi gi alphas g w
 427                     xold p)
 428                    (if
 429                     (<
 430                      (f2cl-lib:fref wt-%data%
 431                                     (np1)
 432                                     ((1 (f2cl-lib:int-add n 1)))
 433                                     wt-%offset%)
 434                      1.0f0)
 435                     (go label140))
 436                    (setf s99 (* 0.5f0 (+ (- s hold) s99)))
 437                    (go label135)
 438                   label140
 439                    (dcopy np1 wt 1 y 1)
 440                    (dcopy np1 yp 1 ypold 1)
 441                    (setf s s99)
 442                    (go label90))))
 443              label150))
 444           (setf iflag 3)
 445           (go end_label)
 446          label160
 447           (setf sa (- s hold))
 448           (setf sb s)
 449           (setf lcode 1)
 450          label170
 451           (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 452               (root sout y1sout sa sb eps eps lcode)
 453             (declare (ignore var-1 var-4 var-5))
 454             (setf sout var-0)
 455             (setf sa var-2)
 456             (setf sb var-3)
 457             (setf lcode var-6))
 458           (if (> lcode 0) (go label190))
 459           (sintrp s y sout wt yp np1 kold phi ivc iv kgi gi alphas g w xold p)
 460           (setf y1sout
 461                   (-
 462                    (f2cl-lib:fref wt-%data%
 463                                   (np1)
 464                                   ((1 (f2cl-lib:int-add n 1)))
 465                                   wt-%offset%)
 466                    1.0f0))
 467           (go label170)
 468          label190
 469           (setf iflag 1)
 470           (if (> lcode 2) (setf iflag 6))
 471           (setf arclen (+ arclen sa))
 472           (sintrp s y sa wt yp np1 kold phi ivc iv kgi gi alphas g w xold p)
 473           (dcopy np1 wt 1 y 1)
 474           (go end_label)
 475          end_label
 476           (return
 477            (values nil
 478                    nil
 479                    iflag
 480                    arctol
 481                    eps
 482                    nil
 483                    nil
 484                    nil
 485                    nfe
 486                    arclen
 487                    nil
 488                    nil
 489                    nil
 490                    lenqr
 491                    nil
 492                    nil
 493                    nil
 494                    nil
 495                    nil
 496                    nil
 497                    nil
 498                    nil)))))))
 499
 500 (in-package #:cl-user)
 501 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 502 (eval-when (:load-toplevel :compile-toplevel :execute)
 503   (setf (gethash 'fortran-to-lisp::fixpds
 504                  fortran-to-lisp::*f2cl-function-info*)
 505           (fortran-to-lisp::make-f2cl-finfo
 506            :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
 507                         (fortran-to-lisp::integer4) (double-float)
 508                         (double-float) (fortran-to-lisp::integer4)
 509                         (array double-float (*)) (fortran-to-lisp::integer4)
 510                         (fortran-to-lisp::integer4) (double-float)
 511                         (array double-float (*)) (array double-float (*))
 512                         (array double-float (*)) (fortran-to-lisp::integer4)
 513                         (array fortran-to-lisp::integer4 (*))
 514                         (array double-float (*)) (array double-float (*))
 515                         (array double-float (*)) (array double-float (*))
 516                         (array double-float (*)) (array double-float (*))
 517                         (array fortran-to-lisp::integer4 (*)))
 518            :return-values '(nil nil fortran-to-lisp::iflag
 519                             fortran-to-lisp::arctol fortran-to-lisp::eps nil
 520                             nil nil fortran-to-lisp::nfe
 521                             fortran-to-lisp::arclen nil nil nil
 522                             fortran-to-lisp::lenqr nil nil nil nil nil nil nil
 523                             nil)
 524            :calls '(fortran-to-lisp::f fortran-to-lisp::rhoa
 525                     fortran-to-lisp::dcopy fortran-to-lisp::root
 526                     fortran-to-lisp::sintrp fortran-to-lisp::stepds))))
 527