Added MACS source

[htsworkflow.git] / htswanalysis / MACS / lib / gsl / gsl-1.11 / doc / gsl-ref.info-4

This is gsl-ref.info, produced by makeinfo version 4.8 from
gsl-ref.texi.

INFO-DIR-SECTION Scientific software
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 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.2 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."

\1f
File: gsl-ref.info,  Node: Providing a function to minimize,  Next: Multimin Iteration,  Prev: Initializing the Multidimensional Minimizer,  Up: Multidimensional Minimization

35.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;

\1f
File: gsl-ref.info,  Node: Multimin Iteration,  Next: Multimin Stopping Criteria,  Prev: Providing a function to minimize,  Up: Multidimensional Minimization

35.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 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.

\1f
File: gsl-ref.info,  Node: Multimin Stopping Criteria,  Next: Multimin Algorithms,  Prev: Multimin Iteration,  Up: Multidimensional Minimization

35.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.

\1f
File: gsl-ref.info,  Node: Multimin Algorithms,  Next: Multimin Examples,  Prev: Multimin Stopping Criteria,  Up: Multidimensional Minimization

35.7 Algorithms
===============

There are several minimization methods available. The best choice of
algorithm depends on the problem.  All of the algorithms use the value
of the function and its gradient at each evaluation point, except for
the Simplex algorithm which uses function values only.

 -- 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 supercedes 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.

 -- Minimizer: gsl_multimin_fminimizer_nmsimplex
     This is the Simplex algorithm of Nelder and Mead. It constructs n
     vectors p_i from the starting vector X and the vector STEP_SIZE as
     follows:

          p_0 = (x_0, x_1, ... , x_n)
          p_1 = (x_0 + step_size_0, x_1, ... , x_n)
          p_2 = (x_0, x_1 + step_size_1, ... , x_n)
          ... = ...
          p_n = (x_0, x_1, ... , x_n+step_size_n)

     These vectors form the n+1 vertices of a simplex in n dimensions.
     On each iteration the algorithm tries to improve the parameter
     vector p_i corresponding to the highest function value by simple
     geometrical transformations.  These are reflection, reflection
     followed by expansion, contraction and multiple contraction. Using
     these transformations the simplex moves through the parameter
     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 routine calculates the minimizer specific characteristic size
     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'.

\1f
File: gsl-ref.info,  Node: Multimin Examples,  Next: Multimin References and Further Reading,  Prev: Multimin Algorithms,  Up: Multidimensional Minimization

35.8 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_nmsimplex;
       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 = (void *)&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.3ef 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.082
         2  5.250e+00  4.000e+00 f() = 290.625 size = 1.372
         3  5.250e+00  4.000e+00 f() = 290.625 size = 1.372
         4  5.500e+00  1.000e+00 f() = 252.500 size = 1.372
         5  2.625e+00  3.500e+00 f() = 101.406 size = 1.823
         6  3.469e+00  1.375e+00 f() = 98.760  size = 1.526
         7  1.820e+00  3.156e+00 f() = 63.467  size = 1.105
         8  1.820e+00  3.156e+00 f() = 63.467  size = 1.105
         9  1.016e+00  2.812e+00 f() = 43.206  size = 1.105
        10  2.041e+00  2.008e+00 f() = 40.838  size = 0.645
        11  1.236e+00  1.664e+00 f() = 32.816  size = 0.645
        12  1.236e+00  1.664e+00 f() = 32.816  size = 0.447
        13  5.225e-01  1.980e+00 f() = 32.288  size = 0.447
        14  1.103e+00  2.073e+00 f() = 30.214  size = 0.345
        15  1.103e+00  2.073e+00 f() = 30.214  size = 0.264
        16  1.103e+00  2.073e+00 f() = 30.214  size = 0.160
        17  9.864e-01  1.934e+00 f() = 30.090  size = 0.132
        18  9.190e-01  1.987e+00 f() = 30.069  size = 0.092
        19  1.028e+00  2.017e+00 f() = 30.013  size = 0.056
        20  1.028e+00  2.017e+00 f() = 30.013  size = 0.046
        21  1.028e+00  2.017e+00 f() = 30.013  size = 0.033
        22  9.874e-01  1.985e+00 f() = 30.006  size = 0.028
        23  9.846e-01  1.995e+00 f() = 30.003  size = 0.023
        24  1.007e+00  2.003e+00 f() = 30.001  size = 0.012
     converged to minimum at
        25  1.007e+00  2.003e+00 f() = 30.001  size = 0.010

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.

\1f
File: gsl-ref.info,  Node: Multimin References and Further Reading,  Prev: Multimin Examples,  Up: Multidimensional Minimization

35.9 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-315.


\1f
File: gsl-ref.info,  Node: Least-Squares Fitting,  Next: Nonlinear Least-Squares Fitting,  Prev: Multidimensional Minimization,  Up: Top

36 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::

\1f
File: gsl-ref.info,  Node: Fitting Overview,  Next: Linear regression,  Up: Least-Squares Fitting

36.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}}.

\1f
File: gsl-ref.info,  Node: Linear regression,  Next: Linear fitting without a constant term,  Prev: Fitting Overview,  Up: Least-Squares Fitting

36.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.

\1f
File: gsl-ref.info,  Node: Linear fitting without a constant term,  Next: Multi-parameter fitting,  Prev: Linear regression,  Up: Least-Squares Fitting

36.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.

\1f
File: gsl-ref.info,  Node: Multi-parameter fitting,  Next: Fitting Examples,  Prev: Linear fitting without a constant term,  Up: Least-Squares Fitting

36.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)
 -- 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)
     These functions compute the best-fit parameters C of the model y =
     X c for the observations Y and the matrix of predictor variables
     X.  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 preallocated workspace provided in WORK. The modified
     Golub-Reinsch SVD algorithm is used, 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.  In the second form of the function the components 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_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)
 -- 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)
     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.  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'.

     The best-fit is found by singular value decomposition of the matrix
     X using the preallocated workspace provided in WORK. Any
     components which have zero singular value (to machine precision)
     are discarded from the fit.  In the second form of the function the
     components 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_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.

\1f
File: gsl-ref.info,  Node: Fitting Examples,  Next: Fitting References and Further Reading,  Prev: Multi-parameter fitting,  Up: Least-Squares Fitting

36.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 <stdio.h>
     #include <gsl/gsl_fit.h>

     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 <stdio.h>
     #include <gsl/gsl_multifit.h>

     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 <stdio.h>
     #include <math.h>
     #include <gsl/gsl_randist.h>

     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.

\1f
File: gsl-ref.info,  Node: Fitting References and Further Reading,  Prev: Fitting Examples,  Up: Least-Squares Fitting

36.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'.

\1f
File: gsl-ref.info,  Node: Nonlinear Least-Squares Fitting,  Next: Basis Splines,  Prev: Least-Squares Fitting,  Up: Top

37 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::

\1f
File: gsl-ref.info,  Node: Overview of Nonlinear Least-Squares Fitting,  Next: Initializing the Nonlinear Least-Squares Solver,  Up: Nonlinear Least-Squares Fitting

37.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).

\1f
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

37.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, 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, 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'.

\1f
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

37.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.

\1f
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

37.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.

\1f
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

37.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.

\1f
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

37.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

     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.

\1f
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

37.7 Minimization Algorithms without Derivatives
================================================

There are no algorithms implemented in this section at the moment.

\1f
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

37.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::.

\1f
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

37.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 <stdlib.h>
     #include <stdio.h>
     #include <gsl/gsl_rng.h>
     #include <gsl/gsl_randist.h>
     #include <gsl/gsl_vector.h>
     #include <gsl/gsl_blas.h>
     #include <gsl/gsl_multifit_nlin.h>

     #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.

\1f
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

37.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.

\1f
File: gsl-ref.info,  Node: Basis Splines,  Next: Physical Constants,  Prev: Nonlinear Least-Squares Fitting,  Up: Top

38 Basis Splines
****************

This chapter describes functions for the computation of smoothing basis
splines (B-splines). The header file `gsl_bspline.h' contains
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::
* Example programs for B-splines::
* References and Further Reading::

\1f
File: gsl-ref.info,  Node: Overview of B-splines,  Next: Initializing the B-splines solver,  Up: Basis Splines

38.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 c_i can be readily
obtained from a least-squares fit.

\1f
File: gsl-ref.info,  Node: Initializing the B-splines solver,  Next: Constructing the knots vector,  Prev: Overview of B-splines,  Up: Basis Splines

38.2 Initializing the B-splines solver
======================================

 -- 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.

\1f
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

38.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'.

\1f
File: gsl-ref.info,  Node: Evaluation of B-spline basis functions,  Next: Example programs for B-splines,  Prev: Constructing the knots vector,  Up: Basis Splines

38.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 B, so that the ith element of B is
     B_i(x). B must be of length n = nbreak + k - 2. This value may
     also be obtained by calling `gsl_bspline_ncoeffs'.  It is far more
     efficient to compute all of the basis functions at once than to
     compute them individually, due to the nature of the defining
     recurrence relation.

 -- 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.

\1f
File: gsl-ref.info,  Node: Example programs for B-splines,  Next: References and Further Reading,  Prev: Evaluation of B-spline basis functions,  Up: Basis Splines

38.5 Example programs for B-splines
===================================

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 [0, 15] with
gaussian noise added.

     #include <stdio.h>
     #include <stdlib.h>
     #include <math.h>
     #include <gsl/gsl_bspline.h>
     #include <gsl/gsl_multifit.h>
     #include <gsl/gsl_rng.h>
     #include <gsl/gsl_randist.h>
     #include <gsl/gsl_statistics.h>

     /* 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;
       double Rsq;
       double dof;

       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;
       Rsq = 1.0 - chisq / gsl_stats_wtss(w->data, 1, y->data, 1, y->size);

       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

\1f
File: gsl-ref.info,  Node: References and Further Reading,  Prev: Example programs for B-splines,  Up: Basis Splines

38.6 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.

A large collection of B-spline routines is available in the PPPACK
library available at `http://www.netlib.org/pppack'.

\1f
File: gsl-ref.info,  Node: Physical Constants,  Next: IEEE floating-point arithmetic,  Prev: Basis Splines,  Up: Top

39 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.

\1f
File: gsl-ref.info,  Node: Fundamental Constants,  Next: Astronomy and Astrophysics,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Astronomy and Astrophysics,  Next: Atomic and Nuclear Physics,  Prev: Fundamental Constants,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Atomic and Nuclear Physics,  Next: Measurement of Time,  Prev: Astronomy and Astrophysics,  Up: Physical Constants

39.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 by Ry = h c R.

`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.

\1f
File: gsl-ref.info,  Node: Measurement of Time,  Next: Imperial Units,  Prev: Atomic and Nuclear Physics,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Imperial Units,  Next: Speed and Nautical Units,  Prev: Measurement of Time,  Up: Physical Constants

39.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).

\1f
File: gsl-ref.info,  Node: Speed and Nautical Units,  Next: Printers Units,  Prev: Imperial Units,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Printers Units,  Next: Volume Area and Length,  Prev: Speed and Nautical Units,  Up: Physical Constants

39.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).

\1f
File: gsl-ref.info,  Node: Volume Area and Length,  Next: Mass and Weight,  Prev: Printers Units,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Mass and Weight,  Next: Thermal Energy and Power,  Prev: Volume Area and Length,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Thermal Energy and Power,  Next: Pressure,  Prev: Mass and Weight,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Pressure,  Next: Viscosity,  Prev: Thermal Energy and Power,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Viscosity,  Next: Light and Illumination,  Prev: Pressure,  Up: Physical Constants

39.12 Viscosity
===============

`GSL_CONST_MKSA_POISE'
     The dynamic viscosity of 1 poise.

`GSL_CONST_MKSA_STOKES'
     The kinematic viscosity of 1 stokes.

\1f
File: gsl-ref.info,  Node: Light and Illumination,  Next: Radioactivity,  Prev: Viscosity,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Radioactivity,  Next: Force and Energy,  Prev: Light and Illumination,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Force and Energy,  Next: Prefixes,  Prev: Radioactivity,  Up: Physical Constants

39.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.

\1f
File: gsl-ref.info,  Node: Prefixes,  Next: Physical Constant Examples,  Prev: Force and Energy,  Up: Physical Constants

39.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

\1f
File: gsl-ref.info,  Node: Physical Constant Examples,  Next: Physical Constant References and Further Reading,  Prev: Prefixes,  Up: Physical Constants

39.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 <stdio.h>
     #include <gsl/gsl_const_mksa.h>

     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

\1f
File: gsl-ref.info,  Node: Physical Constant References and Further Reading,  Prev: Physical Constant Examples,  Up: Physical Constants

39.18 References and Further Reading
====================================

The authoritative sources for physical constants are the 2002 CODATA
recommended values, published in the articles below. Further information
on the values of physical constants is also available from the cited
articles and the NIST website.

     Journal of Physical and Chemical Reference Data, 28(6), 1713-1852,
     1999

     Reviews of Modern Physics, 72(2), 351-495, 2000

     `http://www.physics.nist.gov/cuu/Constants/index.html'

     `http://physics.nist.gov/Pubs/SP811/appenB9.html'

\1f
File: gsl-ref.info,  Node: IEEE floating-point arithmetic,  Next: Debugging Numerical Programs,  Prev: Physical Constants,  Up: Top

40 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::

\1f
File: gsl-ref.info,  Node: Representation of floating point numbers,  Next: Setting up your IEEE environment,  Up: IEEE floating-point arithmetic

40.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 <stdio.h>
     #include <gsl/gsl_ieee_utils.h>

     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.

\1f
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

40.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 <stdio.h>
     #include <gsl/gsl_math.h>
     #include <gsl/gsl_ieee_utils.h>

     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.

\1f
File: gsl-ref.info,  Node: IEEE References and Further Reading,  Prev: Setting up your IEEE environment,  Up: IEEE floating-point arithmetic

40.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.


\1f
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::

\1f
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'.

\1f
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

\1f
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

\1f
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
       -Wtraditional -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.

`-Wtraditional'
     This warns about certain constructs that behave differently in
     traditional and ANSI C. Whether the traditional or ANSI
     interpretation is used might be unpredictable on other compilers.

`-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.

\1f
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

\1f
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.

*Lowell Johnson*
     Implementation of the Mathieu functions.

*Patrick Alken*
     Implementation of non-symmetric and generalized eigensystems and
     B-splines.


   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.

\1f
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
<config.h>' to substitute `gsl_hypot' for each occurrence of `hypot'
when `hypot' is not available.

\1f
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::

\1f
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)

\1f
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)

\1f
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, ...)

\1f
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 <stdio.h>
     #include <gsl/gsl_cblas.h>

     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 ]

\1f
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 <licensing@gnu.org>.

   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'

\1f
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'.

\1f
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.2, November 2002
     Copyright (C) 2000,2001,2002 Free Software Foundation, Inc.
     51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA

     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.

     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 for under this License.  Any other
     attempt to copy, modify, sublicense or distribute the Document is
     void, and will automatically terminate your rights under this
     License.  However, parties who have received copies, or rights,
     from you under this License will not have their licenses
     terminated so long as such parties remain in full compliance.

 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.

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.2 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.