share/hompack/lisp/rootnf.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* ((limit 20))
  21   (declare (type (f2cl-lib:integer4 20 20) limit) (ignorable limit))
  22   (defun rootnf
  23          (n nfe iflag relerr abserr y yp yold ypold a qr alpha tz pivot w wp
  24           par ipar)
  25     (declare (type (array f2cl-lib:integer4 (*)) ipar)
  26              (type (array double-float (*)) par)
  27              (type (array f2cl-lib:integer4 (*)) pivot)
  28              (type (array double-float (*)) wp w tz alpha qr a ypold yold yp y)
  29              (type (double-float) abserr relerr)
  30              (type (f2cl-lib:integer4) iflag nfe n))
  31     (f2cl-lib:with-multi-array-data
  32         ((y double-float y-%data% y-%offset%)
  33          (yp double-float yp-%data% yp-%offset%)
  34          (yold double-float yold-%data% yold-%offset%)
  35          (ypold double-float ypold-%data% ypold-%offset%)
  36          (a double-float a-%data% a-%offset%)
  37          (qr double-float qr-%data% qr-%offset%)
  38          (alpha double-float alpha-%data% alpha-%offset%)
  39          (tz double-float tz-%data% tz-%offset%)
  40          (w double-float w-%data% w-%offset%)
  41          (wp double-float wp-%data% wp-%offset%)
  42          (pivot f2cl-lib:integer4 pivot-%data% pivot-%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 ((judy 0) (jw 0) (lcode 0) (np1 0) (aerr 0.0) (dels 0.0) (f0 0.0)
  86                (f1 0.0) (fp0 0.0) (fp1 0.0) (qsout 0.0) (rerr 0.0) (s 0.0)
  87                (sa 0.0) (sb 0.0) (sout 0.0) (u 0.0))
  88           (declare (type (f2cl-lib:integer4) judy jw lcode np1)
  89                    (type (double-float) aerr dels f0 f1 fp0 fp1 qsout rerr s sa
  90                                         sb sout u))
  91           (setf u (f2cl-lib:d1mach 4))
  92           (setf rerr (max relerr u))
  93           (setf aerr (max abserr 0.0))
  94           (setf np1 (f2cl-lib:int-add n 1))
  95          label100
  96           (f2cl-lib:fdo (judy 1 (f2cl-lib:int-add judy 1))
  97                         ((> judy limit) nil)
  98             (tagbody
  99               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 100                             ((> jw np1) nil)
 101                 (tagbody
 102                   (setf (f2cl-lib:fref tz-%data%
 103                                        (jw)
 104                                        ((1 (f2cl-lib:int-add n 1)))
 105                                        tz-%offset%)
 106                           (-
 107                            (f2cl-lib:fref y-%data%
 108                                           (jw)
 109                                           ((1 (f2cl-lib:int-add n 1)))
 110                                           y-%offset%)
 111                            (f2cl-lib:fref yold-%data%
 112                                           (jw)
 113                                           ((1 (f2cl-lib:int-add n 1)))
 114                                           yold-%offset%)))
 115                  label110))
 116               (setf dels (dnrm2 np1 tz 1))
 117               (setf sa (coerce 0.0f0 'double-float))
 118               (setf sb dels)
 119               (setf lcode 1)
 120              label130
 121               (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6)
 122                   (root sout qsout sa sb rerr aerr lcode)
 123                 (declare (ignore var-1 var-4 var-5))
 124                 (setf sout var-0)
 125                 (setf sa var-2)
 126                 (setf sb var-3)
 127                 (setf lcode var-6))
 128               (if (> lcode 0) (go label140))
 129               (setf qsout
 130                       (-
 131                        (qofs
 132                         (f2cl-lib:fref yold-%data%
 133                                        (1)
 134                                        ((1 (f2cl-lib:int-add n 1)))
 135                                        yold-%offset%)
 136                         (f2cl-lib:fref ypold-%data%
 137                                        (1)
 138                                        ((1 (f2cl-lib:int-add n 1)))
 139                                        ypold-%offset%)
 140                         (f2cl-lib:fref y-%data%
 141                                        (1)
 142                                        ((1 (f2cl-lib:int-add n 1)))
 143                                        y-%offset%)
 144                         (f2cl-lib:fref yp-%data%
 145                                        (1)
 146                                        ((1 (f2cl-lib:int-add n 1)))
 147                                        yp-%offset%)
 148                         dels sout)
 149                        1.0f0))
 150               (go label130)
 151              label140
 152               (cond
 153                 ((> lcode 2)
 154                  (setf iflag 6)
 155                  (go end_label)))
 156               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 157                             ((> jw np1) nil)
 158                 (tagbody
 159                   (setf (f2cl-lib:fref w-%data%
 160                                        (jw)
 161                                        ((1 (f2cl-lib:int-add n 1)))
 162                                        w-%offset%)
 163                           (qofs
 164                            (f2cl-lib:fref yold-%data%
 165                                           (jw)
 166                                           ((1 (f2cl-lib:int-add n 1)))
 167                                           yold-%offset%)
 168                            (f2cl-lib:fref ypold-%data%
 169                                           (jw)
 170                                           ((1 (f2cl-lib:int-add n 1)))
 171                                           ypold-%offset%)
 172                            (f2cl-lib:fref y-%data%
 173                                           (jw)
 174                                           ((1 (f2cl-lib:int-add n 1)))
 175                                           y-%offset%)
 176                            (f2cl-lib:fref yp-%data%
 177                                           (jw)
 178                                           ((1 (f2cl-lib:int-add n 1)))
 179                                           yp-%offset%)
 180                            dels sa))
 181                  label150))
 182               (multiple-value-bind
 183                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 184                      var-9 var-10 var-11 var-12 var-13)
 185                   (tangnf sa w wp ypold a qr alpha tz pivot nfe n iflag par
 186                    ipar)
 187                 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7
 188                                  var-8 var-10 var-12 var-13))
 189                 (setf sa var-0)
 190                 (setf nfe var-9)
 191                 (setf iflag var-11))
 192               (if (> iflag 0) (go end_label))
 193               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 194                             ((> jw np1) nil)
 195                 (tagbody
 196                   (setf (f2cl-lib:fref w-%data%
 197                                        (jw)
 198                                        ((1 (f2cl-lib:int-add n 1)))
 199                                        w-%offset%)
 200                           (+
 201                            (f2cl-lib:fref w-%data%
 202                                           (jw)
 203                                           ((1 (f2cl-lib:int-add n 1)))
 204                                           w-%offset%)
 205                            (f2cl-lib:fref tz-%data%
 206                                           (jw)
 207                                           ((1 (f2cl-lib:int-add n 1)))
 208                                           tz-%offset%)))
 209                  label160))
 210               (multiple-value-bind
 211                     (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8
 212                      var-9 var-10 var-11 var-12 var-13)
 213                   (tangnf sa w wp ypold a qr alpha tz pivot nfe n iflag par
 214                    ipar)
 215                 (declare (ignore var-1 var-2 var-3 var-4 var-5 var-6 var-7
 216                                  var-8 var-10 var-12 var-13))
 217                 (setf sa var-0)
 218                 (setf nfe var-9)
 219                 (setf iflag var-11))
 220               (if (> iflag 0) (go end_label))
 221               (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 222                             ((> jw np1) nil)
 223                 (tagbody
 224                   (setf (f2cl-lib:fref w-%data%
 225                                        (jw)
 226                                        ((1 (f2cl-lib:int-add n 1)))
 227                                        w-%offset%)
 228                           (+
 229                            (f2cl-lib:fref w-%data%
 230                                           (jw)
 231                                           ((1 (f2cl-lib:int-add n 1)))
 232                                           w-%offset%)
 233                            (f2cl-lib:fref tz-%data%
 234                                           (jw)
 235                                           ((1 (f2cl-lib:int-add n 1)))
 236                                           tz-%offset%)))
 237                  label170))
 238               (cond
 239                 ((and
 240                   (<=
 241                    (abs
 242                     (+ (f2cl-lib:fref w (1) ((1 (f2cl-lib:int-add n 1))))
 243                        (- 1.0f0)))
 244                    (+ rerr aerr))
 245                   (<=
 246                    (dnrm2 n
 247                     (f2cl-lib:array-slice tz
 248                                           double-float
 249                                           (2)
 250                                           ((1 (f2cl-lib:int-add n 1))))
 251                     1)
 252                    (+
 253                     (* rerr
 254                        (dnrm2 n
 255                         (f2cl-lib:array-slice w
 256                                               double-float
 257                                               (2)
 258                                               ((1 (f2cl-lib:int-add n 1))))
 259                         1))
 260                     aerr)))
 261                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 262                                ((> jw np1) nil)
 263                    (tagbody
 264                      (setf (f2cl-lib:fref y-%data%
 265                                           (jw)
 266                                           ((1 (f2cl-lib:int-add n 1)))
 267                                           y-%offset%)
 268                              (f2cl-lib:fref w-%data%
 269                                             (jw)
 270                                             ((1 (f2cl-lib:int-add n 1)))
 271                                             w-%offset%))
 272                     label180))
 273                  (go end_label)))
 274               (cond
 275                 ((>
 276                   (*
 277                    (+ (f2cl-lib:fref yold (1) ((1 (f2cl-lib:int-add n 1))))
 278                       (- 1.0f0))
 279                    (+ (f2cl-lib:fref w (1) ((1 (f2cl-lib:int-add n 1))))
 280                       (- 1.0f0)))
 281                   0.0f0)
 282                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 283                                ((> jw np1) nil)
 284                    (tagbody
 285                      (setf (f2cl-lib:fref yold-%data%
 286                                           (jw)
 287                                           ((1 (f2cl-lib:int-add n 1)))
 288                                           yold-%offset%)
 289                              (f2cl-lib:fref w-%data%
 290                                             (jw)
 291                                             ((1 (f2cl-lib:int-add n 1)))
 292                                             w-%offset%))
 293                      (setf (f2cl-lib:fref ypold-%data%
 294                                           (jw)
 295                                           ((1 (f2cl-lib:int-add n 1)))
 296                                           ypold-%offset%)
 297                              (f2cl-lib:fref wp-%data%
 298                                             (jw)
 299                                             ((1 (f2cl-lib:int-add n 1)))
 300                                             wp-%offset%))
 301                     label200)))
 302                 (t
 303                  (f2cl-lib:fdo (jw 1 (f2cl-lib:int-add jw 1))
 304                                ((> jw np1) nil)
 305                    (tagbody
 306                      (setf (f2cl-lib:fref y-%data%
 307                                           (jw)
 308                                           ((1 (f2cl-lib:int-add n 1)))
 309                                           y-%offset%)
 310                              (f2cl-lib:fref w-%data%
 311                                             (jw)
 312                                             ((1 (f2cl-lib:int-add n 1)))
 313                                             w-%offset%))
 314                      (setf (f2cl-lib:fref yp-%data%
 315                                           (jw)
 316                                           ((1 (f2cl-lib:int-add n 1)))
 317                                           yp-%offset%)
 318                              (f2cl-lib:fref wp-%data%
 319                                             (jw)
 320                                             ((1 (f2cl-lib:int-add n 1)))
 321                                             wp-%offset%))
 322                     label210))))
 323              label300))
 324           (setf iflag 6)
 325           (go end_label)
 326          end_label
 327           (return
 328            (values nil
 329                    nfe
 330                    iflag
 331                    nil
 332                    nil
 333                    nil
 334                    nil
 335                    nil
 336                    nil
 337                    nil
 338                    nil
 339                    nil
 340                    nil
 341                    nil
 342                    nil
 343                    nil
 344                    nil
 345                    nil)))))))
 346
 347 (in-package #:cl-user)
 348 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 349 (eval-when (:load-toplevel :compile-toplevel :execute)
 350   (setf (gethash 'fortran-to-lisp::rootnf
 351                  fortran-to-lisp::*f2cl-function-info*)
 352           (fortran-to-lisp::make-f2cl-finfo
 353            :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 354                         (fortran-to-lisp::integer4) (double-float)
 355                         (double-float) (array double-float (*))
 356                         (array double-float (*)) (array double-float (*))
 357                         (array double-float (*)) (array double-float (*))
 358                         (array double-float (*)) (array double-float (*))
 359                         (array double-float (*))
 360                         (array fortran-to-lisp::integer4 (*))
 361                         (array double-float (*)) (array double-float (*))
 362                         (array double-float (*))
 363                         (array fortran-to-lisp::integer4 (*)))
 364            :return-values '(nil fortran-to-lisp::nfe fortran-to-lisp::iflag nil
 365                             nil nil nil nil nil nil nil nil nil nil nil nil nil
 366                             nil)
 367            :calls '(fortran-to-lisp::dnrm2 fortran-to-lisp::tangnf
 368                     fortran-to-lisp::root fortran-to-lisp::d1mach))))
 369