roots-multi.lisp 18.5 KB
 liam committed Jan 14, 2008 1 2 ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08 lhealy committed Jul 24, 2008 3 ;;; Time-stamp: <2008-03-09 19:29:20EDT roots-multi.lisp> lhealy committed Jul 24, 2008 4 ;;; $Id$ liam committed Jan 14, 2008 5 6 7 (in-package :gsl) liam committed Jan 16, 2008 8 9 10 11 12 ;;; I don't like using make-data-from-pointer, but it's the only way ;;; to have access to the GSL functions when given a pointer. ;;; Alternatively, I could provide the GSL pointer and then the only ;;; thing the user could use is vref, or, of course ;;; make-data-from-pointer. liam committed Jan 14, 2008 13 14 15 16 17 18 19 20 21 ;;;;**************************************************************************** ;;;; 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 committed Feb 18, 2008 22 (dimensions size) liam committed Jan 14, 2008 23 24 25 26 27 (parameters :pointer)) (export 'def-mfunction) (defmacro def-mfunction (name dimensions) "Define a function for multivariate root solving." liam committed Jan 21, 2008 28 (def-single-function ,name :success-failure :pointer gsl-mfunction liam committed Jan 14, 2008 29 ((dimensions ,dimensions)) liam committed Jan 21, 2008 30 (gsl-vector-c))) liam committed Jan 14, 2008 31 liam committed Jan 16, 2008 32 33 34 35 36 37 38 (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 committed Feb 18, 2008 39 (dimensions size) liam committed Jan 16, 2008 40 41 (parameters :pointer)) liam committed Jan 14, 2008 42 43 44 45 ;;;;**************************************************************************** ;;;; Initialization ;;;;**************************************************************************** liam committed Feb 03, 2008 46 (defgo mfsolver (type function-derivative initial) liam committed Feb 03, 2008 47 48 49 50 51 52 (list (allocate-mfsolver ,type (dim0 ,initial)) 'free-mfsolver (lambda (symb) (set-mfsolver ,symb ,function-derivative ,initial)))) liam committed Feb 03, 2008 53 (defgo mfdfsolver (type function-derivative initial) liam committed Feb 03, 2008 54 55 56 57 58 (list (allocate-mfdfsolver ,type (dim0 ,initial)) 'free-mfdfsolver (lambda (symb) `(set-mfdfsolver ,symb ,function-derivative ,initial)))) liam committed Jan 28, 2008 59 liam committed Feb 18, 2008 60 (defmfun allocate-mfsolver (type dimension) liam committed Jan 14, 2008 61 "gsl_multiroot_fsolver_alloc" liam committed Feb 18, 2008 62 ((type :pointer) (dimension size)) liam committed Jan 14, 2008 63 :c-return :pointer liam committed Jan 28, 2008 64 :export nil liam committed Feb 03, 2008 65 :index (letm mfsolver) liam committed Feb 04, 2008 66 :documentation ; FDL liam committed Jan 14, 2008 67 68 69 "Allocate an instance of a solver of the type specified for a system of the specified number of dimensions.") liam committed Feb 18, 2008 70 (defmfun allocate-mfdfsolver (type dimension) liam committed Jan 14, 2008 71 "gsl_multiroot_fdfsolver_alloc" liam committed Feb 18, 2008 72 ((type :pointer) (dimension size)) liam committed Jan 14, 2008 73 :c-return :pointer liam committed Jan 28, 2008 74 :export nil liam committed Feb 03, 2008 75 :index (letm mfdfsolver) liam committed Feb 04, 2008 76 :documentation ; FDL liam committed Jan 14, 2008 77 78 79 "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 80 (defmfun set-mfsolver (solver function initial) liam committed Jan 14, 2008 81 82 "gsl_multiroot_fsolver_set" ((solver :pointer) (function :pointer) ((pointer initial) :pointer)) liam committed Jan 28, 2008 83 :export nil liam committed Feb 03, 2008 84 :index (letm mfsolver) liam committed Feb 04, 2008 85 :documentation ; FDL liam committed Jan 14, 2008 86 87 88 "Set or reset an existing solver to use the function and the initial guess gsl-vector.") liam committed Feb 18, 2008 89 (defmfun set-mfdfsolver (solver function-derivative initial) liam committed Jan 14, 2008 90 91 "gsl_multiroot_fdfsolver_set" ((solver :pointer) (function-derivative :pointer) liam committed Jan 16, 2008 92 ((pointer initial) :pointer)) liam committed Jan 28, 2008 93 :export nil liam committed Feb 03, 2008 94 :index (letm mfdfsolver) liam committed Feb 04, 2008 95 :documentation ; FDL liam committed Jan 14, 2008 96 97 98 "Set or reset an existing solver to use the function and derivative (fdf) and the initial guess.") liam committed Feb 18, 2008 99 (defmfun free-mfsolver (solver) liam committed Jan 14, 2008 100 101 102 "gsl_multiroot_fsolver_free" ((solver :pointer)) :c-return :void liam committed Jan 28, 2008 103 :export nil liam committed Feb 03, 2008 104 :index (letm mfsolver) liam committed Feb 04, 2008 105 :documentation ; FDL liam committed Jan 14, 2008 106 107 "Free all the memory associated with the solver.") liam committed Feb 18, 2008 108 (defmfun free-mfdfsolver (solver) liam committed Jan 14, 2008 109 110 111 "gsl_multiroot_fdfsolver_free" ((solver :pointer)) :c-return :void liam committed Jan 28, 2008 112 :export nil liam committed Feb 03, 2008 113 :index (letm mfdfsolver) liam committed Feb 04, 2008 114 :documentation ; FDL liam committed Jan 14, 2008 115 116 "Free all the memory associated with the solver.") liam committed Feb 18, 2008 117 (defmfun mfsolver-name (solver) liam committed Jan 14, 2008 118 119 120 "gsl_multiroot_fsolver_name" ((solver :pointer)) :c-return :string liam committed Feb 04, 2008 121 :documentation ; FDL liam committed Jan 14, 2008 122 123 "The name of the solver.") liam committed Feb 18, 2008 124 (defmfun mfdfsolver-name (solver) liam committed Jan 14, 2008 125 126 127 "gsl_multiroot_fdfsolver_name" ((solver :pointer)) :c-return :string liam committed Feb 04, 2008 128 :documentation ; FDL liam committed Jan 14, 2008 129 130 131 132 133 134 "The name of the solver.") ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** liam committed Feb 18, 2008 135 (defmfun iterate-mfsolver (solver) liam committed Jan 14, 2008 136 137 "gsl_multiroot_fsolver_iterate" ((solver :pointer)) liam committed Feb 04, 2008 138 :documentation ; FDL liam committed Jan 14, 2008 139 140 141 142 143 144 145 "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 146 (defmfun iterate-mfdfsolver (solver) liam committed Jan 14, 2008 147 148 "gsl_multiroot_fdfsolver_iterate" ((solver :pointer)) liam committed Feb 04, 2008 149 :documentation ; FDL liam committed Jan 14, 2008 150 151 152 153 154 155 156 "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 157 (defmfun mfsolver-root (solver) liam committed Jan 14, 2008 158 159 160 161 "gsl_multiroot_fsolver_root" ((solver :pointer)) :c-return (canswer :pointer) :return ((make-data-from-pointer canswer)) liam committed Feb 04, 2008 162 :documentation ; FDL liam committed Jan 14, 2008 163 164 "The current estimate of the root for the solver.") liam committed Feb 18, 2008 165 (defmfun mfdfsolver-root (solver) liam committed Jan 14, 2008 166 167 168 169 170 171 172 "gsl_multiroot_fdfsolver_root" ((solver :pointer)) :c-return (canswer gsl-vector-c) :return ((make-data-from-pointer canswer)) :documentation "The current estimate of the root for the solver.") liam committed Feb 18, 2008 173 (defmfun mfsolver-f (solver) liam committed Jan 14, 2008 174 175 176 177 "gsl_multiroot_fsolver_f" ((solver :pointer)) :c-return (canswer gsl-vector-c) :return ((make-data-from-pointer canswer)) liam committed Feb 04, 2008 178 :documentation ; FDL liam committed Jan 14, 2008 179 180 "The function value f(x) at the current estimate x of the root for the solver.") liam committed Feb 18, 2008 181 (defmfun mfdfsolver-f (solver) liam committed Jan 14, 2008 182 183 184 185 "gsl_multiroot_fdfsolver_f" ((solver :pointer)) :c-return (canswer gsl-vector-c) :return ((make-data-from-pointer canswer)) liam committed Feb 04, 2008 186 :documentation ; FDL liam committed Jan 14, 2008 187 188 "The function value f(x) at the current estimate x of the root for the solver.") liam committed Feb 18, 2008 189 (defmfun mfsolver-dx (solver) liam committed Jan 14, 2008 190 191 192 193 "gsl_multiroot_fsolver_dx" ((solver :pointer)) :c-return (canswer gsl-vector-c) :return ((make-data-from-pointer canswer)) liam committed Feb 04, 2008 194 :documentation ; FDL liam committed Jan 14, 2008 195 196 "The last step dx taken by the solver.") liam committed Feb 18, 2008 197 (defmfun mfdfsolver-dx (solver) liam committed Jan 14, 2008 198 199 200 201 "gsl_multiroot_fsolver_dx" ((solver :pointer)) :c-return (canswer gsl-vector-c) :return ((make-data-from-pointer canswer)) liam committed Feb 04, 2008 202 :documentation ; FDL liam committed Jan 14, 2008 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 "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 224 (defmfun multiroot-test-delta (solver absolute-error relative-error) liam committed Jan 14, 2008 225 226 227 228 229 "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 230 :documentation ; FDL liam committed Jan 14, 2008 231 232 233 234 235 236 237 "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 238 (defmfun multiroot-test-residual (solver absolute-error) liam committed Jan 14, 2008 239 240 241 "gsl_multiroot_test_residual" (((multiroot-slot solver 'f) :pointer) (absolute-error :double)) :c-return :success-failure liam committed Feb 04, 2008 242 :documentation ; FDL liam committed Jan 14, 2008 243 244 245 246 247 248 249 250 251 252 253 254 "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 255 (defmpar *powells-hybrid* "gsl_multiroot_fdfsolver_hybridsj" liam committed Feb 04, 2008 256 ;; FDL liam committed Jan 14, 2008 257 "This is a modified version of Powell's Hybrid method as implemented in liam committed Feb 18, 2008 258 the hybrj algorithm in minpack. Minpack was written by Jorge liam committed Jan 14, 2008 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 295 296 297 298 299 300 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 301 (defmpar *powells-hybrid-unscaled* "gsl_multiroot_fdfsolver_hybridj" liam committed Feb 04, 2008 302 ;; FDL liam committed Jan 14, 2008 303 304 305 306 307 "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 308 (defmpar *newton-mfdfsolver* "gsl_multiroot_fdfsolver_newton" liam committed Feb 04, 2008 309 ;; FDL liam committed Jan 14, 2008 310 311 312 313 314 315 316 317 318 319 320 321 "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 322 (defmpar *gnewton-mfdfsolver* "gsl_multiroot_fdfsolver_gnewton" liam committed Feb 04, 2008 323 ;; FDL liam committed Jan 14, 2008 324 325 326 327 328 329 330 331 332 333 334 335 336 "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 337 (defmpar *hybrid-scaled* "gsl_multiroot_fsolver_hybrids" liam committed Feb 04, 2008 338 339 ;; FDL "This is a version of the Hybrid algorithm which replaces calls to the liam committed Jan 14, 2008 340 341 342 343 344 345 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 346 (defmpar *hybrid-unscaled* "gsl_multiroot_fsolver_hybrid" liam committed Feb 04, 2008 347 ;; FDL liam committed Jan 14, 2008 348 349 350 "A finite difference version of the Hybrid algorithm without internal scaling.") liam committed Feb 18, 2008 351 (defmpar *discrete-newton* "gsl_multiroot_fsolver_dnewton" liam committed Feb 04, 2008 352 ;; FDL liam committed Jan 14, 2008 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 "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 368 (defmpar *broyden* "gsl_multiroot_fsolver_broyden" liam committed Feb 04, 2008 369 ;; FDL liam committed Jan 14, 2008 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 "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 395 396 (setf (maref return 0) (- (* *powell-A* (maref argument 0) (maref argument 1)) liam committed Jan 14, 2008 397 1) liam committed Feb 24, 2008 398 399 (maref return 1) (+ (exp (- (maref argument 0))) (exp (- (maref argument 1))) liam committed Jan 14, 2008 400 401 402 403 404 405 406 407 408 (- (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 409 #| liam committed Jan 16, 2008 410 ;;; One alternative way of writing the function, not recommended. liam committed Jan 14, 2008 411 412 (defun rosenbrock (argument return) "Rosenbrock test function." liam committed Jan 16, 2008 413 414 (with-c-doubles (((vector-data argument) x0 x1) ((vector-data return) f0 f1)) liam committed Jan 14, 2008 415 416 (setf f0 (* *rosenbrock-a* (- 1 x0)) f1 (* *rosenbrock-b* (- x1 (expt x0 2)))))) liam committed Jan 17, 2008 417 |# liam committed Jan 14, 2008 418 liam committed Jan 16, 2008 419 420 421 422 423 424 425 426 ;;; The recommended alternative (defun rosenbrock (argument return) "Rosenbrock test function." (setf (vref return 0) (* *rosenbrock-a* (- 1 (vref argument 0))) (vref return 1) (* *rosenbrock-b* (- (vref argument 1) (expt (vref argument 0) 2))))) liam committed Jan 14, 2008 427 428 429 430 431 (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)) lhealy committed Jul 24, 2008 432 (letm ((vect (vector-double-float #(-10.0d0 -5.0d0)))) liam committed Feb 03, 2008 433 (letm ((solver (mfsolver *hybrid-scaled* rosenbrock vect))) liam committed Jan 14, 2008 434 435 436 437 438 439 440 441 442 (let ((fnval (mfsolver-f solver)) (argval (mfsolver-root solver))) (loop for iter from 0 while (and (< iter max-iter) (not (multiroot-test-residual solver 1.0d-7))) do (iterate-mfsolver solver) (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g" iter liam committed Feb 24, 2008 443 444 445 446 (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1)) liam committed Jan 14, 2008 447 finally (return liam committed Feb 24, 2008 448 449 450 451 (values (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1))))))))) liam committed Jan 16, 2008 452 liam committed Jan 28, 2008 453 liam committed Jan 16, 2008 454 455 456 457 458 459 460 461 462 463 464 (defun rosenbrock-df (argument jacobian) "The partial derivatives of the Rosenbrock functions." (setf (mref jacobian 0 0) (- *rosenbrock-a*) (mref jacobian 0 1) 0.0d0 (mref jacobian 1 0) (* -2 *rosenbrock-b* (vref argument 0)) (mref jacobian 1 1) *rosenbrock-b*)) (defun rosenbrock-fdf (argument value jacobian) (rosenbrock argument value) (rosenbrock-df argument jacobian)) liam committed Jan 21, 2008 465 ;;; Because def-solver-functions and def-single-function bind a symbol liam committed Jan 16, 2008 466 467 468 469 470 471 472 473 474 475 476 477 478 ;;; 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 479 480 481 482 (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1)))) liam committed Jan 16, 2008 483 (let ((max-iter 1000)) lhealy committed Jul 24, 2008 484 (letm ((vect (vector-double-float #(-10.0d0 -5.0d0)))) liam committed Feb 03, 2008 485 486 (letm ((solver (mfdfsolver *gnewton-mfdfsolver* rosenbrock-f vect))) liam committed Jan 16, 2008 487 488 489 490 491 492 493 494 495 496 (let ((fnval (mfdfsolver-f solver)) (argval (mfdfsolver-root solver))) (loop for iter from 0 while (and (< iter max-iter) (not (multiroot-test-residual solver 1.0d-7))) initially (print-state iter argval fnval) do (iterate-mfdfsolver solver) (print-state iter argval fnval) finally (return liam committed Feb 24, 2008 497 498 499 500 (values (maref argval 0) (maref argval 1) (maref fnval 0) (maref fnval 1))))))))))