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