share/hompack/lisp/stepns.lisp

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