path: root/gsl-1.9/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 freedom to copy and modify this
GNU Manual, like GNU software."


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 C11,
          double * Y, double * Y_ERR)
     This function uses the best-fit linear regression coefficient C1
     and its covariance COV11 to compute the fitted function Y and its
     standard deviation Y_ERR for the model Y = c_1 X at the point X.


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

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.

     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.

     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.


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.


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


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


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


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


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.


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.


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.


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.


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.


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


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.


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.


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


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, \dots,
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, \dots, 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.


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.


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


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 is also
     stored in `w->n'.  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.


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(-0.1 x) 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>

     /* number of data points to fit */
     #define N        200

     /* number of fit coefficients */
     #define NCOEFFS  8

     /* 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;

       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;
           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);

       /* 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
     $ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps


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

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


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.


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.


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.


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.


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.


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


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.


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


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.


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.


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.


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.


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.


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.


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.


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.


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


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


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'


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


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.


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.

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 series no
longer converges (the reason is that when we add each term to the sum
the final result is always rounded up.  This is guaranteed to increase
the sum by at least one tick on each iteration).  To avoid this problem
we would need to use a safer converge criterion, such as `while
(fabs(sum - oldsum) > epsilon)', with a suitably chosen value of
epsilon.

   Finally we can see the effect of computing the sum using
single-precision rounding, in the default `round-to-nearest' mode.  In
this case the program thinks it is still using double precision numbers
but the CPU rounds the result of each floating point operation to
single-precision accuracy.  This simulates the effect of writing the
program using single-precision `float' variables instead of `double'
variables.  The iteration stops after about half the number of
iterations and the final result is much less accurate,

     $ GSL_IEEE_MODE="single-precision" ./a.out
     ....
     i=12 sum=2.718281984329223633 error=1.5587e-07

with an error of O(10^-7), which corresponds to single precision
accuracy (about 1 part in 10^7).  Continuing the iterations further
does not decrease the error because all the subsequent results are
rounded to the same value.


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

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.



File: gsl-ref.info,  Node: Debugging Numerical Programs,  Next: Contributors to GSL,  Prev: IEEE floating-point arithmetic,  Up: Top

Appendix A Debugging Numerical Programs
***************************************

This chapter describes some tips and tricks for debugging numerical
programs which use GSL.

* Menu:

* Using gdb::
* Examining floating point registers::
* Handling floating point exceptions::
* GCC warning options for numerical programs::
* Debugging References::


File: gsl-ref.info,  Node: Using gdb,  Next: Examining floating point registers,  Up: Debugging Numerical Programs

A.1 Using gdb
=============

Any errors reported by the library are passed to the function
`gsl_error'.  By running your programs under gdb and setting a
breakpoint in this function you can automatically catch any library
errors.  You can add a breakpoint for every session by putting

     break gsl_error

into your `.gdbinit' file in the directory where your program is
started.

   If the breakpoint catches an error then you can use a backtrace
(`bt') to see the call-tree, and the arguments which possibly caused
the error.  By moving up into the calling function you can investigate
the values of variables at that point.  Here is an example from the
program `fft/test_trap', which contains the following line,

     status = gsl_fft_complex_wavetable_alloc (0, &complex_wavetable);

The function `gsl_fft_complex_wavetable_alloc' takes the length of an
FFT as its first argument.  When this line is executed an error will be
generated because the length of an FFT is not allowed to be zero.

   To debug this problem we start `gdb', using the file `.gdbinit' to
define a breakpoint in `gsl_error',

     $ gdb test_trap

     GDB is free software and you are welcome to distribute copies
     of it under certain conditions; type "show copying" to see
     the conditions.  There is absolutely no warranty for GDB;
     type "show warranty" for details.  GDB 4.16 (i586-debian-linux),
     Copyright 1996 Free Software Foundation, Inc.

     Breakpoint 1 at 0x8050b1e: file error.c, line 14.

When we run the program this breakpoint catches the error and shows the
reason for it.

     (gdb) run
     Starting program: test_trap

     Breakpoint 1, gsl_error (reason=0x8052b0d
         "length n must be positive integer",
         file=0x8052b04 "c_init.c", line=108, gsl_errno=1)
         at error.c:14
     14        if (gsl_error_handler)

The first argument of `gsl_error' is always a string describing the
error.  Now we can look at the backtrace to see what caused the problem,

     (gdb) bt
     #0  gsl_error (reason=0x8052b0d
         "length n must be positive integer",
         file=0x8052b04 "c_init.c", line=108, gsl_errno=1)
         at error.c:14
     #1  0x8049376 in gsl_fft_complex_wavetable_alloc (n=0,
         wavetable=0xbffff778) at c_init.c:108
     #2  0x8048a00 in main (argc=1, argv=0xbffff9bc)
         at test_trap.c:94
     #3  0x80488be in ___crt_dummy__ ()

We can see that the error was generated in the function
`gsl_fft_complex_wavetable_alloc' when it was called with an argument
of N=0.  The original call came from line 94 in the file `test_trap.c'.

   By moving up to the level of the original call we can find the line
that caused the error,

     (gdb) up
     #1  0x8049376 in gsl_fft_complex_wavetable_alloc (n=0,
         wavetable=0xbffff778) at c_init.c:108
     108   GSL_ERROR ("length n must be positive integer", GSL_EDOM);
     (gdb) up
     #2  0x8048a00 in main (argc=1, argv=0xbffff9bc)
         at test_trap.c:94
     94    status = gsl_fft_complex_wavetable_alloc (0,
             &complex_wavetable);

Thus we have found the line that caused the problem.  From this point we
could also print out the values of other variables such as
`complex_wavetable'.


File: gsl-ref.info,  Node: Examining floating point registers,  Next: Handling floating point exceptions,  Prev: Using gdb,  Up: Debugging Numerical Programs

A.2 Examining floating point registers
======================================

The contents of floating point registers can be examined using the
command `info float' (on supported platforms).

     (gdb) info float
          st0: 0xc4018b895aa17a945000  Valid Normal -7.838871e+308
          st1: 0x3ff9ea3f50e4d7275000  Valid Normal 0.0285946
          st2: 0x3fe790c64ce27dad4800  Valid Normal 6.7415931e-08
          st3: 0x3ffaa3ef0df6607d7800  Spec  Normal 0.0400229
          st4: 0x3c028000000000000000  Valid Normal 4.4501477e-308
          st5: 0x3ffef5412c22219d9000  Zero  Normal 0.9580257
          st6: 0x3fff8000000000000000  Valid Normal 1
          st7: 0xc4028b65a1f6d243c800  Valid Normal -1.566206e+309
        fctrl: 0x0272 53 bit; NEAR; mask DENOR UNDER LOS;
        fstat: 0xb9ba flags 0001; top 7; excep DENOR OVERF UNDER LOS
         ftag: 0x3fff
          fip: 0x08048b5c
          fcs: 0x051a0023
       fopoff: 0x08086820
       fopsel: 0x002b

Individual registers can be examined using the variables $REG, where
REG is the register name.

     (gdb) p $st1
     $1 = 0.02859464454261210347719


File: gsl-ref.info,  Node: Handling floating point exceptions,  Next: GCC warning options for numerical programs,  Prev: Examining floating point registers,  Up: Debugging Numerical Programs

A.3 Handling floating point exceptions
======================================

It is possible to stop the program whenever a `SIGFPE' floating point
exception occurs.  This can be useful for finding the cause of an
unexpected infinity or `NaN'.  The current handler settings can be
shown with the command `info signal SIGFPE'.

     (gdb) info signal SIGFPE
     Signal  Stop  Print  Pass to program Description
     SIGFPE  Yes   Yes    Yes             Arithmetic exception

Unless the program uses a signal handler the default setting should be
changed so that SIGFPE is not passed to the program, as this would cause
it to exit.  The command `handle SIGFPE stop nopass' prevents this.

     (gdb) handle SIGFPE stop nopass
     Signal  Stop  Print  Pass to program Description
     SIGFPE  Yes   Yes    No              Arithmetic exception

Depending on the platform it may be necessary to instruct the kernel to
generate signals for floating point exceptions.  For programs using GSL
this can be achieved using the `GSL_IEEE_MODE' environment variable in
conjunction with the function `gsl_ieee_env_setup' as described in
*note IEEE floating-point arithmetic::.

     (gdb) set env GSL_IEEE_MODE=double-precision


File: gsl-ref.info,  Node: GCC warning options for numerical programs,  Next: Debugging References,  Prev: Handling floating point exceptions,  Up: Debugging Numerical Programs

A.4 GCC warning options for numerical programs
==============================================

Writing reliable numerical programs in C requires great care.  The
following GCC warning options are recommended when compiling numerical
programs:

     gcc -ansi -pedantic -Werror -Wall -W
       -Wmissing-prototypes -Wstrict-prototypes
       -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.


File: gsl-ref.info,  Node: Debugging References,  Prev: GCC warning options for numerical programs,  Up: Debugging Numerical Programs

A.5 References and Further Reading
==================================

The following books are essential reading for anyone writing and
debugging numerical programs with GCC and GDB.

     R.M. Stallman, `Using and Porting GNU CC', Free Software
     Foundation, ISBN 1882114388

     R.M. Stallman, R.H. Pesch, `Debugging with GDB: The GNU
     Source-Level Debugger', Free Software Foundation, ISBN 1882114779

For a tutorial introduction to the GNU C Compiler and related programs,
see

     B.J. Gough, `An Introduction to GCC', Network Theory Ltd, ISBN
     0954161793


File: gsl-ref.info,  Node: Contributors to GSL,  Next: Autoconf Macros,  Prev: Debugging Numerical Programs,  Up: Top

Appendix B Contributors to GSL
******************************

(See the AUTHORS file in the distribution for up-to-date information.)

*Mark Galassi*
     Conceived GSL (with James Theiler) and wrote the design document.
     Wrote the simulated annealing package and the relevant chapter in
     the manual.

*James Theiler*
     Conceived GSL (with Mark Galassi).  Wrote the random number
     generators and the relevant chapter in this manual.

*Jim Davies*
     Wrote the statistical routines and the relevant chapter in this
     manual.

*Brian Gough*
     FFTs, numerical integration, random number generators and
     distributions, root finding, minimization and fitting, polynomial
     solvers, complex numbers, physical constants, permutations, vector
     and matrix functions, histograms, statistics, ieee-utils, revised
     CBLAS Level 2 & 3, matrix decompositions, eigensystems, cumulative
     distribution functions, testing, documentation and releases.

*Reid Priedhorsky*
     Wrote and documented the initial version of the root finding
     routines while at Los Alamos National Laboratory, Mathematical
     Modeling and Analysis Group.

*Gerard Jungman*
     Special Functions, Series acceleration, ODEs, BLAS, Linear Algebra,
     Eigensystems, Hankel Transforms.

*Mike Booth*
     Wrote the Monte Carlo library.

*Jorma Olavi Ta"htinen*
     Wrote the initial complex arithmetic functions.

*Thomas Walter*
     Wrote the initial heapsort routines and cholesky decomposition.

*Fabrice Rossi*
     Multidimensional minimization.

*Carlo Perassi*
     Implementation of the random number generators in Knuth's
     `Seminumerical Algorithms', 3rd Ed.

*Szymon Jaroszewicz*
     Wrote the routines for generating combinations.

*Nicolas Darnis*
     Wrote the initial routines 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 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.


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,main)
     AC_CHECK_LIB(gslcblas,main)
     AC_CHECK_LIB(gsl,main)

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 main in -lm... yes
     checking for main in -lgslcblas... yes
     checking for main 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:

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


File: gsl-ref.info,  Node: GSL CBLAS Library,  Next: Free Software Needs Free Documentation,  Prev: Autoconf Macros,  Up: Top

Appendix D GSL CBLAS Library
****************************

The prototypes for the low-level CBLAS functions are declared in the
file `gsl_cblas.h'.  For the definition of the functions consult the
documentation available from Netlib (*note BLAS References and Further
Reading::).

* Menu:

* Level 1 CBLAS Functions::
* Level 2 CBLAS Functions::
* Level 3 CBLAS Functions::
* GSL CBLAS Examples::


File: gsl-ref.info,  Node: Level 1 CBLAS Functions,  Next: Level 2 CBLAS Functions,  Up: GSL CBLAS Library

D.1 Level 1
===========

 -- Function: float cblas_sdsdot (const int N, const float ALPHA, const
          float * X, const int INCX, const float * Y, const int INCY)

 -- Function: double cblas_dsdot (const int N, const float * X, const
          int INCX, const float * Y, const int INCY)

 -- Function: float cblas_sdot (const int N, const float * X, const int
          INCX, const float * Y, const int INCY)

 -- Function: double cblas_ddot (const int N, const double * X, const
          int INCX, const double * Y, const int INCY)

 -- Function: void cblas_cdotu_sub (const int N, const void * X, const
          int INCX, const void * Y, const int INCY, void * DOTU)

 -- Function: void cblas_cdotc_sub (const int N, const void * X, const
          int INCX, const void * Y, const int INCY, void * DOTC)

 -- Function: void cblas_zdotu_sub (const int N, const void * X, const
          int INCX, const void * Y, const int INCY, void * DOTU)

 -- Function: void cblas_zdotc_sub (const int N, const void * X, const
          int INCX, const void * Y, const int INCY, void * DOTC)

 -- Function: float cblas_snrm2 (const int N, const float * X, const
          int INCX)

 -- Function: float cblas_sasum (const int N, const float * X, const
          int INCX)

 -- Function: double cblas_dnrm2 (const int N, const double * X, const
          int INCX)

 -- Function: double cblas_dasum (const int N, const double * X, const
          int INCX)

 -- Function: float cblas_scnrm2 (const int N, const void * X, const
          int INCX)

 -- Function: float cblas_scasum (const int N, const void * X, const
          int INCX)

 -- Function: double cblas_dznrm2 (const int N, const void * X, const
          int INCX)

 -- Function: double cblas_dzasum (const int N, const void * X, const
          int INCX)

 -- Function: CBLAS_INDEX cblas_isamax (const int N, const float * X,
          const int INCX)

 -- Function: CBLAS_INDEX cblas_idamax (const int N, const double * X,
          const int INCX)

 -- Function: CBLAS_INDEX cblas_icamax (const int N, const void * X,
          const int INCX)

 -- Function: CBLAS_INDEX cblas_izamax (const int N, const void * X,
          const int INCX)

 -- Function: void cblas_sswap (const int N, float * X, const int INCX,
          float * Y, const int INCY)

 -- Function: void cblas_scopy (const int N, const float * X, const int
          INCX, float * Y, const int INCY)

 -- Function: void cblas_saxpy (const int N, const float ALPHA, const
          float * X, const int INCX, float * Y, const int INCY)

 -- Function: void cblas_dswap (const int N, double * X, const int
          INCX, double * Y, const int INCY)

 -- Function: void cblas_dcopy (const int N, const double * X, const
          int INCX, double * Y, const int INCY)

 -- Function: void cblas_daxpy (const int N, const double ALPHA, const
          double * X, const int INCX, double * Y, const int INCY)

 -- Function: void cblas_cswap (const int N, void * X, const int INCX,
          void * Y, const int INCY)

 -- Function: void cblas_ccopy (const int N, const void * X, const int
          INCX, void * Y, const int INCY)

 -- Function: void cblas_caxpy (const int N, const void * ALPHA, const
          void * X, const int INCX, void * Y, const int INCY)

 -- Function: void cblas_zswap (const int N, void * X, const int INCX,
          void * Y, const int INCY)

 -- Function: void cblas_zcopy (const int N, const void * X, const int
          INCX, void * Y, const int INCY)

 -- Function: void cblas_zaxpy (const int N, const void * ALPHA, const
          void * X, const int INCX, void * Y, const int INCY)

 -- Function: void cblas_srotg (float * A, float * B, float * C, float
          * S)

 -- Function: void cblas_srotmg (float * D1, float * D2, float * B1,
          const float B2, float * P)

 -- Function: void cblas_srot (const int N, float * X, const int INCX,
          float * Y, const int INCY, const float C, const float S)

 -- Function: void cblas_srotm (const int N, float * X, const int INCX,
          float * Y, const int INCY, const float * P)

 -- Function: void cblas_drotg (double * A, double * B, double * C,
          double * S)

 -- Function: void cblas_drotmg (double * D1, double * D2, double * B1,
          const double B2, double * P)

 -- Function: void cblas_drot (const int N, double * X, const int INCX,
          double * Y, const int INCY, const double C, const double S)

 -- Function: void cblas_drotm (const int N, double * X, const int
          INCX, double * Y, const int INCY, const double * P)

 -- Function: void cblas_sscal (const int N, const float ALPHA, float *
          X, const int INCX)

 -- Function: void cblas_dscal (const int N, const double ALPHA, double
          * X, const int INCX)

 -- Function: void cblas_cscal (const int N, const void * ALPHA, void *
          X, const int INCX)

 -- Function: void cblas_zscal (const int N, const void * ALPHA, void *
          X, const int INCX)

 -- Function: void cblas_csscal (const int N, const float ALPHA, void *
          X, const int INCX)

 -- Function: void cblas_zdscal (const int N, const double ALPHA, void
          * X, const int INCX)


File: gsl-ref.info,  Node: Level 2 CBLAS Functions,  Next: Level 3 CBLAS Functions,  Prev: Level 1 CBLAS Functions,  Up: GSL CBLAS Library

D.2 Level 2
===========

 -- Function: void cblas_sgemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          float ALPHA, const float * A, const int LDA, const float * X,
          const int INCX, const float BETA, float * Y, const int INCY)

 -- Function: void cblas_sgbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          int KL, const int KU, const float ALPHA, const float * A,
          const int LDA, const float * X, const int INCX, const float
          BETA, float * Y, const int INCY)

 -- Function: void cblas_strmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const float * A,
          const int LDA, float * X, const int INCX)

 -- Function: void cblas_stbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          float * A, const int LDA, float * X, const int INCX)

 -- Function: void cblas_stpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const float * AP,
          float * X, const int INCX)

 -- Function: void cblas_strsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const float * A,
          const int LDA, float * X, const int INCX)

 -- Function: void cblas_stbsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          float * A, const int LDA, float * X, const int INCX)

 -- Function: void cblas_stpsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const float * AP,
          float * X, const int INCX)

 -- Function: void cblas_dgemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          double ALPHA, const double * A, const int LDA, const double *
          X, const int INCX, const double BETA, double * Y, const int
          INCY)

 -- Function: void cblas_dgbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          int KL, const int KU, const double ALPHA, const double * A,
          const int LDA, const double * X, const int INCX, const double
          BETA, double * Y, const int INCY)

 -- Function: void cblas_dtrmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const double * A,
          const int LDA, double * X, const int INCX)

 -- Function: void cblas_dtbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          double * A, const int LDA, double * X, const int INCX)

 -- Function: void cblas_dtpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const double * AP,
          double * X, const int INCX)

 -- Function: void cblas_dtrsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const double * A,
          const int LDA, double * X, const int INCX)

 -- Function: void cblas_dtbsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          double * A, const int LDA, double * X, const int INCX)

 -- Function: void cblas_dtpsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const double * AP,
          double * X, const int INCX)

 -- Function: void cblas_cgemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          void * ALPHA, const void * A, const int LDA, const void * X,
          const int INCX, const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_cgbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          int KL, const int KU, const void * ALPHA, const void * A,
          const int LDA, const void * X, const int INCX, const void *
          BETA, void * Y, const int INCY)

 -- Function: void cblas_ctrmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * A,
          const int LDA, void * X, const int INCX)

 -- Function: void cblas_ctbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          void * A, const int LDA, void * X, const int INCX)

 -- Function: void cblas_ctpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * AP,
          void * X, const int INCX)

 -- Function: void cblas_ctrsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * A,
          const int LDA, void * X, const int INCX)

 -- Function: void cblas_ctbsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          void * A, const int LDA, void * X, const int INCX)

 -- Function: void cblas_ctpsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * AP,
          void * X, const int INCX)

 -- Function: void cblas_zgemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          void * ALPHA, const void * A, const int LDA, const void * X,
          const int INCX, const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_zgbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const int M, const int N, const
          int KL, const int KU, const void * ALPHA, const void * A,
          const int LDA, const void * X, const int INCX, const void *
          BETA, void * Y, const int INCY)

 -- Function: void cblas_ztrmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * A,
          const int LDA, void * X, const int INCX)

 -- Function: void cblas_ztbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          void * A, const int LDA, void * X, const int INCX)

 -- Function: void cblas_ztpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * AP,
          void * X, const int INCX)

 -- Function: void cblas_ztrsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * A,
          const int LDA, void * X, const int INCX)

 -- Function: void cblas_ztbsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const int K, const
          void * A, const int LDA, void * X, const int INCX)

 -- Function: void cblas_ztpsv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANSA,
          const enum CBLAS_DIAG DIAG, const int N, const void * AP,
          void * X, const int INCX)

 -- Function: void cblas_ssymv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const
          float * A, const int LDA, const float * X, const int INCX,
          const float BETA, float * Y, const int INCY)

 -- Function: void cblas_ssbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const int K, const float
          ALPHA, const float * A, const int LDA, const float * X, const
          int INCX, const float BETA, float * Y, const int INCY)

 -- Function: void cblas_sspmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const
          float * AP, const float * X, const int INCX, const float
          BETA, float * Y, const int INCY)

 -- Function: void cblas_sger (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const float ALPHA, const float * X, const int
          INCX, const float * Y, const int INCY, float * A, const int
          LDA)

 -- Function: void cblas_ssyr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const float ALPHA, const float
          * X, const int INCX, float * A, const int LDA)

 -- Function: void cblas_sspr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const float ALPHA, const float
          * X, const int INCX, float * AP)

 -- Function: void cblas_ssyr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const
          float * X, const int INCX, const float * Y, const int INCY,
          float * A, const int LDA)

 -- Function: void cblas_sspr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const float ALPHA, const
          float * X, const int INCX, const float * Y, const int INCY,
          float * A)

 -- Function: void cblas_dsymv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * A, const int LDA, const double * X, const int INCX,
          const double BETA, double * Y, const int INCY)

 -- Function: void cblas_dsbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const int K, const double
          ALPHA, const double * A, const int LDA, const double * X,
          const int INCX, const double BETA, double * Y, const int INCY)

 -- Function: void cblas_dspmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * AP, const double * X, const int INCX, const double
          BETA, double * Y, const int INCY)

 -- Function: void cblas_dger (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const double ALPHA, const double * X, const
          int INCX, const double * Y, const int INCY, double * A, const
          int LDA)

 -- Function: void cblas_dsyr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * X, const int INCX, double * A, const int LDA)

 -- Function: void cblas_dspr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * X, const int INCX, double * AP)

 -- Function: void cblas_dsyr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * X, const int INCX, const double * Y, const int INCY,
          double * A, const int LDA)

 -- Function: void cblas_dspr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const double ALPHA, const
          double * X, const int INCX, const double * Y, const int INCY,
          double * A)

 -- Function: void cblas_chemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * A, const int LDA, const void * X, const int INCX,
          const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_chbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const int K, const void *
          ALPHA, const void * A, const int LDA, const void * X, const
          int INCX, const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_chpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * AP, const void * X, const int INCX, const void * BETA,
          void * Y, const int INCY)

 -- Function: void cblas_cgeru (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const void * ALPHA, const void * X, const int
          INCX, const void * Y, const int INCY, void * A, const int LDA)

 -- Function: void cblas_cgerc (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const void * ALPHA, const void * X, const int
          INCX, const void * Y, const int INCY, void * A, const int LDA)

 -- Function: void cblas_cher (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const float ALPHA, const void *
          X, const int INCX, void * A, const int LDA)

 -- Function: void cblas_chpr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const float ALPHA, const void *
          X, const int INCX, void * A)

 -- Function: void cblas_cher2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * X, const int INCX, const void * Y, const int INCY,
          void * A, const int LDA)

 -- Function: void cblas_chpr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * X, const int INCX, const void * Y, const int INCY,
          void * AP)

 -- Function: void cblas_zhemv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * A, const int LDA, const void * X, const int INCX,
          const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_zhbmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const int K, const void *
          ALPHA, const void * A, const int LDA, const void * X, const
          int INCX, const void * BETA, void * Y, const int INCY)

 -- Function: void cblas_zhpmv (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * AP, const void * X, const int INCX, const void * BETA,
          void * Y, const int INCY)

 -- Function: void cblas_zgeru (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const void * ALPHA, const void * X, const int
          INCX, const void * Y, const int INCY, void * A, const int LDA)

 -- Function: void cblas_zgerc (const enum CBLAS_ORDER ORDER, const int
          M, const int N, const void * ALPHA, const void * X, const int
          INCX, const void * Y, const int INCY, void * A, const int LDA)

 -- Function: void cblas_zher (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const double ALPHA, const void
          * X, const int INCX, void * A, const int LDA)

 -- Function: void cblas_zhpr (const enum CBLAS_ORDER ORDER, const enum
          CBLAS_UPLO UPLO, const int N, const double ALPHA, const void
          * X, const int INCX, void * A)

 -- Function: void cblas_zher2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * X, const int INCX, const void * Y, const int INCY,
          void * A, const int LDA)

 -- Function: void cblas_zhpr2 (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const int N, const void * ALPHA, const
          void * X, const int INCX, const void * Y, const int INCY,
          void * AP)


File: gsl-ref.info,  Node: Level 3 CBLAS Functions,  Next: GSL CBLAS Examples,  Prev: Level 2 CBLAS Functions,  Up: GSL CBLAS Library

D.3 Level 3
===========

 -- Function: void cblas_sgemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE
          TRANSB, const int M, const int N, const int K, const float
          ALPHA, const float * A, const int LDA, const float * B, const
          int LDB, const float BETA, float * C, const int LDC)

 -- Function: void cblas_ssymm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const float ALPHA, const float * A, const int
          LDA, const float * B, const int LDB, const float BETA, float
          * C, const int LDC)

 -- Function: void cblas_ssyrk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const float ALPHA, const float * A, const
          int LDA, const float BETA, float * C, const int LDC)

 -- Function: void cblas_ssyr2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const float ALPHA, const float * A, const
          int LDA, const float * B, const int LDB, const float BETA,
          float * C, const int LDC)

 -- Function: void cblas_strmm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const float ALPHA, const float * A, const int
          LDA, float * B, const int LDB)

 -- Function: void cblas_strsm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const float ALPHA, const float * A, const int
          LDA, float * B, const int LDB)

 -- Function: void cblas_dgemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE
          TRANSB, const int M, const int N, const int K, const double
          ALPHA, const double * A, const int LDA, const double * B,
          const int LDB, const double BETA, double * C, const int LDC)

 -- Function: void cblas_dsymm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const double ALPHA, const double * A, const
          int LDA, const double * B, const int LDB, const double BETA,
          double * C, const int LDC)

 -- Function: void cblas_dsyrk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const double ALPHA, const double * A,
          const int LDA, const double BETA, double * C, const int LDC)

 -- Function: void cblas_dsyr2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const double ALPHA, const double * A,
          const int LDA, const double * B, const int LDB, const double
          BETA, double * C, const int LDC)

 -- Function: void cblas_dtrmm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const double ALPHA, const double * A, const
          int LDA, double * B, const int LDB)

 -- Function: void cblas_dtrsm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const double ALPHA, const double * A, const
          int LDA, double * B, const int LDB)

 -- Function: void cblas_cgemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE
          TRANSB, const int M, const int N, const int K, const void *
          ALPHA, const void * A, const int LDA, const void * B, const
          int LDB, const void * BETA, void * C, const int LDC)

 -- Function: void cblas_csymm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, const void * B, const int LDB, const void * BETA, void *
          C, const int LDC)

 -- Function: void cblas_csyrk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * BETA, void * C, const int LDC)

 -- Function: void cblas_csyr2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * B, const int LDB, const void * BETA,
          void * C, const int LDC)

 -- Function: void cblas_ctrmm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, void * B, const int LDB)

 -- Function: void cblas_ctrsm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, void * B, const int LDB)

 -- Function: void cblas_zgemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_TRANSPOSE TRANSA, const enum CBLAS_TRANSPOSE
          TRANSB, const int M, const int N, const int K, const void *
          ALPHA, const void * A, const int LDA, const void * B, const
          int LDB, const void * BETA, void * C, const int LDC)

 -- Function: void cblas_zsymm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, const void * B, const int LDB, const void * BETA, void *
          C, const int LDC)

 -- Function: void cblas_zsyrk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * BETA, void * C, const int LDC)

 -- Function: void cblas_zsyr2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * B, const int LDB, const void * BETA,
          void * C, const int LDC)

 -- Function: void cblas_ztrmm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, void * B, const int LDB)

 -- Function: void cblas_ztrsm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const enum
          CBLAS_TRANSPOSE TRANSA, const enum CBLAS_DIAG DIAG, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, void * B, const int LDB)

 -- Function: void cblas_chemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, const void * B, const int LDB, const void * BETA, void *
          C, const int LDC)

 -- Function: void cblas_cherk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const float ALPHA, const void * A, const
          int LDA, const float BETA, void * C, const int LDC)

 -- Function: void cblas_cher2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * B, const int LDB, const float BETA,
          void * C, const int LDC)

 -- Function: void cblas_zhemm (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_SIDE SIDE, const enum CBLAS_UPLO UPLO, const int
          M, const int N, const void * ALPHA, const void * A, const int
          LDA, const void * B, const int LDB, const void * BETA, void *
          C, const int LDC)

 -- Function: void cblas_zherk (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const double ALPHA, const void * A, const
          int LDA, const double BETA, void * C, const int LDC)

 -- Function: void cblas_zher2k (const enum CBLAS_ORDER ORDER, const
          enum CBLAS_UPLO UPLO, const enum CBLAS_TRANSPOSE TRANS, const
          int N, const int K, const void * ALPHA, const void * A, const
          int LDA, const void * B, const int LDB, const double BETA,
          void * C, const int LDC)

 -- Function: void cblas_xerbla (int P, const char * ROUT, const char *
          FORM, ...)


File: gsl-ref.info,  Node: GSL CBLAS Examples,  Prev: Level 3 CBLAS Functions,  Up: GSL CBLAS Library

D.4 Examples
============

The following program computes the product of two matrices using the
Level-3 BLAS function SGEMM,

     [ 0.11 0.12 0.13 ]  [ 1011 1012 ]     [ 367.76 368.12 ]
     [ 0.21 0.22 0.23 ]  [ 1021 1022 ]  =  [ 674.06 674.72 ]
                         [ 1031 1032 ]

The matrices are stored in row major order but could be stored in column
major order if the first argument of the call to `cblas_sgemm' was
changed to `CblasColMajor'.

     #include <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 ]


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'


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 2, June 1991
     Copyright (C) 1989, 1991 Free Software Foundation, Inc.
     59 Temple Place - Suite 330, Boston, MA  02111-1307, USA

     Everyone is permitted to copy and distribute verbatim copies of this
     license document, but changing it is not allowed.

Preamble
========

The licenses for most software are designed to take away your freedom
to share and change it.  By contrast, the GNU General Public License is
intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users.  This
General Public License applies to most of the Free Software
Foundation's software and to any other program whose authors commit to
using it.  (Some other Free Software Foundation software is covered by
the GNU Library General Public License instead.)  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
this service 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 make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.

   For example, if you distribute copies of such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have.  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.

   We protect your rights with two steps: (1) copyright the software,
and (2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.

   Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software.  If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.

   Finally, any free program is threatened constantly by software
patents.  We wish to avoid the danger that redistributors of a free
program will individually obtain patent licenses, in effect making the
program proprietary.  To prevent this, we have made it clear that any
patent must be licensed for everyone's free use or not licensed at all.

   The precise terms and conditions for copying, distribution and
modification follow.

    TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
  0. This License applies to any program or other work which contains a
     notice placed by the copyright holder saying it may be distributed
     under the terms of this General Public License.  The "Program",
     below, refers to any such program or work, and a "work based on
     the Program" means either the Program or any derivative work under
     copyright law: that is to say, a work containing the Program or a
     portion of it, either verbatim or with modifications and/or
     translated into another language.  (Hereinafter, translation is
     included without limitation in the term "modification".)  Each
     licensee is addressed as "you".

     Activities other than copying, distribution and modification are
     not covered by this License; they are outside its scope.  The act
     of running the Program is not restricted, and the output from the
     Program is covered only if its contents constitute a work based on
     the Program (independent of having been made by running the
     Program).  Whether that is true depends on what the Program does.

  1. You may copy and distribute 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 and disclaimer of warranty; keep intact all the
     notices that refer to this License and to the absence of any
     warranty; and give any other recipients of the Program a copy of
     this License along with the Program.

     You may charge a fee for the physical act of transferring a copy,
     and you may at your option offer warranty protection in exchange
     for a fee.

  2. You may modify your copy or copies of the Program or any portion
     of it, thus forming a work based on the Program, and copy and
     distribute such modifications or work under the terms of Section 1
     above, provided that you also meet all of these conditions:

       a. You must cause the modified files to carry prominent notices
          stating that you changed the files and the date of any change.

       b. You must cause any work that you distribute or publish, that
          in whole or in part contains or is derived from the Program
          or any part thereof, to be licensed as a whole at no charge
          to all third parties under the terms of this License.

       c. If the modified program normally reads commands interactively
          when run, you must cause it, when started running for such
          interactive use in the most ordinary way, to print or display
          an announcement including an appropriate copyright notice and
          a notice that there is no warranty (or else, saying that you
          provide a warranty) and that users may redistribute the
          program under these conditions, and telling the user how to
          view a copy of this License.  (Exception: if the Program
          itself is interactive but does not normally print such an
          announcement, your work based on the Program is not required
          to print an announcement.)

     These requirements apply to the modified work as a whole.  If
     identifiable sections of that work are not derived from the
     Program, and can be reasonably considered independent and separate
     works in themselves, then this License, and its terms, do not
     apply to those sections when you distribute them as separate
     works.  But when you distribute the same sections as part of a
     whole which is a work based on the Program, the distribution of
     the whole must be on the terms of this License, whose permissions
     for other licensees extend to the entire whole, and thus to each
     and every part regardless of who wrote it.

     Thus, it is not the intent of this section to claim rights or
     contest your rights to work written entirely by you; rather, the
     intent is to exercise the right to control the distribution of
     derivative or collective works based on the Program.

     In addition, mere aggregation of another work not based on the
     Program with the Program (or with a work based on the Program) on
     a volume of a storage or distribution medium does not bring the
     other work under the scope of this License.

  3. You may copy and distribute the Program (or a work based on it,
     under Section 2) in object code or executable form under the terms
     of Sections 1 and 2 above provided that you also do one of the
     following:

       a. Accompany it with the complete corresponding machine-readable
          source code, which must be distributed under the terms of
          Sections 1 and 2 above on a medium customarily used for
          software interchange; or,

       b. Accompany it with a written offer, valid for at least three
          years, to give any third party, for a charge no more than your
          cost of physically performing source distribution, a complete
          machine-readable copy of the corresponding source code, to be
          distributed under the terms of Sections 1 and 2 above on a
          medium customarily used for software interchange; or,

       c. Accompany it with the information you received as to the offer
          to distribute corresponding source code.  (This alternative is
          allowed only for noncommercial distribution and only if you
          received the program in object code or executable form with
          such an offer, in accord with Subsection b above.)

     The source code for a work means the preferred form of the work for
     making modifications to it.  For an executable work, complete
     source code means all the source code for all modules it contains,
     plus any associated interface definition files, plus the scripts
     used to control compilation and installation of the executable.
     However, as a special exception, the source code distributed need
     not include anything that is normally distributed (in either
     source or binary form) with the major components (compiler,
     kernel, and so on) of the operating system on which the executable
     runs, unless that component itself accompanies the executable.

     If distribution of executable or object code is made by offering
     access to copy from a designated place, then offering equivalent
     access to copy the source code from the same place counts as
     distribution of the source code, even though third parties are not
     compelled to copy the source along with the object code.

  4. You may not copy, modify, sublicense, or distribute the Program
     except as expressly provided under this License.  Any attempt
     otherwise to copy, modify, sublicense or distribute the Program 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.

  5. You are not required to accept this License, since you have not
     signed it.  However, nothing else grants you permission to modify
     or distribute the Program or its derivative works.  These actions
     are prohibited by law if you do not accept this License.
     Therefore, by modifying or distributing the Program (or any work
     based on the Program), you indicate your acceptance of this
     License to do so, and all its terms and conditions for copying,
     distributing or modifying the Program or works based on it.

  6. Each time you redistribute the Program (or any work based on the
     Program), the recipient automatically receives a license from the
     original licensor to copy, distribute or modify the Program
     subject to these terms and conditions.  You may not impose any
     further restrictions on the recipients' exercise of the rights
     granted herein.  You are not responsible for enforcing compliance
     by third parties to this License.

  7. If, as a consequence of a court judgment or allegation of patent
     infringement or for any other reason (not limited to patent
     issues), 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 distribute so as to satisfy simultaneously
     your obligations under this License and any other pertinent
     obligations, then as a consequence you may not distribute the
     Program at all.  For example, if a patent license would not permit
     royalty-free redistribution of the Program by all those who
     receive copies directly or indirectly through you, then the only
     way you could satisfy both it and this License would be to refrain
     entirely from distribution of the Program.

     If any portion of this section is held invalid or unenforceable
     under any particular circumstance, the balance of the section is
     intended to apply and the section as a whole is intended to apply
     in other circumstances.

     It is not the purpose of this section to induce you to infringe any
     patents or other property right claims or to contest validity of
     any such claims; this section has the sole purpose of protecting
     the integrity of the free software distribution system, which is
     implemented by public license practices.  Many people have made
     generous contributions to the wide range of software distributed
     through that system in reliance on consistent application of that
     system; it is up to the author/donor to decide if he or she is
     willing to distribute software through any other system and a
     licensee cannot impose that choice.

     This section is intended to make thoroughly clear what is believed
     to be a consequence of the rest of this License.

  8. If the distribution and/or use of the Program is restricted in
     certain countries either by patents or by copyrighted interfaces,
     the original copyright holder who places the Program under this
     License may add an explicit geographical distribution limitation
     excluding those countries, so that distribution is permitted only
     in or among countries not thus excluded.  In such case, this
     License incorporates the limitation as if written in the body of
     this License.

  9. The Free Software Foundation may publish revised and/or new
     versions of the 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 a version number of this License which applies
     to it and "any later version", you have the option of following
     the terms and conditions either of that version or of any later
     version published by the Free Software Foundation.  If the Program
     does not specify a version number of this License, you may choose
     any version ever published by the Free Software Foundation.

 10. If you wish to incorporate parts of the Program into other free
     programs whose distribution conditions are different, write to the
     author to ask for permission.  For software which is copyrighted
     by the Free Software Foundation, write to the Free Software
     Foundation; we sometimes make exceptions for this.  Our decision
     will be guided by the two goals of preserving the free status of
     all derivatives of our free software and of promoting the sharing
     and reuse of software generally.

                                NO WARRANTY
 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, 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.

 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
     WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY
     MODIFY AND/OR REDISTRIBUTE 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.

                      END OF TERMS AND CONDITIONS
Appendix: 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
convey 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) YYYY  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 2 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, write to the Free
     Software Foundation, Inc., 59 Temple Place - Suite 330,
     Boston, MA 02111-1307, USA.

   Also add information on how to contact you by electronic and paper
mail.

   If the program is interactive, make it output a short notice like
this when it starts in an interactive mode:

     Gnomovision version 69, Copyright (C) 19YY NAME OF AUTHOR
     Gnomovision 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, the
commands you use may be called something other than `show w' and `show
c'; they could even be mouse-clicks or menu items--whatever suits your
program.

   You should also get your employer (if you work as a programmer) or
your school, if any, to sign a "copyright disclaimer" for the program,
if necessary.  Here is a sample; alter the names:

     Yoyodyne, Inc., hereby disclaims all copyright interest in
     the program `Gnomovision' (which makes passes at compilers)
     written by James Hacker.

     SIGNATURE OF TY COON, 1 April 1989
     Ty Coon, President of Vice

   This 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 Library General Public License instead of this License.


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.