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