share/hompack/lisp/fixpns.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) (cursw (coerce 10.0f0 'double-float)))
  21   (declare (type (f2cl-lib:integer4 1000 1000) limitd)
  22            (type (double-float) cursw)
  23            (ignorable limitd cursw))
  24   (let ((s 0.0)
  25         (relerr 0.0)
  26         (hold 0.0)
  27         (h 0.0)
  28         (curtol 0.0)
  29         (abserr 0.0)
  30         (np1 0)
  31         (nfec 0)
  32         (nc 0)
  33         (limit 0)
  34         (jw 0)
  35         (iz1 0)
  36         (iz0 0)
  37         (iwp 0)
  38         (iw 0)
  39         (itz 0)
  40         (iter 0)
  41         (itangw 0)
  42         (irho 0)
  43         (ipp 0)
  44         (iflagc 0)
  45         (crash nil)
  46         (start nil))
  47     (declare (type (double-float) s relerr hold h curtol abserr)
  48              (type (f2cl-lib:integer4) np1 nfec nc limit jw iz1 iz0 iwp iw itz
  49                                        iter itangw irho ipp iflagc)
  50              (type f2cl-lib:logical crash start))
  51     (defun fixpns
  52            (n y iflag arcre arcae ansre ansae trace$ a nfe arclen yp yold ypold
  53             qr lenqr pivot work sspar par ipar)
  54       (declare (type (array f2cl-lib:integer4 (*)) ipar)
  55                (type (array double-float (*)) par)
  56                (type (array double-float (*)) sspar)
  57                (type (array f2cl-lib:integer4 (*)) pivot)
  58                (type (double-float) arclen ansae ansre arcae arcre)
  59                (type (array double-float (*)) work qr ypold yold yp a y)
  60                (type (f2cl-lib:integer4) lenqr nfe trace$ iflag n))
  61       (f2cl-lib:with-multi-array-data
  62           ((y double-float y-%data% y-%offset%)
  63            (a double-float a-%data% a-%offset%)
  64            (yp double-float yp-%data% yp-%offset%)
  65            (yold double-float yold-%data% yold-%offset%)
  66            (ypold double-float ypold-%data% ypold-%offset%)
  67            (qr double-float qr-%data% qr-%offset%)
  68            (work double-float work-%data% work-%offset%)
  69            (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  70            (sspar double-float sspar-%data% sspar-%offset%)
  71            (par double-float par-%data% par-%offset%)
  72            (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
  73         (prog ()
  74           (declare)
  75           (if (or (<= n 0) (<= ansre 0.0f0) (< ansae 0.0f0)) (setf iflag 7))
  76           (if (and (>= iflag -2) (<= iflag 0)) (go label20))
  77           (if (= iflag 2) (go label120))
  78           (if (= iflag 3) (go label90))
  79           (setf iflag 7)
  80           (go end_label)
  81          label20
  82           (setf arclen (coerce 0.0f0 'double-float))
  83           (if (<= arcre 0.0f0) (setf arcre (* 0.5f0 (f2cl-lib:fsqrt ansre))))
  84           (if (<= arcae 0.0f0) (setf arcae (* 0.5f0 (f2cl-lib:fsqrt ansae))))
  85           (setf nc n)
  86           (setf nfec 0)
  87           (setf iflagc iflag)
  88           (setf np1 (f2cl-lib:int-add n 1))
  89           (setf ipp 1)
  90           (setf irho (f2cl-lib:int-add n 1))
  91           (setf iw (f2cl-lib:int-add irho n))
  92           (setf iwp (f2cl-lib:int-add iw np1))
  93           (setf itz (f2cl-lib:int-add iwp np1))
  94           (setf iz0 (f2cl-lib:int-add itz np1))
  95           (setf iz1 (f2cl-lib:int-add iz0 np1))
  96           (setf itangw (f2cl-lib:int-add iz1 np1))
  97           (setf start f2cl-lib:%true%)
  98           (setf crash f2cl-lib:%false%)
  99           (setf hold (coerce 1.0f0 'double-float))
 100           (setf h (coerce 0.1f0 'double-float))
 101           (setf s (coerce 0.0f0 'double-float))
 102           (setf (f2cl-lib:fref ypold-%data%
 103                                (np1)
 104                                ((1 (f2cl-lib:int-add n 1)))
 105                                ypold-%offset%)
 106                   (coerce 1.0f0 'double-float))
 107           (setf (f2cl-lib:fref yp-%data%
 108                                (np1)
 109                                ((1 (f2cl-lib:int-add n 1)))
 110                                yp-%offset%)
 111                   (coerce 1.0f0 'double-float))
 112           (setf (f2cl-lib:fref y-%data%
 113                                (np1)
 114                                ((1 (f2cl-lib:int-add n 1)))
 115                                y-%offset%)
 116                   (coerce 0.0f0 'double-float))
 117           (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 118                         ((> jw n) nil)
 119             (tagbody
 120               (setf (f2cl-lib:fref ypold-%data%
 121                                    (jw)
 122                                    ((1 (f2cl-lib:int-add n 1)))
 123                                    ypold-%offset%)
 124                       (coerce 0.0f0 'double-float))
 125               (setf (f2cl-lib:fref yp-%data%
 126                                    (jw)
 127                                    ((1 (f2cl-lib:int-add n 1)))
 128                                    yp-%offset%)
 129                       (coerce 0.0f0 'double-float))
 130              label40))
 131           (f2cl-lib:fdo (jw itangw (f2cl-lib:int-add jw 1))
 132                         ((> jw (f2cl-lib:int-add itangw np1 n)) nil)
 133             (tagbody
 134               (setf (f2cl-lib:fref work-%data%
 135                                    (jw)
 136                                    ((1
 137                                      (f2cl-lib:int-add
 138                                       (f2cl-lib:int-mul 13
 139                                                         (f2cl-lib:int-add n 1))
 140                                       (f2cl-lib:int-mul 2 n)
 141                                       lenqr)))
 142                                    work-%offset%)
 143                       (coerce 0.0f0 'double-float))
 144              label50))
 145           (if
 146            (<= (f2cl-lib:fref sspar-%data% (1) ((1 8)) sspar-%offset%) 0.0f0)
 147            (setf (f2cl-lib:fref sspar-%data% (1) ((1 8)) sspar-%offset%)
 148                    (coerce 0.5f0 'double-float)))
 149           (if
 150            (<= (f2cl-lib:fref sspar-%data% (2) ((1 8)) sspar-%offset%) 0.0f0)
 151            (setf (f2cl-lib:fref sspar-%data% (2) ((1 8)) sspar-%offset%)
 152                    (coerce 0.01f0 'double-float)))
 153           (if
 154            (<= (f2cl-lib:fref sspar-%data% (3) ((1 8)) sspar-%offset%) 0.0f0)
 155            (setf (f2cl-lib:fref sspar-%data% (3) ((1 8)) sspar-%offset%)
 156                    (coerce 0.5f0 'double-float)))
 157           (if
 158            (<= (f2cl-lib:fref sspar-%data% (4) ((1 8)) sspar-%offset%) 0.0f0)
 159            (setf (f2cl-lib:fref sspar-%data% (4) ((1 8)) sspar-%offset%)
 160                    (* (+ (f2cl-lib:fsqrt (+ n 1.0f0)) 4.0f0)
 161                       (f2cl-lib:d1mach 4))))
 162           (if
 163            (<= (f2cl-lib:fref sspar-%data% (5) ((1 8)) sspar-%offset%) 0.0f0)
 164            (setf (f2cl-lib:fref sspar-%data% (5) ((1 8)) sspar-%offset%)
 165                    (coerce 1.0f0 'double-float)))
 166           (if
 167            (<= (f2cl-lib:fref sspar-%data% (6) ((1 8)) sspar-%offset%) 0.0f0)
 168            (setf (f2cl-lib:fref sspar-%data% (6) ((1 8)) sspar-%offset%)
 169                    (coerce 0.1f0 'double-float)))
 170           (if
 171            (<= (f2cl-lib:fref sspar-%data% (7) ((1 8)) sspar-%offset%) 0.0f0)
 172            (setf (f2cl-lib:fref sspar-%data% (7) ((1 8)) sspar-%offset%)
 173                    (coerce 3.0f0 'double-float)))
 174           (if
 175            (<= (f2cl-lib:fref sspar-%data% (8) ((1 8)) sspar-%offset%) 0.0f0)
 176            (setf (f2cl-lib:fref sspar-%data% (8) ((1 8)) sspar-%offset%)
 177                    (coerce 2.0f0 'double-float)))
 178           (cond
 179             ((>= iflagc (f2cl-lib:int-sub 1))
 180              (dcopy n y 1 a 1)))
 181          label90
 182           (setf limit limitd)
 183          label120
 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                  (setf arclen s)
 190                  (setf iflag 5)
 191                  (go end_label)))
 192              label140
 193               (setf curtol (* cursw hold))
 194               (setf relerr arcre)
 195               (setf abserr arcae)
 196               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 197                             ((> jw np1) nil)
 198                 (tagbody
 199                   (cond
 200                     ((>
 201                       (abs
 202                        (+ (f2cl-lib:fref yp (jw) ((1 (f2cl-lib:int-add n 1))))
 203                           (-
 204                            (f2cl-lib:fref ypold
 205                                           (jw)
 206                                           ((1 (f2cl-lib:int-add n 1)))))))
 207                       curtol)
 208                      (setf relerr ansre)
 209                      (setf abserr ansae)
 210                      (go label200)))
 211                  label160))
 212              label200
 213               (multiple-value-bind
 214                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 215                      var-9 var-10 var-11 var-12 var-13 var-14 var-15 var-16
 216                      var-17 var-18 var-19 var-20 var-21)
 217                   (stepns nc nfec iflagc start crash hold h relerr abserr s y
 218                    yp yold ypold a qr lenqr pivot work sspar par ipar)
 219                 (declare (ignore var-0 var-10 var-11 var-12 var-13 var-14
 220                                  var-15 var-16 var-17 var-18 var-19 var-20
 221                                  var-21))
 222                 (setf nfec var-1)
 223                 (setf iflagc var-2)
 224                 (setf start var-3)
 225                 (setf crash var-4)
 226                 (setf hold var-5)
 227                 (setf h var-6)
 228                 (setf relerr var-7)
 229                 (setf abserr var-8)
 230                 (setf s var-9))
 231               (cond
 232                 ((> trace$ 0)
 233                  (f2cl-lib:fformat trace
 234                                    ("~%" " STEP" 1 (("~5D")) "~3@T" "NFE =" 1
 235                                     (("~5D")) "~3@T" "ARC LENGTH =" 1
 236                                     (("~9,4,0,'*,F")) "~3@T" "LAMBDA =" 1
 237                                     (("~7,4,0,'*,F")) "~5@T" "X vector:" "~%" t
 238                                     ("~1@T" 6 (("~12,4,2,0,'*,,'EE"))) "~%")
 239                                    iter
 240                                    nfec
 241                                    s
 242                                    (f2cl-lib:fref y-%data%
 243                                                   (np1)
 244                                                   ((1 (f2cl-lib:int-add n 1)))
 245                                                   y-%offset%)
 246                                    (do ((jw 1 (f2cl-lib:int-add jw 1))
 247                                         (%ret nil))
 248                                        ((> jw nc) (nreverse %ret))
 249                                      (declare (type f2cl-lib:integer4 jw))
 250                                      (push
 251                                       (f2cl-lib:fref y-%data%
 252                                                      (jw)
 253                                                      ((1
 254                                                        (f2cl-lib:int-add n 1)))
 255                                                      y-%offset%)
 256                                       %ret)))))
 257               (setf nfe nfec)
 258               (cond
 259                 ((> iflagc 0)
 260                  (setf arclen s)
 261                  (setf iflag iflagc)
 262                  (go end_label)))
 263               (cond
 264                 (crash
 265                  (setf iflag 2)
 266                  (if (< arcre relerr) (setf arcre relerr))
 267                  (if (< ansre relerr) (setf ansre relerr))
 268                  (if (< arcae abserr) (setf arcae abserr))
 269                  (if (< ansae abserr) (setf ansae abserr))
 270                  (setf limit (f2cl-lib:int-sub limit iter))
 271                  (go end_label)))
 272               (cond
 273                 ((>= (f2cl-lib:fref y (np1) ((1 (f2cl-lib:int-add n 1)))) 1.0f0)
 274                  (dcopy np1 yold 1
 275                   (f2cl-lib:array-slice work-%data%
 276                                         double-float
 277                                         (iz0)
 278                                         ((1
 279                                           (f2cl-lib:int-add
 280                                            (f2cl-lib:int-mul 13
 281                                                              (f2cl-lib:int-add
 282                                                               n
 283                                                               1))
 284                                            (f2cl-lib:int-mul 2 n)
 285                                            lenqr)))
 286                                         work-%offset%)
 287                   1)
 288                  (multiple-value-bind
 289                        (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 290                         var-9 var-10 var-11 var-12 var-13 var-14 var-15)
 291                      (rootns nc nfec iflagc ansre ansae y yp yold ypold a qr
 292                       lenqr pivot work par ipar)
 293                    (declare (ignore var-0 var-3 var-4 var-5 var-6 var-7 var-8
 294                                     var-9 var-10 var-11 var-12 var-13 var-14
 295                                     var-15))
 296                    (setf nfec var-1)
 297                    (setf iflagc var-2))
 298                  (setf nfe nfec)
 299                  (setf iflag 1)
 300                  (if (> iflagc 0) (setf iflag iflagc))
 301                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 302                                ((> jw np1) nil)
 303                    (tagbody
 304                      (setf (f2cl-lib:fref work-%data%
 305                                           (jw)
 306                                           ((1
 307                                             (f2cl-lib:int-add
 308                                              (f2cl-lib:int-mul 13
 309                                                                (f2cl-lib:int-add
 310                                                                 n
 311                                                                 1))
 312                                              (f2cl-lib:int-mul 2 n)
 313                                              lenqr)))
 314                                           work-%offset%)
 315                              (-
 316                               (f2cl-lib:fref y-%data%
 317                                              (jw)
 318                                              ((1 (f2cl-lib:int-add n 1)))
 319                                              y-%offset%)
 320                               (f2cl-lib:fref work-%data%
 321                                              ((f2cl-lib:int-sub
 322                                                (f2cl-lib:int-add iz0 jw)
 323                                                1))
 324                                              ((1
 325                                                (f2cl-lib:int-add
 326                                                 (f2cl-lib:int-mul 13
 327                                                                   (f2cl-lib:int-add
 328                                                                    n
 329                                                                    1))
 330                                                 (f2cl-lib:int-mul 2 n)
 331                                                 lenqr)))
 332                                              work-%offset%)))
 333                     label290))
 334                  (setf arclen (+ (- s hold) (dnrm2 np1 work 1)))
 335                  (go end_label)))
 336              label400))
 337           (setf iflag 3)
 338           (setf arclen s)
 339           (go end_label)
 340          end_label
 341           (return
 342            (values nil
 343                    nil
 344                    iflag
 345                    arcre
 346                    arcae
 347                    ansre
 348                    ansae
 349                    nil
 350                    nil
 351                    nfe
 352                    arclen
 353                    nil
 354                    nil
 355                    nil
 356                    nil
 357                    nil
 358                    nil
 359                    nil
 360                    nil
 361                    nil
 362                    nil)))))))
 363
 364 (in-package #:cl-user)
 365 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 366 (eval-when (:load-toplevel :compile-toplevel :execute)
 367   (setf (gethash 'fortran-to-lisp::fixpns
 368                  fortran-to-lisp::*f2cl-function-info*)
 369           (fortran-to-lisp::make-f2cl-finfo
 370            :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
 371                         (fortran-to-lisp::integer4) (double-float)
 372                         (double-float) (double-float) (double-float)
 373                         (fortran-to-lisp::integer4) (array double-float (*))
 374                         (fortran-to-lisp::integer4) (double-float)
 375                         (array double-float (*)) (array double-float (*))
 376                         (array double-float (*)) (array double-float (*))
 377                         (fortran-to-lisp::integer4)
 378                         (array fortran-to-lisp::integer4 (*))
 379                         (array double-float (*)) (array double-float (*))
 380                         (array double-float (*))
 381                         (array fortran-to-lisp::integer4 (*)))
 382            :return-values '(nil nil fortran-to-lisp::iflag
 383                             fortran-to-lisp::arcre fortran-to-lisp::arcae
 384                             fortran-to-lisp::ansre fortran-to-lisp::ansae nil
 385                             nil fortran-to-lisp::nfe fortran-to-lisp::arclen
 386                             nil nil nil nil nil nil nil nil nil nil)
 387            :calls '(fortran-to-lisp::dnrm2 fortran-to-lisp::dcopy
 388                     fortran-to-lisp::rootns fortran-to-lisp::stepns
 389                     fortran-to-lisp::d1mach))))
 390