roots-multi.lisp 18.2 KB
 liam committed Jan 14, 2008 1 2 ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08  Liam Healy committed Aug 31, 2008 3 ;;; Time-stamp: <2008-08-23 23:04:50EDT roots-multi.lisp>  lhealy committed Jul 25, 2008 4 ;;; $Id$  liam committed Jan 14, 2008 5 6 7 8 9 10 11 12 13 14 15  (in-package :gsl) ;;;;**************************************************************************** ;;;; Function definition ;;;;**************************************************************************** (cffi:defcstruct gsl-mfunction ;; See /usr/include/gsl/gsl_multiroots.h "The definition of a function for multiroot finding in GSL." (function :pointer)  Liam Healy committed Aug 31, 2008 16  (dimensions sizet)  liam committed Jan 14, 2008 17 18 19 20 21  (parameters :pointer)) (export 'def-mfunction) (defmacro def-mfunction (name dimensions) "Define a function for multivariate root solving."  liam committed Jan 21, 2008 22  (def-single-function ,name :success-failure :pointer gsl-mfunction  liam committed Jan 14, 2008 23  ((dimensions ,dimensions))  liam committed Jan 21, 2008 24  (gsl-vector-c)))  liam committed Jan 14, 2008 25   liam committed Jan 16, 2008 26 27 28 29 30 31 32 (cffi:defcstruct gsl-mfunction-fdf ;; See /usr/include/gsl/gsl_multiroots.h "The definition of a function and its derivatives for multiroot finding in GSL." (function :pointer) (df :pointer) (fdf :pointer)  Liam Healy committed Aug 31, 2008 33  (dimensions sizet)  liam committed Jan 16, 2008 34 35  (parameters :pointer))  liam committed Jan 14, 2008 36 37 38 39 ;;;;**************************************************************************** ;;;; Initialization ;;;;****************************************************************************  liam committed Feb 03, 2008 40 (defgo mfsolver (type function-derivative initial)  liam committed Feb 03, 2008 41 42 43 44 45 46  (list (allocate-mfsolver ,type (dim0 ,initial)) 'free-mfsolver (lambda (symb) (set-mfsolver ,symb ,function-derivative ,initial))))  liam committed Feb 03, 2008 47 (defgo mfdfsolver (type function-derivative initial)  liam committed Feb 03, 2008 48 49 50 51 52  (list (allocate-mfdfsolver ,type (dim0 ,initial)) 'free-mfdfsolver (lambda (symb) (set-mfdfsolver ,symb ,function-derivative ,initial))))  liam committed Jan 28, 2008 53   liam committed Feb 18, 2008 54 (defmfun allocate-mfsolver (type dimension)  liam committed Jan 14, 2008 55  "gsl_multiroot_fsolver_alloc"  Liam Healy committed Aug 31, 2008 56  ((type :pointer) (dimension sizet))  liam committed Jan 14, 2008 57  :c-return :pointer  liam committed Jan 28, 2008 58  :export nil  liam committed Feb 03, 2008 59  :index (letm mfsolver)  liam committed Feb 04, 2008 60  :documentation ; FDL  liam committed Jan 14, 2008 61 62 63  "Allocate an instance of a solver of the type specified for a system of the specified number of dimensions.")  liam committed Feb 18, 2008 64 (defmfun allocate-mfdfsolver (type dimension)  liam committed Jan 14, 2008 65  "gsl_multiroot_fdfsolver_alloc"  Liam Healy committed Aug 31, 2008 66  ((type :pointer) (dimension sizet))  liam committed Jan 14, 2008 67  :c-return :pointer  liam committed Jan 28, 2008 68  :export nil  liam committed Feb 03, 2008 69  :index (letm mfdfsolver)  liam committed Feb 04, 2008 70  :documentation ; FDL  liam committed Jan 14, 2008 71 72 73  "Allocate an instance of a derivative solver of the type specified for a system of the specified number of dimensions.")  liam committed Feb 18, 2008 74 (defmfun set-mfsolver (solver function initial)  liam committed Jan 14, 2008 75  "gsl_multiroot_fsolver_set"  Liam Healy committed Aug 31, 2008 76  ((solver :pointer) (function :pointer) ((mpointer initial) :pointer))  liam committed Jan 28, 2008 77  :export nil  liam committed Feb 03, 2008 78  :index (letm mfsolver)  liam committed Feb 04, 2008 79  :documentation ; FDL  liam committed Jan 14, 2008 80 81 82  "Set or reset an existing solver to use the function and the initial guess gsl-vector.")  liam committed Feb 18, 2008 83 (defmfun set-mfdfsolver (solver function-derivative initial)  liam committed Jan 14, 2008 84 85  "gsl_multiroot_fdfsolver_set" ((solver :pointer) (function-derivative :pointer)  Liam Healy committed Aug 31, 2008 86  ((mpointer initial) :pointer))  liam committed Jan 28, 2008 87  :export nil  liam committed Feb 03, 2008 88  :index (letm mfdfsolver)  liam committed Feb 04, 2008 89  :documentation ; FDL  liam committed Jan 14, 2008 90 91 92  "Set or reset an existing solver to use the function and derivative (fdf) and the initial guess.")  liam committed Feb 18, 2008 93 (defmfun free-mfsolver (solver)  liam committed Jan 14, 2008 94 95 96  "gsl_multiroot_fsolver_free" ((solver :pointer)) :c-return :void  liam committed Jan 28, 2008 97  :export nil  liam committed Feb 03, 2008 98  :index (letm mfsolver)  liam committed Feb 04, 2008 99  :documentation ; FDL  liam committed Jan 14, 2008 100 101  "Free all the memory associated with the solver.")  liam committed Feb 18, 2008 102 (defmfun free-mfdfsolver (solver)  liam committed Jan 14, 2008 103 104 105  "gsl_multiroot_fdfsolver_free" ((solver :pointer)) :c-return :void  liam committed Jan 28, 2008 106  :export nil  liam committed Feb 03, 2008 107  :index (letm mfdfsolver)  liam committed Feb 04, 2008 108  :documentation ; FDL  liam committed Jan 14, 2008 109 110  "Free all the memory associated with the solver.")  liam committed Feb 18, 2008 111 (defmfun mfsolver-name (solver)  liam committed Jan 14, 2008 112 113 114  "gsl_multiroot_fsolver_name" ((solver :pointer)) :c-return :string  liam committed Feb 04, 2008 115  :documentation ; FDL  liam committed Jan 14, 2008 116 117  "The name of the solver.")  liam committed Feb 18, 2008 118 (defmfun mfdfsolver-name (solver)  liam committed Jan 14, 2008 119 120 121  "gsl_multiroot_fdfsolver_name" ((solver :pointer)) :c-return :string  liam committed Feb 04, 2008 122  :documentation ; FDL  liam committed Jan 14, 2008 123 124 125 126 127 128  "The name of the solver.") ;;;;**************************************************************************** ;;;; Iteration ;;;;****************************************************************************  liam committed Feb 18, 2008 129 (defmfun iterate-mfsolver (solver)  liam committed Jan 14, 2008 130 131  "gsl_multiroot_fsolver_iterate" ((solver :pointer))  liam committed Feb 04, 2008 132  :documentation ; FDL  liam committed Jan 14, 2008 133 134 135 136 137 138 139  "Perform a single iteration of the solver. The following errors may be signalled: :EBADFUNC, the iteration encountered a singular point where the function or its derivative evaluated to infinity or NaN, or :EZERODIV, the derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.")  liam committed Feb 18, 2008 140 (defmfun iterate-mfdfsolver (solver)  liam committed Jan 14, 2008 141 142  "gsl_multiroot_fdfsolver_iterate" ((solver :pointer))  liam committed Feb 04, 2008 143  :documentation ; FDL  liam committed Jan 14, 2008 144 145 146 147 148 149 150  "Perform a single iteration of the solver. The following errors may be signalled: :EBADFUNC, the iteration encountered a singular point where the function or its derivative evaluated to infinity or NaN, or :EZERODIV, the derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.")  liam committed Feb 18, 2008 151 (defmfun mfsolver-root (solver)  liam committed Jan 14, 2008 152 153  "gsl_multiroot_fsolver_root" ((solver :pointer))  Liam Healy committed Aug 31, 2008 154 155  :c-return :pointer :return (:c-return)  liam committed Feb 04, 2008 156  :documentation ; FDL  liam committed Jan 14, 2008 157 158  "The current estimate of the root for the solver.")  liam committed Feb 18, 2008 159 (defmfun mfdfsolver-root (solver)  liam committed Jan 14, 2008 160 161  "gsl_multiroot_fdfsolver_root" ((solver :pointer))  Liam Healy committed Aug 31, 2008 162 163  :c-return :pointer :return (:c-return)  liam committed Jan 14, 2008 164 165 166  :documentation "The current estimate of the root for the solver.")  liam committed Feb 18, 2008 167 (defmfun mfsolver-f (solver)  liam committed Jan 14, 2008 168 169  "gsl_multiroot_fsolver_f" ((solver :pointer))  Liam Healy committed Aug 31, 2008 170 171  :c-return :pointer :return (:c-return)  liam committed Feb 04, 2008 172  :documentation ; FDL  liam committed Jan 14, 2008 173 174  "The function value f(x) at the current estimate x of the root for the solver.")  liam committed Feb 18, 2008 175 (defmfun mfdfsolver-f (solver)  liam committed Jan 14, 2008 176 177  "gsl_multiroot_fdfsolver_f" ((solver :pointer))  Liam Healy committed Aug 31, 2008 178 179  :c-return :pointer :return (:c-return)  liam committed Feb 04, 2008 180  :documentation ; FDL  liam committed Jan 14, 2008 181 182  "The function value f(x) at the current estimate x of the root for the solver.")  liam committed Feb 18, 2008 183 (defmfun mfsolver-dx (solver)  liam committed Jan 14, 2008 184 185  "gsl_multiroot_fsolver_dx" ((solver :pointer))  Liam Healy committed Aug 31, 2008 186 187  :c-return :pointer :return (:c-return)  liam committed Feb 04, 2008 188  :documentation ; FDL  liam committed Jan 14, 2008 189 190  "The last step dx taken by the solver.")  liam committed Feb 18, 2008 191 (defmfun mfdfsolver-dx (solver)  liam committed Jan 14, 2008 192 193  "gsl_multiroot_fsolver_dx" ((solver :pointer))  Liam Healy committed Aug 31, 2008 194 195  :c-return :pointer :return (:c-return)  liam committed Feb 04, 2008 196  :documentation ; FDL  liam committed Jan 14, 2008 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217  "The last step dx taken by the solver.") ;;;;**************************************************************************** ;;;; Search stopping conditions ;;;;**************************************************************************** ;;; The only place we need to pick apart the gsl_multiroot_fsolver ;;; struct is here. We could use mfsolver-dx etc., but then we'd have ;;; to discriminate on mfsolver vs. mfdfsolver. (cffi:defcstruct gsl-multiroot-fsolver ;; See /usr/include/gsl/gsl_multiroots.h (type :pointer) (function :pointer) (x :pointer) (f :pointer) (dx :pointer) (state :pointer)) (defun multiroot-slot (solver slot) (cffi:foreign-slot-value solver 'gsl-multiroot-fsolver slot))  liam committed Feb 18, 2008 218 (defmfun multiroot-test-delta (solver absolute-error relative-error)  liam committed Jan 14, 2008 219 220 221 222 223  "gsl_multiroot_test_delta" (((multiroot-slot solver 'dx) :pointer) ((multiroot-slot solver 'x) :pointer) (absolute-error :double) (relative-error :double)) :c-return :success-continue  liam committed Feb 04, 2008 224  :documentation ; FDL  liam committed Jan 14, 2008 225 226 227 228 229 230 231  "Test for the convergence of the sequence by comparing the last step dx with the absolute error and relative errors given to the current position x. The test returns T if the following condition is achieved: |dx_i| < epsabs + epsrel |x_i| for each component of x and returns NIL otherwise.")  liam committed Feb 18, 2008 232 (defmfun multiroot-test-residual (solver absolute-error)  liam committed Jan 14, 2008 233 234 235  "gsl_multiroot_test_residual" (((multiroot-slot solver 'f) :pointer) (absolute-error :double)) :c-return :success-failure  liam committed Feb 04, 2008 236  :documentation ; FDL  liam committed Jan 14, 2008 237 238 239 240 241 242 243 244 245 246 247 248  "Test the residual value f against the absolute error, returning T if the following condition is achieved: \sum_i |f_i| < absolute_error and returns NIL otherwise. This criterion is suitable for situations where the precise location of the root x is unimportant provided a value can be found where the residual is small enough.") ;;;;**************************************************************************** ;;;; Algorithms using derivatives ;;;;****************************************************************************  liam committed Feb 18, 2008 249 (defmpar *powells-hybrid* "gsl_multiroot_fdfsolver_hybridsj"  liam committed Feb 04, 2008 250  ;; FDL  liam committed Jan 14, 2008 251  "This is a modified version of Powell's Hybrid method as implemented in  liam committed Feb 18, 2008 252  the hybrj algorithm in minpack. Minpack was written by Jorge  liam committed Jan 14, 2008 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294  J. More, Burton S. Garbow and Kenneth E. Hillstrom. The Hybrid algorithm retains the fast convergence of Newton's method but will also reduce the residual when Newton's method is unreliable. The algorithm uses a generalized trust region to keep each step under control. In order to be accepted a proposed new position x' must satisfy the condition |D (x' - x)| < \delta, where D is a diagonal scaling matrix and \delta is the size of the trust region. The components of D are computed internally, using the column norms of the Jacobian to estimate the sensitivity of the residual to each component of x. This improves the behavior of the algorithm for badly scaled functions. On each iteration the algorithm first determines the standard Newton step by solving the system J dx = - f. If this step falls inside the trust region it is used as a trial step in the next stage. If not, the algorithm uses the linear combination of the Newton and gradient directions which is predicted to minimize the norm of the function while staying inside the trust region, dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2. This combination of Newton and gradient directions is referred to as a dogleg step. The proposed step is now tested by evaluating the function at the resulting point, x'. If the step reduces the norm of the function sufficiently then it is accepted and size of the trust region is increased. If the proposed step fails to improve the solution then the size of the trust region is decreased and another trial step is computed. The speed of the algorithm is increased by computing the changes to the Jacobian approximately, using a rank-1 update. If two successive attempts fail to reduce the residual then the full Jacobian is recomputed. The algorithm also monitors the progress of the solution and returns an error if several steps fail to make any improvement, :ENOPROG the iteration is not making any progress, preventing the algorithm from continuing. :ENOPROGJ re-evaluations of the Jacobian indicate that the iteration is not making any progress, preventing the algorithm from continuing.")  liam committed Feb 18, 2008 295 (defmpar *powells-hybrid-unscaled* "gsl_multiroot_fdfsolver_hybridj"  liam committed Feb 04, 2008 296  ;; FDL  liam committed Jan 14, 2008 297 298 299 300 301  "This algorithm is an unscaled version of *powells-hybrid*. The steps are controlled by a spherical trust region |x' - x| < \delta, instead of a generalized region. This can be useful if the generalized region estimated by *powells-hybrid* is inappropriate.")  liam committed Feb 18, 2008 302 (defmpar *newton-mfdfsolver* "gsl_multiroot_fdfsolver_newton"  liam committed Feb 04, 2008 303  ;; FDL  liam committed Jan 14, 2008 304 305 306 307 308 309 310 311 312 313 314 315  "Newton's Method is the standard root-polishing algorithm. The algorithm begins with an initial guess for the location of the solution. On each iteration a linear approximation to the function F is used to estimate the step which will zero all the components of the residual. The iteration is defined by the following sequence, x -> x' = x - J{-1} f(x) where the Jacobian matrix J is computed from the derivative functions provided by f. The step dx is obtained by solving the linear system, J dx = - f(x) using LU decomposition.")  liam committed Feb 18, 2008 316 (defmpar *gnewton-mfdfsolver* "gsl_multiroot_fdfsolver_gnewton"  liam committed Feb 04, 2008 317  ;; FDL  liam committed Jan 14, 2008 318 319 320 321 322 323 324 325 326 327 328 329 330  "A modified version of Newton's method which attempts to improve global convergence by requiring every step to reduce the Euclidean norm of the residual, |f(x)|. If the Newton step leads to an increase in the norm then a reduced step of relative size, t = (\sqrt(1 + 6 r) - 1) / (3 r) is proposed, with r being the ratio of norms |f(x')|^2/|f(x)|^2. This procedure is repeated until a suitable step size is found.") ;;;;**************************************************************************** ;;;; Algorithms without derivatives ;;;;****************************************************************************  liam committed Feb 18, 2008 331 (defmpar *hybrid-scaled* "gsl_multiroot_fsolver_hybrids"  liam committed Feb 04, 2008 332 333  ;; FDL "This is a version of the Hybrid algorithm which replaces calls to the  liam committed Jan 14, 2008 334 335 336 337 338 339  Jacobian function by its finite difference approximation. The finite difference approximation is computed using gsl_multiroots_fdjac with a relative step size of GSL_SQRT_DBL_EPSILON.") ;; Where is this function and parameter? Only thing that shows in the ;; library is gsl_multiroot_fdjacobian.  liam committed Feb 18, 2008 340 (defmpar *hybrid-unscaled* "gsl_multiroot_fsolver_hybrid"  liam committed Feb 04, 2008 341  ;; FDL  liam committed Jan 14, 2008 342 343 344  "A finite difference version of the Hybrid algorithm without internal scaling.")  liam committed Feb 18, 2008 345 (defmpar *discrete-newton* "gsl_multiroot_fsolver_dnewton"  liam committed Feb 04, 2008 346  ;; FDL  liam committed Jan 14, 2008 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361  "The discrete Newton algorithm is the simplest method of solving a multidimensional system. It uses the Newton iteration x -> x - J^{-1} f(x) where the Jacobian matrix J is approximated by taking finite differences of the function f. The approximation scheme used by this implementation is J_{ij} = (f_i(x + \delta_j) - f_i(x)) / \delta_j where \delta_j is a step of size \sqrt\epsilon |x_j| with \epsilon being the machine precision (\epsilon \approx 2.22 \times 10^-16}). The order of convergence of Newton's algorithm is quadratic, but the finite differences require n^2 function evaluations on each iteration. The algorithm may become unstable if the finite differences are not a good approximation to the true derivatives.")  liam committed Feb 18, 2008 362 (defmpar *broyden* "gsl_multiroot_fsolver_broyden"  liam committed Feb 04, 2008 363  ;; FDL  liam committed Jan 14, 2008 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388  "The Broyden algorithm is a version of the discrete Newton algorithm which attempts to avoids the expensive update of the Jacobian matrix on each iteration. The changes to the Jacobian are also approximated, using a rank-1 update, J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df where the vectors dx and df are the changes in x and f. On the first iteration the inverse Jacobian is estimated using finite differences, as in the discrete Newton algorithm. This approximation gives a fast update but is unreliable if the changes are not small, and the estimate of the inverse Jacobian becomes worse as time passes. The algorithm has a tendency to become unstable unless it starts close to the root. The Jacobian is refreshed if this instability is detected (consult the source for details). This algorithm is included only for demonstration purposes, and is not recommended for serious use.") ;;;;**************************************************************************** ;;;; Examples ;;;;**************************************************************************** (defparameter *powell-A* 1.0d4) (defun powell (argument return) "Powell's test function."  liam committed Feb 24, 2008 389 390  (setf (maref return 0) (- (* *powell-A* (maref argument 0) (maref argument 1))  liam committed Jan 14, 2008 391  1)  liam committed Feb 24, 2008 392 393  (maref return 1) (+ (exp (- (maref argument 0))) (exp (- (maref argument 1)))  liam committed Jan 14, 2008 394 395 396 397 398 399 400 401 402  (- (1+ (/ *powell-A*)))))) ;;; (def-mfunction powell 2) ;;; This is the example given in Sec. 34.8. (defparameter *rosenbrock-a* 1.0d0) (defparameter *rosenbrock-b* 10.0d0)  liam committed Jan 17, 2008 403 #|  liam committed Jan 16, 2008 404 ;;; One alternative way of writing the function, not recommended.  liam committed Jan 14, 2008 405 406 (defun rosenbrock (argument return) "Rosenbrock test function."  liam committed Jan 16, 2008 407 408  (with-c-doubles (((vector-data argument) x0 x1) ((vector-data return) f0 f1))  liam committed Jan 14, 2008 409 410  (setf f0 (* *rosenbrock-a* (- 1 x0)) f1 (* *rosenbrock-b* (- x1 (expt x0 2))))))  liam committed Jan 17, 2008 411 |#  liam committed Jan 14, 2008 412   liam committed Jan 16, 2008 413 414 415 ;;; The recommended alternative (defun rosenbrock (argument return) "Rosenbrock test function."  Liam Healy committed Aug 31, 2008 416 417 418 419  (setf (maref return 0) (* *rosenbrock-a* (- 1 (maref argument 0))) (maref return 1) (* *rosenbrock-b* (- (maref argument 1) (expt (maref argument 0) 2)))))  liam committed Jan 16, 2008 420   liam committed Jan 14, 2008 421 422 423 424 425 (def-mfunction rosenbrock 2) (defun roots-multi-example () "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual." (let ((max-iter 1000))  Liam Healy committed Aug 31, 2008 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447  (letm ((vect (vector-double-float (a -10.0d0 -5.0d0))) (solver (mfsolver *hybrid-scaled* rosenbrock vect))) (loop for iter from 0 with fnval and argval while (and (< iter max-iter) (not (multiroot-test-residual solver 1.0d-7))) do (iterate-mfsolver solver) (setf fnval (cl-array (mfsolver-f solver)) argval (cl-array (mfsolver-root solver))) (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g" iter (aref argval 0) (aref argval 1) (aref fnval 0) (aref fnval 1)) finally (return (values (aref argval 0) (aref argval 1) (aref fnval 0) (aref fnval 1)))))))  liam committed Jan 28, 2008 448   liam committed Jan 16, 2008 449 450 (defun rosenbrock-df (argument jacobian) "The partial derivatives of the Rosenbrock functions."  Liam Healy committed Aug 31, 2008 451 452 453 454  (setf (maref jacobian 0 0) (- *rosenbrock-a*) (maref jacobian 0 1) 0.0d0 (maref jacobian 1 0) (* -2 *rosenbrock-b* (maref argument 0)) (maref jacobian 1 1) *rosenbrock-b*))  liam committed Jan 16, 2008 455 456 457 458 459  (defun rosenbrock-fdf (argument value jacobian) (rosenbrock argument value) (rosenbrock-df argument jacobian))  liam committed Jan 21, 2008 460 ;;; Because def-solver-functions and def-single-function bind a symbol  liam committed Jan 16, 2008 461 462 463 464 465 466 467 468 469 470 471 472 473 ;;; of the same name as the first function, and we want both to run, ;;; we'll make an alias function so we can use both. (eval-when (:load-toplevel :execute) (setf (fdefinition 'rosenbrock-f) #'rosenbrock)) (def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2) (defun roots-multi-example-df () "Solving Rosenbrock with derivatives, the example given in Sec. 34.8 of the GSL manual." (flet ((print-state (iter argval fnval) (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g" iter  liam committed Feb 24, 2008 474 475 476 477  (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1))))  liam committed Jan 16, 2008 478  (let ((max-iter 1000))  Liam Healy committed Aug 31, 2008 479  (letm ((vect (vector-double-float (a -10.0d0 -5.0d0))))  liam committed Feb 03, 2008 480 481  (letm ((solver (mfdfsolver *gnewton-mfdfsolver* rosenbrock-f vect)))  Liam Healy committed Aug 31, 2008 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497  (loop for iter from 0 with fnval = (mfdfsolver-f solver) and argval = (mfdfsolver-root solver) while (and (< iter max-iter) (not (multiroot-test-residual solver 1.0d-7))) initially (print-state iter argval fnval) do (iterate-mfdfsolver solver) (setf fnval (mfdfsolver-f solver) argval (mfdfsolver-root solver)) (print-state iter argval fnval) finally (return (values (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1)))))))))`