share/minpack/lisp/hybrd.lisp

   1 ;;; Compiled by f2cl version:
   2 ;;; ("f2cl1.l,v 1.215 2009/04/07 22:05:21 rtoy Exp $"
   3 ;;;  "f2cl2.l,v 1.37 2008/02/22 22:19:33 rtoy Exp $"
   4 ;;;  "f2cl3.l,v 1.6 2008/02/22 22:19:33 rtoy Exp $"
   5 ;;;  "f2cl4.l,v 1.7 2008/02/22 22:19:34 rtoy Exp $"
   6 ;;;  "f2cl5.l,v 1.200 2009/01/19 02:38:17 rtoy Exp $"
   7 ;;;  "f2cl6.l,v 1.48 2008/08/24 00:56:27 rtoy Exp $"
   8 ;;;  "macros.l,v 1.112 2009/01/08 12:57:19 rtoy Exp $")
   9
  10 ;;; Using Lisp CMU Common Lisp 19f (19F)
  11 ;;;
  12 ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls nil)
  13 ;;;           (:coerce-assigns :as-needed) (:array-type ':array)
  14 ;;;           (:array-slicing t) (:declare-common nil)
  15 ;;;           (:float-format double-float))
  16
  17 (in-package :minpack)
  18
  19
  20 (let ((one 1.0) (p1 0.1) (p5 0.5) (p001 0.001) (p0001 1.0e-4) (zero 0.0))
  21   (declare (type (double-float) one p1 p5 p001 p0001 zero))
  22   (defun hybrd
  23          (fcn n x fvec xtol maxfev ml mu epsfcn diag mode factor nprint info
  24           nfev fjac ldfjac r lr qtf wa1 wa2 wa3 wa4)
  25     (declare (type (double-float) factor epsfcn xtol)
  26              (type (array double-float (*)) wa4 wa3 wa2 wa1 qtf r fjac diag
  27                                             fvec x)
  28              (type (f2cl-lib:integer4) lr ldfjac nfev info nprint mode mu ml
  29                                        maxfev n))
  30     (f2cl-lib:with-multi-array-data
  31         ((x double-float x-%data% x-%offset%)
  32          (fvec double-float fvec-%data% fvec-%offset%)
  33          (diag double-float diag-%data% diag-%offset%)
  34          (fjac double-float fjac-%data% fjac-%offset%)
  35          (r double-float r-%data% r-%offset%)
  36          (qtf double-float qtf-%data% qtf-%offset%)
  37          (wa1 double-float wa1-%data% wa1-%offset%)
  38          (wa2 double-float wa2-%data% wa2-%offset%)
  39          (wa3 double-float wa3-%data% wa3-%offset%)
  40          (wa4 double-float wa4-%data% wa4-%offset%))
  41       (prog ((actred 0.0) (delta 0.0) (epsmch 0.0) (fnorm 0.0) (fnorm1 0.0)
  42              (pnorm 0.0) (prered 0.0) (ratio 0.0) (sum 0.0) (temp 0.0)
  43              (xnorm 0.0) (jeval nil) (sing nil)
  44              (iwa (make-array 1 :element-type 'f2cl-lib:integer4)) (i 0)
  45              (iflag 0) (iter 0) (j 0) (jm1 0) (l 0) (msum 0) (ncfail 0)
  46              (ncsuc 0) (nslow1 0) (nslow2 0))
  47         (declare (type (f2cl-lib:integer4) nslow2 nslow1 ncsuc ncfail msum l
  48                                            jm1 j iter iflag i)
  49                  (type (array f2cl-lib:integer4 (1)) iwa)
  50                  (type f2cl-lib:logical sing jeval)
  51                  (type (double-float) xnorm temp sum ratio prered pnorm fnorm1
  52                                       fnorm epsmch delta actred))
  53         '"     **********"
  54         '""
  55         '"     subroutine hybrd"
  56         '""
  57         '"     the purpose of hybrd is to find a zero of a system of"
  58         '"     n nonlinear functions in n variables by a modification"
  59         '"     of the powell hybrid method. the user must provide a"
  60         '"     subroutine which calculates the functions. the jacobian is"
  61         '"     then calculated by a forward-difference approximation."
  62         '""
  63         '"     the subroutine statement is"
  64         '""
  65         '"       subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,"
  66         '"                        diag,mode,factor,nprint,info,nfev,fjac,"
  67         '"                        ldfjac,r,lr,qtf,wa1,wa2,wa3,wa4)"
  68         '""
  69         '"     where"
  70         '""
  71         '"       fcn is the name of the user-supplied subroutine which"
  72         '"         calculates the functions. fcn must be declared"
  73         '"         in an external statement in the user calling"
  74         '"         program, and should be written as follows."
  75         '""
  76         '"         subroutine fcn(n,x,fvec,iflag)"
  77         '"         integer n,iflag"
  78         '"         double precision x(n),fvec(n)"
  79         '"         ----------"
  80         '"         calculate the functions at x and"
  81         '"         return this vector in fvec."
  82         '"         ---------"
  83         '"         return"
  84         '"         end"
  85         '""
  86         '"         the value of iflag should not be changed by fcn unless"
  87         '"         the user wants to terminate execution of hybrd."
  88         '"         in this case set iflag to a negative integer."
  89         '""
  90         '"       n is a positive integer input variable set to the number"
  91         '"         of functions and variables."
  92         '""
  93         '"       x is an array of length n. on input x must contain"
  94         '"         an initial estimate of the solution vector. on output x"
  95         '"         contains the final estimate of the solution vector."
  96         '""
  97         '"       fvec is an output array of length n which contains"
  98         '"         the functions evaluated at the output x."
  99         '""
 100         '"       xtol is a nonnegative input variable. termination"
 101         '"         occurs when the relative error between two consecutive"
 102         '"         iterates is at most xtol."
 103         '""
 104         '"       maxfev is a positive integer input variable. termination"
 105         '"         occurs when the number of calls to fcn is at least maxfev"
 106         '"         by the end of an iteration."
 107         '""
 108         '"       ml is a nonnegative integer input variable which specifies"
 109         '"         the number of subdiagonals within the band of the"
 110         '"         jacobian matrix. if the jacobian is not banded, set"
 111         '"         ml to at least n - 1."
 112         '""
 113         '"       mu is a nonnegative integer input variable which specifies"
 114         '"         the number of superdiagonals within the band of the"
 115         '"         jacobian matrix. if the jacobian is not banded, set"
 116         '"         mu to at least n - 1."
 117         '""
 118         '"       epsfcn is an input variable used in determining a suitable"
 119         '"         step length for the forward-difference approximation. this"
 120         '"         approximation assumes that the relative errors in the"
 121         '"         functions are of the order of epsfcn. if epsfcn is less"
 122         '"         than the machine precision, it is assumed that the relative"
 123         '"         errors in the functions are of the order of the machine"
 124         '"         precision."
 125         '""
 126         '"       diag is an array of length n. if mode = 1 (see"
 127         '"         below), diag is internally set. if mode = 2, diag"
 128         '"         must contain positive entries that serve as"
 129         '"         multiplicative scale factors for the variables."
 130         '""
 131         '"       mode is an integer input variable. if mode = 1, the"
 132         '"         variables will be scaled internally. if mode = 2,"
 133         '"         the scaling is specified by the input diag. other"
 134         '"         values of mode are equivalent to mode = 1."
 135         '""
 136         '"       factor is a positive input variable used in determining the"
 137         '"         initial step bound. this bound is set to the product of"
 138         '"         factor and the euclidean norm of diag*x if nonzero, or else"
 139         '"         to factor itself. in most cases factor should lie in the"
 140         '"         interval (.1,100.). 100. is a generally recommended value."
 141         '""
 142         '"       nprint is an integer input variable that enables controlled"
 143         '"         printing of iterates if it is positive. in this case,"
 144         '"         fcn is called with iflag = 0 at the beginning of the first"
 145         '"         iteration and every nprint iterations thereafter and"
 146         '"         immediately prior to return, with x and fvec available"
 147         '"         for printing. if nprint is not positive, no special calls"
 148         '"         of fcn with iflag = 0 are made."
 149         '""
 150         '"       info is an integer output variable. if the user has"
 151         '"         terminated execution, info is set to the (negative)"
 152         '"         value of iflag. see description of fcn. otherwise,"
 153         '"         info is set as follows."
 154         '""
 155         '"         info = 0   improper input parameters."
 156         '""
 157         '"         info = 1   relative error between two consecutive iterates"
 158         '"                    is at most xtol."
 159         '""
 160         '"         info = 2   number of calls to fcn has reached or exceeded"
 161         '"                    maxfev."
 162         '""
 163         '"         info = 3   xtol is too small. no further improvement in"
 164         '"                    the approximate solution x is possible."
 165         '""
 166         '"         info = 4   iteration is not making good progress, as"
 167         '"                    measured by the improvement from the last"
 168         '"                    five jacobian evaluations."
 169         '""
 170         '"         info = 5   iteration is not making good progress, as"
 171         '"                    measured by the improvement from the last"
 172         '"                    ten iterations."
 173         '""
 174         '"       nfev is an integer output variable set to the number of"
 175         '"         calls to fcn."
 176         '""
 177         '"       fjac is an output n by n array which contains the"
 178         '"         orthogonal matrix q produced by the qr factorization"
 179         '"         of the final approximate jacobian."
 180         '""
 181         '"       ldfjac is a positive integer input variable not less than n"
 182         '"         which specifies the leading dimension of the array fjac."
 183         '""
 184         '"       r is an output array of length lr which contains the"
 185         '"         upper triangular matrix produced by the qr factorization"
 186         '"         of the final approximate jacobian, stored rowwise."
 187         '""
 188         '"       lr is a positive integer input variable not less than"
 189         '"         (n*(n+1))/2."
 190         '""
 191         '"       qtf is an output array of length n which contains"
 192         '"         the vector (q transpose)*fvec."
 193         '""
 194         '"       wa1, wa2, wa3, and wa4 are work arrays of length n."
 195         '""
 196         '"     subprograms called"
 197         '""
 198         '"       user-supplied ...... fcn"
 199         '""
 200         '"       minpack-supplied ... dogleg,dpmpar,enorm,fdjac1,"
 201         '"                            qform,qrfac,r1mpyq,r1updt"
 202         '""
 203         '"       fortran-supplied ... dabs,dmax1,dmin1,min0,mod"
 204         '""
 205         '"     argonne national laboratory. minpack project. march 1980."
 206         '"     burton s. garbow, kenneth e. hillstrom, jorge j. more"
 207         '""
 208         '"     **********"
 209         '""
 210         '"     epsmch is the machine precision."
 211         '""
 212         (setf epsmch (dpmpar 1))
 213         '""
 214         (setf info 0)
 215         (setf iflag 0)
 216         (setf nfev 0)
 217         '""
 218         '"     check the input parameters for errors."
 219         '""
 220         (if
 221          (or (<= n 0)
 222              (< xtol zero)
 223              (<= maxfev 0)
 224              (< ml 0)
 225              (< mu 0)
 226              (<= factor zero)
 227              (< ldfjac n)
 228              (< lr (the f2cl-lib:integer4 (truncate (* n (+ n 1)) 2))))
 229          (go label300))
 230         (if (/= mode 2) (go label20))
 231         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 232                       ((> j n) nil)
 233           (tagbody
 234             (if (<= (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%) zero)
 235                 (go label300))
 236            label10))
 237        label20
 238         '""
 239         '"     evaluate the function at the starting point"
 240         '"     and calculate its norm."
 241         '""
 242         (setf iflag 1)
 243         (multiple-value-bind (var-0 var-1 var-2 var-3)
 244             (funcall fcn n x fvec iflag)
 245           (declare (ignore var-1 var-2))
 246           (when var-0
 247             (setf n var-0))
 248           (when var-3
 249             (setf iflag var-3)))
 250         (setf nfev 1)
 251         (if (< iflag 0) (go label300))
 252         (setf fnorm (enorm n fvec))
 253         '""
 254         '"     determine the number of calls to fcn needed to compute"
 255         '"     the jacobian matrix."
 256         '""
 257         (setf msum (f2cl-lib:min0 (f2cl-lib:int-add ml mu 1) n))
 258         '""
 259         '"     initialize iteration counter and monitors."
 260         '""
 261         (setf iter 1)
 262         (setf ncsuc 0)
 263         (setf ncfail 0)
 264         (setf nslow1 0)
 265         (setf nslow2 0)
 266         '""
 267         '"     beginning of the outer loop."
 268         '""
 269        label30
 270         (setf jeval f2cl-lib:%true%)
 271         '""
 272         '"        calculate the jacobian matrix."
 273         '""
 274         (setf iflag 2)
 275         (multiple-value-bind
 276               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
 277                var-10 var-11)
 278             (fdjac1 fcn n x fvec fjac ldfjac iflag ml mu epsfcn wa1 wa2)
 279           (declare (ignore var-0 var-2 var-3 var-4 var-5 var-7 var-8 var-9
 280                            var-10 var-11))
 281           (setf n var-1)
 282           (setf iflag var-6))
 283         (setf nfev (f2cl-lib:int-add nfev msum))
 284         (if (< iflag 0) (go label300))
 285         '""
 286         '"        compute the qr factorization of the jacobian."
 287         '""
 288         (qrfac n n fjac ldfjac f2cl-lib:%false% iwa 1 wa1 wa2 wa3)
 289         '""
 290         '"        on the first iteration and if mode is 1, scale according"
 291         '"        to the norms of the columns of the initial jacobian."
 292         '""
 293         (if (/= iter 1) (go label70))
 294         (if (= mode 2) (go label50))
 295         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 296                       ((> j n) nil)
 297           (tagbody
 298             (setf (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 299                     (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%))
 300             (if (= (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%) zero)
 301                 (setf (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 302                         one))
 303            label40))
 304        label50
 305         '""
 306         '"        on the first iteration, calculate the norm of the scaled x"
 307         '"        and initialize the step bound delta."
 308         '""
 309         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 310                       ((> j n) nil)
 311           (tagbody
 312             (setf (f2cl-lib:fref wa3-%data% (j) ((1 n)) wa3-%offset%)
 313                     (* (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 314                        (f2cl-lib:fref x-%data% (j) ((1 n)) x-%offset%)))
 315            label60))
 316         (setf xnorm (enorm n wa3))
 317         (setf delta (* factor xnorm))
 318         (if (= delta zero) (setf delta factor))
 319        label70
 320         '""
 321         '"        form (q transpose)*fvec and store in qtf."
 322         '""
 323         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 324                       ((> i n) nil)
 325           (tagbody
 326             (setf (f2cl-lib:fref qtf-%data% (i) ((1 n)) qtf-%offset%)
 327                     (f2cl-lib:fref fvec-%data% (i) ((1 n)) fvec-%offset%))
 328            label80))
 329         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 330                       ((> j n) nil)
 331           (tagbody
 332             (if
 333              (=
 334               (f2cl-lib:fref fjac-%data%
 335                              (j j)
 336                              ((1 ldfjac) (1 n))
 337                              fjac-%offset%)
 338               zero)
 339              (go label110))
 340             (setf sum zero)
 341             (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 342                           ((> i n) nil)
 343               (tagbody
 344                 (setf sum
 345                         (+ sum
 346                            (*
 347                             (f2cl-lib:fref fjac-%data%
 348                                            (i j)
 349                                            ((1 ldfjac) (1 n))
 350                                            fjac-%offset%)
 351                             (f2cl-lib:fref qtf-%data%
 352                                            (i)
 353                                            ((1 n))
 354                                            qtf-%offset%))))
 355                label90))
 356             (setf temp
 357                     (/ (- sum)
 358                        (f2cl-lib:fref fjac-%data%
 359                                       (j j)
 360                                       ((1 ldfjac) (1 n))
 361                                       fjac-%offset%)))
 362             (f2cl-lib:fdo (i j (f2cl-lib:int-add i 1))
 363                           ((> i n) nil)
 364               (tagbody
 365                 (setf (f2cl-lib:fref qtf-%data% (i) ((1 n)) qtf-%offset%)
 366                         (+ (f2cl-lib:fref qtf-%data% (i) ((1 n)) qtf-%offset%)
 367                            (*
 368                             (f2cl-lib:fref fjac-%data%
 369                                            (i j)
 370                                            ((1 ldfjac) (1 n))
 371                                            fjac-%offset%)
 372                             temp)))
 373                label100))
 374            label110
 375            label120))
 376         '""
 377         '"        copy the triangular factor of the qr factorization into r."
 378         '""
 379         (setf sing f2cl-lib:%false%)
 380         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 381                       ((> j n) nil)
 382           (tagbody
 383             (setf l j)
 384             (setf jm1 (f2cl-lib:int-sub j 1))
 385             (if (< jm1 1) (go label140))
 386             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 387                           ((> i jm1) nil)
 388               (tagbody
 389                 (setf (f2cl-lib:fref r-%data% (l) ((1 lr)) r-%offset%)
 390                         (f2cl-lib:fref fjac-%data%
 391                                        (i j)
 392                                        ((1 ldfjac) (1 n))
 393                                        fjac-%offset%))
 394                 (setf l (f2cl-lib:int-sub (f2cl-lib:int-add l n) i))
 395                label130))
 396            label140
 397             (setf (f2cl-lib:fref r-%data% (l) ((1 lr)) r-%offset%)
 398                     (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%))
 399             (if (= (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%) zero)
 400                 (setf sing f2cl-lib:%true%))
 401            label150))
 402         '""
 403         '"        accumulate the orthogonal factor in fjac."
 404         '""
 405         (qform n n fjac ldfjac wa1)
 406         '""
 407         '"        rescale if necessary."
 408         '""
 409         (if (= mode 2) (go label170))
 410         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 411                       ((> j n) nil)
 412           (tagbody
 413             (setf (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 414                     (f2cl-lib:dmax1
 415                      (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 416                      (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%)))
 417            label160))
 418        label170
 419         '""
 420         '"        beginning of the inner loop."
 421         '""
 422        label180
 423         '""
 424         '"           if requested, call fcn to enable printing of iterates."
 425         '""
 426         (if (<= nprint 0) (go label190))
 427         (setf iflag 0)
 428         (if (= (mod (f2cl-lib:int-sub iter 1) nprint) 0)
 429             (multiple-value-bind (var-0 var-1 var-2 var-3)
 430                 (funcall fcn n x fvec iflag)
 431               (declare (ignore var-1 var-2))
 432               (when var-0
 433                 (setf n var-0))
 434               (when var-3
 435                 (setf iflag var-3))))
 436         (if (< iflag 0) (go label300))
 437        label190
 438         '""
 439         '"           determine the direction p."
 440         '""
 441         (dogleg n r lr diag qtf delta wa1 wa2 wa3)
 442         '""
 443         '"           store the direction p and x + p. calculate the norm of p."
 444         '""
 445         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 446                       ((> j n) nil)
 447           (tagbody
 448             (setf (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%)
 449                     (- (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%)))
 450             (setf (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%)
 451                     (+ (f2cl-lib:fref x-%data% (j) ((1 n)) x-%offset%)
 452                        (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%)))
 453             (setf (f2cl-lib:fref wa3-%data% (j) ((1 n)) wa3-%offset%)
 454                     (* (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 455                        (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%)))
 456            label200))
 457         (setf pnorm (enorm n wa3))
 458         '""
 459         '"           on the first iteration, adjust the initial step bound."
 460         '""
 461         (if (= iter 1) (setf delta (f2cl-lib:dmin1 delta pnorm)))
 462         '""
 463         '"           evaluate the function at x + p and calculate its norm."
 464         '""
 465         (setf iflag 1)
 466         (multiple-value-bind (var-0 var-1 var-2 var-3)
 467             (funcall fcn n wa2 wa4 iflag)
 468           (declare (ignore var-1 var-2))
 469           (when var-0
 470             (setf n var-0))
 471           (when var-3
 472             (setf iflag var-3)))
 473         (setf nfev (f2cl-lib:int-add nfev 1))
 474         (if (< iflag 0) (go label300))
 475         (setf fnorm1 (enorm n wa4))
 476         '""
 477         '"           compute the scaled actual reduction."
 478         '""
 479         (setf actred (- one))
 480         (if (< fnorm1 fnorm) (setf actred (- one (expt (/ fnorm1 fnorm) 2))))
 481         '""
 482         '"           compute the scaled predicted reduction."
 483         '""
 484         (setf l 1)
 485         (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 486                       ((> i n) nil)
 487           (tagbody
 488             (setf sum zero)
 489             (f2cl-lib:fdo (j i (f2cl-lib:int-add j 1))
 490                           ((> j n) nil)
 491               (tagbody
 492                 (setf sum
 493                         (+ sum
 494                            (* (f2cl-lib:fref r-%data% (l) ((1 lr)) r-%offset%)
 495                               (f2cl-lib:fref wa1-%data%
 496                                              (j)
 497                                              ((1 n))
 498                                              wa1-%offset%))))
 499                 (setf l (f2cl-lib:int-add l 1))
 500                label210))
 501             (setf (f2cl-lib:fref wa3-%data% (i) ((1 n)) wa3-%offset%)
 502                     (+ (f2cl-lib:fref qtf-%data% (i) ((1 n)) qtf-%offset%)
 503                        sum))
 504            label220))
 505         (setf temp (enorm n wa3))
 506         (setf prered zero)
 507         (if (< temp fnorm) (setf prered (- one (expt (/ temp fnorm) 2))))
 508         '""
 509         '"           compute the ratio of the actual to the predicted"
 510         '"           reduction."
 511         '""
 512         (setf ratio zero)
 513         (if (> prered zero) (setf ratio (/ actred prered)))
 514         '""
 515         '"           update the step bound."
 516         '""
 517         (if (>= ratio p1) (go label230))
 518         (setf ncsuc 0)
 519         (setf ncfail (f2cl-lib:int-add ncfail 1))
 520         (setf delta (* p5 delta))
 521         (go label240)
 522        label230
 523         (setf ncfail 0)
 524         (setf ncsuc (f2cl-lib:int-add ncsuc 1))
 525         (if (or (>= ratio p5) (> ncsuc 1))
 526             (setf delta (f2cl-lib:dmax1 delta (/ pnorm p5))))
 527         (if (<= (f2cl-lib:dabs (- ratio one)) p1) (setf delta (/ pnorm p5)))
 528        label240
 529         '""
 530         '"           test for successful iteration."
 531         '""
 532         (if (< ratio p0001) (go label260))
 533         '""
 534         '"           successful iteration. update x, fvec, and their norms."
 535         '""
 536         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 537                       ((> j n) nil)
 538           (tagbody
 539             (setf (f2cl-lib:fref x-%data% (j) ((1 n)) x-%offset%)
 540                     (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%))
 541             (setf (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%)
 542                     (* (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 543                        (f2cl-lib:fref x-%data% (j) ((1 n)) x-%offset%)))
 544             (setf (f2cl-lib:fref fvec-%data% (j) ((1 n)) fvec-%offset%)
 545                     (f2cl-lib:fref wa4-%data% (j) ((1 n)) wa4-%offset%))
 546            label250))
 547         (setf xnorm (enorm n wa2))
 548         (setf fnorm fnorm1)
 549         (setf iter (f2cl-lib:int-add iter 1))
 550        label260
 551         '""
 552         '"           determine the progress of the iteration."
 553         '""
 554         (setf nslow1 (f2cl-lib:int-add nslow1 1))
 555         (if (>= actred p001) (setf nslow1 0))
 556         (if jeval (setf nslow2 (f2cl-lib:int-add nslow2 1)))
 557         (if (>= actred p1) (setf nslow2 0))
 558         '""
 559         '"           test for convergence."
 560         '""
 561         (if (or (<= delta (* xtol xnorm)) (= fnorm zero)) (setf info 1))
 562         (if (/= info 0) (go label300))
 563         '""
 564         '"           tests for termination and stringent tolerances."
 565         '""
 566         (if (>= nfev maxfev) (setf info 2))
 567         (if (<= (* p1 (f2cl-lib:dmax1 (* p1 delta) pnorm)) (* epsmch xnorm))
 568             (setf info 3))
 569         (if (= nslow2 5) (setf info 4))
 570         (if (= nslow1 10) (setf info 5))
 571         (if (/= info 0) (go label300))
 572         '""
 573         '"           criterion for recalculating jacobian approximation"
 574         '"           by forward differences."
 575         '""
 576         (if (= ncfail 2) (go label290))
 577         '""
 578         '"           calculate the rank one modification to the jacobian"
 579         '"           and update qtf if necessary."
 580         '""
 581         (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
 582                       ((> j n) nil)
 583           (tagbody
 584             (setf sum zero)
 585             (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
 586                           ((> i n) nil)
 587               (tagbody
 588                 (setf sum
 589                         (+ sum
 590                            (*
 591                             (f2cl-lib:fref fjac-%data%
 592                                            (i j)
 593                                            ((1 ldfjac) (1 n))
 594                                            fjac-%offset%)
 595                             (f2cl-lib:fref wa4-%data%
 596                                            (i)
 597                                            ((1 n))
 598                                            wa4-%offset%))))
 599                label270))
 600             (setf (f2cl-lib:fref wa2-%data% (j) ((1 n)) wa2-%offset%)
 601                     (/
 602                      (- sum
 603                         (f2cl-lib:fref wa3-%data% (j) ((1 n)) wa3-%offset%))
 604                      pnorm))
 605             (setf (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%)
 606                     (* (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 607                        (/
 608                         (*
 609                          (f2cl-lib:fref diag-%data% (j) ((1 n)) diag-%offset%)
 610                          (f2cl-lib:fref wa1-%data% (j) ((1 n)) wa1-%offset%))
 611                         pnorm)))
 612             (if (>= ratio p0001)
 613                 (setf (f2cl-lib:fref qtf-%data% (j) ((1 n)) qtf-%offset%) sum))
 614            label280))
 615         '""
 616         '"           compute the qr factorization of the updated jacobian."
 617         '""
 618         (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
 619             (r1updt n n r lr wa1 wa2 wa3 sing)
 620           (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6))
 621           (setf sing var-7))
 622         (r1mpyq n n fjac ldfjac wa2 wa3)
 623         (r1mpyq 1 n qtf 1 wa2 wa3)
 624         '""
 625         '"           end of the inner loop."
 626         '""
 627         (setf jeval f2cl-lib:%false%)
 628         (go label180)
 629        label290
 630         '""
 631         '"        end of the outer loop."
 632         '""
 633         (go label30)
 634        label300
 635         '""
 636         '"     termination, either normal or user imposed."
 637         '""
 638         (if (< iflag 0) (setf info iflag))
 639         (setf iflag 0)
 640         (if (> nprint 0)
 641             (multiple-value-bind (var-0 var-1 var-2 var-3)
 642                 (funcall fcn n x fvec iflag)
 643               (declare (ignore var-1 var-2))
 644               (when var-0
 645                 (setf n var-0))
 646               (when var-3
 647                 (setf iflag var-3))))
 648         (go end_label)
 649         '""
 650         '"     last card of subroutine hybrd."
 651         '""
 652        end_label
 653         (return
 654          (values nil
 655                  n
 656                  nil
 657                  nil
 658                  nil
 659                  nil
 660                  nil
 661                  nil
 662                  nil
 663                  nil
 664                  nil
 665                  nil
 666                  nil
 667                  info
 668                  nfev
 669                  nil
 670                  nil
 671                  nil
 672                  nil
 673                  nil
 674                  nil
 675                  nil
 676                  nil
 677                  nil))))))
 678
 679 (in-package #:cl-user)
 680 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
 681 (eval-when (:load-toplevel :compile-toplevel :execute)
 682   (setf (gethash 'fortran-to-lisp::hybrd fortran-to-lisp::*f2cl-function-info*)
 683           (fortran-to-lisp::make-f2cl-finfo
 684            :arg-types '(t (fortran-to-lisp::integer4) (array double-float (*))
 685                         (array double-float (*)) (double-float)
 686                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 687                         (fortran-to-lisp::integer4) (double-float)
 688                         (array double-float (*)) (fortran-to-lisp::integer4)
 689                         (double-float) (fortran-to-lisp::integer4)
 690                         (fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
 691                         (array double-float (*)) (fortran-to-lisp::integer4)
 692                         (array double-float (*)) (fortran-to-lisp::integer4)
 693                         (array double-float (*)) (array double-float (*))
 694                         (array double-float (*)) (array double-float (*))
 695                         (array double-float (*)))
 696            :return-values '(nil fortran-to-lisp::n nil nil nil nil nil nil nil
 697                             nil nil nil nil fortran-to-lisp::info
 698                             fortran-to-lisp::nfev nil nil nil nil nil nil nil
 699                             nil nil)
 700            :calls '(fortran-to-lisp::r1mpyq fortran-to-lisp::r1updt
 701                     fortran-to-lisp::dogleg fortran-to-lisp::qform
 702                     fortran-to-lisp::qrfac fortran-to-lisp::fdjac1
 703                     fortran-to-lisp::enorm fortran-to-lisp::dpmpar))))
 704