roots-multi.lisp 20 KB
 liam committed Jan 14, 2008 1 2 ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08 Liam M. Healy committed Jan 13, 2012 3 ;;; Time-stamp: <2012-01-13 12:01:17EST roots-multi.lisp> Liam Healy committed Dec 27, 2009 4 ;; Liam M. Healy committed Jan 10, 2011 5 ;; Copyright 2008, 2009, 2011 Liam M. Healy Liam Healy committed Dec 27, 2009 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ;; Distributed under the terms of the GNU General Public License ;; ;; This program is free software: you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by ;; the Free Software Foundation, either version 3 of the License, or ;; (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public License ;; along with this program. If not, see . liam committed Jan 14, 2008 20 21 22 (in-package :gsl) Liam Healy committed Dec 27, 2008 23 24 ;;; /usr/include/gsl/gsl_multiroots.h Liam Healy committed Jan 19, 2009 25 26 ;;; Currently, functions defined for root solving will be passed ;;; scalars and should return scalars as multiple values. A possible Liam Healy committed Jun 30, 2010 27 28 29 ;;; future enhancement is to optionally pass grid:foreign-arrays and return ;;; grid:foreign-arrays instead. This would allow directly manipulation of ;;; grid:foreign-arrays by the user function. Notes Mon Jan 19 2009. Liam Healy committed Jan 19, 2009 30 liam committed Jan 14, 2008 31 32 33 34 ;;;;**************************************************************************** ;;;; Initialization ;;;;**************************************************************************** Liam Healy committed Dec 25, 2008 35 (defmobject multi-dimensional-root-solver-f "gsl_multiroot_fsolver" Liam M. Healy committed Jan 13, 2012 36 ((type :pointer) ((first dimensions) :sizet)) Liam Healy committed Jan 25, 2009 37 38 "multi-dimensional root solver with function only" :documentation ; FDL Liam Healy committed Dec 25, 2008 39 40 41 "Make an instance of a solver of the type specified for a system of the specified number of dimensions. Optionally set or reset an existing solver to use the function and the Liam Healy committed Mar 20, 2009 42 initial guess gsl-vector. If scalarsp is T, the functions will Liam Healy committed Jun 10, 2009 43 be supplied scalars, and should return scalars." Liam Healy committed Jan 25, 2009 44 :initialize-suffix "set" Liam Healy committed Feb 04, 2009 45 :initialize-args ((callback :pointer) ((mpointer initial) :pointer)) Liam Healy committed Mar 19, 2009 46 :callbacks Liam Healy committed Oct 30, 2011 47 (callback (:struct fnstruct-dimension) (dimension) Liam Healy committed Mar 29, 2009 48 49 (function :success-failure Liam Healy committed Jun 30, 2010 50 51 (:input :double :foreign-array dim0) :slug (:output :double :foreign-array dim0))) Liam Healy committed Jan 25, 2009 52 :arglists-function Liam Healy committed Dec 25, 2008 53 (lambda (set) Liam Healy committed Mar 19, 2009 54 ((type &optional function-or-dimension (initial nil ,set) (scalarsp t)) Liam Healy committed Jan 25, 2009 55 (:type type Liam Healy committed Feb 07, 2009 56 :dimensions Liam M. Healy committed Jan 10, 2011 57 (if ,set (grid:dimensions initial) function-or-dimension)) Liam Healy committed Mar 21, 2009 58 59 (:functions (list function-or-dimension) :initial initial :scalarsp scalarsp))) Liam Healy committed Jan 25, 2009 60 :inputs (initial)) Liam Healy committed Dec 25, 2008 61 62 (defmobject multi-dimensional-root-solver-fdf "gsl_multiroot_fdfsolver" Liam M. Healy committed Jan 13, 2012 63 ((type :pointer) ((first dimensions) :sizet)) Liam Healy committed Jan 25, 2009 64 65 "multi-dimensional root solver with function and derivative" :documentation ; FDL Liam Healy committed Dec 25, 2008 66 67 68 "Make an instance of a derivative solver of the type specified for a system of the specified number of dimensions. Optionally set or reset an existing solver to use the function and derivative Liam Healy committed Mar 20, 2009 69 70 (fdf) and the initial guess. If scalarsp is T, the functions will be supplied, and should return scalars." Liam Healy committed Jan 25, 2009 71 :initialize-suffix "set" Liam Healy committed Feb 07, 2009 72 :initialize-args ((callback :pointer) ((mpointer initial) :pointer)) Liam Healy committed Mar 19, 2009 73 :callbacks Liam Healy committed Oct 30, 2011 74 (callback (:struct fnstruct-dimension-fdf) (dimension) Liam Healy committed Mar 19, 2009 75 (function :success-failure Liam Healy committed Jun 30, 2010 76 (:input :double :foreign-array dim0) Liam Healy committed Mar 29, 2009 77 :slug Liam Healy committed Jun 30, 2010 78 (:output :double :foreign-array dim0)) Liam Healy committed Mar 19, 2009 79 (df :success-failure Liam Healy committed Jun 30, 2010 80 (:input :double :foreign-array dim0) Liam Healy committed Mar 29, 2009 81 :slug Liam Healy committed Jun 30, 2010 82 (:output :double :foreign-array dim0 dim0)) Liam Healy committed Mar 20, 2009 83 (fdf :success-failure Liam Healy committed Jun 30, 2010 84 (:input :double :foreign-array dim0) Liam Healy committed Mar 29, 2009 85 :slug Liam Healy committed Jun 30, 2010 86 87 (:output :double :foreign-array dim0) (:output :double :foreign-array dim0 dim0))) Liam Healy committed Jan 25, 2009 88 :arglists-function Liam Healy committed Dec 25, 2008 89 (lambda (set) Liam Healy committed Mar 20, 2009 90 91 ((type &optional function-or-dimension (initial nil ,set) (scalarsp t)) Liam Healy committed Feb 07, 2009 92 93 (:type type :dimensions Liam M. Healy committed Jan 10, 2011 94 (if ,set (grid:dimensions initial) function-or-dimension)) Liam Healy committed Mar 20, 2009 95 (:functions function-or-dimension :initial initial :scalarsp scalarsp))) Liam Healy committed Jan 25, 2009 96 :inputs (initial)) liam committed Jan 14, 2008 97 Liam Healy committed Dec 26, 2008 98 (defmfun name ((solver multi-dimensional-root-solver-f)) liam committed Jan 14, 2008 99 "gsl_multiroot_fsolver_name" Liam Healy committed Dec 26, 2008 100 101 (((mpointer solver) :pointer)) :definition :method liam committed Jan 14, 2008 102 :c-return :string liam committed Feb 04, 2008 103 :documentation ; FDL liam committed Jan 14, 2008 104 105 "The name of the solver.") Liam Healy committed Dec 26, 2008 106 (defmfun name ((solver multi-dimensional-root-solver-fdf)) liam committed Jan 14, 2008 107 "gsl_multiroot_fdfsolver_name" Liam Healy committed Dec 26, 2008 108 109 (((mpointer solver) :pointer)) :definition :method liam committed Jan 14, 2008 110 :c-return :string liam committed Feb 04, 2008 111 :documentation ; FDL liam committed Jan 14, 2008 112 113 114 115 116 117 "The name of the solver.") ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** Liam Healy committed Jan 03, 2009 118 (defmfun iterate ((solver multi-dimensional-root-solver-f)) liam committed Jan 14, 2008 119 "gsl_multiroot_fsolver_iterate" Liam Healy committed Dec 26, 2008 120 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 121 :definition :method Liam Healy committed Mar 19, 2009 122 :callback-object solver liam committed Feb 04, 2008 123 :documentation ; FDL Liam Healy committed Nov 16, 2008 124 125 126 127 128 129 "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a singular point where the function or its derivative evaluated to infinity or NaN, or 'gsl-division-by-zero, the derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") liam committed Jan 14, 2008 130 Liam Healy committed Jan 03, 2009 131 (defmfun iterate ((solver multi-dimensional-root-solver-fdf)) liam committed Jan 14, 2008 132 "gsl_multiroot_fdfsolver_iterate" Liam Healy committed Dec 26, 2008 133 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 134 :definition :method Liam Healy committed Mar 19, 2009 135 :callback-object solver liam committed Feb 04, 2008 136 :documentation ; FDL Liam Healy committed Nov 16, 2008 137 138 139 140 141 142 "Perform a single iteration of the solver. The following errors may be signalled: 'bad-function-supplied, the iteration encountered a singular point where the function or its derivative evaluated to infinity or NaN, or 'gsl-division-by-zero, the derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.") liam committed Jan 14, 2008 143 Liam Healy committed Jan 03, 2009 144 (defmfun solution ((solver multi-dimensional-root-solver-f)) liam committed Jan 14, 2008 145 "gsl_multiroot_fsolver_root" Liam Healy committed Dec 26, 2008 146 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 147 :definition :method Liam Healy committed Mar 19, 2009 148 :callback-object solver Liam Healy committed Dec 29, 2008 149 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 150 :return ((make-foreign-array-from-mpointer crtn)) liam committed Feb 04, 2008 151 :documentation ; FDL liam committed Jan 14, 2008 152 153 "The current estimate of the root for the solver.") Liam Healy committed Jan 03, 2009 154 (defmfun solution ((solver multi-dimensional-root-solver-fdf)) liam committed Jan 14, 2008 155 "gsl_multiroot_fdfsolver_root" Liam Healy committed Dec 26, 2008 156 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 157 :definition :method Liam Healy committed Mar 19, 2009 158 :callback-object solver Liam Healy committed Dec 29, 2008 159 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 160 :return ((make-foreign-array-from-mpointer crtn)) liam committed Jan 14, 2008 161 162 163 :documentation "The current estimate of the root for the solver.") Liam Healy committed Jan 03, 2009 164 (defmfun function-value ((solver multi-dimensional-root-solver-f)) liam committed Jan 14, 2008 165 "gsl_multiroot_fsolver_f" Liam Healy committed Dec 26, 2008 166 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 167 :definition :method Liam Healy committed Dec 29, 2008 168 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 169 :return ((make-foreign-array-from-mpointer crtn)) liam committed Feb 04, 2008 170 :documentation ; FDL liam committed Jan 14, 2008 171 172 "The function value f(x) at the current estimate x of the root for the solver.") Liam Healy committed Jan 03, 2009 173 (defmfun function-value ((solver multi-dimensional-root-solver-fdf)) liam committed Jan 14, 2008 174 "gsl_multiroot_fdfsolver_f" Liam Healy committed Dec 26, 2008 175 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 176 :definition :method Liam Healy committed Dec 29, 2008 177 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 178 :return ((make-foreign-array-from-mpointer crtn)) liam committed Feb 04, 2008 179 :documentation ; FDL liam committed Jan 14, 2008 180 181 "The function value f(x) at the current estimate x of the root for the solver.") Liam Healy committed Jan 03, 2009 182 (defmfun last-step ((solver multi-dimensional-root-solver-f)) liam committed Jan 14, 2008 183 "gsl_multiroot_fsolver_dx" Liam Healy committed Dec 26, 2008 184 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 185 :definition :method Liam Healy committed Dec 29, 2008 186 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 187 :return ((make-foreign-array-from-mpointer crtn)) liam committed Feb 04, 2008 188 :documentation ; FDL liam committed Jan 14, 2008 189 190 "The last step dx taken by the solver.") Liam Healy committed Jan 03, 2009 191 (defmfun last-step ((solver multi-dimensional-root-solver-fdf)) liam committed Jan 14, 2008 192 "gsl_multiroot_fsolver_dx" Liam Healy committed Dec 26, 2008 193 (((mpointer solver) :pointer)) Liam Healy committed Jan 03, 2009 194 :definition :method Liam Healy committed Dec 29, 2008 195 :c-return (crtn :pointer) Liam Healy committed Jun 27, 2010 196 :return ((make-foreign-array-from-mpointer crtn)) liam committed Feb 04, 2008 197 :documentation ; FDL liam committed Jan 14, 2008 198 199 200 201 202 203 204 "The last step dx taken by the solver.") ;;;;**************************************************************************** ;;;; Search stopping conditions ;;;;**************************************************************************** ;;; The only place we need to pick apart the gsl_multiroot_fsolver Liam Healy committed Jan 03, 2009 205 ;;; struct is here. We could use last-step etc., but then we'd have liam committed Jan 14, 2008 206 207 208 ;;; to discriminate on mfsolver vs. mfdfsolver. (defun multiroot-slot (solver slot) Liam Healy committed Oct 30, 2011 209 (cffi:foreign-slot-value (mpointer solver) '(:struct gsl-multiroot-fsolver) slot)) liam committed Jan 14, 2008 210 liam committed Feb 18, 2008 211 (defmfun multiroot-test-delta (solver absolute-error relative-error) liam committed Jan 14, 2008 212 213 214 215 216 "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 217 :documentation ; FDL liam committed Jan 14, 2008 218 219 220 221 222 223 224 "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 225 (defmfun multiroot-test-residual (solver absolute-error) liam committed Jan 14, 2008 226 227 "gsl_multiroot_test_residual" (((multiroot-slot solver 'f) :pointer) (absolute-error :double)) Liam Healy committed Mar 23, 2009 228 :c-return :success-continue liam committed Feb 04, 2008 229 :documentation ; FDL liam committed Jan 14, 2008 230 231 232 233 234 235 236 237 238 239 240 241 "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 Healy committed Feb 16, 2009 242 (defmpar +powells-hybrid+ "gsl_multiroot_fdfsolver_hybridsj" liam committed Feb 04, 2008 243 ;; FDL liam committed Jan 14, 2008 244 "This is a modified version of Powell's Hybrid method as implemented in liam committed Feb 18, 2008 245 the hybrj algorithm in minpack. Minpack was written by Jorge liam committed Jan 14, 2008 246 247 248 249 250 251 252 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 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, Liam Healy committed Nov 16, 2008 281 'no-progress liam committed Jan 14, 2008 282 283 the iteration is not making any progress, preventing the algorithm from continuing. Liam Healy committed Nov 16, 2008 284 'jacobian-not-improving liam committed Jan 14, 2008 285 286 287 re-evaluations of the Jacobian indicate that the iteration is not making any progress, preventing the algorithm from continuing.") Liam Healy committed Feb 16, 2009 288 (defmpar +powells-hybrid-unscaled+ "gsl_multiroot_fdfsolver_hybridj" liam committed Feb 04, 2008 289 ;; FDL liam committed Jan 14, 2008 290 291 292 293 294 "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 Healy committed Feb 16, 2009 295 (defmpar +newton-mfdfsolver+ "gsl_multiroot_fdfsolver_newton" liam committed Feb 04, 2008 296 ;; FDL liam committed Jan 14, 2008 297 298 299 300 301 302 303 304 305 306 307 308 "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 Healy committed Feb 16, 2009 309 (defmpar +gnewton-mfdfsolver+ "gsl_multiroot_fdfsolver_gnewton" liam committed Feb 04, 2008 310 ;; FDL liam committed Jan 14, 2008 311 312 313 314 315 316 317 318 319 320 321 322 323 "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 Healy committed Feb 16, 2009 324 (defmpar +hybrid-scaled+ "gsl_multiroot_fsolver_hybrids" liam committed Feb 04, 2008 325 326 ;; FDL "This is a version of the Hybrid algorithm which replaces calls to the liam committed Jan 14, 2008 327 328 329 330 331 332 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 Healy committed Feb 16, 2009 333 (defmpar +hybrid-unscaled+ "gsl_multiroot_fsolver_hybrid" liam committed Feb 04, 2008 334 ;; FDL liam committed Jan 14, 2008 335 336 337 "A finite difference version of the Hybrid algorithm without internal scaling.") Liam Healy committed Feb 16, 2009 338 (defmpar +discrete-newton+ "gsl_multiroot_fsolver_dnewton" liam committed Feb 04, 2008 339 ;; FDL liam committed Jan 14, 2008 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 "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 Healy committed Feb 16, 2009 355 (defmpar +broyden+ "gsl_multiroot_fsolver_broyden" liam committed Feb 04, 2008 356 ;; FDL liam committed Jan 14, 2008 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 "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 ;;;;**************************************************************************** Liam M. Healy committed Jul 13, 2010 379 380 381 382 383 384 385 ;;; This is the example given in GSL manual, Sec. 35.8. ;;; http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Multidimensional-Root-finding.html ;;; These examples use use scalarsp=T in the ;;; multi-dimensional-root-solver argument. To see how vectors would ;;; be used, see minimization-multi. liam committed Jan 14, 2008 386 (defparameter *powell-A* 1.0d4) Liam Healy committed Jan 19, 2009 387 (defun powell (arg0 arg1) liam committed Jan 14, 2008 388 "Powell's test function." Liam Healy committed Jan 19, 2009 389 390 391 (values (- (* *powell-A* arg0 arg1) 1) (+ (exp (- arg0)) (exp (- arg1)) (- (1+ (/ *powell-A*)))))) Liam Healy committed Jan 24, 2009 392 ;; not used? liam committed Jan 14, 2008 393 394 395 396 (defparameter *rosenbrock-a* 1.0d0) (defparameter *rosenbrock-b* 10.0d0) Liam Healy committed Feb 07, 2009 397 (defun rosenbrock (arg0 arg1) liam committed Jan 16, 2008 398 "Rosenbrock test function." Liam Healy committed Jan 19, 2009 399 400 401 (values (* *rosenbrock-a* (- 1 arg0)) (* *rosenbrock-b* (- arg1 (expt arg0 2))))) liam committed Jan 16, 2008 402 Liam Healy committed Feb 04, 2009 403 (defun roots-multi-example-no-derivative Liam Healy committed Feb 16, 2009 404 (&optional (method +hybrid-scaled+) (print-steps t)) Liam Healy committed Feb 04, 2009 405 406 407 "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual." (let ((max-iter 1000) (solver (make-multi-dimensional-root-solver-f Liam Healy committed Aug 17, 2011 408 409 410 method 'rosenbrock (grid:make-foreign-array 'double-float :initial-contents '(-10.0d0 -5.0d0))))) Liam Healy committed Feb 04, 2009 411 (loop for iter from 0 Liam Healy committed Aug 17, 2011 412 413 414 415 416 with fnval and argval while (and (< iter max-iter) (or (zerop iter) (not (multiroot-test-residual solver 1.0d-7)))) do Liam Healy committed Feb 04, 2009 417 418 419 420 421 422 (iterate solver) (setf fnval (function-value solver) argval (solution solver)) (when print-steps (format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&" iter Liam M. Healy committed May 26, 2011 423 424 425 426 (grid:aref argval 0) (grid:aref argval 1) (grid:aref fnval 0) (grid:aref fnval 1))) Liam Healy committed Aug 17, 2011 427 428 429 430 431 finally (return (values (grid:aref argval 0) (grid:aref argval 1) (grid:aref fnval 0) (grid:aref fnval 1)))))) Liam Healy committed Feb 04, 2009 432 Liam Healy committed Feb 07, 2009 433 (defun rosenbrock-df (arg0 arg1) liam committed Jan 16, 2008 434 "The partial derivatives of the Rosenbrock functions." Liam Healy committed Jan 19, 2009 435 436 437 438 439 440 441 (declare (ignore arg1)) (values (- *rosenbrock-a*) 0.0d0 (* -2 *rosenbrock-b* arg0) *rosenbrock-b*)) ;;; Why is it necessary to define a function that calls the two other functions? Liam Healy committed Feb 07, 2009 442 (defun rosenbrock-fdf (arg0 arg1) Liam Healy committed Jan 19, 2009 443 444 445 446 447 (multiple-value-bind (v0 v1) (rosenbrock arg0 arg1) (multiple-value-bind (j0 j1 j2 j3) (rosenbrock-df arg0 arg1) (values v0 v1 j0 j1 j2 j3)))) liam committed Jan 16, 2008 448 Liam Healy committed Feb 02, 2009 449 (defun roots-multi-example-derivative Liam Healy committed Feb 16, 2009 450 (&optional (method +gnewton-mfdfsolver+) (print-steps t)) liam committed Jan 16, 2008 451 452 453 "Solving Rosenbrock with derivatives, the example given in Sec. 34.8 of the GSL manual." (flet ((print-state (iter argval fnval) Liam Healy committed Jan 19, 2009 454 455 456 (when print-steps (format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&" iter Liam M. Healy committed May 26, 2011 457 458 459 460 (grid:aref argval 0) (grid:aref argval 1) (grid:aref fnval 0) (grid:aref fnval 1))))) Liam Healy committed Jan 23, 2009 461 462 (let ((max-iter 1000) (solver (make-multi-dimensional-root-solver-fdf Liam Healy committed Jan 24, 2009 463 method Liam Healy committed Feb 07, 2009 464 '(rosenbrock rosenbrock-df rosenbrock-fdf) Liam Healy committed Aug 17, 2011 465 466 (grid:make-foreign-array 'double-float :initial-contents '(-10.0d0 -5.0d0))))) Liam Healy committed Jan 23, 2009 467 (loop for iter from 0 Liam Healy committed Aug 17, 2011 468 469 470 471 472 473 with fnval = (function-value solver) and argval = (solution solver) while (and (< iter max-iter) (not (multiroot-test-residual solver 1.0d-7))) initially (print-state iter argval fnval) do Liam Healy committed Jan 23, 2009 474 475 476 477 (iterate solver) (setf fnval (function-value solver) argval (solution solver)) (print-state iter argval fnval) Liam Healy committed Aug 17, 2011 478 479 480 481 482 finally (return (values (grid:aref argval 0) (grid:aref argval 1) (grid:aref fnval 0) (grid:aref fnval 1))))))) Liam Healy committed Jan 19, 2009 483 Liam Healy committed Jan 23, 2009 484 485 ;; To see step-by-step information as the solution progresses, make ;; the last argument T. Liam Healy committed Jan 19, 2009 486 (save-test roots-multi Liam Healy committed Feb 16, 2009 487 488 489 490 491 492 493 494 (roots-multi-example-no-derivative +hybrid-unscaled+ nil) (roots-multi-example-no-derivative +hybrid-scaled+ nil) (roots-multi-example-no-derivative +discrete-newton+ nil) (roots-multi-example-no-derivative +broyden+ nil) (roots-multi-example-derivative +newton-mfdfsolver+ nil) (roots-multi-example-derivative +gnewton-mfdfsolver+ nil) (roots-multi-example-derivative +powells-hybrid+ nil) (roots-multi-example-derivative +powells-hybrid-unscaled+ nil))