This is gsl-ref.info, produced by makeinfo version 4.13 from gsl-ref.texi. INFO-DIR-SECTION Software libraries START-INFO-DIR-ENTRY * gsl-ref: (gsl-ref). GNU Scientific Library - Reference END-INFO-DIR-ENTRY Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 The GSL Team. Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with the Invariant Sections being "GNU General Public License" and "Free Software Needs Free Documentation", the Front-Cover text being "A GNU Manual", and with the Back-Cover Text being (a) (see below). A copy of the license is included in the section entitled "GNU Free Documentation License". (a) The Back-Cover Text is: "You have the freedom to copy and modify this GNU Manual."  File: gsl-ref.info, Node: Initializing the Multidimensional Solver, Next: Providing the multidimensional system of equations to solve, Prev: Overview of Multidimensional Root Finding, Up: Multidimensional Root-Finding 35.2 Initializing the Solver ============================ The following functions initialize a multidimensional solver, either with or without derivatives. The solver itself depends only on the dimension of the problem and the algorithm and can be reused for different problems. -- Function: gsl_multiroot_fsolver * gsl_multiroot_fsolver_alloc (const gsl_multiroot_fsolver_type * T, size_t N) This function returns a pointer to a newly allocated instance of a solver of type T for a system of N dimensions. For example, the following code creates an instance of a hybrid solver, to solve a 3-dimensional system of equations. const gsl_multiroot_fsolver_type * T = gsl_multiroot_fsolver_hybrid; gsl_multiroot_fsolver * s = gsl_multiroot_fsolver_alloc (T, 3); If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of `GSL_ENOMEM'. -- Function: gsl_multiroot_fdfsolver * gsl_multiroot_fdfsolver_alloc (const gsl_multiroot_fdfsolver_type * T, size_t N) This function returns a pointer to a newly allocated instance of a derivative solver of type T for a system of N dimensions. For example, the following code creates an instance of a Newton-Raphson solver, for a 2-dimensional system of equations. const gsl_multiroot_fdfsolver_type * T = gsl_multiroot_fdfsolver_newton; gsl_multiroot_fdfsolver * s = gsl_multiroot_fdfsolver_alloc (T, 2); If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of `GSL_ENOMEM'. -- Function: int gsl_multiroot_fsolver_set (gsl_multiroot_fsolver * S, gsl_multiroot_function * F, const gsl_vector * X) -- Function: int gsl_multiroot_fdfsolver_set (gsl_multiroot_fdfsolver * S, gsl_multiroot_function_fdf * FDF, const gsl_vector * X) These functions set, or reset, an existing solver S to use the function F or function and derivative FDF, and the initial guess X. Note that the initial position is copied from X, this argument is not modified by subsequent iterations. -- Function: void gsl_multiroot_fsolver_free (gsl_multiroot_fsolver * S) -- Function: void gsl_multiroot_fdfsolver_free (gsl_multiroot_fdfsolver * S) These functions free all the memory associated with the solver S. -- Function: const char * gsl_multiroot_fsolver_name (const gsl_multiroot_fsolver * S) -- Function: const char * gsl_multiroot_fdfsolver_name (const gsl_multiroot_fdfsolver * S) These functions return a pointer to the name of the solver. For example, printf ("s is a '%s' solver\n", gsl_multiroot_fdfsolver_name (s)); would print something like `s is a 'newton' solver'.  File: gsl-ref.info, Node: Providing the multidimensional system of equations to solve, Next: Iteration of the multidimensional solver, Prev: Initializing the Multidimensional Solver, Up: Multidimensional Root-Finding 35.3 Providing the function to solve ==================================== You must provide n functions of n variables for the root finders to operate on. In order to allow for general parameters the functions are defined by the following data types: -- Data Type: gsl_multiroot_function This data type defines a general system of functions with parameters. `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `size_t n' the dimension of the system, i.e. the number of components of the vectors X and F. `void * params' a pointer to the parameters of the function. Here is an example using Powell's test function, f_1(x) = A x_0 x_1 - 1, f_2(x) = exp(-x_0) + exp(-x_1) - (1 + 1/A) with A = 10^4. The following code defines a `gsl_multiroot_function' system `F' which you could pass to a solver: struct powell_params { double A; }; int powell (gsl_vector * x, void * p, gsl_vector * f) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); gsl_vector_set (f, 0, A * x0 * x1 - 1); gsl_vector_set (f, 1, (exp(-x0) + exp(-x1) - (1.0 + 1.0/A))); return GSL_SUCCESS } gsl_multiroot_function F; struct powell_params params = { 10000.0 }; F.f = &powell; F.n = 2; F.params = ¶ms; -- Data Type: gsl_multiroot_function_fdf This data type defines a general system of functions with parameters and the corresponding Jacobian matrix of derivatives, `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)' this function should store the N-by-N matrix result J_ij = d f_i(x,params) / d x_j in J for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)' This function should set the values of the F and J as above, for arguments X and parameters PARAMS. This function provides an optimization of the separate functions for f(x) and J(x)--it is always faster to compute the function and its derivative at the same time. `size_t n' the dimension of the system, i.e. the number of components of the vectors X and F. `void * params' a pointer to the parameters of the function. The example of Powell's test function defined above can be extended to include analytic derivatives using the following code, int powell_df (gsl_vector * x, void * p, gsl_matrix * J) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); gsl_matrix_set (J, 0, 0, A * x1); gsl_matrix_set (J, 0, 1, A * x0); gsl_matrix_set (J, 1, 0, -exp(-x0)); gsl_matrix_set (J, 1, 1, -exp(-x1)); return GSL_SUCCESS } int powell_fdf (gsl_vector * x, void * p, gsl_matrix * f, gsl_matrix * J) { struct powell_params * params = *(struct powell_params *)p; const double A = (params->A); const double x0 = gsl_vector_get(x,0); const double x1 = gsl_vector_get(x,1); const double u0 = exp(-x0); const double u1 = exp(-x1); gsl_vector_set (f, 0, A * x0 * x1 - 1); gsl_vector_set (f, 1, u0 + u1 - (1 + 1/A)); gsl_matrix_set (J, 0, 0, A * x1); gsl_matrix_set (J, 0, 1, A * x0); gsl_matrix_set (J, 1, 0, -u0); gsl_matrix_set (J, 1, 1, -u1); return GSL_SUCCESS } gsl_multiroot_function_fdf FDF; FDF.f = &powell_f; FDF.df = &powell_df; FDF.fdf = &powell_fdf; FDF.n = 2; FDF.params = 0; Note that the function `powell_fdf' is able to reuse existing terms from the function when calculating the Jacobian, thus saving time.  File: gsl-ref.info, Node: Iteration of the multidimensional solver, Next: Search Stopping Parameters for the multidimensional solver, Prev: Providing the multidimensional system of equations to solve, Up: Multidimensional Root-Finding 35.4 Iteration ============== The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any solver of the corresponding type. The same functions work for all solvers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_multiroot_fsolver_iterate (gsl_multiroot_fsolver * S) -- Function: int gsl_multiroot_fdfsolver_iterate (gsl_multiroot_fdfsolver * S) These functions perform a single iteration of the solver S. If the iteration encounters an unexpected problem then an error code will be returned, `GSL_EBADFUNC' the iteration encountered a singular point where the function or its derivative evaluated to `Inf' or `NaN'. `GSL_ENOPROG' the iteration is not making any progress, preventing the algorithm from continuing. The solver maintains a current best estimate of the root `s->x' and its function value `s->f' at all times. This information can be accessed with the following auxiliary functions, -- Function: gsl_vector * gsl_multiroot_fsolver_root (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_root (const gsl_multiroot_fdfsolver * S) These functions return the current estimate of the root for the solver S, given by `s->x'. -- Function: gsl_vector * gsl_multiroot_fsolver_f (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_f (const gsl_multiroot_fdfsolver * S) These functions return the function value f(x) at the current estimate of the root for the solver S, given by `s->f'. -- Function: gsl_vector * gsl_multiroot_fsolver_dx (const gsl_multiroot_fsolver * S) -- Function: gsl_vector * gsl_multiroot_fdfsolver_dx (const gsl_multiroot_fdfsolver * S) These functions return the last step dx taken by the solver S, given by `s->dx'.  File: gsl-ref.info, Node: Search Stopping Parameters for the multidimensional solver, Next: Algorithms using Derivatives, Prev: Iteration of the multidimensional solver, Up: Multidimensional Root-Finding 35.5 Search Stopping Parameters =============================== A root finding procedure should stop when one of the following conditions is true: * A multidimensional root has been found to within the user-specified precision. * A user-specified maximum number of iterations has been reached. * An error has occurred. The handling of these conditions is under user control. The functions below allow the user to test the precision of the current result in several standard ways. -- Function: int gsl_multiroot_test_delta (const gsl_vector * DX, const gsl_vector * X, double EPSABS, double EPSREL) This function tests for the convergence of the sequence by comparing the last step DX with the absolute error EPSABS and relative error EPSREL to the current position X. The test returns `GSL_SUCCESS' if the following condition is achieved, |dx_i| < epsabs + epsrel |x_i| for each component of X and returns `GSL_CONTINUE' otherwise. -- Function: int gsl_multiroot_test_residual (const gsl_vector * F, double EPSABS) This function tests the residual value F against the absolute error bound EPSABS. The test returns `GSL_SUCCESS' if the following condition is achieved, \sum_i |f_i| < epsabs and returns `GSL_CONTINUE' 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.  File: gsl-ref.info, Node: Algorithms using Derivatives, Next: Algorithms without Derivatives, Prev: Search Stopping Parameters for the multidimensional solver, Up: Multidimensional Root-Finding 35.6 Algorithms using Derivatives ================================= The root finding algorithms described in this section make use of both the function and its derivative. They require an initial guess for the location of the root, but there is no absolute guarantee of convergence--the function must be suitable for this technique and the initial guess must be sufficiently close to the root for it to work. When the conditions are satisfied then convergence is quadratic. -- Derivative Solver: gsl_multiroot_fdfsolver_hybridsj This is a modified version of Powell's Hybrid method as implemented in the HYBRJ algorithm in MINPACK. Minpack was written by Jorge 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, `GSL_ENOPROG' the iteration is not making any progress, preventing the algorithm from continuing. `GSL_ENOPROGJ' re-evaluations of the Jacobian indicate that the iteration is not making any progress, preventing the algorithm from continuing. -- Derivative Solver: gsl_multiroot_fdfsolver_hybridj This algorithm is an unscaled version of `hybridsj'. 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 `hybridsj' is inappropriate. -- Derivative Solver: gsl_multiroot_fdfsolver_newton 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. If the Jacobian matrix is singular, an error code of `GSL_EDOM' is returned. -- Derivative Solver: gsl_multiroot_fdfsolver_gnewton This is 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.  File: gsl-ref.info, Node: Algorithms without Derivatives, Next: Example programs for Multidimensional Root finding, Prev: Algorithms using Derivatives, Up: Multidimensional Root-Finding 35.7 Algorithms without Derivatives =================================== The algorithms described in this section do not require any derivative information to be supplied by the user. Any derivatives needed are approximated by finite differences. Note that if the finite-differencing step size chosen by these routines is inappropriate, an explicit user-supplied numerical derivative can always be used with the algorithms described in the previous section. -- Solver: gsl_multiroot_fsolver_hybrids This is a version of the Hybrid algorithm which replaces calls to the 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'. Note that this step size will not be suitable for all problems. -- Solver: gsl_multiroot_fsolver_hybrid This is a finite difference version of the Hybrid algorithm without internal scaling. -- Solver: gsl_multiroot_fsolver_dnewton 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. -- Solver: gsl_multiroot_fsolver_broyden 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.  File: gsl-ref.info, Node: Example programs for Multidimensional Root finding, Next: References and Further Reading for Multidimensional Root Finding, Prev: Algorithms without Derivatives, Up: Multidimensional Root-Finding 35.8 Examples ============= The multidimensional solvers are used in a similar way to the one-dimensional root finding algorithms. This first example demonstrates the `hybrids' scaled-hybrid algorithm, which does not require derivatives. The program solves the Rosenbrock system of equations, f_1 (x, y) = a (1 - x) f_2 (x, y) = b (y - x^2) with a = 1, b = 10. The solution of this system lies at (x,y) = (1,1) in a narrow valley. The first stage of the program is to define the system of equations, #include #include #include #include struct rparams { double a; double b; }; int rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f) { double a = ((struct rparams *) params)->a; double b = ((struct rparams *) params)->b; const double x0 = gsl_vector_get (x, 0); const double x1 = gsl_vector_get (x, 1); const double y0 = a * (1 - x0); const double y1 = b * (x1 - x0 * x0); gsl_vector_set (f, 0, y0); gsl_vector_set (f, 1, y1); return GSL_SUCCESS; } The main program begins by creating the function object `f', with the arguments `(x,y)' and parameters `(a,b)'. The solver `s' is initialized to use this function, with the `hybrids' method. int main (void) { const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *s; int status; size_t i, iter = 0; const size_t n = 2; struct rparams p = {1.0, 10.0}; gsl_multiroot_function f = {&rosenbrock_f, n, &p}; double x_init[2] = {-10.0, -5.0}; gsl_vector *x = gsl_vector_alloc (n); gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); T = gsl_multiroot_fsolver_hybrids; s = gsl_multiroot_fsolver_alloc (T, 2); gsl_multiroot_fsolver_set (s, &f, x); print_state (iter, s); do { iter++; status = gsl_multiroot_fsolver_iterate (s); print_state (iter, s); if (status) /* check if solver is stuck */ break; status = gsl_multiroot_test_residual (s->f, 1e-7); } while (status == GSL_CONTINUE && iter < 1000); printf ("status = %s\n", gsl_strerror (status)); gsl_multiroot_fsolver_free (s); gsl_vector_free (x); return 0; } Note that it is important to check the return status of each solver step, in case the algorithm becomes stuck. If an error condition is detected, indicating that the algorithm cannot proceed, then the error can be reported to the user, a new starting point chosen or a different algorithm used. The intermediate state of the solution is displayed by the following function. The solver state contains the vector `s->x' which is the current position, and the vector `s->f' with corresponding function values. int print_state (size_t iter, gsl_multiroot_fsolver * s) { printf ("iter = %3u x = % .3f % .3f " "f(x) = % .3e % .3e\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1)); } Here are the results of running the program. The algorithm is started at (-10,-5) far from the solution. Since the solution is hidden in a narrow valley the earliest steps follow the gradient of the function downhill, in an attempt to reduce the large value of the residual. Once the root has been approximately located, on iteration 8, the Newton behavior takes over and convergence is very rapid. iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 1 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 2 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 3 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 4 x = -3.976 24.827 f(x) = 4.976e+00 9.020e+01 iter = 5 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 iter = 6 x = -1.274 -5.680 f(x) = 2.274e+00 -7.302e+01 iter = 7 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 iter = 8 x = 0.249 0.298 f(x) = 7.511e-01 2.359e+00 iter = 9 x = 1.000 0.878 f(x) = 1.268e-10 -1.218e+00 iter = 10 x = 1.000 0.989 f(x) = 1.124e-11 -1.080e-01 iter = 11 x = 1.000 1.000 f(x) = 0.000e+00 0.000e+00 status = success Note that the algorithm does not update the location on every iteration. Some iterations are used to adjust the trust-region parameter, after trying a step which was found to be divergent, or to recompute the Jacobian, when poor convergence behavior is detected. The next example program adds derivative information, in order to accelerate the solution. There are two derivative functions `rosenbrock_df' and `rosenbrock_fdf'. The latter computes both the function and its derivative simultaneously. This allows the optimization of any common terms. For simplicity we substitute calls to the separate `f' and `df' functions at this point in the code below. int rosenbrock_df (const gsl_vector * x, void *params, gsl_matrix * J) { const double a = ((struct rparams *) params)->a; const double b = ((struct rparams *) params)->b; const double x0 = gsl_vector_get (x, 0); const double df00 = -a; const double df01 = 0; const double df10 = -2 * b * x0; const double df11 = b; gsl_matrix_set (J, 0, 0, df00); gsl_matrix_set (J, 0, 1, df01); gsl_matrix_set (J, 1, 0, df10); gsl_matrix_set (J, 1, 1, df11); return GSL_SUCCESS; } int rosenbrock_fdf (const gsl_vector * x, void *params, gsl_vector * f, gsl_matrix * J) { rosenbrock_f (x, params, f); rosenbrock_df (x, params, J); return GSL_SUCCESS; } The main program now makes calls to the corresponding `fdfsolver' versions of the functions, int main (void) { const gsl_multiroot_fdfsolver_type *T; gsl_multiroot_fdfsolver *s; int status; size_t i, iter = 0; const size_t n = 2; struct rparams p = {1.0, 10.0}; gsl_multiroot_function_fdf f = {&rosenbrock_f, &rosenbrock_df, &rosenbrock_fdf, n, &p}; double x_init[2] = {-10.0, -5.0}; gsl_vector *x = gsl_vector_alloc (n); gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); T = gsl_multiroot_fdfsolver_gnewton; s = gsl_multiroot_fdfsolver_alloc (T, n); gsl_multiroot_fdfsolver_set (s, &f, x); print_state (iter, s); do { iter++; status = gsl_multiroot_fdfsolver_iterate (s); print_state (iter, s); if (status) break; status = gsl_multiroot_test_residual (s->f, 1e-7); } while (status == GSL_CONTINUE && iter < 1000); printf ("status = %s\n", gsl_strerror (status)); gsl_multiroot_fdfsolver_free (s); gsl_vector_free (x); return 0; } The addition of derivative information to the `hybrids' solver does not make any significant difference to its behavior, since it able to approximate the Jacobian numerically with sufficient accuracy. To illustrate the behavior of a different derivative solver we switch to `gnewton'. This is a traditional Newton solver with the constraint that it scales back its step if the full step would lead "uphill". Here is the output for the `gnewton' algorithm, iter = 0 x = -10.000 -5.000 f(x) = 1.100e+01 -1.050e+03 iter = 1 x = -4.231 -65.317 f(x) = 5.231e+00 -8.321e+02 iter = 2 x = 1.000 -26.358 f(x) = -8.882e-16 -2.736e+02 iter = 3 x = 1.000 1.000 f(x) = -2.220e-16 -4.441e-15 status = success The convergence is much more rapid, but takes a wide excursion out to the point (-4.23,-65.3). This could cause the algorithm to go astray in a realistic application. The hybrid algorithm follows the downhill path to the solution more reliably.  File: gsl-ref.info, Node: References and Further Reading for Multidimensional Root Finding, Prev: Example programs for Multidimensional Root finding, Up: Multidimensional Root-Finding 35.9 References and Further Reading =================================== The original version of the Hybrid method is described in the following articles by Powell, M.J.D. Powell, "A Hybrid Method for Nonlinear Equations" (Chap 6, p 87-114) and "A Fortran Subroutine for Solving systems of Nonlinear Algebraic Equations" (Chap 7, p 115-161), in `Numerical Methods for Nonlinear Algebraic Equations', P. Rabinowitz, editor. Gordon and Breach, 1970. The following papers are also relevant to the algorithms described in this section, J.J. More', M.Y. Cosnard, "Numerical Solution of Nonlinear Equations", `ACM Transactions on Mathematical Software', Vol 5, No 1, (1979), p 64-85 C.G. Broyden, "A Class of Methods for Solving Nonlinear Simultaneous Equations", `Mathematics of Computation', Vol 19 (1965), p 577-593 J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained Optimization Software", ACM Transactions on Mathematical Software, Vol 7, No 1 (1981), p 17-41  File: gsl-ref.info, Node: Multidimensional Minimization, Next: Least-Squares Fitting, Prev: Multidimensional Root-Finding, Up: Top 36 Multidimensional Minimization ******************************** This chapter describes routines for finding minima of arbitrary multidimensional functions. The library provides low level components for a variety of iterative minimizers and convergence tests. These can be combined by the user to achieve the desired solution, while providing full access to the intermediate steps of the algorithms. Each class of methods uses the same framework, so that you can switch between minimizers at runtime without needing to recompile your program. Each instance of a minimizer keeps track of its own state, allowing the minimizers to be used in multi-threaded programs. The minimization algorithms can be used to maximize a function by inverting its sign. The header file `gsl_multimin.h' contains prototypes for the minimization functions and related declarations. * Menu: * Multimin Overview:: * Multimin Caveats:: * Initializing the Multidimensional Minimizer:: * Providing a function to minimize:: * Multimin Iteration:: * Multimin Stopping Criteria:: * Multimin Algorithms with Derivatives:: * Multimin Algorithms without Derivatives:: * Multimin Examples:: * Multimin References and Further Reading::  File: gsl-ref.info, Node: Multimin Overview, Next: Multimin Caveats, Up: Multidimensional Minimization 36.1 Overview ============= The problem of multidimensional minimization requires finding a point x such that the scalar function, f(x_1, ..., x_n) takes a value which is lower than at any neighboring point. For smooth functions the gradient g = \nabla f vanishes at the minimum. In general there are no bracketing methods available for the minimization of n-dimensional functions. The algorithms proceed from an initial guess using a search algorithm which attempts to move in a downhill direction. Algorithms making use of the gradient of the function perform a one-dimensional line minimisation along this direction until the lowest point is found to a suitable tolerance. The search direction is then updated with local information from the function and its derivatives, and the whole process repeated until the true n-dimensional minimum is found. Algorithms which do not require the gradient of the function use different strategies. For example, the Nelder-Mead Simplex algorithm maintains n+1 trial parameter vectors as the vertices of a n-dimensional simplex. On each iteration it tries to improve the worst vertex of the simplex by geometrical transformations. The iterations are continued until the overall size of the simplex has decreased sufficiently. Both types of algorithms use a standard framework. The user provides a high-level driver for the algorithms, and the library provides the individual functions necessary for each of the steps. There are three main phases of the iteration. The steps are, * initialize minimizer state, S, for algorithm T * update S using the iteration T * test S for convergence, and repeat iteration if necessary Each iteration step consists either of an improvement to the line-minimisation in the current direction or an update to the search direction itself. The state for the minimizers is held in a `gsl_multimin_fdfminimizer' struct or a `gsl_multimin_fminimizer' struct.  File: gsl-ref.info, Node: Multimin Caveats, Next: Initializing the Multidimensional Minimizer, Prev: Multimin Overview, Up: Multidimensional Minimization 36.2 Caveats ============ Note that the minimization algorithms can only search for one local minimum at a time. When there are several local minima in the search area, the first minimum to be found will be returned; however it is difficult to predict which of the minima this will be. In most cases, no error will be reported if you try to find a local minimum in an area where there is more than one. It is also important to note that the minimization algorithms find local minima; there is no way to determine whether a minimum is a global minimum of the function in question.  File: gsl-ref.info, Node: Initializing the Multidimensional Minimizer, Next: Providing a function to minimize, Prev: Multimin Caveats, Up: Multidimensional Minimization 36.3 Initializing the Multidimensional Minimizer ================================================ The following function initializes a multidimensional minimizer. The minimizer itself depends only on the dimension of the problem and the algorithm and can be reused for different problems. -- Function: gsl_multimin_fdfminimizer * gsl_multimin_fdfminimizer_alloc (const gsl_multimin_fdfminimizer_type * T, size_t N) -- Function: gsl_multimin_fminimizer * gsl_multimin_fminimizer_alloc (const gsl_multimin_fminimizer_type * T, size_t N) This function returns a pointer to a newly allocated instance of a minimizer of type T for an N-dimension function. If there is insufficient memory to create the minimizer then the function returns a null pointer and the error handler is invoked with an error code of `GSL_ENOMEM'. -- Function: int gsl_multimin_fdfminimizer_set (gsl_multimin_fdfminimizer * S, gsl_multimin_function_fdf * FDF, const gsl_vector * X, double STEP_SIZE, double TOL) This function initializes the minimizer S to minimize the function FDF starting from the initial point X. The size of the first trial step is given by STEP_SIZE. The accuracy of the line minimization is specified by TOL. The precise meaning of this parameter depends on the method used. Typically the line minimization is considered successful if the gradient of the function g is orthogonal to the current search direction p to a relative accuracy of TOL, where dot(p,g) < tol |p| |g|. A TOL value of 0.1 is suitable for most purposes, since line minimization only needs to be carried out approximately. Note that setting TOL to zero will force the use of "exact" line-searches, which are extremely expensive. -- Function: int gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * S, gsl_multimin_function * F, const gsl_vector * X, const gsl_vector * STEP_SIZE) This function initializes the minimizer S to minimize the function F, starting from the initial point X. The size of the initial trial steps is given in vector STEP_SIZE. The precise meaning of this parameter depends on the method used. -- Function: void gsl_multimin_fdfminimizer_free (gsl_multimin_fdfminimizer * S) -- Function: void gsl_multimin_fminimizer_free (gsl_multimin_fminimizer * S) This function frees all the memory associated with the minimizer S. -- Function: const char * gsl_multimin_fdfminimizer_name (const gsl_multimin_fdfminimizer * S) -- Function: const char * gsl_multimin_fminimizer_name (const gsl_multimin_fminimizer * S) This function returns a pointer to the name of the minimizer. For example, printf ("s is a '%s' minimizer\n", gsl_multimin_fdfminimizer_name (s)); would print something like `s is a 'conjugate_pr' minimizer'.  File: gsl-ref.info, Node: Providing a function to minimize, Next: Multimin Iteration, Prev: Initializing the Multidimensional Minimizer, Up: Multidimensional Minimization 36.4 Providing a function to minimize ===================================== You must provide a parametric function of n variables for the minimizers to operate on. You may also need to provide a routine which calculates the gradient of the function and a third routine which calculates both the function value and the gradient together. In order to allow for general parameters the functions are defined by the following data types: -- Data Type: gsl_multimin_function_fdf This data type defines a general function of n variables with parameters and the corresponding gradient vector of derivatives, `double (* f) (const gsl_vector * X, void * PARAMS)' this function should return the result f(x,params) for argument X and parameters PARAMS. If the function cannot be computed, an error value of `GSL_NAN' should be returned. `void (* df) (const gsl_vector * X, void * PARAMS, gsl_vector * G)' this function should store the N-dimensional gradient g_i = d f(x,params) / d x_i in the vector G for argument X and parameters PARAMS, returning an appropriate error code if the function cannot be computed. `void (* fdf) (const gsl_vector * X, void * PARAMS, double * f, gsl_vector * G)' This function should set the values of the F and G as above, for arguments X and parameters PARAMS. This function provides an optimization of the separate functions for f(x) and g(x)--it is always faster to compute the function and its derivative at the same time. `size_t n' the dimension of the system, i.e. the number of components of the vectors X. `void * params' a pointer to the parameters of the function. -- Data Type: gsl_multimin_function This data type defines a general function of n variables with parameters, `double (* f) (const gsl_vector * X, void * PARAMS)' this function should return the result f(x,params) for argument X and parameters PARAMS. If the function cannot be computed, an error value of `GSL_NAN' should be returned. `size_t n' the dimension of the system, i.e. the number of components of the vectors X. `void * params' a pointer to the parameters of the function. The following example function defines a simple two-dimensional paraboloid with five parameters, /* Paraboloid centered on (p[0],p[1]), with scale factors (p[2],p[3]) and minimum p[4] */ double my_f (const gsl_vector *v, void *params) { double x, y; double *p = (double *)params; x = gsl_vector_get(v, 0); y = gsl_vector_get(v, 1); return p[2] * (x - p[0]) * (x - p[0]) + p[3] * (y - p[1]) * (y - p[1]) + p[4]; } /* The gradient of f, df = (df/dx, df/dy). */ void my_df (const gsl_vector *v, void *params, gsl_vector *df) { double x, y; double *p = (double *)params; x = gsl_vector_get(v, 0); y = gsl_vector_get(v, 1); gsl_vector_set(df, 0, 2.0 * p[2] * (x - p[0])); gsl_vector_set(df, 1, 2.0 * p[3] * (y - p[1])); } /* Compute both f and df together. */ void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *df) { *f = my_f(x, params); my_df(x, params, df); } The function can be initialized using the following code, gsl_multimin_function_fdf my_func; /* Paraboloid center at (1,2), scale factors (10, 20), minimum value 30 */ double p[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 }; my_func.n = 2; /* number of function components */ my_func.f = &my_f; my_func.df = &my_df; my_func.fdf = &my_fdf; my_func.params = (void *)p;  File: gsl-ref.info, Node: Multimin Iteration, Next: Multimin Stopping Criteria, Prev: Providing a function to minimize, Up: Multidimensional Minimization 36.5 Iteration ============== The following function drives the iteration of each algorithm. The function performs one iteration to update the state of the minimizer. The same function works for all minimizers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_multimin_fdfminimizer_iterate (gsl_multimin_fdfminimizer * S) -- Function: int gsl_multimin_fminimizer_iterate (gsl_multimin_fminimizer * S) These functions perform a single iteration of the minimizer S. If the iteration encounters an unexpected problem then an error code will be returned. The error code `GSL_ENOPROG' signifies that the minimizer is unable to improve on its current estimate, either due to numerical difficulty or because a genuine local minimum has been reached. The minimizer maintains a current best estimate of the minimum at all times. This information can be accessed with the following auxiliary functions, -- Function: gsl_vector * gsl_multimin_fdfminimizer_x (const gsl_multimin_fdfminimizer * S) -- Function: gsl_vector * gsl_multimin_fminimizer_x (const gsl_multimin_fminimizer * S) -- Function: double gsl_multimin_fdfminimizer_minimum (const gsl_multimin_fdfminimizer * S) -- Function: double gsl_multimin_fminimizer_minimum (const gsl_multimin_fminimizer * S) -- Function: gsl_vector * gsl_multimin_fdfminimizer_gradient (const gsl_multimin_fdfminimizer * S) -- Function: double gsl_multimin_fminimizer_size (const gsl_multimin_fminimizer * S) These functions return the current best estimate of the location of the minimum, the value of the function at that point, its gradient, and minimizer specific characteristic size for the minimizer S. -- Function: int gsl_multimin_fdfminimizer_restart (gsl_multimin_fdfminimizer * S) This function resets the minimizer S to use the current point as a new starting point.  File: gsl-ref.info, Node: Multimin Stopping Criteria, Next: Multimin Algorithms with Derivatives, Prev: Multimin Iteration, Up: Multidimensional Minimization 36.6 Stopping Criteria ====================== A minimization procedure should stop when one of the following conditions is true: * A minimum has been found to within the user-specified precision. * A user-specified maximum number of iterations has been reached. * An error has occurred. The handling of these conditions is under user control. The functions below allow the user to test the precision of the current result. -- Function: int gsl_multimin_test_gradient (const gsl_vector * G, double EPSABS) This function tests the norm of the gradient G against the absolute tolerance EPSABS. The gradient of a multidimensional function goes to zero at a minimum. The test returns `GSL_SUCCESS' if the following condition is achieved, |g| < epsabs and returns `GSL_CONTINUE' otherwise. A suitable choice of EPSABS can be made from the desired accuracy in the function for small variations in x. The relationship between these quantities is given by \delta f = g \delta x. -- Function: int gsl_multimin_test_size (const double SIZE, double EPSABS) This function tests the minimizer specific characteristic size (if applicable to the used minimizer) against absolute tolerance EPSABS. The test returns `GSL_SUCCESS' if the size is smaller than tolerance, otherwise `GSL_CONTINUE' is returned.  File: gsl-ref.info, Node: Multimin Algorithms with Derivatives, Next: Multimin Algorithms without Derivatives, Prev: Multimin Stopping Criteria, Up: Multidimensional Minimization 36.7 Algorithms with Derivatives ================================ There are several minimization methods available. The best choice of algorithm depends on the problem. The algorithms described in this section use the value of the function and its gradient at each evaluation point. -- Minimizer: gsl_multimin_fdfminimizer_conjugate_fr This is the Fletcher-Reeves conjugate gradient algorithm. The conjugate gradient algorithm proceeds as a succession of line minimizations. The sequence of search directions is used to build up an approximation to the curvature of the function in the neighborhood of the minimum. An initial search direction P is chosen using the gradient, and line minimization is carried out in that direction. The accuracy of the line minimization is specified by the parameter TOL. The minimum along this line occurs when the function gradient G and the search direction P are orthogonal. The line minimization terminates when dot(p,g) < tol |p| |g|. The search direction is updated using the Fletcher-Reeves formula p' = g' - \beta g where \beta=-|g'|^2/|g|^2, and the line minimization is then repeated for the new search direction. -- Minimizer: gsl_multimin_fdfminimizer_conjugate_pr This is the Polak-Ribiere conjugate gradient algorithm. It is similar to the Fletcher-Reeves method, differing only in the choice of the coefficient \beta. Both methods work well when the evaluation point is close enough to the minimum of the objective function that it is well approximated by a quadratic hypersurface. -- Minimizer: gsl_multimin_fdfminimizer_vector_bfgs2 -- Minimizer: gsl_multimin_fdfminimizer_vector_bfgs These methods use the vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. This is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region. The `bfgs2' version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's `Practical Methods of Optimization', Algorithms 2.6.2 and 2.6.4. It supersedes the original `bfgs' routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance TOL corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches). -- Minimizer: gsl_multimin_fdfminimizer_steepest_descent The steepest descent algorithm follows the downhill gradient of the function at each step. When a downhill step is successful the step-size is increased by a factor of two. If the downhill step leads to a higher function value then the algorithm backtracks and the step size is decreased using the parameter TOL. A suitable value of TOL for most applications is 0.1. The steepest descent method is inefficient and is included only for demonstration purposes.  File: gsl-ref.info, Node: Multimin Algorithms without Derivatives, Next: Multimin Examples, Prev: Multimin Algorithms with Derivatives, Up: Multidimensional Minimization 36.8 Algorithms without Derivatives =================================== The algorithms described in this section use only the value of the function at each evaluation point. -- Minimizer: gsl_multimin_fminimizer_nmsimplex2 -- Minimizer: gsl_multimin_fminimizer_nmsimplex These methods use the Simplex algorithm of Nelder and Mead. Starting from the initial vector X = p_0, the algorithm constructs an additional n vectors p_i using the step size vector s = STEP_SIZE as follows: p_0 = (x_0, x_1, ... , x_n) p_1 = (x_0 + s_0, x_1, ... , x_n) p_2 = (x_0, x_1 + s_1, ... , x_n) ... = ... p_n = (x_0, x_1, ... , x_n + s_n) These vectors form the n+1 vertices of a simplex in n dimensions. On each iteration the algorithm uses simple geometrical transformations to update the vector corresponding to the highest function value. The geometric transformations are reflection, reflection followed by expansion, contraction and multiple contraction. Using these transformations the simplex moves through the space towards the minimum, where it contracts itself. After each iteration, the best vertex is returned. Note, that due to the nature of the algorithm not every step improves the current best parameter vector. Usually several iterations are required. The minimizer-specific characteristic size is calculated as the average distance from the geometrical center of the simplex to all its vertices. This size can be used as a stopping criteria, as the simplex contracts itself near the minimum. The size is returned by the function `gsl_multimin_fminimizer_size'. The `nmsimplex2' version of this minimiser is a new O(N) operations implementation of the earlier O(N^2) operations `nmsimplex' minimiser. It uses the same underlying algorithm, but the simplex updates are computed more efficiently for high-dimensional problems. In addition, the size of simplex is calculated as the RMS distance of each vertex from the center rather than the mean distance, allowing a linear update of this quantity on each step. The memory usage is O(N^2) for both algorithms. -- Minimizer: gsl_multimin_fminimizer_nmsimplex2rand This method is a variant of `nmsimplex2' which initialises the simplex around the starting point X using a randomly-oriented set of basis vectors instead of the fixed coordinate axes. The final dimensions of the simplex are scaled along the coordinate axes by the vector STEP_SIZE. The randomisation uses a simple deterministic generator so that repeated calls to `gsl_multimin_fminimizer_set' for a given solver object will vary the orientation in a well-defined way.  File: gsl-ref.info, Node: Multimin Examples, Next: Multimin References and Further Reading, Prev: Multimin Algorithms without Derivatives, Up: Multidimensional Minimization 36.9 Examples ============= This example program finds the minimum of the paraboloid function defined earlier. The location of the minimum is offset from the origin in x and y, and the function value at the minimum is non-zero. The main program is given below, it requires the example function given earlier in this chapter. int main (void) { size_t iter = 0; int status; const gsl_multimin_fdfminimizer_type *T; gsl_multimin_fdfminimizer *s; /* Position of the minimum (1,2), scale factors 10,20, height 30. */ double par[5] = { 1.0, 2.0, 10.0, 20.0, 30.0 }; gsl_vector *x; gsl_multimin_function_fdf my_func; my_func.n = 2; my_func.f = my_f; my_func.df = my_df; my_func.fdf = my_fdf; my_func.params = par; /* Starting point, x = (5,7) */ x = gsl_vector_alloc (2); gsl_vector_set (x, 0, 5.0); gsl_vector_set (x, 1, 7.0); T = gsl_multimin_fdfminimizer_conjugate_fr; s = gsl_multimin_fdfminimizer_alloc (T, 2); gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4); do { iter++; status = gsl_multimin_fdfminimizer_iterate (s); if (status) break; status = gsl_multimin_test_gradient (s->gradient, 1e-3); if (status == GSL_SUCCESS) printf ("Minimum found at:\n"); printf ("%5d %.5f %.5f %10.5f\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), s->f); } while (status == GSL_CONTINUE && iter < 100); gsl_multimin_fdfminimizer_free (s); gsl_vector_free (x); return 0; } The initial step-size is chosen as 0.01, a conservative estimate in this case, and the line minimization parameter is set at 0.0001. The program terminates when the norm of the gradient has been reduced below 0.001. The output of the program is shown below, x y f 1 4.99629 6.99072 687.84780 2 4.98886 6.97215 683.55456 3 4.97400 6.93501 675.01278 4 4.94429 6.86073 658.10798 5 4.88487 6.71217 625.01340 6 4.76602 6.41506 561.68440 7 4.52833 5.82083 446.46694 8 4.05295 4.63238 261.79422 9 3.10219 2.25548 75.49762 10 2.85185 1.62963 67.03704 11 2.19088 1.76182 45.31640 12 0.86892 2.02622 30.18555 Minimum found at: 13 1.00000 2.00000 30.00000 Note that the algorithm gradually increases the step size as it successfully moves downhill, as can be seen by plotting the successive points. The conjugate gradient algorithm finds the minimum on its second direction because the function is purely quadratic. Additional iterations would be needed for a more complicated function. Here is another example using the Nelder-Mead Simplex algorithm to minimize the same example object function, as above. int main(void) { double par[5] = {1.0, 2.0, 10.0, 20.0, 30.0}; const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2; gsl_multimin_fminimizer *s = NULL; gsl_vector *ss, *x; gsl_multimin_function minex_func; size_t iter = 0; int status; double size; /* Starting point */ x = gsl_vector_alloc (2); gsl_vector_set (x, 0, 5.0); gsl_vector_set (x, 1, 7.0); /* Set initial step sizes to 1 */ ss = gsl_vector_alloc (2); gsl_vector_set_all (ss, 1.0); /* Initialize method and iterate */ minex_func.n = 2; minex_func.f = my_f; minex_func.params = par; s = gsl_multimin_fminimizer_alloc (T, 2); gsl_multimin_fminimizer_set (s, &minex_func, x, ss); do { iter++; status = gsl_multimin_fminimizer_iterate(s); if (status) break; size = gsl_multimin_fminimizer_size (s); status = gsl_multimin_test_size (size, 1e-2); if (status == GSL_SUCCESS) { printf ("converged to minimum at\n"); } printf ("%5d %10.3e %10.3e f() = %7.3f size = %.3f\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), s->fval, size); } while (status == GSL_CONTINUE && iter < 100); gsl_vector_free(x); gsl_vector_free(ss); gsl_multimin_fminimizer_free (s); return status; } The minimum search stops when the Simplex size drops to 0.01. The output is shown below. 1 6.500e+00 5.000e+00 f() = 512.500 size = 1.130 2 5.250e+00 4.000e+00 f() = 290.625 size = 1.409 3 5.250e+00 4.000e+00 f() = 290.625 size = 1.409 4 5.500e+00 1.000e+00 f() = 252.500 size = 1.409 5 2.625e+00 3.500e+00 f() = 101.406 size = 1.847 6 2.625e+00 3.500e+00 f() = 101.406 size = 1.847 7 0.000e+00 3.000e+00 f() = 60.000 size = 1.847 8 2.094e+00 1.875e+00 f() = 42.275 size = 1.321 9 2.578e-01 1.906e+00 f() = 35.684 size = 1.069 10 5.879e-01 2.445e+00 f() = 35.664 size = 0.841 11 1.258e+00 2.025e+00 f() = 30.680 size = 0.476 12 1.258e+00 2.025e+00 f() = 30.680 size = 0.367 13 1.093e+00 1.849e+00 f() = 30.539 size = 0.300 14 8.830e-01 2.004e+00 f() = 30.137 size = 0.172 15 8.830e-01 2.004e+00 f() = 30.137 size = 0.126 16 9.582e-01 2.060e+00 f() = 30.090 size = 0.106 17 1.022e+00 2.004e+00 f() = 30.005 size = 0.063 18 1.022e+00 2.004e+00 f() = 30.005 size = 0.043 19 1.022e+00 2.004e+00 f() = 30.005 size = 0.043 20 1.022e+00 2.004e+00 f() = 30.005 size = 0.027 21 1.022e+00 2.004e+00 f() = 30.005 size = 0.022 22 9.920e-01 1.997e+00 f() = 30.001 size = 0.016 23 9.920e-01 1.997e+00 f() = 30.001 size = 0.013 converged to minimum at 24 9.920e-01 1.997e+00 f() = 30.001 size = 0.008 The simplex size first increases, while the simplex moves towards the minimum. After a while the size begins to decrease as the simplex contracts around the minimum.  File: gsl-ref.info, Node: Multimin References and Further Reading, Prev: Multimin Examples, Up: Multidimensional Minimization 36.10 References and Further Reading ==================================== The conjugate gradient and BFGS methods are described in detail in the following book, R. Fletcher, `Practical Methods of Optimization (Second Edition)' Wiley (1987), ISBN 0471915475. A brief description of multidimensional minimization algorithms and more recent references can be found in, C.W. Ueberhuber, `Numerical Computation (Volume 2)', Chapter 14, Section 4.4 "Minimization Methods", p. 325-335, Springer (1997), ISBN 3-540-62057-5. The simplex algorithm is described in the following paper, J.A. Nelder and R. Mead, `A simplex method for function minimization', Computer Journal vol. 7 (1965), 308-313.  File: gsl-ref.info, Node: Least-Squares Fitting, Next: Nonlinear Least-Squares Fitting, Prev: Multidimensional Minimization, Up: Top 37 Least-Squares Fitting ************************ This chapter describes routines for performing least squares fits to experimental data using linear combinations of functions. The data may be weighted or unweighted, i.e. with known or unknown errors. For weighted data the functions compute the best fit parameters and their associated covariance matrix. For unweighted data the covariance matrix is estimated from the scatter of the points, giving a variance-covariance matrix. The functions are divided into separate versions for simple one- or two-parameter regression and multiple-parameter fits. The functions are declared in the header file `gsl_fit.h'. * Menu: * Fitting Overview:: * Linear regression:: * Linear fitting without a constant term:: * Multi-parameter fitting:: * Fitting Examples:: * Fitting References and Further Reading::  File: gsl-ref.info, Node: Fitting Overview, Next: Linear regression, Up: Least-Squares Fitting 37.1 Overview ============= Least-squares fits are found by minimizing \chi^2 (chi-squared), the weighted sum of squared residuals over n experimental datapoints (x_i, y_i) for the model Y(c,x), \chi^2 = \sum_i w_i (y_i - Y(c, x_i))^2 The p parameters of the model are c = {c_0, c_1, ...}. The weight factors w_i are given by w_i = 1/\sigma_i^2, where \sigma_i is the experimental error on the data-point y_i. The errors are assumed to be Gaussian and uncorrelated. For unweighted data the chi-squared sum is computed without any weight factors. The fitting routines return the best-fit parameters c and their p \times p covariance matrix. The covariance matrix measures the statistical errors on the best-fit parameters resulting from the errors on the data, \sigma_i, and is defined as C_{ab} = <\delta c_a \delta c_b> where < > denotes an average over the Gaussian error distributions of the underlying datapoints. The covariance matrix is calculated by error propagation from the data errors \sigma_i. The change in a fitted parameter \delta c_a caused by a small change in the data \delta y_i is given by \delta c_a = \sum_i (dc_a/dy_i) \delta y_i allowing the covariance matrix to be written in terms of the errors on the data, C_{ab} = \sum_{i,j} (dc_a/dy_i) (dc_b/dy_j) <\delta y_i \delta y_j> For uncorrelated data the fluctuations of the underlying datapoints satisfy <\delta y_i \delta y_j> = \sigma_i^2 \delta_{ij}, giving a corresponding parameter covariance matrix of C_{ab} = \sum_i (1/w_i) (dc_a/dy_i) (dc_b/dy_i) When computing the covariance matrix for unweighted data, i.e. data with unknown errors, the weight factors w_i in this sum are replaced by the single estimate w = 1/\sigma^2, where \sigma^2 is the computed variance of the residuals about the best-fit model, \sigma^2 = \sum (y_i - Y(c,x_i))^2 / (n-p). This is referred to as the "variance-covariance matrix". The standard deviations of the best-fit parameters are given by the square root of the corresponding diagonal elements of the covariance matrix, \sigma_{c_a} = \sqrt{C_{aa}}. The correlation coefficient of the fit parameters c_a and c_b is given by \rho_{ab} = C_{ab} / \sqrt{C_{aa} C_{bb}}.  File: gsl-ref.info, Node: Linear regression, Next: Linear fitting without a constant term, Prev: Fitting Overview, Up: Least-Squares Fitting 37.2 Linear regression ====================== The functions described in this section can be used to perform least-squares fits to a straight line model, Y(c,x) = c_0 + c_1 x. -- Function: int gsl_fit_linear (const double * X, const size_t XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * SUMSQ) This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The errors on Y are assumed unknown so the variance-covariance matrix for the parameters (C0, C1) is estimated from the scatter of the points around the best-fit line and returned via the parameters (COV00, COV01, COV11). The sum of squares of the residuals from the best-fit line is returned in SUMSQ. Note: the correlation coefficient of the data can be computed using `gsl_stats_correlation' (*note Correlation::), it does not depend on the fit. -- Function: int gsl_fit_wlinear (const double * X, const size_t XSTRIDE, const double * W, const size_t WSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C0, double * C1, double * COV00, double * COV01, double * COV11, double * CHISQ) This function computes the best-fit linear regression coefficients (C0,C1) of the model Y = c_0 + c_1 X for the weighted dataset (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The vector W, of length N and stride WSTRIDE, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in Y. The covariance matrix for the parameters (C0, C1) is computed using the weights and returned via the parameters (COV00, COV01, COV11). The weighted sum of squares of the residuals from the best-fit line, \chi^2, is returned in CHISQ. -- Function: int gsl_fit_linear_est (double X, double C0, double C1, double COV00, double COV01, double COV11, double * Y, double * Y_ERR) This function uses the best-fit linear regression coefficients C0, C1 and their covariance COV00, COV01, COV11 to compute the fitted function Y and its standard deviation Y_ERR for the model Y = c_0 + c_1 X at the point X.  File: gsl-ref.info, Node: Linear fitting without a constant term, Next: Multi-parameter fitting, Prev: Linear regression, Up: Least-Squares Fitting 37.3 Linear fitting without a constant term =========================================== The functions described in this section can be used to perform least-squares fits to a straight line model without a constant term, Y = c_1 X. -- Function: int gsl_fit_mul (const double * X, const size_t XSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C1, double * COV11, double * SUMSQ) This function computes the best-fit linear regression coefficient C1 of the model Y = c_1 X for the datasets (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The errors on Y are assumed unknown so the variance of the parameter C1 is estimated from the scatter of the points around the best-fit line and returned via the parameter COV11. The sum of squares of the residuals from the best-fit line is returned in SUMSQ. -- Function: int gsl_fit_wmul (const double * X, const size_t XSTRIDE, const double * W, const size_t WSTRIDE, const double * Y, const size_t YSTRIDE, size_t N, double * C1, double * COV11, double * SUMSQ) This function computes the best-fit linear regression coefficient C1 of the model Y = c_1 X for the weighted datasets (X, Y), two vectors of length N with strides XSTRIDE and YSTRIDE. The vector W, of length N and stride WSTRIDE, specifies the weight of each datapoint. The weight is the reciprocal of the variance for each datapoint in Y. The variance of the parameter C1 is computed using the weights and returned via the parameter COV11. The weighted sum of squares of the residuals from the best-fit line, \chi^2, is returned in CHISQ. -- Function: int gsl_fit_mul_est (double X, double C1, double COV11, double * Y, double * Y_ERR) This function uses the best-fit linear regression coefficient C1 and its covariance COV11 to compute the fitted function Y and its standard deviation Y_ERR for the model Y = c_1 X at the point X.  File: gsl-ref.info, Node: Multi-parameter fitting, Next: Fitting Examples, Prev: Linear fitting without a constant term, Up: Least-Squares Fitting 37.4 Multi-parameter fitting ============================ The functions described in this section perform least-squares fits to a general linear model, y = X c where y is a vector of n observations, X is an n by p matrix of predictor variables, and the elements of the vector c are the p unknown best-fit parameters which are to be estimated. The chi-squared value is given by \chi^2 = \sum_i w_i (y_i - \sum_j X_{ij} c_j)^2. This formulation can be used for fits to any number of functions and/or variables by preparing the n-by-p matrix X appropriately. For example, to fit to a p-th order polynomial in X, use the following matrix, X_{ij} = x_i^j where the index i runs over the observations and the index j runs from 0 to p-1. To fit to a set of p sinusoidal functions with fixed frequencies \omega_1, \omega_2, ..., \omega_p, use, X_{ij} = sin(\omega_j x_i) To fit to p independent variables x_1, x_2, ..., x_p, use, X_{ij} = x_j(i) where x_j(i) is the i-th value of the predictor variable x_j. The functions described in this section are declared in the header file `gsl_multifit.h'. The solution of the general linear least-squares system requires an additional working space for intermediate results, such as the singular value decomposition of the matrix X. -- Function: gsl_multifit_linear_workspace * gsl_multifit_linear_alloc (size_t N, size_t P) This function allocates a workspace for fitting a model to N observations using P parameters. -- Function: void gsl_multifit_linear_free (gsl_multifit_linear_workspace * WORK) This function frees the memory associated with the workspace W. -- Function: int gsl_multifit_linear (const gsl_matrix * X, const gsl_vector * Y, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) This function computes the best-fit parameters C of the model y = X c for the observations Y and the matrix of predictor variables X, using the preallocated workspace provided in WORK. The variance-covariance matrix of the model parameters COV is estimated from the scatter of the observations about the best-fit. The sum of squares of the residuals from the best-fit, \chi^2, is returned in CHISQ. If the coefficient of determination is desired, it can be computed from the expression R^2 = 1 - \chi^2 / TSS, where the total sum of squares (TSS) of the observations Y may be computed from `gsl_stats_tss'. The best-fit is found by singular value decomposition of the matrix X using the modified Golub-Reinsch SVD algorithm, with column scaling to improve the accuracy of the singular values. Any components which have zero singular value (to machine precision) are discarded from the fit. -- Function: int gsl_multifit_wlinear (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) This function computes the best-fit parameters C of the weighted model y = X c for the observations Y with weights W and the matrix of predictor variables X, using the preallocated workspace provided in WORK. The covariance matrix of the model parameters COV is computed with the given weights. The weighted sum of squares of the residuals from the best-fit, \chi^2, is returned in CHISQ. If the coefficient of determination is desired, it can be computed from the expression R^2 = 1 - \chi^2 / WTSS, where the weighted total sum of squares (WTSS) of the observations Y may be computed from `gsl_stats_wtss'. -- Function: int gsl_multifit_linear_svd (const gsl_matrix * X, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) -- Function: int gsl_multifit_wlinear_svd (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) In these functions components of the fit are discarded if the ratio of singular values s_i/s_0 falls below the user-specified tolerance TOL, and the effective rank is returned in RANK. -- Function: int gsl_multifit_linear_usvd (const gsl_matrix * X, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) -- Function: int gsl_multifit_wlinear_usvd (const gsl_matrix * X, const gsl_vector * W, const gsl_vector * Y, double TOL, size_t * RANK, gsl_vector * C, gsl_matrix * COV, double * CHISQ, gsl_multifit_linear_workspace * WORK) These functions compute the fit using an SVD without column scaling. -- Function: int gsl_multifit_linear_est (const gsl_vector * X, const gsl_vector * C, const gsl_matrix * COV, double * Y, double * Y_ERR) This function uses the best-fit multilinear regression coefficients C and their covariance matrix COV to compute the fitted function value Y and its standard deviation Y_ERR for the model y = x.c at the point X. -- Function: int gsl_multifit_linear_residuals (const gsl_matrix * X, const gsl_vector * Y, const gsl_vector * C, gsl_vector * R) This function computes the vector of residuals r = y - X c for the observations Y, coefficients C and matrix of predictor variables X.  File: gsl-ref.info, Node: Fitting Examples, Next: Fitting References and Further Reading, Prev: Multi-parameter fitting, Up: Least-Squares Fitting 37.5 Examples ============= The following program computes a least squares straight-line fit to a simple dataset, and outputs the best-fit line and its associated one standard-deviation error bars. #include #include int main (void) { int i, n = 4; double x[4] = { 1970, 1980, 1990, 2000 }; double y[4] = { 12, 11, 14, 13 }; double w[4] = { 0.1, 0.2, 0.3, 0.4 }; double c0, c1, cov00, cov01, cov11, chisq; gsl_fit_wlinear (x, 1, w, 1, y, 1, n, &c0, &c1, &cov00, &cov01, &cov11, &chisq); printf ("# best fit: Y = %g + %g X\n", c0, c1); printf ("# covariance matrix:\n"); printf ("# [ %g, %g\n# %g, %g]\n", cov00, cov01, cov01, cov11); printf ("# chisq = %g\n", chisq); for (i = 0; i < n; i++) printf ("data: %g %g %g\n", x[i], y[i], 1/sqrt(w[i])); printf ("\n"); for (i = -30; i < 130; i++) { double xf = x[0] + (i/100.0) * (x[n-1] - x[0]); double yf, yf_err; gsl_fit_linear_est (xf, c0, c1, cov00, cov01, cov11, &yf, &yf_err); printf ("fit: %g %g\n", xf, yf); printf ("hi : %g %g\n", xf, yf + yf_err); printf ("lo : %g %g\n", xf, yf - yf_err); } return 0; } The following commands extract the data from the output of the program and display it using the GNU plotutils `graph' utility, $ ./demo > tmp $ more tmp # best fit: Y = -106.6 + 0.06 X # covariance matrix: # [ 39602, -19.9 # -19.9, 0.01] # chisq = 0.8 $ for n in data fit hi lo ; do grep "^$n" tmp | cut -d: -f2 > $n ; done $ graph -T X -X x -Y y -y 0 20 -m 0 -S 2 -Ie data -S 0 -I a -m 1 fit -m 2 hi -m 2 lo The next program performs a quadratic fit y = c_0 + c_1 x + c_2 x^2 to a weighted dataset using the generalised linear fitting function `gsl_multifit_wlinear'. The model matrix X for a quadratic fit is given by, X = [ 1 , x_0 , x_0^2 ; 1 , x_1 , x_1^2 ; 1 , x_2 , x_2^2 ; ... , ... , ... ] where the column of ones corresponds to the constant term c_0. The two remaining columns corresponds to the terms c_1 x and c_2 x^2. The program reads N lines of data in the format (X, Y, ERR) where ERR is the error (standard deviation) in the value Y. #include #include int main (int argc, char **argv) { int i, n; double xi, yi, ei, chisq; gsl_matrix *X, *cov; gsl_vector *y, *w, *c; if (argc != 2) { fprintf (stderr,"usage: fit n < data\n"); exit (-1); } n = atoi (argv[1]); X = gsl_matrix_alloc (n, 3); y = gsl_vector_alloc (n); w = gsl_vector_alloc (n); c = gsl_vector_alloc (3); cov = gsl_matrix_alloc (3, 3); for (i = 0; i < n; i++) { int count = fscanf (stdin, "%lg %lg %lg", &xi, &yi, &ei); if (count != 3) { fprintf (stderr, "error reading file\n"); exit (-1); } printf ("%g %g +/- %g\n", xi, yi, ei); gsl_matrix_set (X, i, 0, 1.0); gsl_matrix_set (X, i, 1, xi); gsl_matrix_set (X, i, 2, xi*xi); gsl_vector_set (y, i, yi); gsl_vector_set (w, i, 1.0/(ei*ei)); } { gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (n, 3); gsl_multifit_wlinear (X, w, y, c, cov, &chisq, work); gsl_multifit_linear_free (work); } #define C(i) (gsl_vector_get(c,(i))) #define COV(i,j) (gsl_matrix_get(cov,(i),(j))) { printf ("# best fit: Y = %g + %g X + %g X^2\n", C(0), C(1), C(2)); printf ("# covariance matrix:\n"); printf ("[ %+.5e, %+.5e, %+.5e \n", COV(0,0), COV(0,1), COV(0,2)); printf (" %+.5e, %+.5e, %+.5e \n", COV(1,0), COV(1,1), COV(1,2)); printf (" %+.5e, %+.5e, %+.5e ]\n", COV(2,0), COV(2,1), COV(2,2)); printf ("# chisq = %g\n", chisq); } gsl_matrix_free (X); gsl_vector_free (y); gsl_vector_free (w); gsl_vector_free (c); gsl_matrix_free (cov); return 0; } A suitable set of data for fitting can be generated using the following program. It outputs a set of points with gaussian errors from the curve y = e^x in the region 0 < x < 2. #include #include #include int main (void) { double x; const gsl_rng_type * T; gsl_rng * r; gsl_rng_env_setup (); T = gsl_rng_default; r = gsl_rng_alloc (T); for (x = 0.1; x < 2; x+= 0.1) { double y0 = exp (x); double sigma = 0.1 * y0; double dy = gsl_ran_gaussian (r, sigma); printf ("%g %g %g\n", x, y0 + dy, sigma); } gsl_rng_free(r); return 0; } The data can be prepared by running the resulting executable program, $ ./generate > exp.dat $ more exp.dat 0.1 0.97935 0.110517 0.2 1.3359 0.12214 0.3 1.52573 0.134986 0.4 1.60318 0.149182 0.5 1.81731 0.164872 0.6 1.92475 0.182212 .... To fit the data use the previous program, with the number of data points given as the first argument. In this case there are 19 data points. $ ./fit 19 < exp.dat 0.1 0.97935 +/- 0.110517 0.2 1.3359 +/- 0.12214 ... # best fit: Y = 1.02318 + 0.956201 X + 0.876796 X^2 # covariance matrix: [ +1.25612e-02, -3.64387e-02, +1.94389e-02 -3.64387e-02, +1.42339e-01, -8.48761e-02 +1.94389e-02, -8.48761e-02, +5.60243e-02 ] # chisq = 23.0987 The parameters of the quadratic fit match the coefficients of the expansion of e^x, taking into account the errors on the parameters and the O(x^3) difference between the exponential and quadratic functions for the larger values of x. The errors on the parameters are given by the square-root of the corresponding diagonal elements of the covariance matrix. The chi-squared per degree of freedom is 1.4, indicating a reasonable fit to the data.  File: gsl-ref.info, Node: Fitting References and Further Reading, Prev: Fitting Examples, Up: Least-Squares Fitting 37.6 References and Further Reading =================================== A summary of formulas and techniques for least squares fitting can be found in the "Statistics" chapter of the Annual Review of Particle Physics prepared by the Particle Data Group, `Review of Particle Properties', R.M. Barnett et al., Physical Review D54, 1 (1996) `http://pdg.lbl.gov/' The Review of Particle Physics is available online at the website given above. The tests used to prepare these routines are based on the NIST Statistical Reference Datasets. The datasets and their documentation are available from NIST at the following website, `http://www.nist.gov/itl/div898/strd/index.html'.  File: gsl-ref.info, Node: Nonlinear Least-Squares Fitting, Next: Basis Splines, Prev: Least-Squares Fitting, Up: Top 38 Nonlinear Least-Squares Fitting ********************************** This chapter describes functions for multidimensional nonlinear least-squares fitting. The library provides low level components for a variety of iterative solvers and convergence tests. These can be combined by the user to achieve the desired solution, with full access to the intermediate steps of the iteration. Each class of methods uses the same framework, so that you can switch between solvers at runtime without needing to recompile your program. Each instance of a solver keeps track of its own state, allowing the solvers to be used in multi-threaded programs. The header file `gsl_multifit_nlin.h' contains prototypes for the multidimensional nonlinear fitting functions and related declarations. * Menu: * Overview of Nonlinear Least-Squares Fitting:: * Initializing the Nonlinear Least-Squares Solver:: * Providing the Function to be Minimized:: * Iteration of the Minimization Algorithm:: * Search Stopping Parameters for Minimization Algorithms:: * Minimization Algorithms using Derivatives:: * Minimization Algorithms without Derivatives:: * Computing the covariance matrix of best fit parameters:: * Example programs for Nonlinear Least-Squares Fitting:: * References and Further Reading for Nonlinear Least-Squares Fitting::  File: gsl-ref.info, Node: Overview of Nonlinear Least-Squares Fitting, Next: Initializing the Nonlinear Least-Squares Solver, Up: Nonlinear Least-Squares Fitting 38.1 Overview ============= The problem of multidimensional nonlinear least-squares fitting requires the minimization of the squared residuals of n functions, f_i, in p parameters, x_i, \Phi(x) = (1/2) || F(x) ||^2 = (1/2) \sum_{i=1}^{n} f_i(x_1, ..., x_p)^2 All algorithms proceed from an initial guess using the linearization, \psi(p) = || F(x+p) || ~=~ || F(x) + J p || where x is the initial point, p is the proposed step and J is the Jacobian matrix J_{ij} = d f_i / d x_j. Additional strategies are used to enlarge the region of convergence. These include requiring a decrease in the norm ||F|| on each step or using a trust region to avoid steps which fall outside the linear regime. To perform a weighted least-squares fit of a nonlinear model Y(x,t) to data (t_i, y_i) with independent Gaussian errors \sigma_i, use function components of the following form, f_i = (Y(x, t_i) - y_i) / \sigma_i Note that the model parameters are denoted by x in this chapter since the non-linear least-squares algorithms are described geometrically (i.e. finding the minimum of a surface). The independent variable of any data to be fitted is denoted by t. With the definition above the Jacobian is J_{ij} =(1 / \sigma_i) d Y_i / d x_j, where Y_i = Y(x,t_i).  File: gsl-ref.info, Node: Initializing the Nonlinear Least-Squares Solver, Next: Providing the Function to be Minimized, Prev: Overview of Nonlinear Least-Squares Fitting, Up: Nonlinear Least-Squares Fitting 38.2 Initializing the Solver ============================ -- Function: gsl_multifit_fsolver * gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T, size_t N, size_t P) This function returns a pointer to a newly allocated instance of a solver of type T for N observations and P parameters. The number of observations N must be greater than or equal to parameters P. If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of `GSL_ENOMEM'. -- Function: gsl_multifit_fdfsolver * gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T, size_t N, size_t P) This function returns a pointer to a newly allocated instance of a derivative solver of type T for N observations and P parameters. For example, the following code creates an instance of a Levenberg-Marquardt solver for 100 data points and 3 parameters, const gsl_multifit_fdfsolver_type * T = gsl_multifit_fdfsolver_lmder; gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc (T, 100, 3); The number of observations N must be greater than or equal to parameters P. If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of `GSL_ENOMEM'. -- Function: int gsl_multifit_fsolver_set (gsl_multifit_fsolver * S, gsl_multifit_function * F, const gsl_vector * X) This function initializes, or reinitializes, an existing solver S to use the function F and the initial guess X. -- Function: int gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * S, gsl_multifit_function_fdf * FDF, const gsl_vector * X) This function initializes, or reinitializes, an existing solver S to use the function and derivative FDF and the initial guess X. -- Function: void gsl_multifit_fsolver_free (gsl_multifit_fsolver * S) -- Function: void gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * S) These functions free all the memory associated with the solver S. -- Function: const char * gsl_multifit_fsolver_name (const gsl_multifit_fsolver * S) -- Function: const char * gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * S) These functions return a pointer to the name of the solver. For example, printf ("s is a '%s' solver\n", gsl_multifit_fdfsolver_name (s)); would print something like `s is a 'lmder' solver'.  File: gsl-ref.info, Node: Providing the Function to be Minimized, Next: Iteration of the Minimization Algorithm, Prev: Initializing the Nonlinear Least-Squares Solver, Up: Nonlinear Least-Squares Fitting 38.3 Providing the Function to be Minimized =========================================== You must provide n functions of p variables for the minimization algorithms to operate on. In order to allow for arbitrary parameters the functions are defined by the following data types: -- Data Type: gsl_multifit_function This data type defines a general system of functions with arbitrary parameters. `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and arbitrary parameters PARAMS, returning an appropriate error code if the function cannot be computed. `size_t n' the number of functions, i.e. the number of components of the vector F. `size_t p' the number of independent variables, i.e. the number of components of the vector X. `void * params' a pointer to the arbitrary parameters of the function. -- Data Type: gsl_multifit_function_fdf This data type defines a general system of functions with arbitrary parameters and the corresponding Jacobian matrix of derivatives, `int (* f) (const gsl_vector * X, void * PARAMS, gsl_vector * F)' this function should store the vector result f(x,params) in F for argument X and arbitrary parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* df) (const gsl_vector * X, void * PARAMS, gsl_matrix * J)' this function should store the N-by-P matrix result J_ij = d f_i(x,params) / d x_j in J for argument X and arbitrary parameters PARAMS, returning an appropriate error code if the function cannot be computed. `int (* fdf) (const gsl_vector * X, void * PARAMS, gsl_vector * F, gsl_matrix * J)' This function should set the values of the F and J as above, for arguments X and arbitrary parameters PARAMS. This function provides an optimization of the separate functions for f(x) and J(x)--it is always faster to compute the function and its derivative at the same time. `size_t n' the number of functions, i.e. the number of components of the vector F. `size_t p' the number of independent variables, i.e. the number of components of the vector X. `void * params' a pointer to the arbitrary parameters of the function. Note that when fitting a non-linear model against experimental data, the data is passed to the functions above using the PARAMS argument and the trial best-fit parameters through the X argument.  File: gsl-ref.info, Node: Iteration of the Minimization Algorithm, Next: Search Stopping Parameters for Minimization Algorithms, Prev: Providing the Function to be Minimized, Up: Nonlinear Least-Squares Fitting 38.4 Iteration ============== The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any solver of the corresponding type. The same functions work for all solvers so that different methods can be substituted at runtime without modifications to the code. -- Function: int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * S) -- Function: int gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * S) These functions perform a single iteration of the solver S. If the iteration encounters an unexpected problem then an error code will be returned. The solver maintains a current estimate of the best-fit parameters at all times. The solver struct S contains the following entries, which can be used to track the progress of the solution: `gsl_vector * x' The current position. `gsl_vector * f' The function value at the current position. `gsl_vector * dx' The difference between the current position and the previous position, i.e. the last step, taken as a vector. `gsl_matrix * J' The Jacobian matrix at the current position (for the `gsl_multifit_fdfsolver' struct only) The best-fit information also can be accessed with the following auxiliary functions, -- Function: gsl_vector * gsl_multifit_fsolver_position (const gsl_multifit_fsolver * S) -- Function: gsl_vector * gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * S) These functions return the current position (i.e. best-fit parameters) `s->x' of the solver S.  File: gsl-ref.info, Node: Search Stopping Parameters for Minimization Algorithms, Next: Minimization Algorithms using Derivatives, Prev: Iteration of the Minimization Algorithm, Up: Nonlinear Least-Squares Fitting 38.5 Search Stopping Parameters =============================== A minimization procedure should stop when one of the following conditions is true: * A minimum has been found to within the user-specified precision. * A user-specified maximum number of iterations has been reached. * An error has occurred. The handling of these conditions is under user control. The functions below allow the user to test the current estimate of the best-fit parameters in several standard ways. -- Function: int gsl_multifit_test_delta (const gsl_vector * DX, const gsl_vector * X, double EPSABS, double EPSREL) This function tests for the convergence of the sequence by comparing the last step DX with the absolute error EPSABS and relative error EPSREL to the current position X. The test returns `GSL_SUCCESS' if the following condition is achieved, |dx_i| < epsabs + epsrel |x_i| for each component of X and returns `GSL_CONTINUE' otherwise. -- Function: int gsl_multifit_test_gradient (const gsl_vector * G, double EPSABS) This function tests the residual gradient G against the absolute error bound EPSABS. Mathematically, the gradient should be exactly zero at the minimum. The test returns `GSL_SUCCESS' if the following condition is achieved, \sum_i |g_i| < epsabs and returns `GSL_CONTINUE' otherwise. This criterion is suitable for situations where the precise location of the minimum, x, is unimportant provided a value can be found where the gradient is small enough. -- Function: int gsl_multifit_gradient (const gsl_matrix * J, const gsl_vector * F, gsl_vector * G) This function computes the gradient G of \Phi(x) = (1/2) ||F(x)||^2 from the Jacobian matrix J and the function values F, using the formula g = J^T f.  File: gsl-ref.info, Node: Minimization Algorithms using Derivatives, Next: Minimization Algorithms without Derivatives, Prev: Search Stopping Parameters for Minimization Algorithms, Up: Nonlinear Least-Squares Fitting 38.6 Minimization Algorithms using Derivatives ============================================== The minimization algorithms described in this section make use of both the function and its derivative. They require an initial guess for the location of the minimum. There is no absolute guarantee of convergence--the function must be suitable for this technique and the initial guess must be sufficiently close to the minimum for it to work. -- Derivative Solver: gsl_multifit_fdfsolver_lmsder This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled LMDER routine in MINPACK. Minpack was written by Jorge J. More', Burton S. Garbow and Kenneth E. Hillstrom. 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 attempts to minimize the linear system |F + J p| subject to the constraint |D p| < \Delta. The solution to this constrained linear system is found using the Levenberg-Marquardt method. 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, and follows the predicted behavior of the function within the trust region, then it is accepted and the size of the trust region is increased. If the proposed step fails to improve the solution, or differs significantly from the expected behavior within the trust region, then the size of the trust region is decreased and another trial step is computed. The algorithm also monitors the progress of the solution and returns an error if the changes in the solution are smaller than the machine precision. The possible error codes are, `GSL_ETOLF' the decrease in the function falls below machine precision `GSL_ETOLX' the change in the position vector falls below machine precision `GSL_ETOLG' the norm of the gradient, relative to the norm of the function, falls below machine precision `GSL_ENOPROG' the routine has made 10 or more attempts to find a suitable trial step without success (but subsequent calls can be made to continue the search).(1) These error codes indicate that further iterations will be unlikely to change the solution from its current value. -- Derivative Solver: gsl_multifit_fdfsolver_lmder This is an unscaled version of the LMDER algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of LMDER converges too slowly, or the function is already scaled appropriately. ---------- Footnotes ---------- (1) The return code `GSL_CONTINUE' was used for this case in versions prior to 1.14.  File: gsl-ref.info, Node: Minimization Algorithms without Derivatives, Next: Computing the covariance matrix of best fit parameters, Prev: Minimization Algorithms using Derivatives, Up: Nonlinear Least-Squares Fitting 38.7 Minimization Algorithms without Derivatives ================================================ There are no algorithms implemented in this section at the moment.  File: gsl-ref.info, Node: Computing the covariance matrix of best fit parameters, Next: Example programs for Nonlinear Least-Squares Fitting, Prev: Minimization Algorithms without Derivatives, Up: Nonlinear Least-Squares Fitting 38.8 Computing the covariance matrix of best fit parameters =========================================================== -- Function: int gsl_multifit_covar (const gsl_matrix * J, double EPSREL, gsl_matrix * COVAR) This function uses the Jacobian matrix J to compute the covariance matrix of the best-fit parameters, COVAR. The parameter EPSREL is used to remove linear-dependent columns when J is rank deficient. The covariance matrix is given by, covar = (J^T J)^{-1} and is computed by QR decomposition of J with column-pivoting. Any columns of R which satisfy |R_{kk}| <= epsrel |R_{11}| are considered linearly-dependent and are excluded from the covariance matrix (the corresponding rows and columns of the covariance matrix are set to zero). If the minimisation uses the weighted least-squares function f_i = (Y(x, t_i) - y_i) / \sigma_i then the covariance matrix above gives the statistical error on the best-fit parameters resulting from the Gaussian errors \sigma_i on the underlying data y_i. This can be verified from the relation \delta f = J \delta c and the fact that the fluctuations in f from the data y_i are normalised by \sigma_i and so satisfy <\delta f \delta f^T> = I. For an unweighted least-squares function f_i = (Y(x, t_i) - y_i) the covariance matrix above should be multiplied by the variance of the residuals about the best-fit \sigma^2 = \sum (y_i - Y(x,t_i))^2 / (n-p) to give the variance-covariance matrix \sigma^2 C. This estimates the statistical error on the best-fit parameters from the scatter of the underlying data. For more information about covariance matrices see *note Fitting Overview::.  File: gsl-ref.info, Node: Example programs for Nonlinear Least-Squares Fitting, Next: References and Further Reading for Nonlinear Least-Squares Fitting, Prev: Computing the covariance matrix of best fit parameters, Up: Nonlinear Least-Squares Fitting 38.9 Examples ============= The following example program fits a weighted exponential model with background to experimental data, Y = A \exp(-\lambda t) + b. The first part of the program sets up the functions `expb_f' and `expb_df' to calculate the model and its Jacobian. The appropriate fitting function is given by, f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i where we have chosen t_i = i. The Jacobian matrix J is the derivative of these functions with respect to the three parameters (A, \lambda, b). It is given by, J_{ij} = d f_i / d x_j where x_0 = A, x_1 = \lambda and x_2 = b. /* expfit.c -- model functions for exponential + background */ struct data { size_t n; double * y; double * sigma; }; int expb_f (const gsl_vector * x, void *data, gsl_vector * f) { size_t n = ((struct data *)data)->n; double *y = ((struct data *)data)->y; double *sigma = ((struct data *) data)->sigma; double A = gsl_vector_get (x, 0); double lambda = gsl_vector_get (x, 1); double b = gsl_vector_get (x, 2); size_t i; for (i = 0; i < n; i++) { /* Model Yi = A * exp(-lambda * i) + b */ double t = i; double Yi = A * exp (-lambda * t) + b; gsl_vector_set (f, i, (Yi - y[i])/sigma[i]); } return GSL_SUCCESS; } int expb_df (const gsl_vector * x, void *data, gsl_matrix * J) { size_t n = ((struct data *)data)->n; double *sigma = ((struct data *) data)->sigma; double A = gsl_vector_get (x, 0); double lambda = gsl_vector_get (x, 1); size_t i; for (i = 0; i < n; i++) { /* Jacobian matrix J(i,j) = dfi / dxj, */ /* where fi = (Yi - yi)/sigma[i], */ /* Yi = A * exp(-lambda * i) + b */ /* and the xj are the parameters (A,lambda,b) */ double t = i; double s = sigma[i]; double e = exp(-lambda * t); gsl_matrix_set (J, i, 0, e/s); gsl_matrix_set (J, i, 1, -t * A * e/s); gsl_matrix_set (J, i, 2, 1/s); } return GSL_SUCCESS; } int expb_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) { expb_f (x, data, f); expb_df (x, data, J); return GSL_SUCCESS; } The main part of the program sets up a Levenberg-Marquardt solver and some simulated random data. The data uses the known parameters (1.0,5.0,0.1) combined with Gaussian noise (standard deviation = 0.1) over a range of 40 timesteps. The initial guess for the parameters is chosen as (0.0, 1.0, 0.0). #include #include #include #include #include #include #include #include "expfit.c" #define N 40 void print_state (size_t iter, gsl_multifit_fdfsolver * s); int main (void) { const gsl_multifit_fdfsolver_type *T; gsl_multifit_fdfsolver *s; int status; unsigned int i, iter = 0; const size_t n = N; const size_t p = 3; gsl_matrix *covar = gsl_matrix_alloc (p, p); double y[N], sigma[N]; struct data d = { n, y, sigma}; gsl_multifit_function_fdf f; double x_init[3] = { 1.0, 0.0, 0.0 }; gsl_vector_view x = gsl_vector_view_array (x_init, p); const gsl_rng_type * type; gsl_rng * r; gsl_rng_env_setup(); type = gsl_rng_default; r = gsl_rng_alloc (type); f.f = &expb_f; f.df = &expb_df; f.fdf = &expb_fdf; f.n = n; f.p = p; f.params = &d; /* This is the data to be fitted */ for (i = 0; i < n; i++) { double t = i; y[i] = 1.0 + 5 * exp (-0.1 * t) + gsl_ran_gaussian (r, 0.1); sigma[i] = 0.1; printf ("data: %u %g %g\n", i, y[i], sigma[i]); }; T = gsl_multifit_fdfsolver_lmsder; s = gsl_multifit_fdfsolver_alloc (T, n, p); gsl_multifit_fdfsolver_set (s, &f, &x.vector); print_state (iter, s); do { iter++; status = gsl_multifit_fdfsolver_iterate (s); printf ("status = %s\n", gsl_strerror (status)); print_state (iter, s); if (status) break; status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4); } while (status == GSL_CONTINUE && iter < 500); gsl_multifit_covar (s->J, 0.0, covar); #define FIT(i) gsl_vector_get(s->x, i) #define ERR(i) sqrt(gsl_matrix_get(covar,i,i)) { double chi = gsl_blas_dnrm2(s->f); double dof = n - p; double c = GSL_MAX_DBL(1, chi / sqrt(dof)); printf("chisq/dof = %g\n", pow(chi, 2.0) / dof); printf ("A = %.5f +/- %.5f\n", FIT(0), c*ERR(0)); printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1)); printf ("b = %.5f +/- %.5f\n", FIT(2), c*ERR(2)); } printf ("status = %s\n", gsl_strerror (status)); gsl_multifit_fdfsolver_free (s); gsl_matrix_free (covar); gsl_rng_free (r); return 0; } void print_state (size_t iter, gsl_multifit_fdfsolver * s) { printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f " "|f(x)| = %g\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->x, 2), gsl_blas_dnrm2 (s->f)); } The iteration terminates when the change in x is smaller than 0.0001, as both an absolute and relative change. Here are the results of running the program: iter: 0 x=1.00000000 0.00000000 0.00000000 |f(x)|=117.349 status=success iter: 1 x=1.64659312 0.01814772 0.64659312 |f(x)|=76.4578 status=success iter: 2 x=2.85876037 0.08092095 1.44796363 |f(x)|=37.6838 status=success iter: 3 x=4.94899512 0.11942928 1.09457665 |f(x)|=9.58079 status=success iter: 4 x=5.02175572 0.10287787 1.03388354 |f(x)|=5.63049 status=success iter: 5 x=5.04520433 0.10405523 1.01941607 |f(x)|=5.44398 status=success iter: 6 x=5.04535782 0.10404906 1.01924871 |f(x)|=5.44397 chisq/dof = 0.800996 A = 5.04536 +/- 0.06028 lambda = 0.10405 +/- 0.00316 b = 1.01925 +/- 0.03782 status = success The approximate values of the parameters are found correctly, and the chi-squared value indicates a good fit (the chi-squared per degree of freedom is approximately 1). In this case the errors on the parameters can be estimated from the square roots of the diagonal elements of the covariance matrix. If the chi-squared value shows a poor fit (i.e. chi^2/dof >> 1) then the error estimates obtained from the covariance matrix will be too small. In the example program the error estimates are multiplied by \sqrt{\chi^2/dof} in this case, a common way of increasing the errors for a poor fit. Note that a poor fit will result from the use an inappropriate model, and the scaled error estimates may then be outside the range of validity for Gaussian errors.  File: gsl-ref.info, Node: References and Further Reading for Nonlinear Least-Squares Fitting, Prev: Example programs for Nonlinear Least-Squares Fitting, Up: Nonlinear Least-Squares Fitting 38.10 References and Further Reading ==================================== The MINPACK algorithm is described in the following article, J.J. More', `The Levenberg-Marquardt Algorithm: Implementation and Theory', Lecture Notes in Mathematics, v630 (1978), ed G. Watson. The following paper is also relevant to the algorithms described in this section, J.J. More', B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained Optimization Software", ACM Transactions on Mathematical Software, Vol 7, No 1 (1981), p 17-41.  File: gsl-ref.info, Node: Basis Splines, Next: Physical Constants, Prev: Nonlinear Least-Squares Fitting, Up: Top 39 Basis Splines **************** This chapter describes functions for the computation of smoothing basis splines (B-splines). A smoothing spline differs from an interpolating spline in that the resulting curve is not required to pass through each datapoint. *Note Interpolation::, for information about interpolating splines. The header file `gsl_bspline.h' contains the prototypes for the bspline functions and related declarations. * Menu: * Overview of B-splines:: * Initializing the B-splines solver:: * Constructing the knots vector:: * Evaluation of B-spline basis functions:: * Evaluation of B-spline basis function derivatives:: * Obtaining Greville abscissae for B-spline basis functions:: * Example programs for B-splines:: * References and Further Reading::  File: gsl-ref.info, Node: Overview of B-splines, Next: Initializing the B-splines solver, Up: Basis Splines 39.1 Overview ============= B-splines are commonly used as basis functions to fit smoothing curves to large data sets. To do this, the abscissa axis is broken up into some number of intervals, where the endpoints of each interval are called "breakpoints". These breakpoints are then converted to "knots" by imposing various continuity and smoothness conditions at each interface. Given a nondecreasing knot vector t = {t_0, t_1, ..., t_{n+k-1}}, the n basis splines of order k are defined by B_(i,1)(x) = (1, t_i <= x < t_(i+1) (0, else B_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x) + [(t_(i+k) - x)/(t_(i+k) - t_(i+1))] B_(i+1,k-1)(x) for i = 0, ..., n-1. The common case of cubic B-splines is given by k = 4. The above recurrence relation can be evaluated in a numerically stable way by the de Boor algorithm. If we define appropriate knots on an interval [a,b] then the B-spline basis functions form a complete set on that interval. Therefore we can expand a smoothing function as f(x) = \sum_i c_i B_(i,k)(x) given enough (x_j, f(x_j)) data pairs. The coefficients c_i can be readily obtained from a least-squares fit.  File: gsl-ref.info, Node: Initializing the B-splines solver, Next: Constructing the knots vector, Prev: Overview of B-splines, Up: Basis Splines 39.2 Initializing the B-splines solver ====================================== The computation of B-spline functions requires a preallocated workspace of type `gsl_bspline_workspace'. If B-spline derivatives are also required, an additional `gsl_bspline_deriv_workspace' is needed. -- Function: gsl_bspline_workspace * gsl_bspline_alloc (const size_t K, const size_t NBREAK) This function allocates a workspace for computing B-splines of order K. The number of breakpoints is given by NBREAK. This leads to n = nbreak + k - 2 basis functions. Cubic B-splines are specified by k = 4. The size of the workspace is O(5k + nbreak). -- Function: void gsl_bspline_free (gsl_bspline_workspace * W) This function frees the memory associated with the workspace W. -- Function: gsl_bspline_deriv_workspace * gsl_bspline_deriv_alloc (const size_t K) This function allocates a workspace for computing the derivatives of a B-spline basis function of order K. The size of the workspace is O(2k^2). -- Function: void gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * W) This function frees the memory associated with the derivative workspace W.  File: gsl-ref.info, Node: Constructing the knots vector, Next: Evaluation of B-spline basis functions, Prev: Initializing the B-splines solver, Up: Basis Splines 39.3 Constructing the knots vector ================================== -- Function: int gsl_bspline_knots (const gsl_vector * BREAKPTS, gsl_bspline_workspace * W) This function computes the knots associated with the given breakpoints and stores them internally in `w->knots'. -- Function: int gsl_bspline_knots_uniform (const double A, const double B, gsl_bspline_workspace * W) This function assumes uniformly spaced breakpoints on [a,b] and constructs the corresponding knot vector using the previously specified NBREAK parameter. The knots are stored in `w->knots'.  File: gsl-ref.info, Node: Evaluation of B-spline basis functions, Next: Evaluation of B-spline basis function derivatives, Prev: Constructing the knots vector, Up: Basis Splines 39.4 Evaluation of B-splines ============================ -- Function: int gsl_bspline_eval (const double X, gsl_vector * B, gsl_bspline_workspace * W) This function evaluates all B-spline basis functions at the position X and stores them in the vector B, so that the i-th element is B_i(x). The vector B must be of length n = nbreak + k - 2. This value may also be obtained by calling `gsl_bspline_ncoeffs'. Computing all the basis functions at once is more efficient than computing them individually, due to the nature of the defining recurrence relation. -- Function: int gsl_bspline_eval_nonzero (const double X, gsl_vector * BK, size_t * ISTART, size_t * IEND, gsl_bspline_workspace * W) This function evaluates all potentially nonzero B-spline basis functions at the position X and stores them in the vector BK, so that the i-th element is B_(istart+i)(x). The last element of BK is B_(iend)(x). The vector BK must be of length k. By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the B_i(x) to be computed without unnecessary terms (such linear combinations occur, for example, when evaluating an interpolated function). -- Function: size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * W) This function returns the number of B-spline coefficients given by n = nbreak + k - 2.  File: gsl-ref.info, Node: Evaluation of B-spline basis function derivatives, Next: Obtaining Greville abscissae for B-spline basis functions, Prev: Evaluation of B-spline basis functions, Up: Basis Splines 39.5 Evaluation of B-spline derivatives ======================================= -- Function: int gsl_bspline_deriv_eval (const double X, const size_t NDERIV, gsl_matrix * DB, gsl_bspline_workspace * W, gsl_bspline_deriv_workspace * DW) This function evaluates all B-spline basis function derivatives of orders 0 through nderiv (inclusive) at the position X and stores them in the matrix DB. The (i,j)-th element of DB is d^jB_i(x)/dx^j. The matrix DB must be of size n = nbreak + k - 2 by nderiv + 1. The value n may also be obtained by calling `gsl_bspline_ncoeffs'. Note that function evaluations are included as the zeroth order derivatives in DB. Computing all the basis function derivatives at once is more efficient than computing them individually, due to the nature of the defining recurrence relation. -- Function: int gsl_bspline_deriv_eval_nonzero (const double X, const size_t NDERIV, gsl_matrix * DB, size_t * ISTART, size_t * IEND, gsl_bspline_workspace * W, gsl_bspline_deriv_workspace * DW) This function evaluates all potentially nonzero B-spline basis function derivatives of orders 0 through nderiv (inclusive) at the position X and stores them in the matrix DB. The (i,j)-th element of DB is d^j/dx^j B_(istart+i)(x). The last row of DB contains d^j/dx^j B_(iend)(x). The matrix DB must be of size k by at least nderiv + 1. Note that function evaluations are included as the zeroth order derivatives in DB. By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the B_i(x) and their derivatives to be computed without unnecessary terms.  File: gsl-ref.info, Node: Obtaining Greville abscissae for B-spline basis functions, Next: Example programs for B-splines, Prev: Evaluation of B-spline basis function derivatives, Up: Basis Splines 39.6 Greville abscissae ======================= The Greville abscissae are defined to be the mean location of k-1 consecutive knots in the knot vector for each basis spline function of order k. Note that the first and last knots in the knot vector are excluded when applying this definition; consequently there are `gsl_bspline_ncoeffs' Greville abscissa. They are often used in B-spline collocation applications and may also be called Marsden-Schoenberg points. The above definition is undefined for k=1. The implementation chooses to return interval midpoints in the degenerate k=1 case. -- Function: double gsl_bspline_greville_abscissa (size_t I, gsl_bspline_workspace *W); Returns the location of the i-th Greville abscissa for the given spline basis. Here, i = 0, ..., `gsl_bspline_ncoeffs(w) - 1'.  File: gsl-ref.info, Node: Example programs for B-splines, Next: References and Further Reading, Prev: Obtaining Greville abscissae for B-spline basis functions, Up: Basis Splines 39.7 Examples ============= The following program computes a linear least squares fit to data using cubic B-spline basis functions with uniform breakpoints. The data is generated from the curve y(x) = \cos(x) \exp(-x/10) on the interval [0, 15] with Gaussian noise added. #include #include #include #include #include #include #include #include /* number of data points to fit */ #define N 200 /* number of fit coefficients */ #define NCOEFFS 12 /* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */ #define NBREAK (NCOEFFS - 2) int main (void) { const size_t n = N; const size_t ncoeffs = NCOEFFS; const size_t nbreak = NBREAK; size_t i, j; gsl_bspline_workspace *bw; gsl_vector *B; double dy; gsl_rng *r; gsl_vector *c, *w; gsl_vector *x, *y; gsl_matrix *X, *cov; gsl_multifit_linear_workspace *mw; double chisq, Rsq, dof, tss; gsl_rng_env_setup(); r = gsl_rng_alloc(gsl_rng_default); /* allocate a cubic bspline workspace (k = 4) */ bw = gsl_bspline_alloc(4, nbreak); B = gsl_vector_alloc(ncoeffs); x = gsl_vector_alloc(n); y = gsl_vector_alloc(n); X = gsl_matrix_alloc(n, ncoeffs); c = gsl_vector_alloc(ncoeffs); w = gsl_vector_alloc(n); cov = gsl_matrix_alloc(ncoeffs, ncoeffs); mw = gsl_multifit_linear_alloc(n, ncoeffs); printf("#m=0,S=0\n"); /* this is the data to be fitted */ for (i = 0; i < n; ++i) { double sigma; double xi = (15.0 / (N - 1)) * i; double yi = cos(xi) * exp(-0.1 * xi); sigma = 0.1 * yi; dy = gsl_ran_gaussian(r, sigma); yi += dy; gsl_vector_set(x, i, xi); gsl_vector_set(y, i, yi); gsl_vector_set(w, i, 1.0 / (sigma * sigma)); printf("%f %f\n", xi, yi); } /* use uniform breakpoints on [0, 15] */ gsl_bspline_knots_uniform(0.0, 15.0, bw); /* construct the fit matrix X */ for (i = 0; i < n; ++i) { double xi = gsl_vector_get(x, i); /* compute B_j(xi) for all j */ gsl_bspline_eval(xi, B, bw); /* fill in row i of X */ for (j = 0; j < ncoeffs; ++j) { double Bj = gsl_vector_get(B, j); gsl_matrix_set(X, i, j, Bj); } } /* do the fit */ gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw); dof = n - ncoeffs; tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size); Rsq = 1.0 - chisq / tss; fprintf(stderr, "chisq/dof = %e, Rsq = %f\n", chisq / dof, Rsq); /* output the smoothed curve */ { double xi, yi, yerr; printf("#m=1,S=0\n"); for (xi = 0.0; xi < 15.0; xi += 0.1) { gsl_bspline_eval(xi, B, bw); gsl_multifit_linear_est(B, c, cov, &yi, &yerr); printf("%f %f\n", xi, yi); } } gsl_rng_free(r); gsl_bspline_free(bw); gsl_vector_free(B); gsl_vector_free(x); gsl_vector_free(y); gsl_matrix_free(X); gsl_vector_free(c); gsl_vector_free(w); gsl_matrix_free(cov); gsl_multifit_linear_free(mw); return 0; } /* main() */ The output can be plotted with GNU `graph'. $ ./a.out > bspline.dat chisq/dof = 1.118217e+00, Rsq = 0.989771 $ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps  File: gsl-ref.info, Node: References and Further Reading, Prev: Example programs for B-splines, Up: Basis Splines 39.8 References and Further Reading =================================== Further information on the algorithms described in this section can be found in the following book, C. de Boor, `A Practical Guide to Splines' (1978), Springer-Verlag, ISBN 0-387-90356-9. Further information of Greville abscissae and B-spline collocation can be found in the following paper, Richard W. Johnson, Higher order B-spline collocation at the Greville abscissae. `Applied Numerical Mathematics'. vol. 52, 2005, 63-75. A large collection of B-spline routines is available in the PPPACK library available at `http://www.netlib.org/pppack', which is also part of SLATEC.  File: gsl-ref.info, Node: Physical Constants, Next: IEEE floating-point arithmetic, Prev: Basis Splines, Up: Top 40 Physical Constants ********************* This chapter describes macros for the values of physical constants, such as the speed of light, c, and gravitational constant, G. The values are available in different unit systems, including the standard MKSA system (meters, kilograms, seconds, amperes) and the CGSM system (centimeters, grams, seconds, gauss), which is commonly used in Astronomy. The definitions of constants in the MKSA system are available in the file `gsl_const_mksa.h'. The constants in the CGSM system are defined in `gsl_const_cgsm.h'. Dimensionless constants, such as the fine structure constant, which are pure numbers are defined in `gsl_const_num.h'. * Menu: * Fundamental Constants:: * Astronomy and Astrophysics:: * Atomic and Nuclear Physics:: * Measurement of Time:: * Imperial Units :: * Speed and Nautical Units:: * Printers Units:: * Volume Area and Length:: * Mass and Weight :: * Thermal Energy and Power:: * Pressure:: * Viscosity:: * Light and Illumination:: * Radioactivity:: * Force and Energy:: * Prefixes:: * Physical Constant Examples:: * Physical Constant References and Further Reading:: The full list of constants is described briefly below. Consult the header files themselves for the values of the constants used in the library.  File: gsl-ref.info, Node: Fundamental Constants, Next: Astronomy and Astrophysics, Up: Physical Constants 40.1 Fundamental Constants ========================== `GSL_CONST_MKSA_SPEED_OF_LIGHT' The speed of light in vacuum, c. `GSL_CONST_MKSA_VACUUM_PERMEABILITY' The permeability of free space, \mu_0. This constant is defined in the MKSA system only. `GSL_CONST_MKSA_VACUUM_PERMITTIVITY' The permittivity of free space, \epsilon_0. This constant is defined in the MKSA system only. `GSL_CONST_MKSA_PLANCKS_CONSTANT_H' Planck's constant, h. `GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR' Planck's constant divided by 2\pi, \hbar. `GSL_CONST_NUM_AVOGADRO' Avogadro's number, N_a. `GSL_CONST_MKSA_FARADAY' The molar charge of 1 Faraday. `GSL_CONST_MKSA_BOLTZMANN' The Boltzmann constant, k. `GSL_CONST_MKSA_MOLAR_GAS' The molar gas constant, R_0. `GSL_CONST_MKSA_STANDARD_GAS_VOLUME' The standard gas volume, V_0. `GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT' The Stefan-Boltzmann radiation constant, \sigma. `GSL_CONST_MKSA_GAUSS' The magnetic field of 1 Gauss.  File: gsl-ref.info, Node: Astronomy and Astrophysics, Next: Atomic and Nuclear Physics, Prev: Fundamental Constants, Up: Physical Constants 40.2 Astronomy and Astrophysics =============================== `GSL_CONST_MKSA_ASTRONOMICAL_UNIT' The length of 1 astronomical unit (mean earth-sun distance), au. `GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT' The gravitational constant, G. `GSL_CONST_MKSA_LIGHT_YEAR' The distance of 1 light-year, ly. `GSL_CONST_MKSA_PARSEC' The distance of 1 parsec, pc. `GSL_CONST_MKSA_GRAV_ACCEL' The standard gravitational acceleration on Earth, g. `GSL_CONST_MKSA_SOLAR_MASS' The mass of the Sun.  File: gsl-ref.info, Node: Atomic and Nuclear Physics, Next: Measurement of Time, Prev: Astronomy and Astrophysics, Up: Physical Constants 40.3 Atomic and Nuclear Physics =============================== `GSL_CONST_MKSA_ELECTRON_CHARGE' The charge of the electron, e. `GSL_CONST_MKSA_ELECTRON_VOLT' The energy of 1 electron volt, eV. `GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS' The unified atomic mass, amu. `GSL_CONST_MKSA_MASS_ELECTRON' The mass of the electron, m_e. `GSL_CONST_MKSA_MASS_MUON' The mass of the muon, m_\mu. `GSL_CONST_MKSA_MASS_PROTON' The mass of the proton, m_p. `GSL_CONST_MKSA_MASS_NEUTRON' The mass of the neutron, m_n. `GSL_CONST_NUM_FINE_STRUCTURE' The electromagnetic fine structure constant \alpha. `GSL_CONST_MKSA_RYDBERG' The Rydberg constant, Ry, in units of energy. This is related to the Rydberg inverse wavelength R_\infty by Ry = h c R_\infty. `GSL_CONST_MKSA_BOHR_RADIUS' The Bohr radius, a_0. `GSL_CONST_MKSA_ANGSTROM' The length of 1 angstrom. `GSL_CONST_MKSA_BARN' The area of 1 barn. `GSL_CONST_MKSA_BOHR_MAGNETON' The Bohr Magneton, \mu_B. `GSL_CONST_MKSA_NUCLEAR_MAGNETON' The Nuclear Magneton, \mu_N. `GSL_CONST_MKSA_ELECTRON_MAGNETIC_MOMENT' The absolute value of the magnetic moment of the electron, \mu_e. The physical magnetic moment of the electron is negative. `GSL_CONST_MKSA_PROTON_MAGNETIC_MOMENT' The magnetic moment of the proton, \mu_p. `GSL_CONST_MKSA_THOMSON_CROSS_SECTION' The Thomson cross section, \sigma_T. `GSL_CONST_MKSA_DEBYE' The electric dipole moment of 1 Debye, D.  File: gsl-ref.info, Node: Measurement of Time, Next: Imperial Units, Prev: Atomic and Nuclear Physics, Up: Physical Constants 40.4 Measurement of Time ======================== `GSL_CONST_MKSA_MINUTE' The number of seconds in 1 minute. `GSL_CONST_MKSA_HOUR' The number of seconds in 1 hour. `GSL_CONST_MKSA_DAY' The number of seconds in 1 day. `GSL_CONST_MKSA_WEEK' The number of seconds in 1 week.  File: gsl-ref.info, Node: Imperial Units, Next: Speed and Nautical Units, Prev: Measurement of Time, Up: Physical Constants 40.5 Imperial Units =================== `GSL_CONST_MKSA_INCH' The length of 1 inch. `GSL_CONST_MKSA_FOOT' The length of 1 foot. `GSL_CONST_MKSA_YARD' The length of 1 yard. `GSL_CONST_MKSA_MILE' The length of 1 mile. `GSL_CONST_MKSA_MIL' The length of 1 mil (1/1000th of an inch).  File: gsl-ref.info, Node: Speed and Nautical Units, Next: Printers Units, Prev: Imperial Units, Up: Physical Constants 40.6 Speed and Nautical Units ============================= `GSL_CONST_MKSA_KILOMETERS_PER_HOUR' The speed of 1 kilometer per hour. `GSL_CONST_MKSA_MILES_PER_HOUR' The speed of 1 mile per hour. `GSL_CONST_MKSA_NAUTICAL_MILE' The length of 1 nautical mile. `GSL_CONST_MKSA_FATHOM' The length of 1 fathom. `GSL_CONST_MKSA_KNOT' The speed of 1 knot.  File: gsl-ref.info, Node: Printers Units, Next: Volume Area and Length, Prev: Speed and Nautical Units, Up: Physical Constants 40.7 Printers Units =================== `GSL_CONST_MKSA_POINT' The length of 1 printer's point (1/72 inch). `GSL_CONST_MKSA_TEXPOINT' The length of 1 TeX point (1/72.27 inch).  File: gsl-ref.info, Node: Volume Area and Length, Next: Mass and Weight, Prev: Printers Units, Up: Physical Constants 40.8 Volume, Area and Length ============================ `GSL_CONST_MKSA_MICRON' The length of 1 micron. `GSL_CONST_MKSA_HECTARE' The area of 1 hectare. `GSL_CONST_MKSA_ACRE' The area of 1 acre. `GSL_CONST_MKSA_LITER' The volume of 1 liter. `GSL_CONST_MKSA_US_GALLON' The volume of 1 US gallon. `GSL_CONST_MKSA_CANADIAN_GALLON' The volume of 1 Canadian gallon. `GSL_CONST_MKSA_UK_GALLON' The volume of 1 UK gallon. `GSL_CONST_MKSA_QUART' The volume of 1 quart. `GSL_CONST_MKSA_PINT' The volume of 1 pint.  File: gsl-ref.info, Node: Mass and Weight, Next: Thermal Energy and Power, Prev: Volume Area and Length, Up: Physical Constants 40.9 Mass and Weight ==================== `GSL_CONST_MKSA_POUND_MASS' The mass of 1 pound. `GSL_CONST_MKSA_OUNCE_MASS' The mass of 1 ounce. `GSL_CONST_MKSA_TON' The mass of 1 ton. `GSL_CONST_MKSA_METRIC_TON' The mass of 1 metric ton (1000 kg). `GSL_CONST_MKSA_UK_TON' The mass of 1 UK ton. `GSL_CONST_MKSA_TROY_OUNCE' The mass of 1 troy ounce. `GSL_CONST_MKSA_CARAT' The mass of 1 carat. `GSL_CONST_MKSA_GRAM_FORCE' The force of 1 gram weight. `GSL_CONST_MKSA_POUND_FORCE' The force of 1 pound weight. `GSL_CONST_MKSA_KILOPOUND_FORCE' The force of 1 kilopound weight. `GSL_CONST_MKSA_POUNDAL' The force of 1 poundal.  File: gsl-ref.info, Node: Thermal Energy and Power, Next: Pressure, Prev: Mass and Weight, Up: Physical Constants 40.10 Thermal Energy and Power ============================== `GSL_CONST_MKSA_CALORIE' The energy of 1 calorie. `GSL_CONST_MKSA_BTU' The energy of 1 British Thermal Unit, btu. `GSL_CONST_MKSA_THERM' The energy of 1 Therm. `GSL_CONST_MKSA_HORSEPOWER' The power of 1 horsepower.  File: gsl-ref.info, Node: Pressure, Next: Viscosity, Prev: Thermal Energy and Power, Up: Physical Constants 40.11 Pressure ============== `GSL_CONST_MKSA_BAR' The pressure of 1 bar. `GSL_CONST_MKSA_STD_ATMOSPHERE' The pressure of 1 standard atmosphere. `GSL_CONST_MKSA_TORR' The pressure of 1 torr. `GSL_CONST_MKSA_METER_OF_MERCURY' The pressure of 1 meter of mercury. `GSL_CONST_MKSA_INCH_OF_MERCURY' The pressure of 1 inch of mercury. `GSL_CONST_MKSA_INCH_OF_WATER' The pressure of 1 inch of water. `GSL_CONST_MKSA_PSI' The pressure of 1 pound per square inch.  File: gsl-ref.info, Node: Viscosity, Next: Light and Illumination, Prev: Pressure, Up: Physical Constants 40.12 Viscosity =============== `GSL_CONST_MKSA_POISE' The dynamic viscosity of 1 poise. `GSL_CONST_MKSA_STOKES' The kinematic viscosity of 1 stokes.  File: gsl-ref.info, Node: Light and Illumination, Next: Radioactivity, Prev: Viscosity, Up: Physical Constants 40.13 Light and Illumination ============================ `GSL_CONST_MKSA_STILB' The luminance of 1 stilb. `GSL_CONST_MKSA_LUMEN' The luminous flux of 1 lumen. `GSL_CONST_MKSA_LUX' The illuminance of 1 lux. `GSL_CONST_MKSA_PHOT' The illuminance of 1 phot. `GSL_CONST_MKSA_FOOTCANDLE' The illuminance of 1 footcandle. `GSL_CONST_MKSA_LAMBERT' The luminance of 1 lambert. `GSL_CONST_MKSA_FOOTLAMBERT' The luminance of 1 footlambert.  File: gsl-ref.info, Node: Radioactivity, Next: Force and Energy, Prev: Light and Illumination, Up: Physical Constants 40.14 Radioactivity =================== `GSL_CONST_MKSA_CURIE' The activity of 1 curie. `GSL_CONST_MKSA_ROENTGEN' The exposure of 1 roentgen. `GSL_CONST_MKSA_RAD' The absorbed dose of 1 rad.  File: gsl-ref.info, Node: Force and Energy, Next: Prefixes, Prev: Radioactivity, Up: Physical Constants 40.15 Force and Energy ====================== `GSL_CONST_MKSA_NEWTON' The SI unit of force, 1 Newton. `GSL_CONST_MKSA_DYNE' The force of 1 Dyne = 10^-5 Newton. `GSL_CONST_MKSA_JOULE' The SI unit of energy, 1 Joule. `GSL_CONST_MKSA_ERG' The energy 1 erg = 10^-7 Joule.  File: gsl-ref.info, Node: Prefixes, Next: Physical Constant Examples, Prev: Force and Energy, Up: Physical Constants 40.16 Prefixes ============== These constants are dimensionless scaling factors. `GSL_CONST_NUM_YOTTA' 10^24 `GSL_CONST_NUM_ZETTA' 10^21 `GSL_CONST_NUM_EXA' 10^18 `GSL_CONST_NUM_PETA' 10^15 `GSL_CONST_NUM_TERA' 10^12 `GSL_CONST_NUM_GIGA' 10^9 `GSL_CONST_NUM_MEGA' 10^6 `GSL_CONST_NUM_KILO' 10^3 `GSL_CONST_NUM_MILLI' 10^-3 `GSL_CONST_NUM_MICRO' 10^-6 `GSL_CONST_NUM_NANO' 10^-9 `GSL_CONST_NUM_PICO' 10^-12 `GSL_CONST_NUM_FEMTO' 10^-15 `GSL_CONST_NUM_ATTO' 10^-18 `GSL_CONST_NUM_ZEPTO' 10^-21 `GSL_CONST_NUM_YOCTO' 10^-24  File: gsl-ref.info, Node: Physical Constant Examples, Next: Physical Constant References and Further Reading, Prev: Prefixes, Up: Physical Constants 40.17 Examples ============== The following program demonstrates the use of the physical constants in a calculation. In this case, the goal is to calculate the range of light-travel times from Earth to Mars. The required data is the average distance of each planet from the Sun in astronomical units (the eccentricities and inclinations of the orbits will be neglected for the purposes of this calculation). The average radius of the orbit of Mars is 1.52 astronomical units, and for the orbit of Earth it is 1 astronomical unit (by definition). These values are combined with the MKSA values of the constants for the speed of light and the length of an astronomical unit to produce a result for the shortest and longest light-travel times in seconds. The figures are converted into minutes before being displayed. #include #include int main (void) { double c = GSL_CONST_MKSA_SPEED_OF_LIGHT; double au = GSL_CONST_MKSA_ASTRONOMICAL_UNIT; double minutes = GSL_CONST_MKSA_MINUTE; /* distance stored in meters */ double r_earth = 1.00 * au; double r_mars = 1.52 * au; double t_min, t_max; t_min = (r_mars - r_earth) / c; t_max = (r_mars + r_earth) / c; printf ("light travel time from Earth to Mars:\n"); printf ("minimum = %.1f minutes\n", t_min / minutes); printf ("maximum = %.1f minutes\n", t_max / minutes); return 0; } Here is the output from the program, light travel time from Earth to Mars: minimum = 4.3 minutes maximum = 21.0 minutes  File: gsl-ref.info, Node: Physical Constant References and Further Reading, Prev: Physical Constant Examples, Up: Physical Constants 40.18 References and Further Reading ==================================== The authoritative sources for physical constants are the 2006 CODATA recommended values, published in the article below. Further information on the values of physical constants is also available from the NIST website. P.J. Mohr, B.N. Taylor, D.B. Newell, "CODATA Recommended Values of the Fundamental Physical Constants: 2006", Reviews of Modern Physics, 80(2), pp. 633-730 (2008). `http://www.physics.nist.gov/cuu/Constants/index.html' `http://physics.nist.gov/Pubs/SP811/appenB9.html'  File: gsl-ref.info, Node: IEEE floating-point arithmetic, Next: Debugging Numerical Programs, Prev: Physical Constants, Up: Top 41 IEEE floating-point arithmetic ********************************* This chapter describes functions for examining the representation of floating point numbers and controlling the floating point environment of your program. The functions described in this chapter are declared in the header file `gsl_ieee_utils.h'. * Menu: * Representation of floating point numbers:: * Setting up your IEEE environment:: * IEEE References and Further Reading::  File: gsl-ref.info, Node: Representation of floating point numbers, Next: Setting up your IEEE environment, Up: IEEE floating-point arithmetic 41.1 Representation of floating point numbers ============================================= The IEEE Standard for Binary Floating-Point Arithmetic defines binary formats for single and double precision numbers. Each number is composed of three parts: a "sign bit" (s), an "exponent" (E) and a "fraction" (f). The numerical value of the combination (s,E,f) is given by the following formula, (-1)^s (1.fffff...) 2^E The sign bit is either zero or one. The exponent ranges from a minimum value E_min to a maximum value E_max depending on the precision. The exponent is converted to an unsigned number e, known as the "biased exponent", for storage by adding a "bias" parameter, e = E + bias. The sequence fffff... represents the digits of the binary fraction f. The binary digits are stored in "normalized form", by adjusting the exponent to give a leading digit of 1. Since the leading digit is always 1 for normalized numbers it is assumed implicitly and does not have to be stored. Numbers smaller than 2^(E_min) are be stored in "denormalized form" with a leading zero, (-1)^s (0.fffff...) 2^(E_min) This allows gradual underflow down to 2^(E_min - p) for p bits of precision. A zero is encoded with the special exponent of 2^(E_min - 1) and infinities with the exponent of 2^(E_max + 1). The format for single precision numbers uses 32 bits divided in the following way, seeeeeeeefffffffffffffffffffffff s = sign bit, 1 bit e = exponent, 8 bits (E_min=-126, E_max=127, bias=127) f = fraction, 23 bits The format for double precision numbers uses 64 bits divided in the following way, seeeeeeeeeeeffffffffffffffffffffffffffffffffffffffffffffffffffff s = sign bit, 1 bit e = exponent, 11 bits (E_min=-1022, E_max=1023, bias=1023) f = fraction, 52 bits It is often useful to be able to investigate the behavior of a calculation at the bit-level and the library provides functions for printing the IEEE representations in a human-readable form. -- Function: void gsl_ieee_fprintf_float (FILE * STREAM, const float * X) -- Function: void gsl_ieee_fprintf_double (FILE * STREAM, const double * X) These functions output a formatted version of the IEEE floating-point number pointed to by X to the stream STREAM. A pointer is used to pass the number indirectly, to avoid any undesired promotion from `float' to `double'. The output takes one of the following forms, `NaN' the Not-a-Number symbol `Inf, -Inf' positive or negative infinity `1.fffff...*2^E, -1.fffff...*2^E' a normalized floating point number `0.fffff...*2^E, -0.fffff...*2^E' a denormalized floating point number `0, -0' positive or negative zero The output can be used directly in GNU Emacs Calc mode by preceding it with `2#' to indicate binary. -- Function: void gsl_ieee_printf_float (const float * X) -- Function: void gsl_ieee_printf_double (const double * X) These functions output a formatted version of the IEEE floating-point number pointed to by X to the stream `stdout'. The following program demonstrates the use of the functions by printing the single and double precision representations of the fraction 1/3. For comparison the representation of the value promoted from single to double precision is also printed. #include #include int main (void) { float f = 1.0/3.0; double d = 1.0/3.0; double fd = f; /* promote from float to double */ printf (" f="); gsl_ieee_printf_float(&f); printf ("\n"); printf ("fd="); gsl_ieee_printf_double(&fd); printf ("\n"); printf (" d="); gsl_ieee_printf_double(&d); printf ("\n"); return 0; } The binary representation of 1/3 is 0.01010101... . The output below shows that the IEEE format normalizes this fraction to give a leading digit of 1, f= 1.01010101010101010101011*2^-2 fd= 1.0101010101010101010101100000000000000000000000000000*2^-2 d= 1.0101010101010101010101010101010101010101010101010101*2^-2 The output also shows that a single-precision number is promoted to double-precision by adding zeros in the binary representation.  File: gsl-ref.info, Node: Setting up your IEEE environment, Next: IEEE References and Further Reading, Prev: Representation of floating point numbers, Up: IEEE floating-point arithmetic 41.2 Setting up your IEEE environment ===================================== The IEEE standard defines several "modes" for controlling the behavior of floating point operations. These modes specify the important properties of computer arithmetic: the direction used for rounding (e.g. whether numbers should be rounded up, down or to the nearest number), the rounding precision and how the program should handle arithmetic exceptions, such as division by zero. Many of these features can now be controlled via standard functions such as `fpsetround', which should be used whenever they are available. Unfortunately in the past there has been no universal API for controlling their behavior--each system has had its own low-level way of accessing them. To help you write portable programs GSL allows you to specify modes in a platform-independent way using the environment variable `GSL_IEEE_MODE'. The library then takes care of all the necessary machine-specific initializations for you when you call the function `gsl_ieee_env_setup'. -- Function: void gsl_ieee_env_setup () This function reads the environment variable `GSL_IEEE_MODE' and attempts to set up the corresponding specified IEEE modes. The environment variable should be a list of keywords, separated by commas, like this, `GSL_IEEE_MODE' = "KEYWORD,KEYWORD,..." where KEYWORD is one of the following mode-names, `single-precision' `double-precision' `extended-precision' `round-to-nearest' `round-down' `round-up' `round-to-zero' `mask-all' `mask-invalid' `mask-denormalized' `mask-division-by-zero' `mask-overflow' `mask-underflow' `trap-inexact' `trap-common' If `GSL_IEEE_MODE' is empty or undefined then the function returns immediately and no attempt is made to change the system's IEEE mode. When the modes from `GSL_IEEE_MODE' are turned on the function prints a short message showing the new settings to remind you that the results of the program will be affected. If the requested modes are not supported by the platform being used then the function calls the error handler and returns an error code of `GSL_EUNSUP'. When options are specified using this method, the resulting mode is based on a default setting of the highest available precision (double precision or extended precision, depending on the platform) in round-to-nearest mode, with all exceptions enabled apart from the INEXACT exception. The INEXACT exception is generated whenever rounding occurs, so it must generally be disabled in typical scientific calculations. All other floating-point exceptions are enabled by default, including underflows and the use of denormalized numbers, for safety. They can be disabled with the individual `mask-' settings or together using `mask-all'. The following adjusted combination of modes is convenient for many purposes, GSL_IEEE_MODE="double-precision,"\ "mask-underflow,"\ "mask-denormalized" This choice ignores any errors relating to small numbers (either denormalized, or underflowing to zero) but traps overflows, division by zero and invalid operations. Note that on the x86 series of processors this function sets both the original x87 mode and the newer MXCSR mode, which controls SSE floating-point operations. The SSE floating-point units do not have a precision-control bit, and always work in double-precision. The single-precision and extended-precision keywords have no effect in this case. To demonstrate the effects of different rounding modes consider the following program which computes e, the base of natural logarithms, by summing a rapidly-decreasing series, e = 1 + 1/2! + 1/3! + 1/4! + ... = 2.71828182846... #include #include #include int main (void) { double x = 1, oldsum = 0, sum = 0; int i = 0; gsl_ieee_env_setup (); /* read GSL_IEEE_MODE */ do { i++; oldsum = sum; sum += x; x = x / i; printf ("i=%2d sum=%.18f error=%g\n", i, sum, sum - M_E); if (i > 30) break; } while (sum != oldsum); return 0; } Here are the results of running the program in `round-to-nearest' mode. This is the IEEE default so it isn't really necessary to specify it here, $ GSL_IEEE_MODE="round-to-nearest" ./a.out i= 1 sum=1.000000000000000000 error=-1.71828 i= 2 sum=2.000000000000000000 error=-0.718282 .... i=18 sum=2.718281828459045535 error=4.44089e-16 i=19 sum=2.718281828459045535 error=4.44089e-16 After nineteen terms the sum converges to within 4 \times 10^-16 of the correct value. If we now change the rounding mode to `round-down' the final result is less accurate, $ GSL_IEEE_MODE="round-down" ./a.out i= 1 sum=1.000000000000000000 error=-1.71828 .... i=19 sum=2.718281828459041094 error=-3.9968e-15 The result is about 4 \times 10^-15 below the correct value, an order of magnitude worse than the result obtained in the `round-to-nearest' mode. If we change to rounding mode to `round-up' then the final result is higher than the correct value (when we add each term to the sum the final result is always rounded up, which increases the sum by at least one tick until the added term underflows to zero). To avoid this problem we would need to use a safer converge criterion, such as `while (fabs(sum - oldsum) > epsilon)', with a suitably chosen value of epsilon. Finally we can see the effect of computing the sum using single-precision rounding, in the default `round-to-nearest' mode. In this case the program thinks it is still using double precision numbers but the CPU rounds the result of each floating point operation to single-precision accuracy. This simulates the effect of writing the program using single-precision `float' variables instead of `double' variables. The iteration stops after about half the number of iterations and the final result is much less accurate, $ GSL_IEEE_MODE="single-precision" ./a.out .... i=12 sum=2.718281984329223633 error=1.5587e-07 with an error of O(10^-7), which corresponds to single precision accuracy (about 1 part in 10^7). Continuing the iterations further does not decrease the error because all the subsequent results are rounded to the same value.  File: gsl-ref.info, Node: IEEE References and Further Reading, Prev: Setting up your IEEE environment, Up: IEEE floating-point arithmetic 41.3 References and Further Reading =================================== The reference for the IEEE standard is, ANSI/IEEE Std 754-1985, IEEE Standard for Binary Floating-Point Arithmetic. A more pedagogical introduction to the standard can be found in the following paper, David Goldberg: What Every Computer Scientist Should Know About Floating-Point Arithmetic. `ACM Computing Surveys', Vol. 23, No. 1 (March 1991), pages 5-48. Corrigendum: `ACM Computing Surveys', Vol. 23, No. 3 (September 1991), page 413. and see also the sections by B. A. Wichmann and Charles B. Dunham in Surveyor's Forum: "What Every Computer Scientist Should Know About Floating-Point Arithmetic". `ACM Computing Surveys', Vol. 24, No. 3 (September 1992), page 319. A detailed textbook on IEEE arithmetic and its practical use is available from SIAM Press, Michael L. Overton, `Numerical Computing with IEEE Floating Point Arithmetic', SIAM Press, ISBN 0898715717.  File: gsl-ref.info, Node: Debugging Numerical Programs, Next: Contributors to GSL, Prev: IEEE floating-point arithmetic, Up: Top Appendix A Debugging Numerical Programs *************************************** This chapter describes some tips and tricks for debugging numerical programs which use GSL. * Menu: * Using gdb:: * Examining floating point registers:: * Handling floating point exceptions:: * GCC warning options for numerical programs:: * Debugging References::  File: gsl-ref.info, Node: Using gdb, Next: Examining floating point registers, Up: Debugging Numerical Programs A.1 Using gdb ============= Any errors reported by the library are passed to the function `gsl_error'. By running your programs under gdb and setting a breakpoint in this function you can automatically catch any library errors. You can add a breakpoint for every session by putting break gsl_error into your `.gdbinit' file in the directory where your program is started. If the breakpoint catches an error then you can use a backtrace (`bt') to see the call-tree, and the arguments which possibly caused the error. By moving up into the calling function you can investigate the values of variables at that point. Here is an example from the program `fft/test_trap', which contains the following line, status = gsl_fft_complex_wavetable_alloc (0, &complex_wavetable); The function `gsl_fft_complex_wavetable_alloc' takes the length of an FFT as its first argument. When this line is executed an error will be generated because the length of an FFT is not allowed to be zero. To debug this problem we start `gdb', using the file `.gdbinit' to define a breakpoint in `gsl_error', $ gdb test_trap GDB is free software and you are welcome to distribute copies of it under certain conditions; type "show copying" to see the conditions. There is absolutely no warranty for GDB; type "show warranty" for details. GDB 4.16 (i586-debian-linux), Copyright 1996 Free Software Foundation, Inc. Breakpoint 1 at 0x8050b1e: file error.c, line 14. When we run the program this breakpoint catches the error and shows the reason for it. (gdb) run Starting program: test_trap Breakpoint 1, gsl_error (reason=0x8052b0d "length n must be positive integer", file=0x8052b04 "c_init.c", line=108, gsl_errno=1) at error.c:14 14 if (gsl_error_handler) The first argument of `gsl_error' is always a string describing the error. Now we can look at the backtrace to see what caused the problem, (gdb) bt #0 gsl_error (reason=0x8052b0d "length n must be positive integer", file=0x8052b04 "c_init.c", line=108, gsl_errno=1) at error.c:14 #1 0x8049376 in gsl_fft_complex_wavetable_alloc (n=0, wavetable=0xbffff778) at c_init.c:108 #2 0x8048a00 in main (argc=1, argv=0xbffff9bc) at test_trap.c:94 #3 0x80488be in ___crt_dummy__ () We can see that the error was generated in the function `gsl_fft_complex_wavetable_alloc' when it was called with an argument of N=0. The original call came from line 94 in the file `test_trap.c'. By moving up to the level of the original call we can find the line that caused the error, (gdb) up #1 0x8049376 in gsl_fft_complex_wavetable_alloc (n=0, wavetable=0xbffff778) at c_init.c:108 108 GSL_ERROR ("length n must be positive integer", GSL_EDOM); (gdb) up #2 0x8048a00 in main (argc=1, argv=0xbffff9bc) at test_trap.c:94 94 status = gsl_fft_complex_wavetable_alloc (0, &complex_wavetable); Thus we have found the line that caused the problem. From this point we could also print out the values of other variables such as `complex_wavetable'.  File: gsl-ref.info, Node: Examining floating point registers, Next: Handling floating point exceptions, Prev: Using gdb, Up: Debugging Numerical Programs A.2 Examining floating point registers ====================================== The contents of floating point registers can be examined using the command `info float' (on supported platforms). (gdb) info float st0: 0xc4018b895aa17a945000 Valid Normal -7.838871e+308 st1: 0x3ff9ea3f50e4d7275000 Valid Normal 0.0285946 st2: 0x3fe790c64ce27dad4800 Valid Normal 6.7415931e-08 st3: 0x3ffaa3ef0df6607d7800 Spec Normal 0.0400229 st4: 0x3c028000000000000000 Valid Normal 4.4501477e-308 st5: 0x3ffef5412c22219d9000 Zero Normal 0.9580257 st6: 0x3fff8000000000000000 Valid Normal 1 st7: 0xc4028b65a1f6d243c800 Valid Normal -1.566206e+309 fctrl: 0x0272 53 bit; NEAR; mask DENOR UNDER LOS; fstat: 0xb9ba flags 0001; top 7; excep DENOR OVERF UNDER LOS ftag: 0x3fff fip: 0x08048b5c fcs: 0x051a0023 fopoff: 0x08086820 fopsel: 0x002b Individual registers can be examined using the variables $REG, where REG is the register name. (gdb) p $st1 $1 = 0.02859464454261210347719  File: gsl-ref.info, Node: Handling floating point exceptions, Next: GCC warning options for numerical programs, Prev: Examining floating point registers, Up: Debugging Numerical Programs A.3 Handling floating point exceptions ====================================== It is possible to stop the program whenever a `SIGFPE' floating point exception occurs. This can be useful for finding the cause of an unexpected infinity or `NaN'. The current handler settings can be shown with the command `info signal SIGFPE'. (gdb) info signal SIGFPE Signal Stop Print Pass to program Description SIGFPE Yes Yes Yes Arithmetic exception Unless the program uses a signal handler the default setting should be changed so that SIGFPE is not passed to the program, as this would cause it to exit. The command `handle SIGFPE stop nopass' prevents this. (gdb) handle SIGFPE stop nopass Signal Stop Print Pass to program Description SIGFPE Yes Yes No Arithmetic exception Depending on the platform it may be necessary to instruct the kernel to generate signals for floating point exceptions. For programs using GSL this can be achieved using the `GSL_IEEE_MODE' environment variable in conjunction with the function `gsl_ieee_env_setup' as described in *note IEEE floating-point arithmetic::. (gdb) set env GSL_IEEE_MODE=double-precision  File: gsl-ref.info, Node: GCC warning options for numerical programs, Next: Debugging References, Prev: Handling floating point exceptions, Up: Debugging Numerical Programs A.4 GCC warning options for numerical programs ============================================== Writing reliable numerical programs in C requires great care. The following GCC warning options are recommended when compiling numerical programs: gcc -ansi -pedantic -Werror -Wall -W -Wmissing-prototypes -Wstrict-prototypes -Wconversion -Wshadow -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings -Wnested-externs -fshort-enums -fno-common -Dinline= -g -O2 For details of each option consult the manual `Using and Porting GCC'. The following table gives a brief explanation of what types of errors these options catch. `-ansi -pedantic' Use ANSI C, and reject any non-ANSI extensions. These flags help in writing portable programs that will compile on other systems. `-Werror' Consider warnings to be errors, so that compilation stops. This prevents warnings from scrolling off the top of the screen and being lost. You won't be able to compile the program until it is completely warning-free. `-Wall' This turns on a set of warnings for common programming problems. You need `-Wall', but it is not enough on its own. `-O2' Turn on optimization. The warnings for uninitialized variables in `-Wall' rely on the optimizer to analyze the code. If there is no optimization then these warnings aren't generated. `-W' This turns on some extra warnings not included in `-Wall', such as missing return values and comparisons between signed and unsigned integers. `-Wmissing-prototypes -Wstrict-prototypes' Warn if there are any missing or inconsistent prototypes. Without prototypes it is harder to detect problems with incorrect arguments. `-Wconversion' The main use of this option is to warn about conversions from signed to unsigned integers. For example, `unsigned int x = -1'. If you need to perform such a conversion you can use an explicit cast. `-Wshadow' This warns whenever a local variable shadows another local variable. If two variables have the same name then it is a potential source of confusion. `-Wpointer-arith -Wcast-qual -Wcast-align' These options warn if you try to do pointer arithmetic for types which don't have a size, such as `void', if you remove a `const' cast from a pointer, or if you cast a pointer to a type which has a different size, causing an invalid alignment. `-Wwrite-strings' This option gives string constants a `const' qualifier so that it will be a compile-time error to attempt to overwrite them. `-fshort-enums' This option makes the type of `enum' as short as possible. Normally this makes an `enum' different from an `int'. Consequently any attempts to assign a pointer-to-int to a pointer-to-enum will generate a cast-alignment warning. `-fno-common' This option prevents global variables being simultaneously defined in different object files (you get an error at link time). Such a variable should be defined in one file and referred to in other files with an `extern' declaration. `-Wnested-externs' This warns if an `extern' declaration is encountered within a function. `-Dinline=' The `inline' keyword is not part of ANSI C. Thus if you want to use `-ansi' with a program which uses inline functions you can use this preprocessor definition to remove the `inline' keywords. `-g' It always makes sense to put debugging symbols in the executable so that you can debug it using `gdb'. The only effect of debugging symbols is to increase the size of the file, and you can use the `strip' command to remove them later if necessary.  File: gsl-ref.info, Node: Debugging References, Prev: GCC warning options for numerical programs, Up: Debugging Numerical Programs A.5 References and Further Reading ================================== The following books are essential reading for anyone writing and debugging numerical programs with GCC and GDB. R.M. Stallman, `Using and Porting GNU CC', Free Software Foundation, ISBN 1882114388 R.M. Stallman, R.H. Pesch, `Debugging with GDB: The GNU Source-Level Debugger', Free Software Foundation, ISBN 1882114779 For a tutorial introduction to the GNU C Compiler and related programs, see B.J. Gough, `An Introduction to GCC', Network Theory Ltd, ISBN 0954161793  File: gsl-ref.info, Node: Contributors to GSL, Next: Autoconf Macros, Prev: Debugging Numerical Programs, Up: Top Appendix B Contributors to GSL ****************************** (See the AUTHORS file in the distribution for up-to-date information.) *Mark Galassi* Conceived GSL (with James Theiler) and wrote the design document. Wrote the simulated annealing package and the relevant chapter in the manual. *James Theiler* Conceived GSL (with Mark Galassi). Wrote the random number generators and the relevant chapter in this manual. *Jim Davies* Wrote the statistical routines and the relevant chapter in this manual. *Brian Gough* FFTs, numerical integration, random number generators and distributions, root finding, minimization and fitting, polynomial solvers, complex numbers, physical constants, permutations, vector and matrix functions, histograms, statistics, ieee-utils, revised CBLAS Level 2 & 3, matrix decompositions, eigensystems, cumulative distribution functions, testing, documentation and releases. *Reid Priedhorsky* Wrote and documented the initial version of the root finding routines while at Los Alamos National Laboratory, Mathematical Modeling and Analysis Group. *Gerard Jungman* Special Functions, Series acceleration, ODEs, BLAS, Linear Algebra, Eigensystems, Hankel Transforms. *Mike Booth* Wrote the Monte Carlo library. *Jorma Olavi Ta"htinen* Wrote the initial complex arithmetic functions. *Thomas Walter* Wrote the initial heapsort routines and Cholesky decomposition. *Fabrice Rossi* Multidimensional minimization. *Carlo Perassi* Implementation of the random number generators in Knuth's `Seminumerical Algorithms', 3rd Ed. *Szymon Jaroszewicz* Wrote the routines for generating combinations. *Nicolas Darnis* Wrote the cyclic functions and the initial functions for canonical permutations. *Jason H. Stover* Wrote the major cumulative distribution functions. *Ivo Alxneit* Wrote the routines for wavelet transforms. *Tuomo Keskitalo* Improved the implementation of the ODE solvers and wrote the ode-initval2 routines. *Lowell Johnson* Implementation of the Mathieu functions. *Patrick Alken* Implementation of non-symmetric and generalized eigensystems and B-splines. *Rhys Ulerich* Wrote the multiset routines. *Pavel Holoborodko* Wrote the fixed order Gauss-Legendre quadrature routines. *Pedro Gonnet* Wrote the CQUAD integration routines. Thanks to Nigel Lowry for help in proofreading the manual. The non-symmetric eigensystems routines contain code based on the LAPACK linear algebra library. LAPACK is distributed under the following license: Copyright (c) 1992-2006 The University of Tennessee. All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution. * Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.  File: gsl-ref.info, Node: Autoconf Macros, Next: GSL CBLAS Library, Prev: Contributors to GSL, Up: Top Appendix C Autoconf Macros ************************** For applications using `autoconf' the standard macro `AC_CHECK_LIB' can be used to link with GSL automatically from a `configure' script. The library itself depends on the presence of a CBLAS and math library as well, so these must also be located before linking with the main `libgsl' file. The following commands should be placed in the `configure.ac' file to perform these tests, AC_CHECK_LIB([m],[cos]) AC_CHECK_LIB([gslcblas],[cblas_dgemm]) AC_CHECK_LIB([gsl],[gsl_blas_dgemm]) It is important to check for `libm' and `libgslcblas' before `libgsl', otherwise the tests will fail. Assuming the libraries are found the output during the configure stage looks like this, checking for cos in -lm... yes checking for cblas_dgemm in -lgslcblas... yes checking for gsl_blas_dgemm in -lgsl... yes If the library is found then the tests will define the macros `HAVE_LIBGSL', `HAVE_LIBGSLCBLAS', `HAVE_LIBM' and add the options `-lgsl -lgslcblas -lm' to the variable `LIBS'. The tests above will find any version of the library. They are suitable for general use, where the versions of the functions are not important. An alternative macro is available in the file `gsl.m4' to test for a specific version of the library. To use this macro simply add the following line to your `configure.in' file instead of the tests above: AX_PATH_GSL(GSL_VERSION, [action-if-found], [action-if-not-found]) The argument `GSL_VERSION' should be the two or three digit MAJOR.MINOR or MAJOR.MINOR.MICRO version number of the release you require. A suitable choice for `action-if-not-found' is, AC_MSG_ERROR(could not find required version of GSL) Then you can add the variables `GSL_LIBS' and `GSL_CFLAGS' to your Makefile.am files to obtain the correct compiler flags. `GSL_LIBS' is equal to the output of the `gsl-config --libs' command and `GSL_CFLAGS' is equal to `gsl-config --cflags' command. For example, libfoo_la_LDFLAGS = -lfoo $(GSL_LIBS) -lgslcblas Note that the macro `AX_PATH_GSL' needs to use the C compiler so it should appear in the `configure.in' file before the macro `AC_LANG_CPLUSPLUS' for programs that use C++. To test for `inline' the following test should be placed in your `configure.in' file, AC_C_INLINE if test "$ac_cv_c_inline" != no ; then AC_DEFINE(HAVE_INLINE,1) AC_SUBST(HAVE_INLINE) fi and the macro will then be defined in the compilation flags or by including the file `config.h' before any library headers. The following autoconf test will check for `extern inline', dnl Check for "extern inline", using a modified version dnl of the test for AC_C_INLINE from acspecific.mt dnl AC_CACHE_CHECK([for extern inline], ac_cv_c_extern_inline, [ac_cv_c_extern_inline=no AC_TRY_COMPILE([extern $ac_cv_c_inline double foo(double x); extern $ac_cv_c_inline double foo(double x) { return x+1.0; }; double foo (double x) { return x + 1.0; };], [ foo(1.0) ], [ac_cv_c_extern_inline="yes"]) ]) if test "$ac_cv_c_extern_inline" != no ; then AC_DEFINE(HAVE_INLINE,1) AC_SUBST(HAVE_INLINE) fi The substitution of portability functions can be made automatically if you use `autoconf'. For example, to test whether the BSD function `hypot' is available you can include the following line in the configure file `configure.in' for your application, AC_CHECK_FUNCS(hypot) and place the following macro definitions in the file `config.h.in', /* Substitute gsl_hypot for missing system hypot */ #ifndef HAVE_HYPOT #define hypot gsl_hypot #endif The application source files can then use the include command `#include ' to substitute `gsl_hypot' for each occurrence of `hypot' when `hypot' is not available.  File: gsl-ref.info, Node: GSL CBLAS Library, Next: Free Software Needs Free Documentation, Prev: Autoconf Macros, Up: Top Appendix D GSL CBLAS Library **************************** The prototypes for the low-level CBLAS functions are declared in the file `gsl_cblas.h'. For the definition of the functions consult the documentation available from Netlib (*note BLAS References and Further Reading::). * Menu: * Level 1 CBLAS Functions:: * Level 2 CBLAS Functions:: * Level 3 CBLAS Functions:: * GSL CBLAS Examples::  File: gsl-ref.info, Node: Level 1 CBLAS Functions, Next: Level 2 CBLAS Functions, Up: GSL CBLAS Library D.1 Level 1 =========== -- Function: float cblas_sdsdot (const int N, const float ALPHA, const float * X, const int INCX, const float * Y, const int INCY) -- Function: double cblas_dsdot (const int N, const float * X, const int INCX, const float * Y, const int INCY) -- Function: float cblas_sdot (const int N, const float * X, const int INCX, const float * Y, const int INCY) -- Function: double cblas_ddot (const int N, const double * X, const int INCX, const double * Y, const int INCY) -- Function: void cblas_cdotu_sub (const int N, const void * X, const int INCX, const void * Y, const int INCY, void * DOTU) -- Function: void cblas_cdotc_sub (const int N, const void * X, const int INCX, const void * Y, const int INCY, void * DOTC) -- Function: void cblas_zdotu_sub (const int N, const void * X, const int INCX, const void * Y, const int INCY, void * DOTU) -- Function: void cblas_zdotc_sub (const int N, const void * X, const int INCX, const void * Y, const int INCY, void * DOTC) -- Function: float cblas_snrm2 (const int N, const float * X, const int INCX) -- Function: float cblas_sasum (const int N, const float * X, const int INCX) -- Function: double cblas_dnrm2 (const int N, const double * X, const int INCX) -- Function: double cblas_dasum (const int N, const double * X, const int INCX) -- Function: float cblas_scnrm2 (const int N, const void * X, const int INCX) -- Function: float cblas_scasum (const int N, const void * X, const int INCX) -- Function: double cblas_dznrm2 (const int N, const void * X, const int INCX) -- Function: double cblas_dzasum (const int N, const void * X, const int INCX) -- Function: CBLAS_INDEX cblas_isamax (const int N, const float * X, const int INCX) -- Function: CBLAS_INDEX cblas_idamax (const int N, const double * X, const int INCX) -- Function: CBLAS_INDEX cblas_icamax (const int N, const void * X, const int INCX) -- Function: CBLAS_INDEX cblas_izamax (const int N, const void * X, const int INCX) -- Function: void cblas_sswap (const int N, float * X, const int INCX, float * Y, const int INCY) -- Function: void cblas_scopy (const int N, const float * X, const int INCX, float * Y, const int INCY) -- Function: void cblas_saxpy (const int N, const float ALPHA, const float * X, const int INCX, float * Y, const int INCY) -- Function: void cblas_dswap (const int N, double * X, const int INCX, double * Y, const int INCY) -- Function: void cblas_dcopy (const int N, const double * X, const int INCX, double * Y, const int INCY) -- Function: void cblas_daxpy (const int N, const double ALPHA, const double * X, const int INCX, double * Y, const int INCY) -- Function: void cblas_cswap (const int N, void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_ccopy (const int N, const void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_caxpy (const int N, const void * ALPHA, const void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_zswap (const int N, void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_zcopy (const int N, const void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_zaxpy (const int N, const void * ALPHA, const void * X, const int INCX, void * Y, const int INCY) -- Function: void cblas_srotg (float * A, float * B, float * C, float * S) -- Function: void cblas_srotmg (float * D1, float * D2, float * B1, const float B2, float * P) -- Function: void cblas_srot (const int N, float * X, const int INCX, float * Y, const int INCY, const float C, const float S) -- Function: void cblas_srotm (const int N, float * X, const int INCX, float * Y, const int INCY, const float * P) -- Function: void cblas_drotg (double * A, double * B, double * C, double * S) -- Function: void cblas_drotmg (double * D1, double * D2, double * B1, const double B2, double * P) -- Function: void cblas_drot (const int N, double * X, const int INCX, double * Y, const int INCY, const double C, const double S) -- Function: void cblas_drotm (const int N, double * X, const int INCX, double * Y, const int INCY, const double * P) -- Function: void cblas_sscal (const int N, const float ALPHA, float * X, const int INCX) -- Function: void cblas_dscal (const int N, const double ALPHA, double * X, const int INCX) -- Function: void cblas_cscal (const int N, const void * ALPHA, void * X, const int INCX) -- Function: void cblas_zscal (const int N, const void * ALPHA, void * X, const int INCX) -- Function: void cblas_csscal (const int N, const float ALPHA, void * X, const int INCX) -- Function: void cblas_zdscal (const int N, const double ALPHA, void * X, const int INCX)  File: gsl-ref.info, Node: Level 2 CBLAS Functions, Next: Level 3 CBLAS Functions, Prev: Level 1 CBLAS Functions, Up: GSL CBLAS Library D.2 Level 2 =========== -- Function: void cblas_sgemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const float ALPHA, const float * A, const int LDA, const float * X, const int INCX, const float BETA, float * Y, const int INCY) -- Function: void cblas_sgbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const int KL, const int KU, const float ALPHA, const float * A, const int LDA, const float * X, const int INCX, const float BETA, float * Y, const int INCY) -- Function: void cblas_strmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const float * A, const int LDA, float * X, const int INCX) -- Function: void cblas_stbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const float * A, const int LDA, float * X, const int INCX) -- Function: void cblas_stpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const float * AP, float * X, const int INCX) -- Function: void cblas_strsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const float * A, const int LDA, float * X, const int INCX) -- Function: void cblas_stbsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const float * A, const int LDA, float * X, const int INCX) -- Function: void cblas_stpsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const float * AP, float * X, const int INCX) -- Function: void cblas_dgemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const double ALPHA, const double * A, const int LDA, const double * X, const int INCX, const double BETA, double * Y, const int INCY) -- Function: void cblas_dgbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const int KL, const int KU, const double ALPHA, const double * A, const int LDA, const double * X, const int INCX, const double BETA, double * Y, const int INCY) -- Function: void cblas_dtrmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const double * A, const int LDA, double * X, const int INCX) -- Function: void cblas_dtbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const double * A, const int LDA, double * X, const int INCX) -- Function: void cblas_dtpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const double * AP, double * X, const int INCX) -- Function: void cblas_dtrsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const double * A, const int LDA, double * X, const int INCX) -- Function: void cblas_dtbsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const double * A, const int LDA, double * X, const int INCX) -- Function: void cblas_dtpsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const double * AP, double * X, const int INCX) -- Function: void cblas_cgemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_cgbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const int KL, const int KU, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_ctrmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ctbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ctpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * AP, void * X, const int INCX) -- Function: void cblas_ctrsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ctbsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ctpsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * AP, void * X, const int INCX) -- Function: void cblas_zgemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_zgbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const int KL, const int KU, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_ztrmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ztbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ztpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * AP, void * X, const int INCX) -- Function: void cblas_ztrsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ztbsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const int K, const void * A, const int LDA, void * X, const int INCX) -- Function: void cblas_ztpsv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int N, const void * AP, void * X, const int INCX) -- Function: void cblas_ssymv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * A, const int LDA, const float * X, const int INCX, const float BETA, float * Y, const int INCY) -- Function: void cblas_ssbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const int K, const float ALPHA, const float * A, const int LDA, const float * X, const int INCX, const float BETA, float * Y, const int INCY) -- Function: void cblas_sspmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * AP, const float * X, const int INCX, const float BETA, float * Y, const int INCY) -- Function: void cblas_sger (const enum CBLAS_ORDER ORDER, const int M, const int N, const float ALPHA, const float * X, const int INCX, const float * Y, const int INCY, float * A, const int LDA) -- Function: void cblas_ssyr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * X, const int INCX, float * A, const int LDA) -- Function: void cblas_sspr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * X, const int INCX, float * AP) -- Function: void cblas_ssyr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * X, const int INCX, const float * Y, const int INCY, float * A, const int LDA) -- Function: void cblas_sspr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const float * X, const int INCX, const float * Y, const int INCY, float * A) -- Function: void cblas_dsymv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * A, const int LDA, const double * X, const int INCX, const double BETA, double * Y, const int INCY) -- Function: void cblas_dsbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const int K, const double ALPHA, const double * A, const int LDA, const double * X, const int INCX, const double BETA, double * Y, const int INCY) -- Function: void cblas_dspmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * AP, const double * X, const int INCX, const double BETA, double * Y, const int INCY) -- Function: void cblas_dger (const enum CBLAS_ORDER ORDER, const int M, const int N, const double ALPHA, const double * X, const int INCX, const double * Y, const int INCY, double * A, const int LDA) -- Function: void cblas_dsyr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * X, const int INCX, double * A, const int LDA) -- Function: void cblas_dspr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * X, const int INCX, double * AP) -- Function: void cblas_dsyr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * X, const int INCX, const double * Y, const int INCY, double * A, const int LDA) -- Function: void cblas_dspr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const double * X, const int INCX, const double * Y, const int INCY, double * A) -- Function: void cblas_chemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_chbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_chpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * AP, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_cgeru (const enum CBLAS_ORDER ORDER, const int M, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_cgerc (const enum CBLAS_ORDER ORDER, const int M, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_cher (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const void * X, const int INCX, void * A, const int LDA) -- Function: void cblas_chpr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const void * X, const int INCX, void * A) -- Function: void cblas_cher2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_chpr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * AP) -- Function: void cblas_zhemv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_zhbmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_zhpmv (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * AP, const void * X, const int INCX, const void * BETA, void * Y, const int INCY) -- Function: void cblas_zgeru (const enum CBLAS_ORDER ORDER, const int M, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_zgerc (const enum CBLAS_ORDER ORDER, const int M, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_zher (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const void * X, const int INCX, void * A, const int LDA) -- Function: void cblas_zhpr (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const void * X, const int INCX, void * A) -- Function: void cblas_zher2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * A, const int LDA) -- Function: void cblas_zhpr2 (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const void * X, const int INCX, const void * Y, const int INCY, void * AP)  File: gsl-ref.info, Node: Level 3 CBLAS Functions, Next: GSL CBLAS Examples, Prev: Level 2 CBLAS Functions, Up: GSL CBLAS Library D.3 Level 3 =========== -- Function: void cblas_sgemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE TRANSB, const int M, const int N, const int K, const float ALPHA, const float * A, const int LDA, const float * B, const int LDB, const float BETA, float * C, const int LDC) -- Function: void cblas_ssymm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const float ALPHA, const float * A, const int LDA, const float * B, const int LDB, const float BETA, float * C, const int LDC) -- Function: void cblas_ssyrk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const float ALPHA, const float * A, const int LDA, const float BETA, float * C, const int LDC) -- Function: void cblas_ssyr2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const float ALPHA, const float * A, const int LDA, const float * B, const int LDB, const float BETA, float * C, const int LDC) -- Function: void cblas_strmm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const float ALPHA, const float * A, const int LDA, float * B, const int LDB) -- Function: void cblas_strsm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const float ALPHA, const float * A, const int LDA, float * B, const int LDB) -- Function: void cblas_dgemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE TRANSB, const int M, const int N, const int K, const double ALPHA, const double * A, const int LDA, const double * B, const int LDB, const double BETA, double * C, const int LDC) -- Function: void cblas_dsymm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const double ALPHA, const double * A, const int LDA, const double * B, const int LDB, const double BETA, double * C, const int LDC) -- Function: void cblas_dsyrk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const double ALPHA, const double * A, const int LDA, const double BETA, double * C, const int LDC) -- Function: void cblas_dsyr2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const double ALPHA, const double * A, const int LDA, const double * B, const int LDB, const double BETA, double * C, const int LDC) -- Function: void cblas_dtrmm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const double ALPHA, const double * A, const int LDA, double * B, const int LDB) -- Function: void cblas_dtrsm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const double ALPHA, const double * A, const int LDA, double * B, const int LDB) -- Function: void cblas_cgemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE TRANSB, const int M, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_csymm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_csyrk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * BETA, void * C, const int LDC) -- Function: void cblas_csyr2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_ctrmm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const void * ALPHA, const void * A, const int LDA, void * B, const int LDB) -- Function: void cblas_ctrsm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const void * ALPHA, const void * A, const int LDA, void * B, const int LDB) -- Function: void cblas_zgemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE TRANSB, const int M, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_zsymm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_zsyrk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * BETA, void * C, const int LDC) -- Function: void cblas_zsyr2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_ztrmm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const void * ALPHA, const void * A, const int LDA, void * B, const int LDB) -- Function: void cblas_ztrsm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int M, const int N, const void * ALPHA, const void * A, const int LDA, void * B, const int LDB) -- Function: void cblas_chemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_cherk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const float ALPHA, const void * A, const int LDA, const float BETA, void * C, const int LDC) -- Function: void cblas_cher2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const float BETA, void * C, const int LDC) -- Function: void cblas_zhemm (const enum CBLAS_ORDER ORDER, const enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int M, const int N, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const void * BETA, void * C, const int LDC) -- Function: void cblas_zherk (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const double ALPHA, const void * A, const int LDA, const double BETA, void * C, const int LDC) -- Function: void cblas_zher2k (const enum CBLAS_ORDER ORDER, const enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const int N, const int K, const void * ALPHA, const void * A, const int LDA, const void * B, const int LDB, const double BETA, void * C, const int LDC) -- Function: void cblas_xerbla (int P, const char * ROUT, const char * FORM, ...)  File: gsl-ref.info, Node: GSL CBLAS Examples, Prev: Level 3 CBLAS Functions, Up: GSL CBLAS Library D.4 Examples ============ The following program computes the product of two matrices using the Level-3 BLAS function SGEMM, [ 0.11 0.12 0.13 ] [ 1011 1012 ] [ 367.76 368.12 ] [ 0.21 0.22 0.23 ] [ 1021 1022 ] = [ 674.06 674.72 ] [ 1031 1032 ] The matrices are stored in row major order but could be stored in column major order if the first argument of the call to `cblas_sgemm' was changed to `CblasColMajor'. #include #include int main (void) { int lda = 3; float A[] = { 0.11, 0.12, 0.13, 0.21, 0.22, 0.23 }; int ldb = 2; float B[] = { 1011, 1012, 1021, 1022, 1031, 1032 }; int ldc = 2; float C[] = { 0.00, 0.00, 0.00, 0.00 }; /* Compute C = A B */ cblas_sgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans, 2, 2, 3, 1.0, A, lda, B, ldb, 0.0, C, ldc); printf ("[ %g, %g\n", C[0], C[1]); printf (" %g, %g ]\n", C[2], C[3]); return 0; } To compile the program use the following command line, $ gcc -Wall demo.c -lgslcblas There is no need to link with the main library `-lgsl' in this case as the CBLAS library is an independent unit. Here is the output from the program, $ ./a.out [ 367.76, 368.12 674.06, 674.72 ]  File: gsl-ref.info, Node: Free Software Needs Free Documentation, Next: GNU General Public License, Prev: GSL CBLAS Library, Up: Top Free Software Needs Free Documentation ************************************** The following article was written by Richard Stallman, founder of the GNU Project. The biggest deficiency in the free software community today is not in the software--it is the lack of good free documentation that we can include with the free software. Many of our most important programs do not come with free reference manuals and free introductory texts. Documentation is an essential part of any software package; when an important free software package does not come with a free manual and a free tutorial, that is a major gap. We have many such gaps today. Consider Perl, for instance. The tutorial manuals that people normally use are non-free. How did this come about? Because the authors of those manuals published them with restrictive terms--no copying, no modification, source files not available--which exclude them from the free software world. That wasn't the first time this sort of thing happened, and it was far from the last. Many times we have heard a GNU user eagerly describe a manual that he is writing, his intended contribution to the community, only to learn that he had ruined everything by signing a publication contract to make it non-free. Free documentation, like free software, is a matter of freedom, not price. The problem with the non-free manual is not that publishers charge a price for printed copies--that in itself is fine. (The Free Software Foundation sells printed copies of manuals, too.) The problem is the restrictions on the use of the manual. Free manuals are available in source code form, and give you permission to copy and modify. Non-free manuals do not allow this. The criteria of freedom for a free manual are roughly the same as for free software. Redistribution (including the normal kinds of commercial redistribution) must be permitted, so that the manual can accompany every copy of the program, both on-line and on paper. Permission for modification of the technical content is crucial too. When people modify the software, adding or changing features, if they are conscientious they will change the manual too--so they can provide accurate and clear documentation for the modified program. A manual that leaves you no choice but to write a new manual to document a changed version of the program is not really available to our community. Some kinds of limits on the way modification is handled are acceptable. For example, requirements to preserve the original author's copyright notice, the distribution terms, or the list of authors, are ok. It is also no problem to require modified versions to include notice that they were modified. Even entire sections that may not be deleted or changed are acceptable, as long as they deal with nontechnical topics (like this one). These kinds of restrictions are acceptable because they don't obstruct the community's normal use of the manual. However, it must be possible to modify all the _technical_ content of the manual, and then distribute the result in all the usual media, through all the usual channels. Otherwise, the restrictions obstruct the use of the manual, it is not free, and we need another manual to replace it. Please spread the word about this issue. Our community continues to lose manuals to proprietary publishing. If we spread the word that free software needs free reference manuals and free tutorials, perhaps the next person who wants to contribute by writing documentation will realize, before it is too late, that only free manuals contribute to the free software community. If you are writing documentation, please insist on publishing it under the GNU Free Documentation License or another free documentation license. Remember that this decision requires your approval--you don't have to let the publisher decide. Some commercial publishers will use a free license if you insist, but they will not propose the option; it is up to you to raise the issue and say firmly that this is what you want. If the publisher you are dealing with refuses, please try other publishers. If you're not sure whether a proposed license is free, write to . You can encourage commercial publishers to sell more free, copylefted manuals and tutorials by buying them, and particularly by buying copies from the publishers that paid for their writing or for major improvements. Meanwhile, try to avoid buying non-free documentation at all. Check the distribution terms of a manual before you buy it, and insist that whoever seeks your business must respect your freedom. Check the history of the book, and try reward the publishers that have paid or pay the authors to work on it. The Free Software Foundation maintains a list of free documentation published by other publishers: `http://www.fsf.org/doc/other-free-books.html'  File: gsl-ref.info, Node: GNU General Public License, Next: GNU Free Documentation License, Prev: Free Software Needs Free Documentation, Up: Top GNU General Public License ************************** Version 3, 29 June 2007 Copyright (C) 2007 Free Software Foundation, Inc. `http://fsf.org/' Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. Preamble ======== The GNU General Public License is a free, copyleft license for software and other kinds of works. The licenses for most software and other practical works are designed to take away your freedom to share and change the works. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program-to make sure it remains free software for all its users. We, the Free Software Foundation, use the GNU General Public License for most of our software; it applies also to any other work released this way by its authors. You can apply it to your programs, too. When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for them if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs, and that you know you can do these things. To protect your rights, we need to prevent others from denying you these rights or asking you to surrender the rights. Therefore, you have certain responsibilities if you distribute copies of the software, or if you modify it: responsibilities to respect the freedom of others. For example, if you distribute copies of such a program, whether gratis or for a fee, you must pass on to the recipients the same freedoms that you received. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights. Developers that use the GNU GPL protect your rights with two steps: (1) assert copyright on the software, and (2) offer you this License giving you legal permission to copy, distribute and/or modify it. For the developers' and authors' protection, the GPL clearly explains that there is no warranty for this free software. For both users' and authors' sake, the GPL requires that modified versions be marked as changed, so that their problems will not be attributed erroneously to authors of previous versions. Some devices are designed to deny users access to install or run modified versions of the software inside them, although the manufacturer can do so. This is fundamentally incompatible with the aim of protecting users' freedom to change the software. The systematic pattern of such abuse occurs in the area of products for individuals to use, which is precisely where it is most unacceptable. Therefore, we have designed this version of the GPL to prohibit the practice for those products. If such problems arise substantially in other domains, we stand ready to extend this provision to those domains in future versions of the GPL, as needed to protect the freedom of users. Finally, every program is threatened constantly by software patents. States should not allow patents to restrict development and use of software on general-purpose computers, but in those that do, we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary. To prevent this, the GPL assures that patents cannot be used to render the program non-free. The precise terms and conditions for copying, distribution and modification follow. TERMS AND CONDITIONS ==================== 0. Definitions. "This License" refers to version 3 of the GNU General Public License. "Copyright" also means copyright-like laws that apply to other kinds of works, such as semiconductor masks. "The Program" refers to any copyrightable work licensed under this License. Each licensee is addressed as "you". "Licensees" and "recipients" may be individuals or organizations. To "modify" a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission, other than the making of an exact copy. The resulting work is called a "modified version" of the earlier work or a work "based on" the earlier work. A "covered work" means either the unmodified Program or a work based on the Program. To "propagate" a work means to do anything with it that, without permission, would make you directly or secondarily liable for infringement under applicable copyright law, except executing it on a computer or modifying a private copy. Propagation includes copying, distribution (with or without modification), making available to the public, and in some countries other activities as well. To "convey" a work means any kind of propagation that enables other parties to make or receive copies. Mere interaction with a user through a computer network, with no transfer of a copy, is not conveying. An interactive user interface displays "Appropriate Legal Notices" to the extent that it includes a convenient and prominently visible feature that (1) displays an appropriate copyright notice, and (2) tells the user that there is no warranty for the work (except to the extent that warranties are provided), that licensees may convey the work under this License, and how to view a copy of this License. If the interface presents a list of user commands or options, such as a menu, a prominent item in the list meets this criterion. 1. Source Code. The "source code" for a work means the preferred form of the work for making modifications to it. "Object code" means any non-source form of a work. A "Standard Interface" means an interface that either is an official standard defined by a recognized standards body, or, in the case of interfaces specified for a particular programming language, one that is widely used among developers working in that language. The "System Libraries" of an executable work include anything, other than the work as a whole, that (a) is included in the normal form of packaging a Major Component, but which is not part of that Major Component, and (b) serves only to enable use of the work with that Major Component, or to implement a Standard Interface for which an implementation is available to the public in source code form. A "Major Component", in this context, means a major essential component (kernel, window system, and so on) of the specific operating system (if any) on which the executable work runs, or a compiler used to produce the work, or an object code interpreter used to run it. The "Corresponding Source" for a work in object code form means all the source code needed to generate, install, and (for an executable work) run the object code and to modify the work, including scripts to control those activities. However, it does not include the work's System Libraries, or general-purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work. For example, Corresponding Source includes interface definition files associated with source files for the work, and the source code for shared libraries and dynamically linked subprograms that the work is specifically designed to require, such as by intimate data communication or control flow between those subprograms and other parts of the work. The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source. The Corresponding Source for a work in source code form is that same work. 2. Basic Permissions. All rights granted under this License are granted for the term of copyright on the Program, and are irrevocable provided the stated conditions are met. This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work. This License acknowledges your rights of fair use or other equivalent, as provided by copyright law. You may make, run and propagate covered works that you do not convey, without conditions so long as your license otherwise remains in force. You may convey covered works to others for the sole purpose of having them make modifications exclusively for you, or provide you with facilities for running those works, provided that you comply with the terms of this License in conveying all material for which you do not control copyright. Those thus making or running the covered works for you must do so exclusively on your behalf, under your direction and control, on terms that prohibit them from making any copies of your copyrighted material outside their relationship with you. Conveying under any other circumstances is permitted solely under the conditions stated below. Sublicensing is not allowed; section 10 makes it unnecessary. 3. Protecting Users' Legal Rights From Anti-Circumvention Law. No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention of such measures. When you convey a covered work, you waive any legal power to forbid circumvention of technological measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work, and you disclaim any intention to limit operation or modification of the work as a means of enforcing, against the work's users, your or third parties' legal rights to forbid circumvention of technological measures. 4. Conveying Verbatim Copies. You may convey verbatim copies of the Program's source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice; keep intact all notices stating that this License and any non-permissive terms added in accord with section 7 apply to the code; keep intact all notices of the absence of any warranty; and give all recipients a copy of this License along with the Program. You may charge any price or no price for each copy that you convey, and you may offer support or warranty protection for a fee. 5. Conveying Modified Source Versions. You may convey a work based on the Program, or the modifications to produce it from the Program, in the form of source code under the terms of section 4, provided that you also meet all of these conditions: a. The work must carry prominent notices stating that you modified it, and giving a relevant date. b. The work must carry prominent notices stating that it is released under this License and any conditions added under section 7. This requirement modifies the requirement in section 4 to "keep intact all notices". c. You must license the entire work, as a whole, under this License to anyone who comes into possession of a copy. This License will therefore apply, along with any applicable section 7 additional terms, to the whole of the work, and all its parts, regardless of how they are packaged. This License gives no permission to license the work in any other way, but it does not invalidate such permission if you have separately received it. d. If the work has interactive user interfaces, each must display Appropriate Legal Notices; however, if the Program has interactive interfaces that do not display Appropriate Legal Notices, your work need not make them do so. A compilation of a covered work with other separate and independent works, which are not by their nature extensions of the covered work, and which are not combined with it such as to form a larger program, in or on a volume of a storage or distribution medium, is called an "aggregate" if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation's users beyond what the individual works permit. Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate. 6. Conveying Non-Source Forms. You may convey a covered work in object code form under the terms of sections 4 and 5, provided that you also convey the machine-readable Corresponding Source under the terms of this License, in one of these ways: a. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by the Corresponding Source fixed on a durable physical medium customarily used for software interchange. b. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by a written offer, valid for at least three years and valid for as long as you offer spare parts or customer support for that product model, to give anyone who possesses the object code either (1) a copy of the Corresponding Source for all the software in the product that is covered by this License, on a durable physical medium customarily used for software interchange, for a price no more than your reasonable cost of physically performing this conveying of source, or (2) access to copy the Corresponding Source from a network server at no charge. c. Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source. This alternative is allowed only occasionally and noncommercially, and only if you received the object code with such an offer, in accord with subsection 6b. d. Convey the object code by offering access from a designated place (gratis or for a charge), and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge. You need not require recipients to copy the Corresponding Source along with the object code. If the place to copy the object code is a network server, the Corresponding Source may be on a different server (operated by you or a third party) that supports equivalent copying facilities, provided you maintain clear directions next to the object code saying where to find the Corresponding Source. Regardless of what server hosts the Corresponding Source, you remain obligated to ensure that it is available for as long as needed to satisfy these requirements. e. Convey the object code using peer-to-peer transmission, provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d. A separable portion of the object code, whose source code is excluded from the Corresponding Source as a System Library, need not be included in conveying the object code work. A "User Product" is either (1) a "consumer product", which means any tangible personal property which is normally used for personal, family, or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, "normally used" refers to a typical or common use of that class of product, regardless of the status of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial commercial, industrial or non-consumer uses, unless such uses represent the only significant mode of use of the product. "Installation Information" for a User Product means any methods, procedures, authorization keys, or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source. The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been made. If you convey an object code work under this section in, or with, or specifically for use in, a User Product, and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term (regardless of how the transaction is characterized), the Corresponding Source conveyed under this section must be accompanied by the Installation Information. But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product (for example, the work has been installed in ROM). The requirement to provide Installation Information does not include a requirement to continue to provide support service, warranty, or updates for a work that has been modified or installed by the recipient, or for the User Product in which it has been modified or installed. Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network. Corresponding Source conveyed, and Installation Information provided, in accord with this section must be in a format that is publicly documented (and with an implementation available to the public in source code form), and must require no special password or key for unpacking, reading or copying. 7. Additional Terms. "Additional permissions" are terms that supplement the terms of this License by making exceptions from one or more of its conditions. Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License, to the extent that they are valid under applicable law. If additional permissions apply only to part of the Program, that part may be used separately under those permissions, but the entire Program remains governed by this License without regard to the additional permissions. When you convey a copy of a covered work, you may at your option remove any additional permissions from that copy, or from any part of it. (Additional permissions may be written to require their own removal in certain cases when you modify the work.) You may place additional permissions on material, added by you to a covered work, for which you have or can give appropriate copyright permission. Notwithstanding any other provision of this License, for material you add to a covered work, you may (if authorized by the copyright holders of that material) supplement the terms of this License with terms: a. Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License; or b. Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it; or c. Prohibiting misrepresentation of the origin of that material, or requiring that modified versions of such material be marked in reasonable ways as different from the original version; or d. Limiting the use for publicity purposes of names of licensors or authors of the material; or e. Declining to grant rights under trademark law for use of some trade names, trademarks, or service marks; or f. Requiring indemnification of licensors and authors of that material by anyone who conveys the material (or modified versions of it) with contractual assumptions of liability to the recipient, for any liability that these contractual assumptions directly impose on those licensors and authors. All other non-permissive additional terms are considered "further restrictions" within the meaning of section 10. If the Program as you received it, or any part of it, contains a notice stating that it is governed by this License along with a term that is a further restriction, you may remove that term. If a license document contains a further restriction but permits relicensing or conveying under this License, you may add to a covered work material governed by the terms of that license document, provided that the further restriction does not survive such relicensing or conveying. If you add terms to a covered work in accord with this section, you must place, in the relevant source files, a statement of the additional terms that apply to those files, or a notice indicating where to find the applicable terms. Additional terms, permissive or non-permissive, may be stated in the form of a separately written license, or stated as exceptions; the above requirements apply either way. 8. Termination. You may not propagate or modify a covered work except as expressly provided under this License. Any attempt otherwise to propagate or modify it is void, and will automatically terminate your rights under this License (including any patent licenses granted under the third paragraph of section 11). However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation. Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice. Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, you do not qualify to receive new licenses for the same material under section 10. 9. Acceptance Not Required for Having Copies. You are not required to accept this License in order to receive or run a copy of the Program. Ancillary propagation of a covered work occurring solely as a consequence of using peer-to-peer transmission to receive a copy likewise does not require acceptance. However, nothing other than this License grants you permission to propagate or modify any covered work. These actions infringe copyright if you do not accept this License. Therefore, by modifying or propagating a covered work, you indicate your acceptance of this License to do so. 10. Automatic Licensing of Downstream Recipients. Each time you convey a covered work, the recipient automatically receives a license from the original licensors, to run, modify and propagate that work, subject to this License. You are not responsible for enforcing compliance by third parties with this License. An "entity transaction" is a transaction transferring control of an organization, or substantially all assets of one, or subdividing an organization, or merging organizations. If propagation of a covered work results from an entity transaction, each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party's predecessor in interest had or could give under the previous paragraph, plus a right to possession of the Corresponding Source of the work from the predecessor in interest, if the predecessor has it or can get it with reasonable efforts. You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License. For example, you may not impose a license fee, royalty, or other charge for exercise of rights granted under this License, and you may not initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging that any patent claim is infringed by making, using, selling, offering for sale, or importing the Program or any portion of it. 11. Patents. A "contributor" is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based. The work thus licensed is called the contributor's "contributor version". A contributor's "essential patent claims" are all patent claims owned or controlled by the contributor, whether already acquired or hereafter acquired, that would be infringed by some manner, permitted by this License, of making, using, or selling its contributor version, but do not include claims that would be infringed only as a consequence of further modification of the contributor version. For purposes of this definition, "control" includes the right to grant patent sublicenses in a manner consistent with the requirements of this License. Each contributor grants you a non-exclusive, worldwide, royalty-free patent license under the contributor's essential patent claims, to make, use, sell, offer for sale, import and otherwise run, modify and propagate the contents of its contributor version. In the following three paragraphs, a "patent license" is any express agreement or commitment, however denominated, not to enforce a patent (such as an express permission to practice a patent or covenant not to sue for patent infringement). To "grant" such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party. If you convey a covered work, knowingly relying on a patent license, and the Corresponding Source of the work is not available for anyone to copy, free of charge and under the terms of this License, through a publicly available network server or other readily accessible means, then you must either (1) cause the Corresponding Source to be so available, or (2) arrange to deprive yourself of the benefit of the patent license for this particular work, or (3) arrange, in a manner consistent with the requirements of this License, to extend the patent license to downstream recipients. "Knowingly relying" means you have actual knowledge that, but for the patent license, your conveying the covered work in a country, or your recipient's use of the covered work in a country, would infringe one or more identifiable patents in that country that you have reason to believe are valid. If, pursuant to or in connection with a single transaction or arrangement, you convey, or propagate by procuring conveyance of, a covered work, and grant a patent license to some of the parties receiving the covered work authorizing them to use, propagate, modify or convey a specific copy of the covered work, then the patent license you grant is automatically extended to all recipients of the covered work and works based on it. A patent license is "discriminatory" if it does not include within the scope of its coverage, prohibits the exercise of, or is conditioned on the non-exercise of one or more of the rights that are specifically granted under this License. You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software, under which you make payment to the third party based on the extent of your activity of conveying the work, and under which the third party grants, to any of the parties who would receive the covered work from you, a discriminatory patent license (a) in connection with copies of the covered work conveyed by you (or copies made from those copies), or (b) primarily for and in connection with specific products or compilations that contain the covered work, unless you entered into that arrangement, or that patent license was granted, prior to 28 March 2007. Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you under applicable patent law. 12. No Surrender of Others' Freedom. If conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not convey it at all. For example, if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program, the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program. 13. Use with the GNU Affero General Public License. Notwithstanding any other provision of this License, you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work, and to convey the resulting work. The terms of this License will continue to apply to the part which is the covered work, but the special requirements of the GNU Affero General Public License, section 13, concerning interaction through a network will apply to the combination as such. 14. Revised Versions of this License. The Free Software Foundation may publish revised and/or new versions of the GNU General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. Each version is given a distinguishing version number. If the Program specifies that a certain numbered version of the GNU General Public License "or any later version" applies to it, you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of the GNU General Public License, you may choose any version ever published by the Free Software Foundation. If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used, that proxy's public statement of acceptance of a version permanently authorizes you to choose that version for the Program. Later license versions may give you additional or different permissions. However, no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version. 15. Disclaimer of Warranty. THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. 16. Limitation of Liability. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. 17. Interpretation of Sections 15 and 16. If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms, reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee. END OF TERMS AND CONDITIONS =========================== How to Apply These Terms to Your New Programs ============================================= If you develop a new program, and you want it to be of the greatest possible use to the public, the best way to achieve this is to make it free software which everyone can redistribute and change under these terms. To do so, attach the following notices to the program. It is safest to attach them to the start of each source file to most effectively state the exclusion of warranty; and each file should have at least the "copyright" line and a pointer to where the full notice is found. one line to give the program's name and a brief idea of whAT IT DOES. Copyright (C) YEAR NAME OF AUTHOR 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 `http://www.gnu.org/licenses/'. Also add information on how to contact you by electronic and paper mail. If the program does terminal interaction, make it output a short notice like this when it starts in an interactive mode: PROGRAM Copyright (C) YEAR NAME OF AUTHOR This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. This is free software, and you are welcome to redistribute it under certain conditions; type `show c' for details. The hypothetical commands `show w' and `show c' should show the appropriate parts of the General Public License. Of course, your program's commands might be different; for a GUI interface, you would use an "about box". You should also get your employer (if you work as a programmer) or school, if any, to sign a "copyright disclaimer" for the program, if necessary. For more information on this, and how to apply and follow the GNU GPL, see `http://www.gnu.org/licenses/'. The GNU General Public License does not permit incorporating your program into proprietary programs. If your program is a subroutine library, you may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read `http://www.gnu.org/philosophy/why-not-lgpl.html'.  File: gsl-ref.info, Node: GNU Free Documentation License, Next: Function Index, Prev: GNU General Public License, Up: Top GNU Free Documentation License ****************************** Version 1.3, 3 November 2008 Copyright (C) 2000, 2001, 2002, 2007, 2008 Free Software Foundation, Inc. `http://fsf.org/' Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. 0. PREAMBLE The purpose of this License is to make a manual, textbook, or other functional and useful document "free" in the sense of freedom: to assure everyone the effective freedom to copy and redistribute it, with or without modifying it, either commercially or noncommercially. Secondarily, this License preserves for the author and publisher a way to get credit for their work, while not being considered responsible for modifications made by others. This License is a kind of "copyleft", which means that derivative works of the document must themselves be free in the same sense. It complements the GNU General Public License, which is a copyleft license designed for free software. We have designed this License in order to use it for manuals for free software, because free software needs free documentation: a free program should come with manuals providing the same freedoms that the software does. But this License is not limited to software manuals; it can be used for any textual work, regardless of subject matter or whether it is published as a printed book. We recommend this License principally for works whose purpose is instruction or reference. 1. APPLICABILITY AND DEFINITIONS This License applies to any manual or other work, in any medium, that contains a notice placed by the copyright holder saying it can be distributed under the terms of this License. Such a notice grants a world-wide, royalty-free license, unlimited in duration, to use that work under the conditions stated herein. The "Document", below, refers to any such manual or work. Any member of the public is a licensee, and is addressed as "you". You accept the license if you copy, modify or distribute the work in a way requiring permission under copyright law. A "Modified Version" of the Document means any work containing the Document or a portion of it, either copied verbatim, or with modifications and/or translated into another language. A "Secondary Section" is a named appendix or a front-matter section of the Document that deals exclusively with the relationship of the publishers or authors of the Document to the Document's overall subject (or to related matters) and contains nothing that could fall directly within that overall subject. (Thus, if the Document is in part a textbook of mathematics, a Secondary Section may not explain any mathematics.) The relationship could be a matter of historical connection with the subject or with related matters, or of legal, commercial, philosophical, ethical or political position regarding them. The "Invariant Sections" are certain Secondary Sections whose titles are designated, as being those of Invariant Sections, in the notice that says that the Document is released under this License. If a section does not fit the above definition of Secondary then it is not allowed to be designated as Invariant. The Document may contain zero Invariant Sections. If the Document does not identify any Invariant Sections then there are none. The "Cover Texts" are certain short passages of text that are listed, as Front-Cover Texts or Back-Cover Texts, in the notice that says that the Document is released under this License. A Front-Cover Text may be at most 5 words, and a Back-Cover Text may be at most 25 words. A "Transparent" copy of the Document means a machine-readable copy, represented in a format whose specification is available to the general public, that is suitable for revising the document straightforwardly with generic text editors or (for images composed of pixels) generic paint programs or (for drawings) some widely available drawing editor, and that is suitable for input to text formatters or for automatic translation to a variety of formats suitable for input to text formatters. A copy made in an otherwise Transparent file format whose markup, or absence of markup, has been arranged to thwart or discourage subsequent modification by readers is not Transparent. An image format is not Transparent if used for any substantial amount of text. A copy that is not "Transparent" is called "Opaque". Examples of suitable formats for Transparent copies include plain ASCII without markup, Texinfo input format, LaTeX input format, SGML or XML using a publicly available DTD, and standard-conforming simple HTML, PostScript or PDF designed for human modification. Examples of transparent image formats include PNG, XCF and JPG. Opaque formats include proprietary formats that can be read and edited only by proprietary word processors, SGML or XML for which the DTD and/or processing tools are not generally available, and the machine-generated HTML, PostScript or PDF produced by some word processors for output purposes only. The "Title Page" means, for a printed book, the title page itself, plus such following pages as are needed to hold, legibly, the material this License requires to appear in the title page. For works in formats which do not have any title page as such, "Title Page" means the text near the most prominent appearance of the work's title, preceding the beginning of the body of the text. The "publisher" means any person or entity that distributes copies of the Document to the public. A section "Entitled XYZ" means a named subunit of the Document whose title either is precisely XYZ or contains XYZ in parentheses following text that translates XYZ in another language. (Here XYZ stands for a specific section name mentioned below, such as "Acknowledgements", "Dedications", "Endorsements", or "History".) To "Preserve the Title" of such a section when you modify the Document means that it remains a section "Entitled XYZ" according to this definition. The Document may include Warranty Disclaimers next to the notice which states that this License applies to the Document. These Warranty Disclaimers are considered to be included by reference in this License, but only as regards disclaiming warranties: any other implication that these Warranty Disclaimers may have is void and has no effect on the meaning of this License. 2. VERBATIM COPYING You may copy and distribute the Document in any medium, either commercially or noncommercially, provided that this License, the copyright notices, and the license notice saying this License applies to the Document are reproduced in all copies, and that you add no other conditions whatsoever to those of this License. You may not use technical measures to obstruct or control the reading or further copying of the copies you make or distribute. However, you may accept compensation in exchange for copies. If you distribute a large enough number of copies you must also follow the conditions in section 3. You may also lend copies, under the same conditions stated above, and you may publicly display copies. 3. COPYING IN QUANTITY If you publish printed copies (or copies in media that commonly have printed covers) of the Document, numbering more than 100, and the Document's license notice requires Cover Texts, you must enclose the copies in covers that carry, clearly and legibly, all these Cover Texts: Front-Cover Texts on the front cover, and Back-Cover Texts on the back cover. Both covers must also clearly and legibly identify you as the publisher of these copies. The front cover must present the full title with all words of the title equally prominent and visible. You may add other material on the covers in addition. Copying with changes limited to the covers, as long as they preserve the title of the Document and satisfy these conditions, can be treated as verbatim copying in other respects. If the required texts for either cover are too voluminous to fit legibly, you should put the first ones listed (as many as fit reasonably) on the actual cover, and continue the rest onto adjacent pages. If you publish or distribute Opaque copies of the Document numbering more than 100, you must either include a machine-readable Transparent copy along with each Opaque copy, or state in or with each Opaque copy a computer-network location from which the general network-using public has access to download using public-standard network protocols a complete Transparent copy of the Document, free of added material. If you use the latter option, you must take reasonably prudent steps, when you begin distribution of Opaque copies in quantity, to ensure that this Transparent copy will remain thus accessible at the stated location until at least one year after the last time you distribute an Opaque copy (directly or through your agents or retailers) of that edition to the public. It is requested, but not required, that you contact the authors of the Document well before redistributing any large number of copies, to give them a chance to provide you with an updated version of the Document. 4. MODIFICATIONS You may copy and distribute a Modified Version of the Document under the conditions of sections 2 and 3 above, provided that you release the Modified Version under precisely this License, with the Modified Version filling the role of the Document, thus licensing distribution and modification of the Modified Version to whoever possesses a copy of it. In addition, you must do these things in the Modified Version: A. Use in the Title Page (and on the covers, if any) a title distinct from that of the Document, and from those of previous versions (which should, if there were any, be listed in the History section of the Document). You may use the same title as a previous version if the original publisher of that version gives permission. B. List on the Title Page, as authors, one or more persons or entities responsible for authorship of the modifications in the Modified Version, together with at least five of the principal authors of the Document (all of its principal authors, if it has fewer than five), unless they release you from this requirement. C. State on the Title page the name of the publisher of the Modified Version, as the publisher. D. Preserve all the copyright notices of the Document. E. Add an appropriate copyright notice for your modifications adjacent to the other copyright notices. F. Include, immediately after the copyright notices, a license notice giving the public permission to use the Modified Version under the terms of this License, in the form shown in the Addendum below. G. Preserve in that license notice the full lists of Invariant Sections and required Cover Texts given in the Document's license notice. H. Include an unaltered copy of this License. I. Preserve the section Entitled "History", Preserve its Title, and add to it an item stating at least the title, year, new authors, and publisher of the Modified Version as given on the Title Page. If there is no section Entitled "History" in the Document, create one stating the title, year, authors, and publisher of the Document as given on its Title Page, then add an item describing the Modified Version as stated in the previous sentence. J. Preserve the network location, if any, given in the Document for public access to a Transparent copy of the Document, and likewise the network locations given in the Document for previous versions it was based on. These may be placed in the "History" section. You may omit a network location for a work that was published at least four years before the Document itself, or if the original publisher of the version it refers to gives permission. K. For any section Entitled "Acknowledgements" or "Dedications", Preserve the Title of the section, and preserve in the section all the substance and tone of each of the contributor acknowledgements and/or dedications given therein. L. Preserve all the Invariant Sections of the Document, unaltered in their text and in their titles. Section numbers or the equivalent are not considered part of the section titles. M. Delete any section Entitled "Endorsements". Such a section may not be included in the Modified Version. N. Do not retitle any existing section to be Entitled "Endorsements" or to conflict in title with any Invariant Section. O. Preserve any Warranty Disclaimers. If the Modified Version includes new front-matter sections or appendices that qualify as Secondary Sections and contain no material copied from the Document, you may at your option designate some or all of these sections as invariant. To do this, add their titles to the list of Invariant Sections in the Modified Version's license notice. These titles must be distinct from any other section titles. You may add a section Entitled "Endorsements", provided it contains nothing but endorsements of your Modified Version by various parties--for example, statements of peer review or that the text has been approved by an organization as the authoritative definition of a standard. You may add a passage of up to five words as a Front-Cover Text, and a passage of up to 25 words as a Back-Cover Text, to the end of the list of Cover Texts in the Modified Version. Only one passage of Front-Cover Text and one of Back-Cover Text may be added by (or through arrangements made by) any one entity. If the Document already includes a cover text for the same cover, previously added by you or by arrangement made by the same entity you are acting on behalf of, you may not add another; but you may replace the old one, on explicit permission from the previous publisher that added the old one. The author(s) and publisher(s) of the Document do not by this License give permission to use their names for publicity for or to assert or imply endorsement of any Modified Version. 5. COMBINING DOCUMENTS You may combine the Document with other documents released under this License, under the terms defined in section 4 above for modified versions, provided that you include in the combination all of the Invariant Sections of all of the original documents, unmodified, and list them all as Invariant Sections of your combined work in its license notice, and that you preserve all their Warranty Disclaimers. The combined work need only contain one copy of this License, and multiple identical Invariant Sections may be replaced with a single copy. If there are multiple Invariant Sections with the same name but different contents, make the title of each such section unique by adding at the end of it, in parentheses, the name of the original author or publisher of that section if known, or else a unique number. Make the same adjustment to the section titles in the list of Invariant Sections in the license notice of the combined work. In the combination, you must combine any sections Entitled "History" in the various original documents, forming one section Entitled "History"; likewise combine any sections Entitled "Acknowledgements", and any sections Entitled "Dedications". You must delete all sections Entitled "Endorsements." 6. COLLECTIONS OF DOCUMENTS You may make a collection consisting of the Document and other documents released under this License, and replace the individual copies of this License in the various documents with a single copy that is included in the collection, provided that you follow the rules of this License for verbatim copying of each of the documents in all other respects. You may extract a single document from such a collection, and distribute it individually under this License, provided you insert a copy of this License into the extracted document, and follow this License in all other respects regarding verbatim copying of that document. 7. AGGREGATION WITH INDEPENDENT WORKS A compilation of the Document or its derivatives with other separate and independent documents or works, in or on a volume of a storage or distribution medium, is called an "aggregate" if the copyright resulting from the compilation is not used to limit the legal rights of the compilation's users beyond what the individual works permit. When the Document is included in an aggregate, this License does not apply to the other works in the aggregate which are not themselves derivative works of the Document. If the Cover Text requirement of section 3 is applicable to these copies of the Document, then if the Document is less than one half of the entire aggregate, the Document's Cover Texts may be placed on covers that bracket the Document within the aggregate, or the electronic equivalent of covers if the Document is in electronic form. Otherwise they must appear on printed covers that bracket the whole aggregate. 8. TRANSLATION Translation is considered a kind of modification, so you may distribute translations of the Document under the terms of section 4. Replacing Invariant Sections with translations requires special permission from their copyright holders, but you may include translations of some or all Invariant Sections in addition to the original versions of these Invariant Sections. You may include a translation of this License, and all the license notices in the Document, and any Warranty Disclaimers, provided that you also include the original English version of this License and the original versions of those notices and disclaimers. In case of a disagreement between the translation and the original version of this License or a notice or disclaimer, the original version will prevail. If a section in the Document is Entitled "Acknowledgements", "Dedications", or "History", the requirement (section 4) to Preserve its Title (section 1) will typically require changing the actual title. 9. TERMINATION You may not copy, modify, sublicense, or distribute the Document except as expressly provided under this License. Any attempt otherwise to copy, modify, sublicense, or distribute it is void, and will automatically terminate your rights under this License. However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation. Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice. Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, receipt of a copy of some or all of the same material does not give you any rights to use it. 10. FUTURE REVISIONS OF THIS LICENSE The Free Software Foundation may publish new, revised versions of the GNU Free Documentation License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. See `http://www.gnu.org/copyleft/'. Each version of the License is given a distinguishing version number. If the Document specifies that a particular numbered version of this License "or any later version" applies to it, you have the option of following the terms and conditions either of that specified version or of any later version that has been published (not as a draft) by the Free Software Foundation. If the Document does not specify a version number of this License, you may choose any version ever published (not as a draft) by the Free Software Foundation. If the Document specifies that a proxy can decide which future versions of this License can be used, that proxy's public statement of acceptance of a version permanently authorizes you to choose that version for the Document. 11. RELICENSING "Massive Multiauthor Collaboration Site" (or "MMC Site") means any World Wide Web server that publishes copyrightable works and also provides prominent facilities for anybody to edit those works. A public wiki that anybody can edit is an example of such a server. A "Massive Multiauthor Collaboration" (or "MMC") contained in the site means any set of copyrightable works thus published on the MMC site. "CC-BY-SA" means the Creative Commons Attribution-Share Alike 3.0 license published by Creative Commons Corporation, a not-for-profit corporation with a principal place of business in San Francisco, California, as well as future copyleft versions of that license published by that same organization. "Incorporate" means to publish or republish a Document, in whole or in part, as part of another Document. An MMC is "eligible for relicensing" if it is licensed under this License, and if all works that were first published under this License somewhere other than this MMC, and subsequently incorporated in whole or in part into the MMC, (1) had no cover texts or invariant sections, and (2) were thus incorporated prior to November 1, 2008. The operator of an MMC Site may republish an MMC contained in the site under CC-BY-SA on the same site at any time before August 1, 2009, provided the MMC is eligible for relicensing. ADDENDUM: How to use this License for your documents ==================================================== To use this License in a document you have written, include a copy of the License in the document and put the following copyright and license notices just after the title page: Copyright (C) YEAR YOUR NAME. Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled ``GNU Free Documentation License''. If you have Invariant Sections, Front-Cover Texts and Back-Cover Texts, replace the "with...Texts." line with this: with the Invariant Sections being LIST THEIR TITLES, with the Front-Cover Texts being LIST, and with the Back-Cover Texts being LIST. If you have Invariant Sections without Cover Texts, or some other combination of the three, merge those two alternatives to suit the situation. If your document contains nontrivial examples of program code, we recommend releasing these examples in parallel under your choice of free software license, such as the GNU General Public License, to permit their use in free software.