share/odepack/src/dprjs.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-2013-11 (20E 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 single-float))
  16
  17 (in-package "ODEPACK")
  18
  19
  20 (defun dprjs (neq y yh nyh ewt ftem savf wk iwk f jac)
  21   (declare (type (f2cl-lib:integer4) nyh)
  22            (type (array double-float (*)) wk savf ftem ewt yh y)
  23            (type (array f2cl-lib:integer4 (*)) iwk neq))
  24   (let ()
  25     (symbol-macrolet ((el0 (aref (dls001-part-0 *dls001-common-block*) 210))
  26                       (h (aref (dls001-part-0 *dls001-common-block*) 211))
  27                       (rc (aref (dls001-part-0 *dls001-common-block*) 215))
  28                       (tn (aref (dls001-part-0 *dls001-common-block*) 216))
  29                       (uround (aref (dls001-part-0 *dls001-common-block*) 217))
  30                       (icf (aref (dls001-part-1 *dls001-common-block*) 12))
  31                       (ierpj (aref (dls001-part-1 *dls001-common-block*) 13))
  32                       (jcur (aref (dls001-part-1 *dls001-common-block*) 15))
  33                       (miter (aref (dls001-part-1 *dls001-common-block*) 26))
  34                       (n (aref (dls001-part-1 *dls001-common-block*) 31))
  35                       (nst (aref (dls001-part-1 *dls001-common-block*) 33))
  36                       (nfe (aref (dls001-part-1 *dls001-common-block*) 34))
  37                       (nje (aref (dls001-part-1 *dls001-common-block*) 35))
  38                       (con0 (aref (dlss01-part-0 *dlss01-common-block*) 0))
  39                       (conmin (aref (dlss01-part-0 *dlss01-common-block*) 1))
  40                       (ccmxj (aref (dlss01-part-0 *dlss01-common-block*) 2))
  41                       (psmall (aref (dlss01-part-0 *dlss01-common-block*) 3))
  42                       (rbig (aref (dlss01-part-0 *dlss01-common-block*) 4))
  43                       (iplost (aref (dlss01-part-1 *dlss01-common-block*) 0))
  44                       (iesp (aref (dlss01-part-1 *dlss01-common-block*) 1))
  45                       (iys (aref (dlss01-part-1 *dlss01-common-block*) 3))
  46                       (iba (aref (dlss01-part-1 *dlss01-common-block*) 4))
  47                       (ibian (aref (dlss01-part-1 *dlss01-common-block*) 5))
  48                       (ibjan (aref (dlss01-part-1 *dlss01-common-block*) 6))
  49                       (ibjgp (aref (dlss01-part-1 *dlss01-common-block*) 7))
  50                       (ipian (aref (dlss01-part-1 *dlss01-common-block*) 8))
  51                       (ipjan (aref (dlss01-part-1 *dlss01-common-block*) 9))
  52                       (ipigp (aref (dlss01-part-1 *dlss01-common-block*) 11))
  53                       (ipr (aref (dlss01-part-1 *dlss01-common-block*) 12))
  54                       (ipc (aref (dlss01-part-1 *dlss01-common-block*) 13))
  55                       (ipic (aref (dlss01-part-1 *dlss01-common-block*) 14))
  56                       (ipisp (aref (dlss01-part-1 *dlss01-common-block*) 15))
  57                       (iprsp (aref (dlss01-part-1 *dlss01-common-block*) 16))
  58                       (ipa (aref (dlss01-part-1 *dlss01-common-block*) 17))
  59                       (msbj (aref (dlss01-part-1 *dlss01-common-block*) 26))
  60                       (nslj (aref (dlss01-part-1 *dlss01-common-block*) 27))
  61                       (ngp (aref (dlss01-part-1 *dlss01-common-block*) 28))
  62                       (nlu (aref (dlss01-part-1 *dlss01-common-block*) 29))
  63                       (nsp (aref (dlss01-part-1 *dlss01-common-block*) 31)))
  64       (f2cl-lib:with-multi-array-data
  65           ((neq f2cl-lib:integer4 neq-%data% neq-%offset%)
  66            (iwk f2cl-lib:integer4 iwk-%data% iwk-%offset%)
  67            (y double-float y-%data% y-%offset%)
  68            (yh double-float yh-%data% yh-%offset%)
  69            (ewt double-float ewt-%data% ewt-%offset%)
  70            (ftem double-float ftem-%data% ftem-%offset%)
  71            (savf double-float savf-%data% savf-%offset%)
  72            (wk double-float wk-%data% wk-%offset%))
  73         (prog ((ng 0) (kmin 0) (kmax 0) (k 0) (jmin 0) (jmax 0) (jok 0) (jj 0)
  74                (j 0) (imul 0) (i 0) (srur 0.0d0) (rcont 0.0d0) (rcon 0.0d0)
  75                (r0 0.0d0) (r 0.0d0) (pij 0.0d0) (hl0 0.0d0) (fac 0.0d0)
  76                (di 0.0d0) (con 0.0d0))
  77           (declare (type (double-float) con di fac hl0 pij r r0 rcon rcont
  78                                         srur)
  79                    (type (f2cl-lib:integer4) i imul j jj jok jmax jmin k kmax
  80                                              kmin ng))
  81           (setf hl0 (* h el0))
  82           (setf con (- hl0))
  83           (if (= miter 3) (go label300))
  84           (setf jok 1)
  85           (if (or (= nst 0) (>= nst (f2cl-lib:int-add nslj msbj)))
  86               (setf jok 0))
  87           (if (and (= icf 1) (< (abs (- rc 1.0d0)) ccmxj)) (setf jok 0))
  88           (if (= icf 2) (setf jok 0))
  89           (if (= jok 1) (go label250))
  90          label20
  91           (setf jcur 1)
  92           (setf nje (f2cl-lib:int-add nje 1))
  93           (setf nslj nst)
  94           (setf iplost 0)
  95           (setf conmin (abs con))
  96           (f2cl-lib:computed-goto (label100 label200) miter)
  97          label100
  98           (setf kmin (f2cl-lib:fref iwk-%data% (ipian) ((1 *)) iwk-%offset%))
  99           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 100                         ((> j n) nil)
 101             (tagbody
 102               (setf kmax
 103                       (f2cl-lib:int-sub
 104                        (f2cl-lib:fref iwk-%data%
 105                                       ((f2cl-lib:int-add ipian j))
 106                                       ((1 *))
 107                                       iwk-%offset%)
 108                        1))
 109               (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 110                             ((> i n) nil)
 111                 (tagbody
 112                  label110
 113                   (setf (f2cl-lib:fref ftem-%data% (i) ((1 *)) ftem-%offset%)
 114                           0.0d0)))
 115               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 116                   (funcall jac
 117                            neq
 118                            tn
 119                            y
 120                            j
 121                            (f2cl-lib:array-slice iwk-%data%
 122                                                  f2cl-lib:integer4
 123                                                  (ipian)
 124                                                  ((1 *))
 125                                                  iwk-%offset%)
 126                            (f2cl-lib:array-slice iwk-%data%
 127                                                  f2cl-lib:integer4
 128                                                  (ipjan)
 129                                                  ((1 *))
 130                                                  iwk-%offset%)
 131                            ftem)
 132                 (declare (ignore var-0 var-2 var-4 var-5 var-6))
 133                 (when var-1
 134                   (setf tn var-1))
 135                 (when var-3
 136                   (setf j var-3)))
 137               (f2cl-lib:fdo (k kmin (f2cl-lib:int-add k 1))
 138                             ((> k kmax) nil)
 139                 (tagbody
 140                   (setf i
 141                           (f2cl-lib:fref iwk-%data%
 142                                          ((f2cl-lib:int-add ibjan k))
 143                                          ((1 *))
 144                                          iwk-%offset%))
 145                   (setf (f2cl-lib:fref wk-%data%
 146                                        ((f2cl-lib:int-add iba k))
 147                                        ((1 *))
 148                                        wk-%offset%)
 149                           (*
 150                            (f2cl-lib:fref ftem-%data%
 151                                           (i)
 152                                           ((1 *))
 153                                           ftem-%offset%)
 154                            con))
 155                   (if (= i j)
 156                       (setf (f2cl-lib:fref wk-%data%
 157                                            ((f2cl-lib:int-add iba k))
 158                                            ((1 *))
 159                                            wk-%offset%)
 160                               (+
 161                                (f2cl-lib:fref wk-%data%
 162                                               ((f2cl-lib:int-add iba k))
 163                                               ((1 *))
 164                                               wk-%offset%)
 165                                1.0d0)))
 166                  label120))
 167               (setf kmin (f2cl-lib:int-add kmax 1))
 168              label130))
 169           (go label290)
 170          label200
 171           (setf fac (dvnorm n savf ewt))
 172           (setf r0 (* 1000.0d0 (abs h) uround n fac))
 173           (if (= r0 0.0d0) (setf r0 1.0d0))
 174           (setf srur (f2cl-lib:fref wk-%data% (1) ((1 *)) wk-%offset%))
 175           (setf jmin (f2cl-lib:fref iwk-%data% (ipigp) ((1 *)) iwk-%offset%))
 176           (f2cl-lib:fdo (ng 1 (f2cl-lib:int-add ng 1))
 177                         ((> ng ngp) nil)
 178             (tagbody
 179               (setf jmax
 180                       (f2cl-lib:int-sub
 181                        (f2cl-lib:fref iwk-%data%
 182                                       ((f2cl-lib:int-add ipigp ng))
 183                                       ((1 *))
 184                                       iwk-%offset%)
 185                        1))
 186               (f2cl-lib:fdo (j jmin (f2cl-lib:int-add j 1))
 187                             ((> j jmax) nil)
 188                 (tagbody
 189                   (setf jj
 190                           (f2cl-lib:fref iwk-%data%
 191                                          ((f2cl-lib:int-add ibjgp j))
 192                                          ((1 *))
 193                                          iwk-%offset%))
 194                   (setf r
 195                           (max
 196                            (* srur
 197                               (abs
 198                                (f2cl-lib:fref y-%data%
 199                                               (jj)
 200                                               ((1 *))
 201                                               y-%offset%)))
 202                            (/ r0
 203                               (f2cl-lib:fref ewt-%data%
 204                                              (jj)
 205                                              ((1 *))
 206                                              ewt-%offset%))))
 207                  label210
 208                   (setf (f2cl-lib:fref y-%data% (jj) ((1 *)) y-%offset%)
 209                           (+ (f2cl-lib:fref y-%data% (jj) ((1 *)) y-%offset%)
 210                              r))))
 211               (multiple-value-bind (var-0 var-1 var-2 var-3)
 212                   (funcall f neq tn y ftem)
 213                 (declare (ignore var-0 var-2 var-3))
 214                 (when var-1
 215                   (setf tn var-1)))
 216               (f2cl-lib:fdo (j jmin (f2cl-lib:int-add j 1))
 217                             ((> j jmax) nil)
 218                 (tagbody
 219                   (setf jj
 220                           (f2cl-lib:fref iwk-%data%
 221                                          ((f2cl-lib:int-add ibjgp j))
 222                                          ((1 *))
 223                                          iwk-%offset%))
 224                   (setf (f2cl-lib:fref y-%data% (jj) ((1 *)) y-%offset%)
 225                           (f2cl-lib:fref yh-%data%
 226                                          (jj 1)
 227                                          ((1 nyh) (1 *))
 228                                          yh-%offset%))
 229                   (setf r
 230                           (max
 231                            (* srur
 232                               (abs
 233                                (f2cl-lib:fref y-%data%
 234                                               (jj)
 235                                               ((1 *))
 236                                               y-%offset%)))
 237                            (/ r0
 238                               (f2cl-lib:fref ewt-%data%
 239                                              (jj)
 240                                              ((1 *))
 241                                              ewt-%offset%))))
 242                   (setf fac (/ (- hl0) r))
 243                   (setf kmin
 244                           (f2cl-lib:fref iwk-%data%
 245                                          ((f2cl-lib:int-add ibian jj))
 246                                          ((1 *))
 247                                          iwk-%offset%))
 248                   (setf kmax
 249                           (f2cl-lib:int-sub
 250                            (f2cl-lib:fref iwk-%data%
 251                                           ((f2cl-lib:int-add ibian jj 1))
 252                                           ((1 *))
 253                                           iwk-%offset%)
 254                            1))
 255                   (f2cl-lib:fdo (k kmin (f2cl-lib:int-add k 1))
 256                                 ((> k kmax) nil)
 257                     (tagbody
 258                       (setf i
 259                               (f2cl-lib:fref iwk-%data%
 260                                              ((f2cl-lib:int-add ibjan k))
 261                                              ((1 *))
 262                                              iwk-%offset%))
 263                       (setf (f2cl-lib:fref wk-%data%
 264                                            ((f2cl-lib:int-add iba k))
 265                                            ((1 *))
 266                                            wk-%offset%)
 267                               (*
 268                                (-
 269                                 (f2cl-lib:fref ftem-%data%
 270                                                (i)
 271                                                ((1 *))
 272                                                ftem-%offset%)
 273                                 (f2cl-lib:fref savf-%data%
 274                                                (i)
 275                                                ((1 *))
 276                                                savf-%offset%))
 277                                fac))
 278                       (if (= i jj)
 279                           (setf (f2cl-lib:fref wk-%data%
 280                                                ((f2cl-lib:int-add iba k))
 281                                                ((1 *))
 282                                                wk-%offset%)
 283                                   (+
 284                                    (f2cl-lib:fref wk-%data%
 285                                                   ((f2cl-lib:int-add iba k))
 286                                                   ((1 *))
 287                                                   wk-%offset%)
 288                                    1.0d0)))
 289                      label220))
 290                  label230))
 291               (setf jmin (f2cl-lib:int-add jmax 1))
 292              label240))
 293           (setf nfe (f2cl-lib:int-add nfe ngp))
 294           (go label290)
 295          label250
 296           (setf jcur 0)
 297           (setf rcon (/ con con0))
 298           (setf rcont (/ (abs con) conmin))
 299           (if (and (> rcont rbig) (= iplost 1)) (go label20))
 300           (setf kmin (f2cl-lib:fref iwk-%data% (ipian) ((1 *)) iwk-%offset%))
 301           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 302                         ((> j n) nil)
 303             (tagbody
 304               (setf kmax
 305                       (f2cl-lib:int-sub
 306                        (f2cl-lib:fref iwk-%data%
 307                                       ((f2cl-lib:int-add ipian j))
 308                                       ((1 *))
 309                                       iwk-%offset%)
 310                        1))
 311               (f2cl-lib:fdo (k kmin (f2cl-lib:int-add k 1))
 312                             ((> k kmax) nil)
 313                 (tagbody
 314                   (setf i
 315                           (f2cl-lib:fref iwk-%data%
 316                                          ((f2cl-lib:int-add ibjan k))
 317                                          ((1 *))
 318                                          iwk-%offset%))
 319                   (setf pij
 320                           (f2cl-lib:fref wk-%data%
 321                                          ((f2cl-lib:int-add iba k))
 322                                          ((1 *))
 323                                          wk-%offset%))
 324                   (if (/= i j) (go label260))
 325                   (setf pij (- pij 1.0d0))
 326                   (if (>= (abs pij) psmall) (go label260))
 327                   (setf iplost 1)
 328                   (setf conmin (min (abs con0) conmin))
 329                  label260
 330                   (setf pij (* pij rcon))
 331                   (if (= i j) (setf pij (+ pij 1.0d0)))
 332                   (setf (f2cl-lib:fref wk-%data%
 333                                        ((f2cl-lib:int-add iba k))
 334                                        ((1 *))
 335                                        wk-%offset%)
 336                           pij)
 337                  label270))
 338               (setf kmin (f2cl-lib:int-add kmax 1))
 339              label275))
 340          label290
 341           (setf nlu (f2cl-lib:int-add nlu 1))
 342           (setf con0 con)
 343           (setf ierpj 0)
 344           (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 345                         ((> i n) nil)
 346             (tagbody
 347              label295
 348               (setf (f2cl-lib:fref ftem-%data% (i) ((1 *)) ftem-%offset%)
 349                       0.0d0)))
 350           (multiple-value-bind
 351                 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 352                  var-10 var-11 var-12 var-13 var-14)
 353               (cdrv n
 354                (f2cl-lib:array-slice iwk-%data%
 355                                      f2cl-lib:integer4
 356                                      (ipr)
 357                                      ((1 *))
 358                                      iwk-%offset%)
 359                (f2cl-lib:array-slice iwk-%data%
 360                                      f2cl-lib:integer4
 361                                      (ipc)
 362                                      ((1 *))
 363                                      iwk-%offset%)
 364                (f2cl-lib:array-slice iwk-%data%
 365                                      f2cl-lib:integer4
 366                                      (ipic)
 367                                      ((1 *))
 368                                      iwk-%offset%)
 369                (f2cl-lib:array-slice iwk-%data%
 370                                      f2cl-lib:integer4
 371                                      (ipian)
 372                                      ((1 *))
 373                                      iwk-%offset%)
 374                (f2cl-lib:array-slice iwk-%data%
 375                                      f2cl-lib:integer4
 376                                      (ipjan)
 377                                      ((1 *))
 378                                      iwk-%offset%)
 379                (f2cl-lib:array-slice wk-%data%
 380                                      double-float
 381                                      (ipa)
 382                                      ((1 *))
 383                                      wk-%offset%)
 384                ftem ftem nsp
 385                (f2cl-lib:array-slice iwk-%data%
 386                                      f2cl-lib:integer4
 387                                      (ipisp)
 388                                      ((1 *))
 389                                      iwk-%offset%)
 390                (f2cl-lib:array-slice wk-%data%
 391                                      double-float
 392                                      (iprsp)
 393                                      ((1 *))
 394                                      wk-%offset%)
 395                iesp 2 iys)
 396             (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7
 397                              var-8 var-9 var-10 var-11 var-13))
 398             (setf iesp var-12)
 399             (setf iys var-14))
 400           (if (= iys 0) (go end_label))
 401           (setf imul (the f2cl-lib:integer4 (truncate (- iys 1) n)))
 402           (setf ierpj -2)
 403           (if (= imul 8) (setf ierpj 1))
 404           (if (= imul 10) (setf ierpj -1))
 405           (go end_label)
 406          label300
 407           (setf jcur 1)
 408           (setf nje (f2cl-lib:int-add nje 1))
 409           (setf (f2cl-lib:fref wk-%data% (2) ((1 *)) wk-%offset%) hl0)
 410           (setf ierpj 0)
 411           (setf r (* el0 0.1d0))
 412           (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 413                         ((> i n) nil)
 414             (tagbody
 415              label310
 416               (setf (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
 417                       (+ (f2cl-lib:fref y-%data% (i) ((1 *)) y-%offset%)
 418                          (* r
 419                             (-
 420                              (* h
 421                                 (f2cl-lib:fref savf-%data%
 422                                                (i)
 423                                                ((1 *))
 424                                                savf-%offset%))
 425                              (f2cl-lib:fref yh-%data%
 426                                             (i 2)
 427                                             ((1 nyh) (1 *))
 428                                             yh-%offset%)))))))
 429           (multiple-value-bind (var-0 var-1 var-2 var-3)
 430               (funcall f
 431                        neq
 432                        tn
 433                        y
 434                        (f2cl-lib:array-slice wk-%data%
 435                                              double-float
 436                                              (3)
 437                                              ((1 *))
 438                                              wk-%offset%))
 439             (declare (ignore var-0 var-2 var-3))
 440             (when var-1
 441               (setf tn var-1)))
 442           (setf nfe (f2cl-lib:int-add nfe 1))
 443           (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 444                         ((> i n) nil)
 445             (tagbody
 446               (setf r0
 447                       (-
 448                        (* h
 449                           (f2cl-lib:fref savf-%data%
 450                                          (i)
 451                                          ((1 *))
 452                                          savf-%offset%))
 453                        (f2cl-lib:fref yh-%data%
 454                                       (i 2)
 455                                       ((1 nyh) (1 *))
 456                                       yh-%offset%)))
 457               (setf di
 458                       (- (* 0.1d0 r0)
 459                          (* h
 460                             (-
 461                              (f2cl-lib:fref wk-%data%
 462                                             ((f2cl-lib:int-add i 2))
 463                                             ((1 *))
 464                                             wk-%offset%)
 465                              (f2cl-lib:fref savf-%data%
 466                                             (i)
 467                                             ((1 *))
 468                                             savf-%offset%)))))
 469               (setf (f2cl-lib:fref wk-%data%
 470                                    ((f2cl-lib:int-add i 2))
 471                                    ((1 *))
 472                                    wk-%offset%)
 473                       1.0d0)
 474               (if
 475                (< (abs r0)
 476                   (/ uround
 477                      (f2cl-lib:fref ewt-%data% (i) ((1 *)) ewt-%offset%)))
 478                (go label320))
 479               (if (= (abs di) 0.0d0) (go label330))
 480               (setf (f2cl-lib:fref wk-%data%
 481                                    ((f2cl-lib:int-add i 2))
 482                                    ((1 *))
 483                                    wk-%offset%)
 484                       (/ (* 0.1d0 r0) di))
 485              label320))
 486           (go end_label)
 487          label330
 488           (setf ierpj 2)
 489           (go end_label)
 490          end_label
 491           (return (values nil nil nil nil nil nil nil nil nil nil nil)))))))
 492
 493 (in-package #:cl-user)
 494 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 495 (eval-when (:load-toplevel :compile-toplevel :execute)
 496   (setf (gethash 'fortran-to-lisp::dprjs fortran-to-lisp::*f2cl-function-info*)
 497           (fortran-to-lisp::make-f2cl-finfo
 498            :arg-types '((array fortran-to-lisp::integer4 (*))
 499                         (array double-float (*)) (array double-float (*))
 500                         (fortran-to-lisp::integer4) (array double-float (*))
 501                         (array double-float (*)) (array double-float (*))
 502                         (array double-float (*))
 503                         (array fortran-to-lisp::integer4 (*)) t t)
 504            :return-values '(nil nil nil nil nil nil nil nil nil nil nil)
 505            :calls '(fortran-to-lisp::cdrv fortran-to-lisp::dvnorm))))
 506