share/hompack/lisp/fixpqs.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 ((wk 0.0)
  23         (s 0.0)
  24         (relerr 0.0)
  25         (hold 0.0)
  26         (h 0.0)
  27         (abserr 0.0)
  28         (pcgwk 0)
  29         (np1 0)
  30         (limit 0)
  31         (jw 0)
  32         (iter 0)
  33         (iflagc 0)
  34         (start nil)
  35         (crash nil))
  36     (declare (type (double-float) wk s relerr hold h abserr)
  37              (type (f2cl-lib:integer4) pcgwk np1 limit jw iter iflagc)
  38              (type f2cl-lib:logical start crash))
  39     (defun fixpqs
  40            (n y iflag arcre arcae ansre ansae trace$ a nfe arclen yp yold ypold
  41             qr lenqr pivot pp rhovec z0 dz t$ work sspar par ipar)
  42       (declare (type (array f2cl-lib:integer4 (*)) ipar)
  43                (type (array double-float (*)) par)
  44                (type (array double-float (*)) sspar)
  45                (type (array f2cl-lib:integer4 (*)) pivot)
  46                (type (double-float) arclen ansae ansre arcae arcre)
  47                (type (array double-float (*)) work t$ dz z0 rhovec pp qr ypold
  48                                               yold yp a y)
  49                (type (f2cl-lib:integer4) lenqr nfe trace$ iflag n))
  50       (f2cl-lib:with-multi-array-data
  51           ((y double-float y-%data% y-%offset%)
  52            (a double-float a-%data% a-%offset%)
  53            (yp double-float yp-%data% yp-%offset%)
  54            (yold double-float yold-%data% yold-%offset%)
  55            (ypold double-float ypold-%data% ypold-%offset%)
  56            (qr double-float qr-%data% qr-%offset%)
  57            (pp double-float pp-%data% pp-%offset%)
  58            (rhovec double-float rhovec-%data% rhovec-%offset%)
  59            (z0 double-float z0-%data% z0-%offset%)
  60            (dz double-float dz-%data% dz-%offset%)
  61            (t$ double-float t$-%data% t$-%offset%)
  62            (work double-float work-%data% work-%offset%)
  63            (pivot f2cl-lib:integer4 pivot-%data% pivot-%offset%)
  64            (sspar double-float sspar-%data% sspar-%offset%)
  65            (par double-float par-%data% par-%offset%)
  66            (ipar f2cl-lib:integer4 ipar-%data% ipar-%offset%))
  67         (prog ()
  68           (declare)
  69           (if (or (<= n 0) (<= ansre 0.0f0) (< ansae 0.0f0)) (setf iflag 7))
  70           (if (and (>= iflag -2) (<= iflag 0)) (go label10))
  71           (if (= iflag 2) (go label50))
  72           (if (= iflag 3) (go label40))
  73           (setf iflag 7)
  74           (go end_label)
  75          label10
  76           (setf arclen (coerce 0.0f0 'double-float))
  77           (if (<= arcre 0.0f0) (setf arcre (* 0.5f0 (f2cl-lib:fsqrt ansre))))
  78           (if (<= arcae 0.0f0) (setf arcae (* 0.5f0 (f2cl-lib:fsqrt ansae))))
  79           (setf nfe 0)
  80           (setf iflagc iflag)
  81           (setf np1 (f2cl-lib:int-add n 1))
  82           (setf pcgwk (f2cl-lib:int-add (f2cl-lib:int-mul 2 n) 3))
  83           (setf start f2cl-lib:%true%)
  84           (setf crash f2cl-lib:%false%)
  85           (setf relerr arcre)
  86           (setf abserr arcae)
  87           (setf hold (coerce 1.0f0 'double-float))
  88           (setf h (coerce 0.1f0 'double-float))
  89           (setf s (coerce 0.0f0 'double-float))
  90           (setf (f2cl-lib:fref ypold-%data%
  91                                (np1)
  92                                ((1 (f2cl-lib:int-add n 1)))
  93                                ypold-%offset%)
  94                   (coerce 1.0f0 'double-float))
  95           (setf (f2cl-lib:fref y-%data%
  96                                (np1)
  97                                ((1 (f2cl-lib:int-add n 1)))
  98                                y-%offset%)
  99                   (coerce 0.0f0 'double-float))
 100           (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 101                         ((> jw n) nil)
 102             (tagbody
 103               (setf (f2cl-lib:fref ypold-%data%
 104                                    (jw)
 105                                    ((1 (f2cl-lib:int-add n 1)))
 106                                    ypold-%offset%)
 107                       (coerce 0.0f0 'double-float))
 108              label20))
 109           (if
 110            (<= (f2cl-lib:fref sspar-%data% (1) ((1 4)) sspar-%offset%) 0.0f0)
 111            (setf (f2cl-lib:fref sspar-%data% (1) ((1 4)) sspar-%offset%)
 112                    (* (+ (f2cl-lib:fsqrt (+ n 1.0f0)) 4.0f0)
 113                       (f2cl-lib:d1mach 4))))
 114           (if
 115            (<= (f2cl-lib:fref sspar-%data% (2) ((1 4)) sspar-%offset%) 0.0f0)
 116            (setf (f2cl-lib:fref sspar-%data% (2) ((1 4)) sspar-%offset%)
 117                    (coerce 1.0f0 'double-float)))
 118           (if
 119            (<= (f2cl-lib:fref sspar-%data% (3) ((1 4)) sspar-%offset%) 0.0f0)
 120            (setf (f2cl-lib:fref sspar-%data% (3) ((1 4)) sspar-%offset%)
 121                    (coerce 0.1f0 'double-float)))
 122           (if
 123            (<= (f2cl-lib:fref sspar-%data% (4) ((1 4)) sspar-%offset%) 0.0f0)
 124            (setf (f2cl-lib:fref sspar-%data% (4) ((1 4)) sspar-%offset%)
 125                    (coerce 7.0f0 'double-float)))
 126           (cond
 127             ((>= iflagc (f2cl-lib:int-sub 1))
 128              (dcopy n y 1 a 1)))
 129          label40
 130           (setf limit limitd)
 131          label50
 132           (f2cl-lib:fdo (iter 1 (f2cl-lib:int-add iter 1))
 133                         ((> iter limit) nil)
 134             (tagbody
 135               (cond
 136                 ((< (f2cl-lib:fref y (np1) ((1 (f2cl-lib:int-add n 1)))) 0.0f0)
 137                  (setf arclen s)
 138                  (setf iflag 5)
 139                  (go end_label)))
 140               (multiple-value-bind
 141                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 142                      var-9 var-10 var-11 var-12 var-13 var-14 var-15 var-16
 143                      var-17 var-18 var-19 var-20 var-21 var-22 var-23 var-24
 144                      var-25 var-26 var-27)
 145                   (stepqs n nfe iflagc lenqr start crash hold h wk relerr
 146                    abserr s y yp yold ypold a qr pivot pp rhovec z0 dz t$ work
 147                    sspar par ipar)
 148                 (declare (ignore var-0 var-3 var-12 var-13 var-14 var-15 var-16
 149                                  var-17 var-18 var-19 var-20 var-21 var-22
 150                                  var-23 var-24 var-25 var-26 var-27))
 151                 (setf nfe var-1)
 152                 (setf iflagc var-2)
 153                 (setf start var-4)
 154                 (setf crash var-5)
 155                 (setf hold var-6)
 156                 (setf h var-7)
 157                 (setf wk var-8)
 158                 (setf relerr var-9)
 159                 (setf abserr var-10)
 160                 (setf s var-11))
 161               (cond
 162                 ((> trace$ 0)
 163                  (f2cl-lib:fformat trace
 164                                    ("~%" " STEP" 1 (("~5D")) "~3@T" "NFE =" 1
 165                                     (("~5D")) "~3@T" "ARC LENGTH =" 1
 166                                     (("~9,4,0,'*,F")) "~3@T" "LAMBDA =" 1
 167                                     (("~7,4,0,'*,F")) "~5@T" "X vector:" "~%" t
 168                                     ("~1@T" 6 (("~12,4,2,0,'*,,'EE"))) "~%")
 169                                    iter
 170                                    nfe
 171                                    s
 172                                    (f2cl-lib:fref y-%data%
 173                                                   (np1)
 174                                                   ((1 (f2cl-lib:int-add n 1)))
 175                                                   y-%offset%)
 176                                    (do ((jw 1 (f2cl-lib:int-add jw 1))
 177                                         (%ret nil))
 178                                        ((> jw n) (nreverse %ret))
 179                                      (declare (type f2cl-lib:integer4 jw))
 180                                      (push
 181                                       (f2cl-lib:fref y-%data%
 182                                                      (jw)
 183                                                      ((1
 184                                                        (f2cl-lib:int-add n 1)))
 185                                                      y-%offset%)
 186                                       %ret)))))
 187               (cond
 188                 ((> iflagc 0)
 189                  (setf arclen s)
 190                  (setf iflag iflagc)
 191                  (go end_label)))
 192               (cond
 193                 (crash
 194                  (setf iflag 2)
 195                  (cond
 196                    ((< arcre relerr)
 197                     (setf arcre relerr)
 198                     (setf ansre relerr)))
 199                  (if (< arcae abserr) (setf arcae abserr))
 200                  (setf limit (f2cl-lib:int-sub limit iter))
 201                  (go end_label)))
 202               (if
 203                (>=
 204                 (f2cl-lib:fref y-%data%
 205                                (np1)
 206                                ((1 (f2cl-lib:int-add n 1)))
 207                                y-%offset%)
 208                 1.0f0)
 209                (go label500))
 210              label400))
 211           (setf arclen s)
 212           (setf iflag 3)
 213           (go end_label)
 214          label500
 215           (dcopy np1 yold 1 t$ 1)
 216           (multiple-value-bind
 217                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 218                  var-10 var-11 var-12 var-13 var-14 var-15 var-16 var-17 var-18
 219                  var-19)
 220               (rootqs n nfe iflagc lenqr ansre ansae y yp yold ypold a qr pivot
 221                pp rhovec z0 dz
 222                (f2cl-lib:array-slice work-%data%
 223                                      double-float
 224                                      (pcgwk)
 225                                      ((1
 226                                        (f2cl-lib:int-add
 227                                         (f2cl-lib:int-mul 8
 228                                                           (f2cl-lib:int-add n
 229                                                                             1))
 230                                         lenqr)))
 231                                      work-%offset%)
 232                par ipar)
 233             (declare (ignore var-0 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 234                              var-10 var-11 var-12 var-13 var-14 var-15 var-16
 235                              var-17 var-18 var-19))
 236             (setf nfe var-1)
 237             (setf iflagc var-2))
 238           (setf iflag 1)
 239           (if (> iflagc 0) (setf iflag iflagc))
 240           (dcopy np1 y 1 dz 1)
 241           (setf wk (coerce -1.0f0 'double-float))
 242           (daxpy np1 wk t$ 1 dz 1)
 243           (setf arclen (+ (- s hold) (dnrm2 np1 dz 1)))
 244           (go end_label)
 245          end_label
 246           (return
 247            (values nil
 248                    nil
 249                    iflag
 250                    arcre
 251                    arcae
 252                    ansre
 253                    nil
 254                    nil
 255                    nil
 256                    nfe
 257                    arclen
 258                    nil
 259                    nil
 260                    nil
 261                    nil
 262                    nil
 263                    nil
 264                    nil
 265                    nil
 266                    nil
 267                    nil
 268                    nil
 269                    nil
 270                    nil
 271                    nil
 272                    nil)))))))
 273
 274 (in-package #-gcl #:cl-user #+gcl "CL-USER")
 275 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 276 (eval-when (:load-toplevel :compile-toplevel :execute)
 277   (setf (gethash 'fortran-to-lisp::fixpqs
 278                  fortran-to-lisp::*f2cl-function-info*)
 279           (fortran-to-lisp::make-f2cl-finfo
 280            :arg-types '((fortran-to-lisp::integer4) (array double-float (*))
 281                         (fortran-to-lisp::integer4) (double-float)
 282                         (double-float) (double-float) (double-float)
 283                         (fortran-to-lisp::integer4) (array double-float (*))
 284                         (fortran-to-lisp::integer4) (double-float)
 285                         (array double-float (*)) (array double-float (*))
 286                         (array double-float (*)) (array double-float (*))
 287                         (fortran-to-lisp::integer4)
 288                         (array fortran-to-lisp::integer4 (*))
 289                         (array double-float (*)) (array double-float (*))
 290                         (array double-float (*)) (array double-float (*))
 291                         (array double-float (*)) (array double-float (*))
 292                         (array double-float (*)) (array double-float (*))
 293                         (array fortran-to-lisp::integer4 (*)))
 294            :return-values '(nil nil fortran-to-lisp::iflag
 295                             fortran-to-lisp::arcre fortran-to-lisp::arcae
 296                             fortran-to-lisp::ansre nil nil nil
 297                             fortran-to-lisp::nfe fortran-to-lisp::arclen nil
 298                             nil nil nil nil nil nil nil nil nil nil nil nil nil
 299                             nil)
 300            :calls '(fortran-to-lisp::dnrm2 fortran-to-lisp::daxpy
 301                     fortran-to-lisp::dcopy fortran-to-lisp::rootqs
 302                     fortran-to-lisp::stepqs fortran-to-lisp::d1mach))))
 303