1 This is gsl-ref.info, produced by makeinfo version 4.8 from
4 INFO-DIR-SECTION Scientific software
6 * gsl-ref: (gsl-ref). GNU Scientific Library - Reference
9 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
10 2005, 2006, 2007 The GSL Team.
12 Permission is granted to copy, distribute and/or modify this document
13 under the terms of the GNU Free Documentation License, Version 1.2 or
14 any later version published by the Free Software Foundation; with the
15 Invariant Sections being "GNU General Public License" and "Free Software
16 Needs Free Documentation", the Front-Cover text being "A GNU Manual",
17 and with the Back-Cover Text being (a) (see below). A copy of the
18 license is included in the section entitled "GNU Free Documentation
21 (a) The Back-Cover Text is: "You have the freedom to copy and modify
25 File: gsl-ref.info, Node: Sorting vectors, Next: Selecting the k smallest or largest elements, Prev: Sorting objects, Up: Sorting
30 The following functions will sort the elements of an array or vector,
31 either directly or indirectly. They are defined for all real and
32 integer types using the normal suffix rules. For example, the `float'
33 versions of the array functions are `gsl_sort_float' and
34 `gsl_sort_float_index'. The corresponding vector functions are
35 `gsl_sort_vector_float' and `gsl_sort_vector_float_index'. The
36 prototypes are available in the header files `gsl_sort_float.h'
37 `gsl_sort_vector_float.h'. The complete set of prototypes can be
38 included using the header files `gsl_sort.h' and `gsl_sort_vector.h'.
40 There are no functions for sorting complex arrays or vectors, since
41 the ordering of complex numbers is not uniquely defined. To sort a
42 complex vector by magnitude compute a real vector containing the
43 magnitudes of the complex elements, and sort this vector indirectly.
44 The resulting index gives the appropriate ordering of the original
47 -- Function: void gsl_sort (double * DATA, size_t STRIDE, size_t N)
48 This function sorts the N elements of the array DATA with stride
49 STRIDE into ascending numerical order.
51 -- Function: void gsl_sort_vector (gsl_vector * V)
52 This function sorts the elements of the vector V into ascending
55 -- Function: void gsl_sort_index (size_t * P, const double * DATA,
56 size_t STRIDE, size_t N)
57 This function indirectly sorts the N elements of the array DATA
58 with stride STRIDE into ascending order, storing the resulting
59 permutation in P. The array P must be allocated with a sufficient
60 length to store the N elements of the permutation. The elements
61 of P give the index of the array element which would have been
62 stored in that position if the array had been sorted in place.
63 The array DATA is not changed.
65 -- Function: int gsl_sort_vector_index (gsl_permutation * P, const
67 This function indirectly sorts the elements of the vector V into
68 ascending order, storing the resulting permutation in P. The
69 elements of P give the index of the vector element which would
70 have been stored in that position if the vector had been sorted in
71 place. The first element of P gives the index of the least element
72 in V, and the last element of P gives the index of the greatest
73 element in V. The vector V is not changed.
76 File: gsl-ref.info, Node: Selecting the k smallest or largest elements, Next: Computing the rank, Prev: Sorting vectors, Up: Sorting
78 11.3 Selecting the k smallest or largest elements
79 =================================================
81 The functions described in this section select the k smallest or
82 largest elements of a data set of size N. The routines use an O(kN)
83 direct insertion algorithm which is suited to subsets that are small
84 compared with the total size of the dataset. For example, the routines
85 are useful for selecting the 10 largest values from one million data
86 points, but not for selecting the largest 100,000 values. If the
87 subset is a significant part of the total dataset it may be faster to
88 sort all the elements of the dataset directly with an O(N \log N)
89 algorithm and obtain the smallest or largest values that way.
91 -- Function: int gsl_sort_smallest (double * DEST, size_t K, const
92 double * SRC, size_t STRIDE, size_t N)
93 This function copies the K smallest elements of the array SRC, of
94 size N and stride STRIDE, in ascending numerical order into the
95 array DEST. The size K of the subset must be less than or equal
96 to N. The data SRC is not modified by this operation.
98 -- Function: int gsl_sort_largest (double * DEST, size_t K, const
99 double * SRC, size_t STRIDE, size_t N)
100 This function copies the K largest elements of the array SRC, of
101 size N and stride STRIDE, in descending numerical order into the
102 array DEST. K must be less than or equal to N. The data SRC is not
103 modified by this operation.
105 -- Function: int gsl_sort_vector_smallest (double * DEST, size_t K,
106 const gsl_vector * V)
107 -- Function: int gsl_sort_vector_largest (double * DEST, size_t K,
108 const gsl_vector * V)
109 These functions copy the K smallest or largest elements of the
110 vector V into the array DEST. K must be less than or equal to the
111 length of the vector V.
113 The following functions find the indices of the k smallest or
114 largest elements of a dataset,
116 -- Function: int gsl_sort_smallest_index (size_t * P, size_t K, const
117 double * SRC, size_t STRIDE, size_t N)
118 This function stores the indices of the K smallest elements of the
119 array SRC, of size N and stride STRIDE, in the array P. The
120 indices are chosen so that the corresponding data is in ascending
121 numerical order. K must be less than or equal to N. The data SRC
122 is not modified by this operation.
124 -- Function: int gsl_sort_largest_index (size_t * P, size_t K, const
125 double * SRC, size_t STRIDE, size_t N)
126 This function stores the indices of the K largest elements of the
127 array SRC, of size N and stride STRIDE, in the array P. The
128 indices are chosen so that the corresponding data is in descending
129 numerical order. K must be less than or equal to N. The data SRC
130 is not modified by this operation.
132 -- Function: int gsl_sort_vector_smallest_index (size_t * P, size_t K,
133 const gsl_vector * V)
134 -- Function: int gsl_sort_vector_largest_index (size_t * P, size_t K,
135 const gsl_vector * V)
136 These functions store the indices of the K smallest or largest
137 elements of the vector V in the array P. K must be less than or
138 equal to the length of the vector V.
141 File: gsl-ref.info, Node: Computing the rank, Next: Sorting Examples, Prev: Selecting the k smallest or largest elements, Up: Sorting
143 11.4 Computing the rank
144 =======================
146 The "rank" of an element is its order in the sorted data. The rank is
147 the inverse of the index permutation, P. It can be computed using the
150 for (i = 0; i < p->size; i++)
152 size_t pi = p->data[i];
156 This can be computed directly from the function
157 `gsl_permutation_inverse(rank,p)'.
159 The following function will print the rank of each element of the
163 print_rank (gsl_vector * v)
167 gsl_permutation * perm = gsl_permutation_alloc(n);
168 gsl_permutation * rank = gsl_permutation_alloc(n);
170 gsl_sort_vector_index (perm, v);
171 gsl_permutation_inverse (rank, perm);
173 for (i = 0; i < n; i++)
175 double vi = gsl_vector_get(v, i);
176 printf ("element = %d, value = %g, rank = %d\n",
177 i, vi, rank->data[i]);
180 gsl_permutation_free (perm);
181 gsl_permutation_free (rank);
185 File: gsl-ref.info, Node: Sorting Examples, Next: Sorting References and Further Reading, Prev: Computing the rank, Up: Sorting
190 The following example shows how to use the permutation P to print the
191 elements of the vector V in ascending order,
193 gsl_sort_vector_index (p, v);
195 for (i = 0; i < v->size; i++)
197 double vpi = gsl_vector_get (v, p->data[i]);
198 printf ("order = %d, value = %g\n", i, vpi);
201 The next example uses the function `gsl_sort_smallest' to select the 5
202 smallest numbers from 100000 uniform random variates stored in an array,
204 #include <gsl/gsl_rng.h>
205 #include <gsl/gsl_sort_double.h>
210 const gsl_rng_type * T;
213 size_t i, k = 5, N = 100000;
215 double * x = malloc (N * sizeof(double));
216 double * small = malloc (k * sizeof(double));
221 r = gsl_rng_alloc (T);
223 for (i = 0; i < N; i++)
225 x[i] = gsl_rng_uniform(r);
228 gsl_sort_smallest (small, k, x, 1, N);
230 printf ("%d smallest values from %d\n", k, N);
232 for (i = 0; i < k; i++)
234 printf ("%d: %.18f\n", i, small[i]);
242 The output lists the 5 smallest values, in ascending order,
245 5 smallest values from 100000
246 0: 0.000003489200025797
247 1: 0.000008199829608202
248 2: 0.000008953968062997
249 3: 0.000010712770745158
250 4: 0.000033531803637743
253 File: gsl-ref.info, Node: Sorting References and Further Reading, Prev: Sorting Examples, Up: Sorting
255 11.6 References and Further Reading
256 ===================================
258 The subject of sorting is covered extensively in Knuth's `Sorting and
261 Donald E. Knuth, `The Art of Computer Programming: Sorting and
262 Searching' (Vol 3, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850.
264 The Heapsort algorithm is described in the following book,
266 Robert Sedgewick, `Algorithms in C', Addison-Wesley, ISBN
270 File: gsl-ref.info, Node: BLAS Support, Next: Linear Algebra, Prev: Sorting, Up: Top
275 The Basic Linear Algebra Subprograms (BLAS) define a set of fundamental
276 operations on vectors and matrices which can be used to create optimized
277 higher-level linear algebra functionality.
279 The library provides a low-level layer which corresponds directly to
280 the C-language BLAS standard, referred to here as "CBLAS", and a
281 higher-level interface for operations on GSL vectors and matrices.
282 Users who are interested in simple operations on GSL vector and matrix
283 objects should use the high-level layer, which is declared in the file
284 `gsl_blas.h'. This should satisfy the needs of most users. Note that
285 GSL matrices are implemented using dense-storage so the interface only
286 includes the corresponding dense-storage BLAS functions. The full BLAS
287 functionality for band-format and packed-format matrices is available
288 through the low-level CBLAS interface.
290 The interface for the `gsl_cblas' layer is specified in the file
291 `gsl_cblas.h'. This interface corresponds to the BLAS Technical
292 Forum's draft standard for the C interface to legacy BLAS
293 implementations. Users who have access to other conforming CBLAS
294 implementations can use these in place of the version provided by the
295 library. Note that users who have only a Fortran BLAS library can use
296 a CBLAS conformant wrapper to convert it into a CBLAS library. A
297 reference CBLAS wrapper for legacy Fortran implementations exists as
298 part of the draft CBLAS standard and can be obtained from Netlib. The
299 complete set of CBLAS functions is listed in an appendix (*note GSL
302 There are three levels of BLAS operations,
305 Vector operations, e.g. y = \alpha x + y
308 Matrix-vector operations, e.g. y = \alpha A x + \beta y
311 Matrix-matrix operations, e.g. C = \alpha A B + C
313 Each routine has a name which specifies the operation, the type of
314 matrices involved and their precisions. Some of the most common
315 operations and their names are given below,
318 scalar product, x^T y
321 vector sum, \alpha x + y
324 matrix-vector product, A x
327 matrix-vector solve, inv(A) x
330 matrix-matrix product, A B
333 matrix-matrix solve, inv(A) B
335 The types of matrices are,
370 Each operation is defined for four precisions,
384 Thus, for example, the name SGEMM stands for "single-precision general
385 matrix-matrix multiply" and ZGEMM stands for "double-precision complex
386 matrix-matrix multiply".
390 * GSL BLAS Interface::
392 * BLAS References and Further Reading::
395 File: gsl-ref.info, Node: GSL BLAS Interface, Next: BLAS Examples, Up: BLAS Support
397 12.1 GSL BLAS Interface
398 =======================
400 GSL provides dense vector and matrix objects, based on the relevant
401 built-in types. The library provides an interface to the BLAS
402 operations which apply to these objects. The interface to this
403 functionality is given in the file `gsl_blas.h'.
407 * Level 1 GSL BLAS Interface::
408 * Level 2 GSL BLAS Interface::
409 * Level 3 GSL BLAS Interface::
412 File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS Interface, Up: GSL BLAS Interface
417 -- Function: int gsl_blas_sdsdot (float ALPHA, const gsl_vector_float
418 * X, const gsl_vector_float * Y, float * RESULT)
419 This function computes the sum \alpha + x^T y for the vectors X
420 and Y, returning the result in RESULT.
422 -- Function: int gsl_blas_sdot (const gsl_vector_float * X, const
423 gsl_vector_float * Y, float * RESULT)
424 -- Function: int gsl_blas_dsdot (const gsl_vector_float * X, const
425 gsl_vector_float * Y, double * RESULT)
426 -- Function: int gsl_blas_ddot (const gsl_vector * X, const gsl_vector
427 * Y, double * RESULT)
428 These functions compute the scalar product x^T y for the vectors X
429 and Y, returning the result in RESULT.
431 -- Function: int gsl_blas_cdotu (const gsl_vector_complex_float * X,
432 const gsl_vector_complex_float * Y, gsl_complex_float * DOTU)
433 -- Function: int gsl_blas_zdotu (const gsl_vector_complex * X, const
434 gsl_vector_complex * Y, gsl_complex * DOTU)
435 These functions compute the complex scalar product x^T y for the
436 vectors X and Y, returning the result in RESULT
438 -- Function: int gsl_blas_cdotc (const gsl_vector_complex_float * X,
439 const gsl_vector_complex_float * Y, gsl_complex_float * DOTC)
440 -- Function: int gsl_blas_zdotc (const gsl_vector_complex * X, const
441 gsl_vector_complex * Y, gsl_complex * DOTC)
442 These functions compute the complex conjugate scalar product x^H y
443 for the vectors X and Y, returning the result in RESULT
445 -- Function: float gsl_blas_snrm2 (const gsl_vector_float * X)
446 -- Function: double gsl_blas_dnrm2 (const gsl_vector * X)
447 These functions compute the Euclidean norm ||x||_2 = \sqrt {\sum
448 x_i^2} of the vector X.
450 -- Function: float gsl_blas_scnrm2 (const gsl_vector_complex_float * X)
451 -- Function: double gsl_blas_dznrm2 (const gsl_vector_complex * X)
452 These functions compute the Euclidean norm of the complex vector X,
454 ||x||_2 = \sqrt {\sum (\Re(x_i)^2 + \Im(x_i)^2)}.
456 -- Function: float gsl_blas_sasum (const gsl_vector_float * X)
457 -- Function: double gsl_blas_dasum (const gsl_vector * X)
458 These functions compute the absolute sum \sum |x_i| of the
459 elements of the vector X.
461 -- Function: float gsl_blas_scasum (const gsl_vector_complex_float * X)
462 -- Function: double gsl_blas_dzasum (const gsl_vector_complex * X)
463 These functions compute the sum of the magnitudes of the real and
464 imaginary parts of the complex vector X, \sum |\Re(x_i)| +
467 -- Function: CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * X)
468 -- Function: CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * X)
469 -- Function: CBLAS_INDEX_t gsl_blas_icamax (const
470 gsl_vector_complex_float * X)
471 -- Function: CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex *
473 These functions return the index of the largest element of the
474 vector X. The largest element is determined by its absolute
475 magnitude for real vectors and by the sum of the magnitudes of the
476 real and imaginary parts |\Re(x_i)| + |\Im(x_i)| for complex
477 vectors. If the largest value occurs several times then the index
478 of the first occurrence is returned.
480 -- Function: int gsl_blas_sswap (gsl_vector_float * X,
481 gsl_vector_float * Y)
482 -- Function: int gsl_blas_dswap (gsl_vector * X, gsl_vector * Y)
483 -- Function: int gsl_blas_cswap (gsl_vector_complex_float * X,
484 gsl_vector_complex_float * Y)
485 -- Function: int gsl_blas_zswap (gsl_vector_complex * X,
486 gsl_vector_complex * Y)
487 These functions exchange the elements of the vectors X and Y.
489 -- Function: int gsl_blas_scopy (const gsl_vector_float * X,
490 gsl_vector_float * Y)
491 -- Function: int gsl_blas_dcopy (const gsl_vector * X, gsl_vector * Y)
492 -- Function: int gsl_blas_ccopy (const gsl_vector_complex_float * X,
493 gsl_vector_complex_float * Y)
494 -- Function: int gsl_blas_zcopy (const gsl_vector_complex * X,
495 gsl_vector_complex * Y)
496 These functions copy the elements of the vector X into the vector
499 -- Function: int gsl_blas_saxpy (float ALPHA, const gsl_vector_float *
500 X, gsl_vector_float * Y)
501 -- Function: int gsl_blas_daxpy (double ALPHA, const gsl_vector * X,
503 -- Function: int gsl_blas_caxpy (const gsl_complex_float ALPHA, const
504 gsl_vector_complex_float * X, gsl_vector_complex_float * Y)
505 -- Function: int gsl_blas_zaxpy (const gsl_complex ALPHA, const
506 gsl_vector_complex * X, gsl_vector_complex * Y)
507 These functions compute the sum y = \alpha x + y for the vectors X
510 -- Function: void gsl_blas_sscal (float ALPHA, gsl_vector_float * X)
511 -- Function: void gsl_blas_dscal (double ALPHA, gsl_vector * X)
512 -- Function: void gsl_blas_cscal (const gsl_complex_float ALPHA,
513 gsl_vector_complex_float * X)
514 -- Function: void gsl_blas_zscal (const gsl_complex ALPHA,
515 gsl_vector_complex * X)
516 -- Function: void gsl_blas_csscal (float ALPHA,
517 gsl_vector_complex_float * X)
518 -- Function: void gsl_blas_zdscal (double ALPHA, gsl_vector_complex *
520 These functions rescale the vector X by the multiplicative factor
523 -- Function: int gsl_blas_srotg (float A[], float B[], float C[],
525 -- Function: int gsl_blas_drotg (double A[], double B[], double C[],
527 These functions compute a Givens rotation (c,s) which zeroes the
530 [ c s ] [ a ] = [ r ]
533 The variables A and B are overwritten by the routine.
535 -- Function: int gsl_blas_srot (gsl_vector_float * X, gsl_vector_float
536 * Y, float C, float S)
537 -- Function: int gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const
538 double C, const double S)
539 These functions apply a Givens rotation (x', y') = (c x + s y, -s
540 x + c y) to the vectors X, Y.
542 -- Function: int gsl_blas_srotmg (float D1[], float D2[], float B1[],
544 -- Function: int gsl_blas_drotmg (double D1[], double D2[], double
545 B1[], double B2, double P[])
546 These functions compute a modified Givens transformation. The
547 modified Givens transformation is defined in the original Level-1
548 BLAS specification, given in the references.
550 -- Function: int gsl_blas_srotm (gsl_vector_float * X,
551 gsl_vector_float * Y, const float P[])
552 -- Function: int gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const
554 These functions apply a modified Givens transformation.
557 File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS Interface, Prev: Level 1 GSL BLAS Interface, Up: GSL BLAS Interface
562 -- Function: int gsl_blas_sgemv (CBLAS_TRANSPOSE_t TRANSA, float
563 ALPHA, const gsl_matrix_float * A, const gsl_vector_float *
564 X, float BETA, gsl_vector_float * Y)
565 -- Function: int gsl_blas_dgemv (CBLAS_TRANSPOSE_t TRANSA, double
566 ALPHA, const gsl_matrix * A, const gsl_vector * X, double
567 BETA, gsl_vector * Y)
568 -- Function: int gsl_blas_cgemv (CBLAS_TRANSPOSE_t TRANSA, const
569 gsl_complex_float ALPHA, const gsl_matrix_complex_float * A,
570 const gsl_vector_complex_float * X, const gsl_complex_float
571 BETA, gsl_vector_complex_float * Y)
572 -- Function: int gsl_blas_zgemv (CBLAS_TRANSPOSE_t TRANSA, const
573 gsl_complex ALPHA, const gsl_matrix_complex * A, const
574 gsl_vector_complex * X, const gsl_complex BETA,
575 gsl_vector_complex * Y)
576 These functions compute the matrix-vector product and sum y =
577 \alpha op(A) x + \beta y, where op(A) = A, A^T, A^H for TRANSA =
578 `CblasNoTrans', `CblasTrans', `CblasConjTrans'.
580 -- Function: int gsl_blas_strmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
581 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_float * A,
582 gsl_vector_float * X)
583 -- Function: int gsl_blas_dtrmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
584 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix * A, gsl_vector *
586 -- Function: int gsl_blas_ctrmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
587 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float *
588 A, gsl_vector_complex_float * X)
589 -- Function: int gsl_blas_ztrmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
590 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex * A,
591 gsl_vector_complex * X)
592 These functions compute the matrix-vector product x = op(A) x for
593 the triangular matrix A, where op(A) = A, A^T, A^H for TRANSA =
594 `CblasNoTrans', `CblasTrans', `CblasConjTrans'. When UPLO is
595 `CblasUpper' then the upper triangle of A is used, and when UPLO
596 is `CblasLower' then the lower triangle of A is used. If DIAG is
597 `CblasNonUnit' then the diagonal of the matrix is used, but if
598 DIAG is `CblasUnit' then the diagonal elements of the matrix A are
599 taken as unity and are not referenced.
601 -- Function: int gsl_blas_strsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
602 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_float * A,
603 gsl_vector_float * X)
604 -- Function: int gsl_blas_dtrsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
605 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix * A, gsl_vector *
607 -- Function: int gsl_blas_ctrsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
608 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float *
609 A, gsl_vector_complex_float * X)
610 -- Function: int gsl_blas_ztrsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
611 TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex * A,
612 gsl_vector_complex * X)
613 These functions compute inv(op(A)) x for X, where op(A) = A, A^T,
614 A^H for TRANSA = `CblasNoTrans', `CblasTrans', `CblasConjTrans'.
615 When UPLO is `CblasUpper' then the upper triangle of A is used,
616 and when UPLO is `CblasLower' then the lower triangle of A is
617 used. If DIAG is `CblasNonUnit' then the diagonal of the matrix
618 is used, but if DIAG is `CblasUnit' then the diagonal elements of
619 the matrix A are taken as unity and are not referenced.
621 -- Function: int gsl_blas_ssymv (CBLAS_UPLO_t UPLO, float ALPHA, const
622 gsl_matrix_float * A, const gsl_vector_float * X, float BETA,
623 gsl_vector_float * Y)
624 -- Function: int gsl_blas_dsymv (CBLAS_UPLO_t UPLO, double ALPHA,
625 const gsl_matrix * A, const gsl_vector * X, double BETA,
627 These functions compute the matrix-vector product and sum y =
628 \alpha A x + \beta y for the symmetric matrix A. Since the matrix
629 A is symmetric only its upper half or lower half need to be
630 stored. When UPLO is `CblasUpper' then the upper triangle and
631 diagonal of A are used, and when UPLO is `CblasLower' then the
632 lower triangle and diagonal of A are used.
634 -- Function: int gsl_blas_chemv (CBLAS_UPLO_t UPLO, const
635 gsl_complex_float ALPHA, const gsl_matrix_complex_float * A,
636 const gsl_vector_complex_float * X, const gsl_complex_float
637 BETA, gsl_vector_complex_float * Y)
638 -- Function: int gsl_blas_zhemv (CBLAS_UPLO_t UPLO, const gsl_complex
639 ALPHA, const gsl_matrix_complex * A, const gsl_vector_complex
640 * X, const gsl_complex BETA, gsl_vector_complex * Y)
641 These functions compute the matrix-vector product and sum y =
642 \alpha A x + \beta y for the hermitian matrix A. Since the matrix
643 A is hermitian only its upper half or lower half need to be
644 stored. When UPLO is `CblasUpper' then the upper triangle and
645 diagonal of A are used, and when UPLO is `CblasLower' then the
646 lower triangle and diagonal of A are used. The imaginary elements
647 of the diagonal are automatically assumed to be zero and are not
650 -- Function: int gsl_blas_sger (float ALPHA, const gsl_vector_float *
651 X, const gsl_vector_float * Y, gsl_matrix_float * A)
652 -- Function: int gsl_blas_dger (double ALPHA, const gsl_vector * X,
653 const gsl_vector * Y, gsl_matrix * A)
654 -- Function: int gsl_blas_cgeru (const gsl_complex_float ALPHA, const
655 gsl_vector_complex_float * X, const gsl_vector_complex_float
656 * Y, gsl_matrix_complex_float * A)
657 -- Function: int gsl_blas_zgeru (const gsl_complex ALPHA, const
658 gsl_vector_complex * X, const gsl_vector_complex * Y,
659 gsl_matrix_complex * A)
660 These functions compute the rank-1 update A = \alpha x y^T + A of
663 -- Function: int gsl_blas_cgerc (const gsl_complex_float ALPHA, const
664 gsl_vector_complex_float * X, const gsl_vector_complex_float
665 * Y, gsl_matrix_complex_float * A)
666 -- Function: int gsl_blas_zgerc (const gsl_complex ALPHA, const
667 gsl_vector_complex * X, const gsl_vector_complex * Y,
668 gsl_matrix_complex * A)
669 These functions compute the conjugate rank-1 update A = \alpha x
670 y^H + A of the matrix A.
672 -- Function: int gsl_blas_ssyr (CBLAS_UPLO_t UPLO, float ALPHA, const
673 gsl_vector_float * X, gsl_matrix_float * A)
674 -- Function: int gsl_blas_dsyr (CBLAS_UPLO_t UPLO, double ALPHA, const
675 gsl_vector * X, gsl_matrix * A)
676 These functions compute the symmetric rank-1 update A = \alpha x
677 x^T + A of the symmetric matrix A. Since the matrix A is
678 symmetric only its upper half or lower half need to be stored.
679 When UPLO is `CblasUpper' then the upper triangle and diagonal of
680 A are used, and when UPLO is `CblasLower' then the lower triangle
681 and diagonal of A are used.
683 -- Function: int gsl_blas_cher (CBLAS_UPLO_t UPLO, float ALPHA, const
684 gsl_vector_complex_float * X, gsl_matrix_complex_float * A)
685 -- Function: int gsl_blas_zher (CBLAS_UPLO_t UPLO, double ALPHA, const
686 gsl_vector_complex * X, gsl_matrix_complex * A)
687 These functions compute the hermitian rank-1 update A = \alpha x
688 x^H + A of the hermitian matrix A. Since the matrix A is
689 hermitian only its upper half or lower half need to be stored.
690 When UPLO is `CblasUpper' then the upper triangle and diagonal of
691 A are used, and when UPLO is `CblasLower' then the lower triangle
692 and diagonal of A are used. The imaginary elements of the
693 diagonal are automatically set to zero.
695 -- Function: int gsl_blas_ssyr2 (CBLAS_UPLO_t UPLO, float ALPHA, const
696 gsl_vector_float * X, const gsl_vector_float * Y,
697 gsl_matrix_float * A)
698 -- Function: int gsl_blas_dsyr2 (CBLAS_UPLO_t UPLO, double ALPHA,
699 const gsl_vector * X, const gsl_vector * Y, gsl_matrix * A)
700 These functions compute the symmetric rank-2 update A = \alpha x
701 y^T + \alpha y x^T + A of the symmetric matrix A. Since the
702 matrix A is symmetric only its upper half or lower half need to be
703 stored. When UPLO is `CblasUpper' then the upper triangle and
704 diagonal of A are used, and when UPLO is `CblasLower' then the
705 lower triangle and diagonal of A are used.
707 -- Function: int gsl_blas_cher2 (CBLAS_UPLO_t UPLO, const
708 gsl_complex_float ALPHA, const gsl_vector_complex_float * X,
709 const gsl_vector_complex_float * Y, gsl_matrix_complex_float
711 -- Function: int gsl_blas_zher2 (CBLAS_UPLO_t UPLO, const gsl_complex
712 ALPHA, const gsl_vector_complex * X, const gsl_vector_complex
713 * Y, gsl_matrix_complex * A)
714 These functions compute the hermitian rank-2 update A = \alpha x
715 y^H + \alpha^* y x^H A of the hermitian matrix A. Since the
716 matrix A is hermitian only its upper half or lower half need to be
717 stored. When UPLO is `CblasUpper' then the upper triangle and
718 diagonal of A are used, and when UPLO is `CblasLower' then the
719 lower triangle and diagonal of A are used. The imaginary elements
720 of the diagonal are automatically set to zero.
723 File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS Interface, Up: GSL BLAS Interface
728 -- Function: int gsl_blas_sgemm (CBLAS_TRANSPOSE_t TRANSA,
729 CBLAS_TRANSPOSE_t TRANSB, float ALPHA, const gsl_matrix_float
730 * A, const gsl_matrix_float * B, float BETA, gsl_matrix_float
732 -- Function: int gsl_blas_dgemm (CBLAS_TRANSPOSE_t TRANSA,
733 CBLAS_TRANSPOSE_t TRANSB, double ALPHA, const gsl_matrix * A,
734 const gsl_matrix * B, double BETA, gsl_matrix * C)
735 -- Function: int gsl_blas_cgemm (CBLAS_TRANSPOSE_t TRANSA,
736 CBLAS_TRANSPOSE_t TRANSB, const gsl_complex_float ALPHA,
737 const gsl_matrix_complex_float * A, const
738 gsl_matrix_complex_float * B, const gsl_complex_float BETA,
739 gsl_matrix_complex_float * C)
740 -- Function: int gsl_blas_zgemm (CBLAS_TRANSPOSE_t TRANSA,
741 CBLAS_TRANSPOSE_t TRANSB, const gsl_complex ALPHA, const
742 gsl_matrix_complex * A, const gsl_matrix_complex * B, const
743 gsl_complex BETA, gsl_matrix_complex * C)
744 These functions compute the matrix-matrix product and sum C =
745 \alpha op(A) op(B) + \beta C where op(A) = A, A^T, A^H for TRANSA
746 = `CblasNoTrans', `CblasTrans', `CblasConjTrans' and similarly for
747 the parameter TRANSB.
749 -- Function: int gsl_blas_ssymm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
750 float ALPHA, const gsl_matrix_float * A, const
751 gsl_matrix_float * B, float BETA, gsl_matrix_float * C)
752 -- Function: int gsl_blas_dsymm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
753 double ALPHA, const gsl_matrix * A, const gsl_matrix * B,
754 double BETA, gsl_matrix * C)
755 -- Function: int gsl_blas_csymm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
756 const gsl_complex_float ALPHA, const gsl_matrix_complex_float
757 * A, const gsl_matrix_complex_float * B, const
758 gsl_complex_float BETA, gsl_matrix_complex_float * C)
759 -- Function: int gsl_blas_zsymm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
760 const gsl_complex ALPHA, const gsl_matrix_complex * A, const
761 gsl_matrix_complex * B, const gsl_complex BETA,
762 gsl_matrix_complex * C)
763 These functions compute the matrix-matrix product and sum C =
764 \alpha A B + \beta C for SIDE is `CblasLeft' and C = \alpha B A +
765 \beta C for SIDE is `CblasRight', where the matrix A is symmetric.
766 When UPLO is `CblasUpper' then the upper triangle and diagonal of
767 A are used, and when UPLO is `CblasLower' then the lower triangle
768 and diagonal of A are used.
770 -- Function: int gsl_blas_chemm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
771 const gsl_complex_float ALPHA, const gsl_matrix_complex_float
772 * A, const gsl_matrix_complex_float * B, const
773 gsl_complex_float BETA, gsl_matrix_complex_float * C)
774 -- Function: int gsl_blas_zhemm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
775 const gsl_complex ALPHA, const gsl_matrix_complex * A, const
776 gsl_matrix_complex * B, const gsl_complex BETA,
777 gsl_matrix_complex * C)
778 These functions compute the matrix-matrix product and sum C =
779 \alpha A B + \beta C for SIDE is `CblasLeft' and C = \alpha B A +
780 \beta C for SIDE is `CblasRight', where the matrix A is hermitian.
781 When UPLO is `CblasUpper' then the upper triangle and diagonal of
782 A are used, and when UPLO is `CblasLower' then the lower triangle
783 and diagonal of A are used. The imaginary elements of the
784 diagonal are automatically set to zero.
786 -- Function: int gsl_blas_strmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
787 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, float ALPHA,
788 const gsl_matrix_float * A, gsl_matrix_float * B)
789 -- Function: int gsl_blas_dtrmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
790 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, double ALPHA,
791 const gsl_matrix * A, gsl_matrix * B)
792 -- Function: int gsl_blas_ctrmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
793 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const
794 gsl_complex_float ALPHA, const gsl_matrix_complex_float * A,
795 gsl_matrix_complex_float * B)
796 -- Function: int gsl_blas_ztrmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
797 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const
798 gsl_complex ALPHA, const gsl_matrix_complex * A,
799 gsl_matrix_complex * B)
800 These functions compute the matrix-matrix product B = \alpha op(A)
801 B for SIDE is `CblasLeft' and B = \alpha B op(A) for SIDE is
802 `CblasRight'. The matrix A is triangular and op(A) = A, A^T, A^H
803 for TRANSA = `CblasNoTrans', `CblasTrans', `CblasConjTrans'. When
804 UPLO is `CblasUpper' then the upper triangle of A is used, and
805 when UPLO is `CblasLower' then the lower triangle of A is used.
806 If DIAG is `CblasNonUnit' then the diagonal of A is used, but if
807 DIAG is `CblasUnit' then the diagonal elements of the matrix A are
808 taken as unity and are not referenced.
810 -- Function: int gsl_blas_strsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
811 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, float ALPHA,
812 const gsl_matrix_float * A, gsl_matrix_float * B)
813 -- Function: int gsl_blas_dtrsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
814 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, double ALPHA,
815 const gsl_matrix * A, gsl_matrix * B)
816 -- Function: int gsl_blas_ctrsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
817 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const
818 gsl_complex_float ALPHA, const gsl_matrix_complex_float * A,
819 gsl_matrix_complex_float * B)
820 -- Function: int gsl_blas_ztrsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO,
821 CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const
822 gsl_complex ALPHA, const gsl_matrix_complex * A,
823 gsl_matrix_complex * B)
824 These functions compute the inverse-matrix matrix product B =
825 \alpha op(inv(A))B for SIDE is `CblasLeft' and B = \alpha B
826 op(inv(A)) for SIDE is `CblasRight'. The matrix A is triangular
827 and op(A) = A, A^T, A^H for TRANSA = `CblasNoTrans', `CblasTrans',
828 `CblasConjTrans'. When UPLO is `CblasUpper' then the upper
829 triangle of A is used, and when UPLO is `CblasLower' then the
830 lower triangle of A is used. If DIAG is `CblasNonUnit' then the
831 diagonal of A is used, but if DIAG is `CblasUnit' then the
832 diagonal elements of the matrix A are taken as unity and are not
835 -- Function: int gsl_blas_ssyrk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
836 TRANS, float ALPHA, const gsl_matrix_float * A, float BETA,
837 gsl_matrix_float * C)
838 -- Function: int gsl_blas_dsyrk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
839 TRANS, double ALPHA, const gsl_matrix * A, double BETA,
841 -- Function: int gsl_blas_csyrk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
842 TRANS, const gsl_complex_float ALPHA, const
843 gsl_matrix_complex_float * A, const gsl_complex_float BETA,
844 gsl_matrix_complex_float * C)
845 -- Function: int gsl_blas_zsyrk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
846 TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A,
847 const gsl_complex BETA, gsl_matrix_complex * C)
848 These functions compute a rank-k update of the symmetric matrix C,
849 C = \alpha A A^T + \beta C when TRANS is `CblasNoTrans' and C =
850 \alpha A^T A + \beta C when TRANS is `CblasTrans'. Since the
851 matrix C is symmetric only its upper half or lower half need to be
852 stored. When UPLO is `CblasUpper' then the upper triangle and
853 diagonal of C are used, and when UPLO is `CblasLower' then the
854 lower triangle and diagonal of C are used.
856 -- Function: int gsl_blas_cherk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
857 TRANS, float ALPHA, const gsl_matrix_complex_float * A, float
858 BETA, gsl_matrix_complex_float * C)
859 -- Function: int gsl_blas_zherk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
860 TRANS, double ALPHA, const gsl_matrix_complex * A, double
861 BETA, gsl_matrix_complex * C)
862 These functions compute a rank-k update of the hermitian matrix C,
863 C = \alpha A A^H + \beta C when TRANS is `CblasNoTrans' and C =
864 \alpha A^H A + \beta C when TRANS is `CblasTrans'. Since the
865 matrix C is hermitian only its upper half or lower half need to be
866 stored. When UPLO is `CblasUpper' then the upper triangle and
867 diagonal of C are used, and when UPLO is `CblasLower' then the
868 lower triangle and diagonal of C are used. The imaginary elements
869 of the diagonal are automatically set to zero.
871 -- Function: int gsl_blas_ssyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
872 TRANS, float ALPHA, const gsl_matrix_float * A, const
873 gsl_matrix_float * B, float BETA, gsl_matrix_float * C)
874 -- Function: int gsl_blas_dsyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
875 TRANS, double ALPHA, const gsl_matrix * A, const gsl_matrix *
876 B, double BETA, gsl_matrix * C)
877 -- Function: int gsl_blas_csyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
878 TRANS, const gsl_complex_float ALPHA, const
879 gsl_matrix_complex_float * A, const gsl_matrix_complex_float
880 * B, const gsl_complex_float BETA, gsl_matrix_complex_float *
882 -- Function: int gsl_blas_zsyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
883 TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A,
884 const gsl_matrix_complex * B, const gsl_complex BETA,
885 gsl_matrix_complex * C)
886 These functions compute a rank-2k update of the symmetric matrix C,
887 C = \alpha A B^T + \alpha B A^T + \beta C when TRANS is
888 `CblasNoTrans' and C = \alpha A^T B + \alpha B^T A + \beta C when
889 TRANS is `CblasTrans'. Since the matrix C is symmetric only its
890 upper half or lower half need to be stored. When UPLO is
891 `CblasUpper' then the upper triangle and diagonal of C are used,
892 and when UPLO is `CblasLower' then the lower triangle and diagonal
895 -- Function: int gsl_blas_cher2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
896 TRANS, const gsl_complex_float ALPHA, const
897 gsl_matrix_complex_float * A, const gsl_matrix_complex_float
898 * B, float BETA, gsl_matrix_complex_float * C)
899 -- Function: int gsl_blas_zher2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t
900 TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A,
901 const gsl_matrix_complex * B, double BETA, gsl_matrix_complex
903 These functions compute a rank-2k update of the hermitian matrix C,
904 C = \alpha A B^H + \alpha^* B A^H + \beta C when TRANS is
905 `CblasNoTrans' and C = \alpha A^H B + \alpha^* B^H A + \beta C when
906 TRANS is `CblasConjTrans'. Since the matrix C is hermitian only
907 its upper half or lower half need to be stored. When UPLO is
908 `CblasUpper' then the upper triangle and diagonal of C are used,
909 and when UPLO is `CblasLower' then the lower triangle and diagonal
910 of C are used. The imaginary elements of the diagonal are
911 automatically set to zero.
914 File: gsl-ref.info, Node: BLAS Examples, Next: BLAS References and Further Reading, Prev: GSL BLAS Interface, Up: BLAS Support
919 The following program computes the product of two matrices using the
920 Level-3 BLAS function DGEMM,
922 [ 0.11 0.12 0.13 ] [ 1011 1012 ] [ 367.76 368.12 ]
923 [ 0.21 0.22 0.23 ] [ 1021 1022 ] = [ 674.06 674.72 ]
926 The matrices are stored in row major order, according to the C
927 convention for arrays.
930 #include <gsl/gsl_blas.h>
935 double a[] = { 0.11, 0.12, 0.13,
938 double b[] = { 1011, 1012,
942 double c[] = { 0.00, 0.00,
945 gsl_matrix_view A = gsl_matrix_view_array(a, 2, 3);
946 gsl_matrix_view B = gsl_matrix_view_array(b, 3, 2);
947 gsl_matrix_view C = gsl_matrix_view_array(c, 2, 2);
949 /* Compute C = A B */
951 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
952 1.0, &A.matrix, &B.matrix,
955 printf ("[ %g, %g\n", c[0], c[1]);
956 printf (" %g, %g ]\n", c[2], c[3]);
961 Here is the output from the program,
968 File: gsl-ref.info, Node: BLAS References and Further Reading, Prev: BLAS Examples, Up: BLAS Support
970 12.3 References and Further Reading
971 ===================================
973 Information on the BLAS standards, including both the legacy and draft
974 interface standards, is available online from the BLAS Homepage and
975 BLAS Technical Forum web-site.
978 `http://www.netlib.org/blas/'
980 `BLAS Technical Forum'
981 `http://www.netlib.org/cgi-bin/checkout/blast/blast.pl'
983 The following papers contain the specifications for Level 1, Level 2 and
986 C. Lawson, R. Hanson, D. Kincaid, F. Krogh, "Basic Linear Algebra
987 Subprograms for Fortran Usage", `ACM Transactions on Mathematical
988 Software', Vol. 5 (1979), Pages 308-325.
990 J.J. Dongarra, J. DuCroz, S. Hammarling, R. Hanson, "An Extended
991 Set of Fortran Basic Linear Algebra Subprograms", `ACM
992 Transactions on Mathematical Software', Vol. 14, No. 1 (1988),
995 J.J. Dongarra, I. Duff, J. DuCroz, S. Hammarling, "A Set of Level
996 3 Basic Linear Algebra Subprograms", `ACM Transactions on
997 Mathematical Software', Vol. 16 (1990), Pages 1-28.
999 Postscript versions of the latter two papers are available from
1000 `http://www.netlib.org/blas/'. A CBLAS wrapper for Fortran BLAS
1001 libraries is available from the same location.
1004 File: gsl-ref.info, Node: Linear Algebra, Next: Eigensystems, Prev: BLAS Support, Up: Top
1009 This chapter describes functions for solving linear systems. The
1010 library provides linear algebra operations which operate directly on
1011 the `gsl_vector' and `gsl_matrix' objects. These routines use the
1012 standard algorithms from Golub & Van Loan's `Matrix Computations' with
1013 Level-1 and Level-2 BLAS calls for efficiency.
1015 The functions described in this chapter are declared in the header
1016 file `gsl_linalg.h'.
1020 * LU Decomposition::
1021 * QR Decomposition::
1022 * QR Decomposition with Column Pivoting::
1023 * Singular Value Decomposition::
1024 * Cholesky Decomposition::
1025 * Tridiagonal Decomposition of Real Symmetric Matrices::
1026 * Tridiagonal Decomposition of Hermitian Matrices::
1027 * Hessenberg Decomposition of Real Matrices::
1028 * Hessenberg-Triangular Decomposition of Real Matrices::
1029 * Bidiagonalization::
1030 * Householder Transformations::
1031 * Householder solver for linear systems::
1032 * Tridiagonal Systems::
1034 * Linear Algebra Examples::
1035 * Linear Algebra References and Further Reading::
1038 File: gsl-ref.info, Node: LU Decomposition, Next: QR Decomposition, Up: Linear Algebra
1040 13.1 LU Decomposition
1041 =====================
1043 A general square matrix A has an LU decomposition into upper and lower
1044 triangular matrices,
1048 where P is a permutation matrix, L is unit lower triangular matrix and
1049 U is upper triangular matrix. For square matrices this decomposition
1050 can be used to convert the linear system A x = b into a pair of
1051 triangular systems (L y = P b, U x = y), which can be solved by forward
1052 and back-substitution. Note that the LU decomposition is valid for
1055 -- Function: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation
1057 -- Function: int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A,
1058 gsl_permutation * P, int * SIGNUM)
1059 These functions factorize the square matrix A into the LU
1060 decomposition PA = LU. On output the diagonal and upper
1061 triangular part of the input matrix A contain the matrix U. The
1062 lower triangular part of the input matrix (excluding the diagonal)
1063 contains L. The diagonal elements of L are unity, and are not
1066 The permutation matrix P is encoded in the permutation P. The j-th
1067 column of the matrix P is given by the k-th column of the identity
1068 matrix, where k = p_j the j-th element of the permutation vector.
1069 The sign of the permutation is given by SIGNUM. It has the value
1070 (-1)^n, where n is the number of interchanges in the permutation.
1072 The algorithm used in the decomposition is Gaussian Elimination
1073 with partial pivoting (Golub & Van Loan, `Matrix Computations',
1076 -- Function: int gsl_linalg_LU_solve (const gsl_matrix * LU, const
1077 gsl_permutation * P, const gsl_vector * B, gsl_vector * X)
1078 -- Function: int gsl_linalg_complex_LU_solve (const gsl_matrix_complex
1079 * LU, const gsl_permutation * P, const gsl_vector_complex *
1080 B, gsl_vector_complex * X)
1081 These functions solve the square system A x = b using the LU
1082 decomposition of A into (LU, P) given by `gsl_linalg_LU_decomp' or
1083 `gsl_linalg_complex_LU_decomp'.
1085 -- Function: int gsl_linalg_LU_svx (const gsl_matrix * LU, const
1086 gsl_permutation * P, gsl_vector * X)
1087 -- Function: int gsl_linalg_complex_LU_svx (const gsl_matrix_complex *
1088 LU, const gsl_permutation * P, gsl_vector_complex * X)
1089 These functions solve the square system A x = b in-place using the
1090 LU decomposition of A into (LU,P). On input X should contain the
1091 right-hand side b, which is replaced by the solution on output.
1093 -- Function: int gsl_linalg_LU_refine (const gsl_matrix * A, const
1094 gsl_matrix * LU, const gsl_permutation * P, const gsl_vector
1095 * B, gsl_vector * X, gsl_vector * RESIDUAL)
1096 -- Function: int gsl_linalg_complex_LU_refine (const
1097 gsl_matrix_complex * A, const gsl_matrix_complex * LU, const
1098 gsl_permutation * P, const gsl_vector_complex * B,
1099 gsl_vector_complex * X, gsl_vector_complex * RESIDUAL)
1100 These functions apply an iterative improvement to X, the solution
1101 of A x = b, using the LU decomposition of A into (LU,P). The
1102 initial residual r = A x - b is also computed and stored in
1105 -- Function: int gsl_linalg_LU_invert (const gsl_matrix * LU, const
1106 gsl_permutation * P, gsl_matrix * INVERSE)
1107 -- Function: int gsl_linalg_complex_LU_invert (const
1108 gsl_matrix_complex * LU, const gsl_permutation * P,
1109 gsl_matrix_complex * INVERSE)
1110 These functions compute the inverse of a matrix A from its LU
1111 decomposition (LU,P), storing the result in the matrix INVERSE.
1112 The inverse is computed by solving the system A x = b for each
1113 column of the identity matrix. It is preferable to avoid direct
1114 use of the inverse whenever possible, as the linear solver
1115 functions can obtain the same result more efficiently and reliably
1116 (consult any introductory textbook on numerical linear algebra for
1119 -- Function: double gsl_linalg_LU_det (gsl_matrix * LU, int SIGNUM)
1120 -- Function: gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex
1122 These functions compute the determinant of a matrix A from its LU
1123 decomposition, LU. The determinant is computed as the product of
1124 the diagonal elements of U and the sign of the row permutation
1127 -- Function: double gsl_linalg_LU_lndet (gsl_matrix * LU)
1128 -- Function: double gsl_linalg_complex_LU_lndet (gsl_matrix_complex *
1130 These functions compute the logarithm of the absolute value of the
1131 determinant of a matrix A, \ln|\det(A)|, from its LU
1132 decomposition, LU. This function may be useful if the direct
1133 computation of the determinant would overflow or underflow.
1135 -- Function: int gsl_linalg_LU_sgndet (gsl_matrix * LU, int SIGNUM)
1136 -- Function: gsl_complex gsl_linalg_complex_LU_sgndet
1137 (gsl_matrix_complex * LU, int SIGNUM)
1138 These functions compute the sign or phase factor of the
1139 determinant of a matrix A, \det(A)/|\det(A)|, from its LU
1143 File: gsl-ref.info, Node: QR Decomposition, Next: QR Decomposition with Column Pivoting, Prev: LU Decomposition, Up: Linear Algebra
1145 13.2 QR Decomposition
1146 =====================
1148 A general rectangular M-by-N matrix A has a QR decomposition into the
1149 product of an orthogonal M-by-M square matrix Q (where Q^T Q = I) and
1150 an M-by-N right-triangular matrix R,
1154 This decomposition can be used to convert the linear system A x = b
1155 into the triangular system R x = Q^T b, which can be solved by
1156 back-substitution. Another use of the QR decomposition is to compute an
1157 orthonormal basis for a set of vectors. The first N columns of Q form
1158 an orthonormal basis for the range of A, ran(A), when A has full column
1161 -- Function: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector *
1163 This function factorizes the M-by-N matrix A into the QR
1164 decomposition A = Q R. On output the diagonal and upper
1165 triangular part of the input matrix contain the matrix R. The
1166 vector TAU and the columns of the lower triangular part of the
1167 matrix A contain the Householder coefficients and Householder
1168 vectors which encode the orthogonal matrix Q. The vector TAU must
1169 be of length k=\min(M,N). The matrix Q is related to these
1170 components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i
1171 v_i^T and v_i is the Householder vector v_i =
1172 (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage
1173 scheme as used by LAPACK.
1175 The algorithm used to perform the decomposition is Householder QR
1176 (Golub & Van Loan, `Matrix Computations', Algorithm 5.2.1).
1178 -- Function: int gsl_linalg_QR_solve (const gsl_matrix * QR, const
1179 gsl_vector * TAU, const gsl_vector * B, gsl_vector * X)
1180 This function solves the square system A x = b using the QR
1181 decomposition of A into (QR, TAU) given by `gsl_linalg_QR_decomp'.
1182 The least-squares solution for rectangular systems can be found
1183 using `gsl_linalg_QR_lssolve'.
1185 -- Function: int gsl_linalg_QR_svx (const gsl_matrix * QR, const
1186 gsl_vector * TAU, gsl_vector * X)
1187 This function solves the square system A x = b in-place using the
1188 QR decomposition of A into (QR,TAU) given by
1189 `gsl_linalg_QR_decomp'. On input X should contain the right-hand
1190 side b, which is replaced by the solution on output.
1192 -- Function: int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const
1193 gsl_vector * TAU, const gsl_vector * B, gsl_vector * X,
1194 gsl_vector * RESIDUAL)
1195 This function finds the least squares solution to the
1196 overdetermined system A x = b where the matrix A has more rows than
1197 columns. The least squares solution minimizes the Euclidean norm
1198 of the residual, ||Ax - b||.The routine uses the QR decomposition
1199 of A into (QR, TAU) given by `gsl_linalg_QR_decomp'. The solution
1200 is returned in X. The residual is computed as a by-product and
1203 -- Function: int gsl_linalg_QR_QTvec (const gsl_matrix * QR, const
1204 gsl_vector * TAU, gsl_vector * V)
1205 This function applies the matrix Q^T encoded in the decomposition
1206 (QR,TAU) to the vector V, storing the result Q^T v in V. The
1207 matrix multiplication is carried out directly using the encoding
1208 of the Householder vectors without needing to form the full matrix
1211 -- Function: int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const
1212 gsl_vector * TAU, gsl_vector * V)
1213 This function applies the matrix Q encoded in the decomposition
1214 (QR,TAU) to the vector V, storing the result Q v in V. The matrix
1215 multiplication is carried out directly using the encoding of the
1216 Householder vectors without needing to form the full matrix Q.
1218 -- Function: int gsl_linalg_QR_QTmat (const gsl_matrix * QR, const
1219 gsl_vector * TAU, gsl_matrix * A)
1220 This function applies the matrix Q^T encoded in the decomposition
1221 (QR,TAU) to the matrix A, storing the result Q^T A in A. The
1222 matrix multiplication is carried out directly using the encoding
1223 of the Householder vectors without needing to form the full matrix
1226 -- Function: int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const
1227 gsl_vector * B, gsl_vector * X)
1228 This function solves the triangular system R x = b for X. It may
1229 be useful if the product b' = Q^T b has already been computed
1230 using `gsl_linalg_QR_QTvec'.
1232 -- Function: int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector
1234 This function solves the triangular system R x = b for X in-place.
1235 On input X should contain the right-hand side b and is replaced by
1236 the solution on output. This function may be useful if the product
1237 b' = Q^T b has already been computed using `gsl_linalg_QR_QTvec'.
1239 -- Function: int gsl_linalg_QR_unpack (const gsl_matrix * QR, const
1240 gsl_vector * TAU, gsl_matrix * Q, gsl_matrix * R)
1241 This function unpacks the encoded QR decomposition (QR,TAU) into
1242 the matrices Q and R, where Q is M-by-M and R is M-by-N.
1244 -- Function: int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix *
1245 R, const gsl_vector * B, gsl_vector * X)
1246 This function solves the system R x = Q^T b for X. It can be used
1247 when the QR decomposition of a matrix is available in unpacked
1250 -- Function: int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R,
1251 gsl_vector * W, const gsl_vector * V)
1252 This function performs a rank-1 update w v^T of the QR
1253 decomposition (Q, R). The update is given by Q'R' = Q R + w v^T
1254 where the output matrices Q' and R' are also orthogonal and right
1255 triangular. Note that W is destroyed by the update.
1257 -- Function: int gsl_linalg_R_solve (const gsl_matrix * R, const
1258 gsl_vector * B, gsl_vector * X)
1259 This function solves the triangular system R x = b for the N-by-N
1262 -- Function: int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector *
1264 This function solves the triangular system R x = b in-place. On
1265 input X should contain the right-hand side b, which is replaced by
1266 the solution on output.
1269 File: gsl-ref.info, Node: QR Decomposition with Column Pivoting, Next: Singular Value Decomposition, Prev: QR Decomposition, Up: Linear Algebra
1271 13.3 QR Decomposition with Column Pivoting
1272 ==========================================
1274 The QR decomposition can be extended to the rank deficient case by
1275 introducing a column permutation P,
1279 The first r columns of Q form an orthonormal basis for the range of A
1280 for a matrix with column rank r. This decomposition can also be used
1281 to convert the linear system A x = b into the triangular system R y =
1282 Q^T b, x = P y, which can be solved by back-substitution and
1283 permutation. We denote the QR decomposition with column pivoting by
1284 QRP^T since A = Q R P^T.
1286 -- Function: int gsl_linalg_QRPT_decomp (gsl_matrix * A, gsl_vector *
1287 TAU, gsl_permutation * P, int * SIGNUM, gsl_vector * NORM)
1288 This function factorizes the M-by-N matrix A into the QRP^T
1289 decomposition A = Q R P^T. On output the diagonal and upper
1290 triangular part of the input matrix contain the matrix R. The
1291 permutation matrix P is stored in the permutation P. The sign of
1292 the permutation is given by SIGNUM. It has the value (-1)^n, where
1293 n is the number of interchanges in the permutation. The vector TAU
1294 and the columns of the lower triangular part of the matrix A
1295 contain the Householder coefficients and vectors which encode the
1296 orthogonal matrix Q. The vector TAU must be of length
1297 k=\min(M,N). The matrix Q is related to these components by, Q =
1298 Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the
1299 Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)).
1300 This is the same storage scheme as used by LAPACK. The vector
1301 NORM is a workspace of length N used for column pivoting.
1303 The algorithm used to perform the decomposition is Householder QR
1304 with column pivoting (Golub & Van Loan, `Matrix Computations',
1307 -- Function: int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A,
1308 gsl_matrix * Q, gsl_matrix * R, gsl_vector * TAU,
1309 gsl_permutation * P, int * SIGNUM, gsl_vector * NORM)
1310 This function factorizes the matrix A into the decomposition A = Q
1311 R P^T without modifying A itself and storing the output in the
1312 separate matrices Q and R.
1314 -- Function: int gsl_linalg_QRPT_solve (const gsl_matrix * QR, const
1315 gsl_vector * TAU, const gsl_permutation * P, const gsl_vector
1316 * B, gsl_vector * X)
1317 This function solves the square system A x = b using the QRP^T
1318 decomposition of A into (QR, TAU, P) given by
1319 `gsl_linalg_QRPT_decomp'.
1321 -- Function: int gsl_linalg_QRPT_svx (const gsl_matrix * QR, const
1322 gsl_vector * TAU, const gsl_permutation * P, gsl_vector * X)
1323 This function solves the square system A x = b in-place using the
1324 QRP^T decomposition of A into (QR,TAU,P). On input X should
1325 contain the right-hand side b, which is replaced by the solution
1328 -- Function: int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, const
1329 gsl_matrix * R, const gsl_permutation * P, const gsl_vector *
1331 This function solves the square system R P^T x = Q^T b for X. It
1332 can be used when the QR decomposition of a matrix is available in
1333 unpacked form as (Q, R).
1335 -- Function: int gsl_linalg_QRPT_update (gsl_matrix * Q, gsl_matrix *
1336 R, const gsl_permutation * P, gsl_vector * U, const
1338 This function performs a rank-1 update w v^T of the QRP^T
1339 decomposition (Q, R, P). The update is given by Q'R' = Q R + w v^T
1340 where the output matrices Q' and R' are also orthogonal and right
1341 triangular. Note that W is destroyed by the update. The
1342 permutation P is not changed.
1344 -- Function: int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR, const
1345 gsl_permutation * P, const gsl_vector * B, gsl_vector * X)
1346 This function solves the triangular system R P^T x = b for the
1347 N-by-N matrix R contained in QR.
1349 -- Function: int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR, const
1350 gsl_permutation * P, gsl_vector * X)
1351 This function solves the triangular system R P^T x = b in-place
1352 for the N-by-N matrix R contained in QR. On input X should contain
1353 the right-hand side b, which is replaced by the solution on output.
1356 File: gsl-ref.info, Node: Singular Value Decomposition, Next: Cholesky Decomposition, Prev: QR Decomposition with Column Pivoting, Up: Linear Algebra
1358 13.4 Singular Value Decomposition
1359 =================================
1361 A general rectangular M-by-N matrix A has a singular value
1362 decomposition (SVD) into the product of an M-by-N orthogonal matrix U,
1363 an N-by-N diagonal matrix of singular values S and the transpose of an
1364 N-by-N orthogonal square matrix V,
1368 The singular values \sigma_i = S_{ii} are all non-negative and are
1369 generally chosen to form a non-increasing sequence \sigma_1 >= \sigma_2
1370 >= ... >= \sigma_N >= 0.
1372 The singular value decomposition of a matrix has many practical uses.
1373 The condition number of the matrix is given by the ratio of the largest
1374 singular value to the smallest singular value. The presence of a zero
1375 singular value indicates that the matrix is singular. The number of
1376 non-zero singular values indicates the rank of the matrix. In practice
1377 singular value decomposition of a rank-deficient matrix will not produce
1378 exact zeroes for singular values, due to finite numerical precision.
1379 Small singular values should be edited by choosing a suitable tolerance.
1381 For a rank-deficient matrix, the null space of A is given by the
1382 columns of V corresponding to the zero singular values. Similarly, the
1383 range of A is given by columns of U corresponding to the non-zero
1386 -- Function: int gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V,
1387 gsl_vector * S, gsl_vector * WORK)
1388 This function factorizes the M-by-N matrix A into the singular
1389 value decomposition A = U S V^T for M >= N. On output the matrix
1390 A is replaced by U. The diagonal elements of the singular value
1391 matrix S are stored in the vector S. The singular values are
1392 non-negative and form a non-increasing sequence from S_1 to S_N.
1393 The matrix V contains the elements of V in untransposed form. To
1394 form the product U S V^T it is necessary to take the transpose of
1395 V. A workspace of length N is required in WORK.
1397 This routine uses the Golub-Reinsch SVD algorithm.
1399 -- Function: int gsl_linalg_SV_decomp_mod (gsl_matrix * A, gsl_matrix
1400 * X, gsl_matrix * V, gsl_vector * S, gsl_vector * WORK)
1401 This function computes the SVD using the modified Golub-Reinsch
1402 algorithm, which is faster for M>>N. It requires the vector WORK
1403 of length N and the N-by-N matrix X as additional working space.
1405 -- Function: int gsl_linalg_SV_decomp_jacobi (gsl_matrix * A,
1406 gsl_matrix * V, gsl_vector * S)
1407 This function computes the SVD of the M-by-N matrix A using
1408 one-sided Jacobi orthogonalization for M >= N. The Jacobi method
1409 can compute singular values to higher relative accuracy than
1410 Golub-Reinsch algorithms (see references for details).
1412 -- Function: int gsl_linalg_SV_solve (gsl_matrix * U, gsl_matrix * V,
1413 gsl_vector * S, const gsl_vector * B, gsl_vector * X)
1414 This function solves the system A x = b using the singular value
1415 decomposition (U, S, V) of A given by `gsl_linalg_SV_decomp'.
1417 Only non-zero singular values are used in computing the solution.
1418 The parts of the solution corresponding to singular values of zero
1419 are ignored. Other singular values can be edited out by setting
1420 them to zero before calling this function.
1422 In the over-determined case where A has more rows than columns the
1423 system is solved in the least squares sense, returning the solution
1424 X which minimizes ||A x - b||_2.
1427 File: gsl-ref.info, Node: Cholesky Decomposition, Next: Tridiagonal Decomposition of Real Symmetric Matrices, Prev: Singular Value Decomposition, Up: Linear Algebra
1429 13.5 Cholesky Decomposition
1430 ===========================
1432 A symmetric, positive definite square matrix A has a Cholesky
1433 decomposition into a product of a lower triangular matrix L and its
1438 This is sometimes referred to as taking the square-root of a matrix. The
1439 Cholesky decomposition can only be carried out when all the eigenvalues
1440 of the matrix are positive. This decomposition can be used to convert
1441 the linear system A x = b into a pair of triangular systems (L y = b,
1442 L^T x = y), which can be solved by forward and back-substitution.
1444 -- Function: int gsl_linalg_cholesky_decomp (gsl_matrix * A)
1445 -- Function: int gsl_linalg_complex_cholesky_decomp
1446 (gsl_matrix_complex * A)
1447 These functions factorize the symmetric, positive-definite square
1448 matrix A into the Cholesky decomposition A = L L^T (or A = L L^H
1449 for the complex case). On input, the values from the diagonal and
1450 lower-triangular part of the matrix A are used (the upper
1451 triangular part is ignored). On output the diagonal and lower
1452 triangular part of the input matrix A contain the matrix L, while
1453 the upper triangular part of the input matrix is overwritten with
1454 L^T (the diagonal terms being identical for both L and L^T). If
1455 the matrix is not positive-definite then the decomposition will
1456 fail, returning the error code `GSL_EDOM'.
1458 When testing whether a matrix is positive-definite, disable the
1459 error handler first to avoid triggering an error.
1461 -- Function: int gsl_linalg_cholesky_solve (const gsl_matrix *
1462 CHOLESKY, const gsl_vector * B, gsl_vector * X)
1463 -- Function: int gsl_linalg_complex_cholesky_solve (const
1464 gsl_matrix_complex * CHOLESKY, const gsl_vector_complex * B,
1465 gsl_vector_complex * X)
1466 These functions solve the system A x = b using the Cholesky
1467 decomposition of A into the matrix CHOLESKY given by
1468 `gsl_linalg_cholesky_decomp' or
1469 `gsl_linalg_complex_cholesky_decomp'.
1471 -- Function: int gsl_linalg_cholesky_svx (const gsl_matrix * CHOLESKY,
1473 -- Function: int gsl_linalg_complex_cholesky_svx (const
1474 gsl_matrix_complex * CHOLESKY, gsl_vector_complex * X)
1475 These functions solve the system A x = b in-place using the
1476 Cholesky decomposition of A into the matrix CHOLESKY given by
1477 `gsl_linalg_cholesky_decomp' or
1478 `gsl_linalg_complex_cholesky_decomp'. On input X should contain
1479 the right-hand side b, which is replaced by the solution on output.
1482 File: gsl-ref.info, Node: Tridiagonal Decomposition of Real Symmetric Matrices, Next: Tridiagonal Decomposition of Hermitian Matrices, Prev: Cholesky Decomposition, Up: Linear Algebra
1484 13.6 Tridiagonal Decomposition of Real Symmetric Matrices
1485 =========================================================
1487 A symmetric matrix A can be factorized by similarity transformations
1492 where Q is an orthogonal matrix and T is a symmetric tridiagonal matrix.
1494 -- Function: int gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector
1496 This function factorizes the symmetric square matrix A into the
1497 symmetric tridiagonal decomposition Q T Q^T. On output the
1498 diagonal and subdiagonal part of the input matrix A contain the
1499 tridiagonal matrix T. The remaining lower triangular part of the
1500 input matrix contains the Householder vectors which, together with
1501 the Householder coefficients TAU, encode the orthogonal matrix Q.
1502 This storage scheme is the same as used by LAPACK. The upper
1503 triangular part of A is not referenced.
1505 -- Function: int gsl_linalg_symmtd_unpack (const gsl_matrix * A, const
1506 gsl_vector * TAU, gsl_matrix * Q, gsl_vector * DIAG,
1507 gsl_vector * SUBDIAG)
1508 This function unpacks the encoded symmetric tridiagonal
1509 decomposition (A, TAU) obtained from `gsl_linalg_symmtd_decomp'
1510 into the orthogonal matrix Q, the vector of diagonal elements DIAG
1511 and the vector of subdiagonal elements SUBDIAG.
1513 -- Function: int gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
1514 gsl_vector * DIAG, gsl_vector * SUBDIAG)
1515 This function unpacks the diagonal and subdiagonal of the encoded
1516 symmetric tridiagonal decomposition (A, TAU) obtained from
1517 `gsl_linalg_symmtd_decomp' into the vectors DIAG and SUBDIAG.
1520 File: gsl-ref.info, Node: Tridiagonal Decomposition of Hermitian Matrices, Next: Hessenberg Decomposition of Real Matrices, Prev: Tridiagonal Decomposition of Real Symmetric Matrices, Up: Linear Algebra
1522 13.7 Tridiagonal Decomposition of Hermitian Matrices
1523 ====================================================
1525 A hermitian matrix A can be factorized by similarity transformations
1530 where U is a unitary matrix and T is a real symmetric tridiagonal
1533 -- Function: int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A,
1534 gsl_vector_complex * TAU)
1535 This function factorizes the hermitian matrix A into the symmetric
1536 tridiagonal decomposition U T U^T. On output the real parts of
1537 the diagonal and subdiagonal part of the input matrix A contain
1538 the tridiagonal matrix T. The remaining lower triangular part of
1539 the input matrix contains the Householder vectors which, together
1540 with the Householder coefficients TAU, encode the orthogonal matrix
1541 Q. This storage scheme is the same as used by LAPACK. The upper
1542 triangular part of A and imaginary parts of the diagonal are not
1545 -- Function: int gsl_linalg_hermtd_unpack (const gsl_matrix_complex *
1546 A, const gsl_vector_complex * TAU, gsl_matrix_complex * Q,
1547 gsl_vector * DIAG, gsl_vector * SUBDIAG)
1548 This function unpacks the encoded tridiagonal decomposition (A,
1549 TAU) obtained from `gsl_linalg_hermtd_decomp' into the unitary
1550 matrix U, the real vector of diagonal elements DIAG and the real
1551 vector of subdiagonal elements SUBDIAG.
1553 -- Function: int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex
1554 * A, gsl_vector * DIAG, gsl_vector * SUBDIAG)
1555 This function unpacks the diagonal and subdiagonal of the encoded
1556 tridiagonal decomposition (A, TAU) obtained from the
1557 `gsl_linalg_hermtd_decomp' into the real vectors DIAG and SUBDIAG.
1560 File: gsl-ref.info, Node: Hessenberg Decomposition of Real Matrices, Next: Hessenberg-Triangular Decomposition of Real Matrices, Prev: Tridiagonal Decomposition of Hermitian Matrices, Up: Linear Algebra
1562 13.8 Hessenberg Decomposition of Real Matrices
1563 ==============================================
1565 A general real matrix A can be decomposed by orthogonal similarity
1566 transformations into the form
1570 where U is orthogonal and H is an upper Hessenberg matrix, meaning
1571 that it has zeros below the first subdiagonal. The Hessenberg reduction
1572 is the first step in the Schur decomposition for the nonsymmetric
1573 eigenvalue problem, but has applications in other areas as well.
1575 -- Function: int gsl_linalg_hessenberg_decomp (gsl_matrix * A,
1577 This function computes the Hessenberg decomposition of the matrix
1578 A by applying the similarity transformation H = U^T A U. On
1579 output, H is stored in the upper portion of A. The information
1580 required to construct the matrix U is stored in the lower
1581 triangular portion of A. U is a product of N - 2 Householder
1582 matrices. The Householder vectors are stored in the lower portion
1583 of A (below the subdiagonal) and the Householder coefficients are
1584 stored in the vector TAU. TAU must be of length N.
1586 -- Function: int gsl_linalg_hessenberg_unpack (gsl_matrix * H,
1587 gsl_vector * TAU, gsl_matrix * U)
1588 This function constructs the orthogonal matrix U from the
1589 information stored in the Hessenberg matrix H along with the
1590 vector TAU. H and TAU are outputs from
1591 `gsl_linalg_hessenberg_decomp'.
1593 -- Function: int gsl_linalg_hessenberg_unpack_accum (gsl_matrix * H,
1594 gsl_vector * TAU, gsl_matrix * V)
1595 This function is similar to `gsl_linalg_hessenberg_unpack', except
1596 it accumulates the matrix U into V, so that V' = VU. The matrix V
1597 must be initialized prior to calling this function. Setting V to
1598 the identity matrix provides the same result as
1599 `gsl_linalg_hessenberg_unpack'. If H is order N, then V must have
1600 N columns but may have any number of rows.
1602 -- Function: int gsl_linalg_hessenberg_set_zero (gsl_matrix * H)
1603 This function sets the lower triangular portion of H, below the
1604 subdiagonal, to zero. It is useful for clearing out the
1605 Householder vectors after calling `gsl_linalg_hessenberg_decomp'.
1608 File: gsl-ref.info, Node: Hessenberg-Triangular Decomposition of Real Matrices, Next: Bidiagonalization, Prev: Hessenberg Decomposition of Real Matrices, Up: Linear Algebra
1610 13.9 Hessenberg-Triangular Decomposition of Real Matrices
1611 =========================================================
1613 A general real matrix pair (A, B) can be decomposed by orthogonal
1614 similarity transformations into the form
1619 where U and V are orthogonal, H is an upper Hessenberg matrix, and R
1620 is upper triangular. The Hessenberg-Triangular reduction is the first
1621 step in the generalized Schur decomposition for the generalized
1624 -- Function: int gsl_linalg_hesstri_decomp (gsl_matrix * A, gsl_matrix
1625 * B, gsl_matrix * U, gsl_matrix * V, gsl_vector * WORK)
1626 This function computes the Hessenberg-Triangular decomposition of
1627 the matrix pair (A, B). On output, H is stored in A, and R is
1628 stored in B. If U and V are provided (they may be null), the
1629 similarity transformations are stored in them. Additional
1630 workspace of length N is needed in WORK.
1633 File: gsl-ref.info, Node: Bidiagonalization, Next: Householder Transformations, Prev: Hessenberg-Triangular Decomposition of Real Matrices, Up: Linear Algebra
1635 13.10 Bidiagonalization
1636 =======================
1638 A general matrix A can be factorized by similarity transformations into
1643 where U and V are orthogonal matrices and B is a N-by-N bidiagonal
1644 matrix with non-zero entries only on the diagonal and superdiagonal.
1645 The size of U is M-by-N and the size of V is N-by-N.
1647 -- Function: int gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector
1648 * TAU_U, gsl_vector * TAU_V)
1649 This function factorizes the M-by-N matrix A into bidiagonal form
1650 U B V^T. The diagonal and superdiagonal of the matrix B are
1651 stored in the diagonal and superdiagonal of A. The orthogonal
1652 matrices U and V are stored as compressed Householder vectors in
1653 the remaining elements of A. The Householder coefficients are
1654 stored in the vectors TAU_U and TAU_V. The length of TAU_U must
1655 equal the number of elements in the diagonal of A and the length
1656 of TAU_V should be one element shorter.
1658 -- Function: int gsl_linalg_bidiag_unpack (const gsl_matrix * A, const
1659 gsl_vector * TAU_U, gsl_matrix * U, const gsl_vector * TAU_V,
1660 gsl_matrix * V, gsl_vector * DIAG, gsl_vector * SUPERDIAG)
1661 This function unpacks the bidiagonal decomposition of A given by
1662 `gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate
1663 orthogonal matrices U, V and the diagonal vector DIAG and
1664 superdiagonal SUPERDIAG. Note that U is stored as a compact
1665 M-by-N orthogonal matrix satisfying U^T U = I for efficiency.
1667 -- Function: int gsl_linalg_bidiag_unpack2 (gsl_matrix * A, gsl_vector
1668 * TAU_U, gsl_vector * TAU_V, gsl_matrix * V)
1669 This function unpacks the bidiagonal decomposition of A given by
1670 `gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate
1671 orthogonal matrices U, V and the diagonal vector DIAG and
1672 superdiagonal SUPERDIAG. The matrix U is stored in-place in A.
1674 -- Function: int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A,
1675 gsl_vector * DIAG, gsl_vector * SUPERDIAG)
1676 This function unpacks the diagonal and superdiagonal of the
1677 bidiagonal decomposition of A given by `gsl_linalg_bidiag_decomp',
1678 into the diagonal vector DIAG and superdiagonal vector SUPERDIAG.
1681 File: gsl-ref.info, Node: Householder Transformations, Next: Householder solver for linear systems, Prev: Bidiagonalization, Up: Linear Algebra
1683 13.11 Householder Transformations
1684 =================================
1686 A Householder transformation is a rank-1 modification of the identity
1687 matrix which can be used to zero out selected elements of a vector. A
1688 Householder matrix P takes the form,
1692 where v is a vector (called the "Householder vector") and \tau = 2/(v^T
1693 v). The functions described in this section use the rank-1 structure
1694 of the Householder matrix to create and apply Householder
1695 transformations efficiently.
1697 -- Function: double gsl_linalg_householder_transform (gsl_vector * V)
1698 -- Function: gsl_complex gsl_linalg_complex_householder_transform
1699 (gsl_vector_complex * V)
1700 This function prepares a Householder transformation P = I - \tau v
1701 v^T which can be used to zero all the elements of the input vector
1702 except the first. On output the transformation is stored in the
1703 vector V and the scalar \tau is returned.
1705 -- Function: int gsl_linalg_householder_hm (double tau, const
1706 gsl_vector * v, gsl_matrix * A)
1707 -- Function: int gsl_linalg_complex_householder_hm (gsl_complex tau,
1708 const gsl_vector_complex * v, gsl_matrix_complex * A)
1709 This function applies the Householder matrix P defined by the
1710 scalar TAU and the vector V to the left-hand side of the matrix A.
1711 On output the result P A is stored in A.
1713 -- Function: int gsl_linalg_householder_mh (double tau, const
1714 gsl_vector * v, gsl_matrix * A)
1715 -- Function: int gsl_linalg_complex_householder_mh (gsl_complex tau,
1716 const gsl_vector_complex * v, gsl_matrix_complex * A)
1717 This function applies the Householder matrix P defined by the
1718 scalar TAU and the vector V to the right-hand side of the matrix
1719 A. On output the result A P is stored in A.
1721 -- Function: int gsl_linalg_householder_hv (double tau, const
1722 gsl_vector * v, gsl_vector * w)
1723 -- Function: int gsl_linalg_complex_householder_hv (gsl_complex tau,
1724 const gsl_vector_complex * v, gsl_vector_complex * w)
1725 This function applies the Householder transformation P defined by
1726 the scalar TAU and the vector V to the vector W. On output the
1727 result P w is stored in W.
1730 File: gsl-ref.info, Node: Householder solver for linear systems, Next: Tridiagonal Systems, Prev: Householder Transformations, Up: Linear Algebra
1732 13.12 Householder solver for linear systems
1733 ===========================================
1735 -- Function: int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector
1736 * B, gsl_vector * X)
1737 This function solves the system A x = b directly using Householder
1738 transformations. On output the solution is stored in X and B is
1739 not modified. The matrix A is destroyed by the Householder
1742 -- Function: int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * X)
1743 This function solves the system A x = b in-place using Householder
1744 transformations. On input X should contain the right-hand side b,
1745 which is replaced by the solution on output. The matrix A is
1746 destroyed by the Householder transformations.
1749 File: gsl-ref.info, Node: Tridiagonal Systems, Next: Balancing, Prev: Householder solver for linear systems, Up: Linear Algebra
1751 13.13 Tridiagonal Systems
1752 =========================
1754 The functions described in this section efficiently solve symmetric,
1755 non-symmetric and cyclic tridiagonal systems with minimal storage.
1756 Note that the current implementations of these functions use a variant
1757 of Cholesky decomposition, so the tridiagonal matrix must be positive
1758 definite. For non-positive definite matrices, the functions return the
1759 error code `GSL_ESING'.
1761 -- Function: int gsl_linalg_solve_tridiag (const gsl_vector * DIAG,
1762 const gsl_vector * E, const gsl_vector * F, const gsl_vector
1763 * B, gsl_vector * X)
1764 This function solves the general N-by-N system A x = b where A is
1765 tridiagonal (N >= 2). The super-diagonal and sub-diagonal vectors
1766 E and F must be one element shorter than the diagonal vector DIAG.
1767 The form of A for the 4-by-4 case is shown below,
1775 -- Function: int gsl_linalg_solve_symm_tridiag (const gsl_vector *
1776 DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector
1778 This function solves the general N-by-N system A x = b where A is
1779 symmetric tridiagonal (N >= 2). The off-diagonal vector E must be
1780 one element shorter than the diagonal vector DIAG. The form of A
1781 for the 4-by-4 case is shown below,
1788 -- Function: int gsl_linalg_solve_cyc_tridiag (const gsl_vector *
1789 DIAG, const gsl_vector * E, const gsl_vector * F, const
1790 gsl_vector * B, gsl_vector * X)
1791 This function solves the general N-by-N system A x = b where A is
1792 cyclic tridiagonal (N >= 3). The cyclic super-diagonal and
1793 sub-diagonal vectors E and F must have the same number of elements
1794 as the diagonal vector DIAG. The form of A for the 4-by-4 case is
1797 A = ( d_0 e_0 0 f_3 )
1802 -- Function: int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector *
1803 DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector
1805 This function solves the general N-by-N system A x = b where A is
1806 symmetric cyclic tridiagonal (N >= 3). The cyclic off-diagonal
1807 vector E must have the same number of elements as the diagonal
1808 vector DIAG. The form of A for the 4-by-4 case is shown below,
1810 A = ( d_0 e_0 0 e_3 )
1816 File: gsl-ref.info, Node: Balancing, Next: Linear Algebra Examples, Prev: Tridiagonal Systems, Up: Linear Algebra
1821 The process of balancing a matrix applies similarity transformations to
1822 make the rows and columns have comparable norms. This is useful, for
1823 example, to reduce roundoff errors in the solution of eigenvalue
1824 problems. Balancing a matrix A consists of replacing A with a similar
1829 where D is a diagonal matrix whose entries are powers of the
1830 floating point radix.
1832 -- Function: int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector
1834 This function replaces the matrix A with its balanced counterpart
1835 and stores the diagonal elements of the similarity transformation
1839 File: gsl-ref.info, Node: Linear Algebra Examples, Next: Linear Algebra References and Further Reading, Prev: Balancing, Up: Linear Algebra
1844 The following program solves the linear system A x = b. The system to
1847 [ 0.18 0.60 0.57 0.96 ] [x0] [1.0]
1848 [ 0.41 0.24 0.99 0.58 ] [x1] = [2.0]
1849 [ 0.14 0.30 0.97 0.66 ] [x2] [3.0]
1850 [ 0.51 0.13 0.19 0.85 ] [x3] [4.0]
1852 and the solution is found using LU decomposition of the matrix A.
1855 #include <gsl/gsl_linalg.h>
1860 double a_data[] = { 0.18, 0.60, 0.57, 0.96,
1861 0.41, 0.24, 0.99, 0.58,
1862 0.14, 0.30, 0.97, 0.66,
1863 0.51, 0.13, 0.19, 0.85 };
1865 double b_data[] = { 1.0, 2.0, 3.0, 4.0 };
1868 = gsl_matrix_view_array (a_data, 4, 4);
1871 = gsl_vector_view_array (b_data, 4);
1873 gsl_vector *x = gsl_vector_alloc (4);
1877 gsl_permutation * p = gsl_permutation_alloc (4);
1879 gsl_linalg_LU_decomp (&m.matrix, p, &s);
1881 gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
1884 gsl_vector_fprintf (stdout, x, "%g");
1886 gsl_permutation_free (p);
1887 gsl_vector_free (x);
1891 Here is the output from the program,
1898 This can be verified by multiplying the solution x by the original
1899 matrix A using GNU OCTAVE,
1901 octave> A = [ 0.18, 0.60, 0.57, 0.96;
1902 0.41, 0.24, 0.99, 0.58;
1903 0.14, 0.30, 0.97, 0.66;
1904 0.51, 0.13, 0.19, 0.85 ];
1906 octave> x = [ -4.05205; -12.6056; 1.66091; 8.69377];
1915 This reproduces the original right-hand side vector, b, in accordance
1916 with the equation A x = b.
1919 File: gsl-ref.info, Node: Linear Algebra References and Further Reading, Prev: Linear Algebra Examples, Up: Linear Algebra
1921 13.16 References and Further Reading
1922 ====================================
1924 Further information on the algorithms described in this section can be
1925 found in the following book,
1927 G. H. Golub, C. F. Van Loan, `Matrix Computations' (3rd Ed, 1996),
1928 Johns Hopkins University Press, ISBN 0-8018-5414-8.
1930 The LAPACK library is described in the following manual,
1932 `LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM,
1935 `http://www.netlib.org/lapack'
1937 The LAPACK source code can be found at the website above, along with an
1938 online copy of the users guide.
1940 The Modified Golub-Reinsch algorithm is described in the following
1943 T.F. Chan, "An Improved Algorithm for Computing the Singular Value
1944 Decomposition", `ACM Transactions on Mathematical Software', 8
1947 The Jacobi algorithm for singular value decomposition is described in
1948 the following papers,
1950 J.C. Nash, "A one-sided transformation method for the singular
1951 value decomposition and algebraic eigenproblem", `Computer
1952 Journal', Volume 18, Number 1 (1973), p 74-76
1954 James Demmel, Kresimir Veselic, "Jacobi's Method is more accurate
1955 than QR", `Lapack Working Note 15' (LAWN-15), October 1989.
1956 Available from netlib, `http://www.netlib.org/lapack/' in the
1957 `lawns' or `lawnspdf' directories.
1960 File: gsl-ref.info, Node: Eigensystems, Next: Fast Fourier Transforms, Prev: Linear Algebra, Up: Top
1965 This chapter describes functions for computing eigenvalues and
1966 eigenvectors of matrices. There are routines for real symmetric, real
1967 nonsymmetric, complex hermitian, real generalized symmetric-definite,
1968 complex generalized hermitian-definite, and real generalized
1969 nonsymmetric eigensystems. Eigenvalues can be computed with or without
1970 eigenvectors. The hermitian and real symmetric matrix algorithms are
1971 symmetric bidiagonalization followed by QR reduction. The nonsymmetric
1972 algorithm is the Francis QR double-shift. The generalized nonsymmetric
1973 algorithm is the QZ method due to Moler and Stewart.
1975 The functions described in this chapter are declared in the header
1980 * Real Symmetric Matrices::
1981 * Complex Hermitian Matrices::
1982 * Real Nonsymmetric Matrices::
1983 * Real Generalized Symmetric-Definite Eigensystems::
1984 * Complex Generalized Hermitian-Definite Eigensystems::
1985 * Real Generalized Nonsymmetric Eigensystems::
1986 * Sorting Eigenvalues and Eigenvectors::
1987 * Eigenvalue and Eigenvector Examples::
1988 * Eigenvalue and Eigenvector References::
1991 File: gsl-ref.info, Node: Real Symmetric Matrices, Next: Complex Hermitian Matrices, Up: Eigensystems
1993 14.1 Real Symmetric Matrices
1994 ============================
1996 For real symmetric matrices, the library uses the symmetric
1997 bidiagonalization and QR reduction method. This is described in Golub
1998 & van Loan, section 8.3. The computed eigenvalues are accurate to an
1999 absolute accuracy of \epsilon ||A||_2, where \epsilon is the machine
2002 -- Function: gsl_eigen_symm_workspace * gsl_eigen_symm_alloc (const
2004 This function allocates a workspace for computing eigenvalues of
2005 N-by-N real symmetric matrices. The size of the workspace is
2008 -- Function: void gsl_eigen_symm_free (gsl_eigen_symm_workspace * W)
2009 This function frees the memory associated with the workspace W.
2011 -- Function: int gsl_eigen_symm (gsl_matrix * A, gsl_vector * EVAL,
2012 gsl_eigen_symm_workspace * W)
2013 This function computes the eigenvalues of the real symmetric matrix
2014 A. Additional workspace of the appropriate size must be provided
2015 in W. The diagonal and lower triangular part of A are destroyed
2016 during the computation, but the strict upper triangular part is
2017 not referenced. The eigenvalues are stored in the vector EVAL and
2020 -- Function: gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const
2022 This function allocates a workspace for computing eigenvalues and
2023 eigenvectors of N-by-N real symmetric matrices. The size of the
2026 -- Function: void gsl_eigen_symmv_free (gsl_eigen_symmv_workspace * W)
2027 This function frees the memory associated with the workspace W.
2029 -- Function: int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * EVAL,
2030 gsl_matrix * EVEC, gsl_eigen_symmv_workspace * W)
2031 This function computes the eigenvalues and eigenvectors of the real
2032 symmetric matrix A. Additional workspace of the appropriate size
2033 must be provided in W. The diagonal and lower triangular part of
2034 A are destroyed during the computation, but the strict upper
2035 triangular part is not referenced. The eigenvalues are stored in
2036 the vector EVAL and are unordered. The corresponding eigenvectors
2037 are stored in the columns of the matrix EVEC. For example, the
2038 eigenvector in the first column corresponds to the first
2039 eigenvalue. The eigenvectors are guaranteed to be mutually
2040 orthogonal and normalised to unit magnitude.
2043 File: gsl-ref.info, Node: Complex Hermitian Matrices, Next: Real Nonsymmetric Matrices, Prev: Real Symmetric Matrices, Up: Eigensystems
2045 14.2 Complex Hermitian Matrices
2046 ===============================
2048 -- Function: gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const
2050 This function allocates a workspace for computing eigenvalues of
2051 N-by-N complex hermitian matrices. The size of the workspace is
2054 -- Function: void gsl_eigen_herm_free (gsl_eigen_herm_workspace * W)
2055 This function frees the memory associated with the workspace W.
2057 -- Function: int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector *
2058 EVAL, gsl_eigen_herm_workspace * W)
2059 This function computes the eigenvalues of the complex hermitian
2060 matrix A. Additional workspace of the appropriate size must be
2061 provided in W. The diagonal and lower triangular part of A are
2062 destroyed during the computation, but the strict upper triangular
2063 part is not referenced. The imaginary parts of the diagonal are
2064 assumed to be zero and are not referenced. The eigenvalues are
2065 stored in the vector EVAL and are unordered.
2067 -- Function: gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const
2069 This function allocates a workspace for computing eigenvalues and
2070 eigenvectors of N-by-N complex hermitian matrices. The size of
2071 the workspace is O(5n).
2073 -- Function: void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * W)
2074 This function frees the memory associated with the workspace W.
2076 -- Function: int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector *
2077 EVAL, gsl_matrix_complex * EVEC, gsl_eigen_hermv_workspace *
2079 This function computes the eigenvalues and eigenvectors of the
2080 complex hermitian matrix A. Additional workspace of the
2081 appropriate size must be provided in W. The diagonal and lower
2082 triangular part of A are destroyed during the computation, but the
2083 strict upper triangular part is not referenced. The imaginary
2084 parts of the diagonal are assumed to be zero and are not
2085 referenced. The eigenvalues are stored in the vector EVAL and are
2086 unordered. The corresponding complex eigenvectors are stored in
2087 the columns of the matrix EVEC. For example, the eigenvector in
2088 the first column corresponds to the first eigenvalue. The
2089 eigenvectors are guaranteed to be mutually orthogonal and
2090 normalised to unit magnitude.
2093 File: gsl-ref.info, Node: Real Nonsymmetric Matrices, Next: Real Generalized Symmetric-Definite Eigensystems, Prev: Complex Hermitian Matrices, Up: Eigensystems
2095 14.3 Real Nonsymmetric Matrices
2096 ===============================
2098 The solution of the real nonsymmetric eigensystem problem for a matrix
2099 A involves computing the Schur decomposition
2103 where Z is an orthogonal matrix of Schur vectors and T, the Schur
2104 form, is quasi upper triangular with diagonal 1-by-1 blocks which are
2105 real eigenvalues of A, and diagonal 2-by-2 blocks whose eigenvalues are
2106 complex conjugate eigenvalues of A. The algorithm used is the double
2107 shift Francis method.
2109 -- Function: gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc
2111 This function allocates a workspace for computing eigenvalues of
2112 N-by-N real nonsymmetric matrices. The size of the workspace is
2115 -- Function: void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace
2117 This function frees the memory associated with the workspace W.
2119 -- Function: void gsl_eigen_nonsymm_params (const int COMPUTE_T, const
2120 int BALANCE, gsl_eigen_nonsymm_workspace * W)
2121 This function sets some parameters which determine how the
2122 eigenvalue problem is solved in subsequent calls to
2123 `gsl_eigen_nonsymm'.
2125 If COMPUTE_T is set to 1, the full Schur form T will be computed
2126 by `gsl_eigen_nonsymm'. If it is set to 0, T will not be computed
2127 (this is the default setting). Computing the full Schur form T
2128 requires approximately 1.5-2 times the number of flops.
2130 If BALANCE is set to 1, a balancing transformation is applied to
2131 the matrix prior to computing eigenvalues. This transformation is
2132 designed to make the rows and columns of the matrix have comparable
2133 norms, and can result in more accurate eigenvalues for matrices
2134 whose entries vary widely in magnitude. See *Note Balancing:: for
2135 more information. Note that the balancing transformation does not
2136 preserve the orthogonality of the Schur vectors, so if you wish to
2137 compute the Schur vectors with `gsl_eigen_nonsymm_Z' you will
2138 obtain the Schur vectors of the balanced matrix instead of the
2139 original matrix. The relationship will be
2141 T = Q^t D^(-1) A D Q
2143 where Q is the matrix of Schur vectors for the balanced matrix, and
2144 D is the balancing transformation. Then `gsl_eigen_nonsymm_Z' will
2145 compute a matrix Z which satisfies
2149 with Z = D Q. Note that Z will not be orthogonal. For this reason,
2150 balancing is not performed by default.
2152 -- Function: int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex
2153 * EVAL, gsl_eigen_nonsymm_workspace * W)
2154 This function computes the eigenvalues of the real nonsymmetric
2155 matrix A and stores them in the vector EVAL. If T is desired, it
2156 is stored in the upper portion of A on output. Otherwise, on
2157 output, the diagonal of A will contain the 1-by-1 real eigenvalues
2158 and 2-by-2 complex conjugate eigenvalue systems, and the rest of A
2159 is destroyed. In rare cases, this function may fail to find all
2160 eigenvalues. If this happens, an error code is returned and the
2161 number of converged eigenvalues is stored in `w->n_evals'. The
2162 converged eigenvalues are stored in the beginning of EVAL.
2164 -- Function: int gsl_eigen_nonsymm_Z (gsl_matrix * A,
2165 gsl_vector_complex * EVAL, gsl_matrix * Z,
2166 gsl_eigen_nonsymm_workspace * W)
2167 This function is identical to `gsl_eigen_nonsymm' except it also
2168 computes the Schur vectors and stores them into Z.
2170 -- Function: gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc
2172 This function allocates a workspace for computing eigenvalues and
2173 eigenvectors of N-by-N real nonsymmetric matrices. The size of the
2176 -- Function: void gsl_eigen_nonsymmv_free
2177 (gsl_eigen_nonsymmv_workspace * W)
2178 This function frees the memory associated with the workspace W.
2180 -- Function: int gsl_eigen_nonsymmv (gsl_matrix * A,
2181 gsl_vector_complex * EVAL, gsl_matrix_complex * EVEC,
2182 gsl_eigen_nonsymmv_workspace * W)
2183 This function computes eigenvalues and right eigenvectors of the
2184 N-by-N real nonsymmetric matrix A. It first calls
2185 `gsl_eigen_nonsymm' to compute the eigenvalues, Schur form T, and
2186 Schur vectors. Then it finds eigenvectors of T and backtransforms
2187 them using the Schur vectors. The Schur vectors are destroyed in
2188 the process, but can be saved by using `gsl_eigen_nonsymmv_Z'. The
2189 computed eigenvectors are normalized to have unit magnitude. On
2190 output, the upper portion of A contains the Schur form T. If
2191 `gsl_eigen_nonsymm' fails, no eigenvectors are computed, and an
2192 error code is returned.
2194 -- Function: int gsl_eigen_nonsymmv_Z (gsl_matrix * A,
2195 gsl_vector_complex * EVAL, gsl_matrix_complex * EVEC,
2196 gsl_matrix * Z, gsl_eigen_nonsymmv_workspace * W)
2197 This function is identical to `gsl_eigen_nonsymmv' except it also
2198 saves the Schur vectors into Z.
2201 File: gsl-ref.info, Node: Real Generalized Symmetric-Definite Eigensystems, Next: Complex Generalized Hermitian-Definite Eigensystems, Prev: Real Nonsymmetric Matrices, Up: Eigensystems
2203 14.4 Real Generalized Symmetric-Definite Eigensystems
2204 =====================================================
2206 The real generalized symmetric-definite eigenvalue problem is to find
2207 eigenvalues \lambda and eigenvectors x such that
2209 where A and B are symmetric matrices, and B is positive-definite.
2210 This problem reduces to the standard symmetric eigenvalue problem by
2211 applying the Cholesky decomposition to B:
2213 A x = \lambda L L^t x
2214 ( L^{-1} A L^{-t} ) L^t x = \lambda L^t x
2215 Therefore, the problem becomes C y = \lambda y where C = L^{-1} A
2216 L^{-t} is symmetric, and y = L^t x. The standard symmetric eigensolver
2217 can be applied to the matrix C. The resulting eigenvectors are
2218 backtransformed to find the vectors of the original problem. The
2219 eigenvalues and eigenvectors of the generalized symmetric-definite
2220 eigenproblem are always real.
2222 -- Function: gsl_eigen_gensymm_workspace * gsl_eigen_gensymm_alloc
2224 This function allocates a workspace for computing eigenvalues of
2225 N-by-N real generalized symmetric-definite eigensystems. The size
2226 of the workspace is O(2n).
2228 -- Function: void gsl_eigen_gensymm_free (gsl_eigen_gensymm_workspace
2230 This function frees the memory associated with the workspace W.
2232 -- Function: int gsl_eigen_gensymm (gsl_matrix * A, gsl_matrix * B,
2233 gsl_vector * EVAL, gsl_eigen_gensymm_workspace * W)
2234 This function computes the eigenvalues of the real generalized
2235 symmetric-definite matrix pair (A, B), and stores them in EVAL,
2236 using the method outlined above. On output, B contains its
2237 Cholesky decomposition and A is destroyed.
2239 -- Function: gsl_eigen_gensymmv_workspace * gsl_eigen_gensymmv_alloc
2241 This function allocates a workspace for computing eigenvalues and
2242 eigenvectors of N-by-N real generalized symmetric-definite
2243 eigensystems. The size of the workspace is O(4n).
2245 -- Function: void gsl_eigen_gensymmv_free
2246 (gsl_eigen_gensymmv_workspace * W)
2247 This function frees the memory associated with the workspace W.
2249 -- Function: int gsl_eigen_gensymmv (gsl_matrix * A, gsl_matrix * B,
2250 gsl_vector * EVAL, gsl_matrix * EVEC,
2251 gsl_eigen_gensymmv_workspace * W)
2252 This function computes the eigenvalues and eigenvectors of the real
2253 generalized symmetric-definite matrix pair (A, B), and stores them
2254 in EVAL and EVEC respectively. The computed eigenvectors are
2255 normalized to have unit magnitude. On output, B contains its
2256 Cholesky decomposition and A is destroyed.
2259 File: gsl-ref.info, Node: Complex Generalized Hermitian-Definite Eigensystems, Next: Real Generalized Nonsymmetric Eigensystems, Prev: Real Generalized Symmetric-Definite Eigensystems, Up: Eigensystems
2261 14.5 Complex Generalized Hermitian-Definite Eigensystems
2262 ========================================================
2264 The complex generalized hermitian-definite eigenvalue problem is to find
2265 eigenvalues \lambda and eigenvectors x such that
2267 where A and B are hermitian matrices, and B is positive-definite.
2268 Similarly to the real case, this can be reduced to C y = \lambda y where
2269 C = L^{-1} A L^{-H} is hermitian, and y = L^H x. The standard hermitian
2270 eigensolver can be applied to the matrix C. The resulting eigenvectors
2271 are backtransformed to find the vectors of the original problem. The
2272 eigenvalues of the generalized hermitian-definite eigenproblem are
2275 -- Function: gsl_eigen_genherm_workspace * gsl_eigen_genherm_alloc
2277 This function allocates a workspace for computing eigenvalues of
2278 N-by-N complex generalized hermitian-definite eigensystems. The
2279 size of the workspace is O(3n).
2281 -- Function: void gsl_eigen_genherm_free (gsl_eigen_genherm_workspace
2283 This function frees the memory associated with the workspace W.
2285 -- Function: int gsl_eigen_genherm (gsl_matrix_complex * A,
2286 gsl_matrix_complex * B, gsl_vector * EVAL,
2287 gsl_eigen_genherm_workspace * W)
2288 This function computes the eigenvalues of the complex generalized
2289 hermitian-definite matrix pair (A, B), and stores them in EVAL,
2290 using the method outlined above. On output, B contains its
2291 Cholesky decomposition and A is destroyed.
2293 -- Function: gsl_eigen_genhermv_workspace * gsl_eigen_genhermv_alloc
2295 This function allocates a workspace for computing eigenvalues and
2296 eigenvectors of N-by-N complex generalized hermitian-definite
2297 eigensystems. The size of the workspace is O(5n).
2299 -- Function: void gsl_eigen_genhermv_free
2300 (gsl_eigen_genhermv_workspace * W)
2301 This function frees the memory associated with the workspace W.
2303 -- Function: int gsl_eigen_genhermv (gsl_matrix_complex * A,
2304 gsl_matrix_complex * B, gsl_vector * EVAL, gsl_matrix_complex
2305 * EVEC, gsl_eigen_genhermv_workspace * W)
2306 This function computes the eigenvalues and eigenvectors of the
2307 complex generalized hermitian-definite matrix pair (A, B), and
2308 stores them in EVAL and EVEC respectively. The computed
2309 eigenvectors are normalized to have unit magnitude. On output, B
2310 contains its Cholesky decomposition and A is destroyed.
2313 File: gsl-ref.info, Node: Real Generalized Nonsymmetric Eigensystems, Next: Sorting Eigenvalues and Eigenvectors, Prev: Complex Generalized Hermitian-Definite Eigensystems, Up: Eigensystems
2315 14.6 Real Generalized Nonsymmetric Eigensystems
2316 ===============================================
2318 Given two square matrices (A, B), the generalized nonsymmetric
2319 eigenvalue problem is to find eigenvalues \lambda and eigenvectors x
2322 We may also define the problem as finding eigenvalues \mu and
2323 eigenvectors y such that
2325 Note that these two problems are equivalent (with \lambda = 1/\mu)
2326 if neither \lambda nor \mu is zero. If say, \lambda is zero, then it is
2327 still a well defined eigenproblem, but its alternate problem involving
2328 \mu is not. Therefore, to allow for zero (and infinite) eigenvalues,
2329 the problem which is actually solved is
2330 \beta A x = \alpha B x
2331 The eigensolver routines below will return two values \alpha and
2332 \beta and leave it to the user to perform the divisions \lambda =
2333 \alpha / \beta and \mu = \beta / \alpha.
2335 If the determinant of the matrix pencil A - \lambda B is zero for
2336 all \lambda, the problem is said to be singular; otherwise it is called
2337 regular. Singularity normally leads to some \alpha = \beta = 0 which
2338 means the eigenproblem is ill-conditioned and generally does not have
2339 well defined eigenvalue solutions. The routines below are intended for
2340 regular matrix pencils and could yield unpredictable results when
2341 applied to singular pencils.
2343 The solution of the real generalized nonsymmetric eigensystem
2344 problem for a matrix pair (A, B) involves computing the generalized
2348 where Q and Z are orthogonal matrices of left and right Schur
2349 vectors respectively, and (S, T) is the generalized Schur form whose
2350 diagonal elements give the \alpha and \beta values. The algorithm used
2351 is the QZ method due to Moler and Stewart (see references).
2353 -- Function: gsl_eigen_gen_workspace * gsl_eigen_gen_alloc (const
2355 This function allocates a workspace for computing eigenvalues of
2356 N-by-N real generalized nonsymmetric eigensystems. The size of the
2359 -- Function: void gsl_eigen_gen_free (gsl_eigen_gen_workspace * W)
2360 This function frees the memory associated with the workspace W.
2362 -- Function: void gsl_eigen_gen_params (const int COMPUTE_S, const int
2363 COMPUTE_T, const int BALANCE, gsl_eigen_gen_workspace * W)
2364 This function sets some parameters which determine how the
2365 eigenvalue problem is solved in subsequent calls to
2368 If COMPUTE_S is set to 1, the full Schur form S will be computed
2369 by `gsl_eigen_gen'. If it is set to 0, S will not be computed
2370 (this is the default setting). S is a quasi upper triangular
2371 matrix with 1-by-1 and 2-by-2 blocks on its diagonal. 1-by-1
2372 blocks correspond to real eigenvalues, and 2-by-2 blocks
2373 correspond to complex eigenvalues.
2375 If COMPUTE_T is set to 1, the full Schur form T will be computed
2376 by `gsl_eigen_gen'. If it is set to 0, T will not be computed
2377 (this is the default setting). T is an upper triangular matrix
2378 with non-negative elements on its diagonal. Any 2-by-2 blocks in
2379 S will correspond to a 2-by-2 diagonal block in T.
2381 The BALANCE parameter is currently ignored, since generalized
2382 balancing is not yet implemented.
2384 -- Function: int gsl_eigen_gen (gsl_matrix * A, gsl_matrix * B,
2385 gsl_vector_complex * ALPHA, gsl_vector * BETA,
2386 gsl_eigen_gen_workspace * W)
2387 This function computes the eigenvalues of the real generalized
2388 nonsymmetric matrix pair (A, B), and stores them as pairs in
2389 (ALPHA, BETA), where ALPHA is complex and BETA is real. If \beta_i
2390 is non-zero, then \lambda = \alpha_i / \beta_i is an eigenvalue.
2391 Likewise, if \alpha_i is non-zero, then \mu = \beta_i / \alpha_i
2392 is an eigenvalue of the alternate problem \mu A y = B y. The
2393 elements of BETA are normalized to be non-negative.
2395 If S is desired, it is stored in A on output. If T is desired, it
2396 is stored in B on output. The ordering of eigenvalues in (ALPHA,
2397 BETA) follows the ordering of the diagonal blocks in the Schur
2398 forms S and T. In rare cases, this function may fail to find all
2399 eigenvalues. If this occurs, an error code is returned.
2401 -- Function: int gsl_eigen_gen_QZ (gsl_matrix * A, gsl_matrix * B,
2402 gsl_vector_complex * ALPHA, gsl_vector * BETA, gsl_matrix *
2403 Q, gsl_matrix * Z, gsl_eigen_gen_workspace * W)
2404 This function is identical to `gsl_eigen_gen' except it also
2405 computes the left and right Schur vectors and stores them into Q
2408 -- Function: gsl_eigen_genv_workspace * gsl_eigen_genv_alloc (const
2410 This function allocates a workspace for computing eigenvalues and
2411 eigenvectors of N-by-N real generalized nonsymmetric eigensystems.
2412 The size of the workspace is O(7n).
2414 -- Function: void gsl_eigen_genv_free (gsl_eigen_genv_workspace * W)
2415 This function frees the memory associated with the workspace W.
2417 -- Function: int gsl_eigen_genv (gsl_matrix * A, gsl_matrix * B,
2418 gsl_vector_complex * ALPHA, gsl_vector * BETA,
2419 gsl_matrix_complex * EVEC, gsl_eigen_genv_workspace * W)
2420 This function computes eigenvalues and right eigenvectors of the
2421 N-by-N real generalized nonsymmetric matrix pair (A, B). The
2422 eigenvalues are stored in (ALPHA, BETA) and the eigenvectors are
2423 stored in EVEC. It first calls `gsl_eigen_gen' to compute the
2424 eigenvalues, Schur forms, and Schur vectors. Then it finds
2425 eigenvectors of the Schur forms and backtransforms them using the
2426 Schur vectors. The Schur vectors are destroyed in the process, but
2427 can be saved by using `gsl_eigen_genv_QZ'. The computed
2428 eigenvectors are normalized to have unit magnitude. On output, (A,
2429 B) contains the generalized Schur form (S, T). If `gsl_eigen_gen'
2430 fails, no eigenvectors are computed, and an error code is returned.
2432 -- Function: int gsl_eigen_genv_QZ (gsl_matrix * A, gsl_matrix * B,
2433 gsl_vector_complex * ALPHA, gsl_vector * BETA,
2434 gsl_matrix_complex * EVEC, gsl_matrix * Q, gsl_matrix * Z,
2435 gsl_eigen_genv_workspace * W)
2436 This function is identical to `gsl_eigen_genv' except it also
2437 computes the left and right Schur vectors and stores them into Q
2441 File: gsl-ref.info, Node: Sorting Eigenvalues and Eigenvectors, Next: Eigenvalue and Eigenvector Examples, Prev: Real Generalized Nonsymmetric Eigensystems, Up: Eigensystems
2443 14.7 Sorting Eigenvalues and Eigenvectors
2444 =========================================
2446 -- Function: int gsl_eigen_symmv_sort (gsl_vector * EVAL, gsl_matrix *
2447 EVEC, gsl_eigen_sort_t SORT_TYPE)
2448 This function simultaneously sorts the eigenvalues stored in the
2449 vector EVAL and the corresponding real eigenvectors stored in the
2450 columns of the matrix EVEC into ascending or descending order
2451 according to the value of the parameter SORT_TYPE,
2453 `GSL_EIGEN_SORT_VAL_ASC'
2454 ascending order in numerical value
2456 `GSL_EIGEN_SORT_VAL_DESC'
2457 descending order in numerical value
2459 `GSL_EIGEN_SORT_ABS_ASC'
2460 ascending order in magnitude
2462 `GSL_EIGEN_SORT_ABS_DESC'
2463 descending order in magnitude
2466 -- Function: int gsl_eigen_hermv_sort (gsl_vector * EVAL,
2467 gsl_matrix_complex * EVEC, gsl_eigen_sort_t SORT_TYPE)
2468 This function simultaneously sorts the eigenvalues stored in the
2469 vector EVAL and the corresponding complex eigenvectors stored in
2470 the columns of the matrix EVEC into ascending or descending order
2471 according to the value of the parameter SORT_TYPE as shown above.
2473 -- Function: int gsl_eigen_nonsymmv_sort (gsl_vector_complex * EVAL,
2474 gsl_matrix_complex * EVEC, gsl_eigen_sort_t SORT_TYPE)
2475 This function simultaneously sorts the eigenvalues stored in the
2476 vector EVAL and the corresponding complex eigenvectors stored in
2477 the columns of the matrix EVEC into ascending or descending order
2478 according to the value of the parameter SORT_TYPE as shown above.
2479 Only `GSL_EIGEN_SORT_ABS_ASC' and `GSL_EIGEN_SORT_ABS_DESC' are
2480 supported due to the eigenvalues being complex.
2482 -- Function: int gsl_eigen_gensymmv_sort (gsl_vector * EVAL,
2483 gsl_matrix * EVEC, gsl_eigen_sort_t SORT_TYPE)
2484 This function simultaneously sorts the eigenvalues stored in the
2485 vector EVAL and the corresponding real eigenvectors stored in the
2486 columns of the matrix EVEC into ascending or descending order
2487 according to the value of the parameter SORT_TYPE as shown above.
2489 -- Function: int gsl_eigen_genhermv_sort (gsl_vector * EVAL,
2490 gsl_matrix_complex * EVEC, gsl_eigen_sort_t SORT_TYPE)
2491 This function simultaneously sorts the eigenvalues stored in the
2492 vector EVAL and the corresponding complex eigenvectors stored in
2493 the columns of the matrix EVEC into ascending or descending order
2494 according to the value of the parameter SORT_TYPE as shown above.
2496 -- Function: int gsl_eigen_genv_sort (gsl_vector_complex * ALPHA,
2497 gsl_vector * BETA, gsl_matrix_complex * EVEC,
2498 gsl_eigen_sort_t SORT_TYPE)
2499 This function simultaneously sorts the eigenvalues stored in the
2500 vectors (ALPHA, BETA) and the corresponding complex eigenvectors
2501 stored in the columns of the matrix EVEC into ascending or
2502 descending order according to the value of the parameter SORT_TYPE
2503 as shown above. Only `GSL_EIGEN_SORT_ABS_ASC' and
2504 `GSL_EIGEN_SORT_ABS_DESC' are supported due to the eigenvalues
2508 File: gsl-ref.info, Node: Eigenvalue and Eigenvector Examples, Next: Eigenvalue and Eigenvector References, Prev: Sorting Eigenvalues and Eigenvectors, Up: Eigensystems
2513 The following program computes the eigenvalues and eigenvectors of the
2514 4-th order Hilbert matrix, H(i,j) = 1/(i + j + 1).
2517 #include <gsl/gsl_math.h>
2518 #include <gsl/gsl_eigen.h>
2523 double data[] = { 1.0 , 1/2.0, 1/3.0, 1/4.0,
2524 1/2.0, 1/3.0, 1/4.0, 1/5.0,
2525 1/3.0, 1/4.0, 1/5.0, 1/6.0,
2526 1/4.0, 1/5.0, 1/6.0, 1/7.0 };
2529 = gsl_matrix_view_array (data, 4, 4);
2531 gsl_vector *eval = gsl_vector_alloc (4);
2532 gsl_matrix *evec = gsl_matrix_alloc (4, 4);
2534 gsl_eigen_symmv_workspace * w =
2535 gsl_eigen_symmv_alloc (4);
2537 gsl_eigen_symmv (&m.matrix, eval, evec, w);
2539 gsl_eigen_symmv_free (w);
2541 gsl_eigen_symmv_sort (eval, evec,
2542 GSL_EIGEN_SORT_ABS_ASC);
2547 for (i = 0; i < 4; i++)
2550 = gsl_vector_get (eval, i);
2551 gsl_vector_view evec_i
2552 = gsl_matrix_column (evec, i);
2554 printf ("eigenvalue = %g\n", eval_i);
2555 printf ("eigenvector = \n");
2556 gsl_vector_fprintf (stdout,
2557 &evec_i.vector, "%g");
2561 gsl_vector_free (eval);
2562 gsl_matrix_free (evec);
2567 Here is the beginning of the output from the program,
2570 eigenvalue = 9.67023e-05
2578 This can be compared with the corresponding output from GNU OCTAVE,
2580 octave> [v,d] = eig(hilb(4));
2592 0.029193 0.179186 -0.582076 0.792608
2593 -0.328712 -0.741918 0.370502 0.451923
2594 0.791411 0.100228 0.509579 0.322416
2595 -0.514553 0.638283 0.514048 0.252161
2597 Note that the eigenvectors can differ by a change of sign, since the
2598 sign of an eigenvector is arbitrary.
2600 The following program illustrates the use of the nonsymmetric
2601 eigensolver, by computing the eigenvalues and eigenvectors of the
2602 Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4).
2605 #include <gsl/gsl_math.h>
2606 #include <gsl/gsl_eigen.h>
2611 double data[] = { -1.0, 1.0, -1.0, 1.0,
2612 -8.0, 4.0, -2.0, 1.0,
2613 27.0, 9.0, 3.0, 1.0,
2614 64.0, 16.0, 4.0, 1.0 };
2617 = gsl_matrix_view_array (data, 4, 4);
2619 gsl_vector_complex *eval = gsl_vector_complex_alloc (4);
2620 gsl_matrix_complex *evec = gsl_matrix_complex_alloc (4, 4);
2622 gsl_eigen_nonsymmv_workspace * w =
2623 gsl_eigen_nonsymmv_alloc (4);
2625 gsl_eigen_nonsymmv (&m.matrix, eval, evec, w);
2627 gsl_eigen_nonsymmv_free (w);
2629 gsl_eigen_nonsymmv_sort (eval, evec,
2630 GSL_EIGEN_SORT_ABS_DESC);
2635 for (i = 0; i < 4; i++)
2638 = gsl_vector_complex_get (eval, i);
2639 gsl_vector_complex_view evec_i
2640 = gsl_matrix_complex_column (evec, i);
2642 printf ("eigenvalue = %g + %gi\n",
2643 GSL_REAL(eval_i), GSL_IMAG(eval_i));
2644 printf ("eigenvector = \n");
2645 for (j = 0; j < 4; ++j)
2647 gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
2648 printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));
2653 gsl_vector_complex_free(eval);
2654 gsl_matrix_complex_free(evec);
2659 Here is the beginning of the output from the program,
2662 eigenvalue = -6.41391 + 0i
2668 eigenvalue = 5.54555 + 3.08545i
2670 -0.043487 + -0.0076308i
2671 0.0642377 + -0.142127i
2672 -0.515253 + 0.0405118i
2673 -0.840592 + -0.00148565i
2676 This can be compared with the corresponding output from GNU OCTAVE,
2678 octave> [v,d] = eig(vander([-1 -2 3 4]));
2690 Columns 1 through 3:
2692 -0.09988 + 0.00000i -0.04350 - 0.00755i -0.04350 + 0.00755i
2693 -0.11125 + 0.00000i 0.06399 - 0.14224i 0.06399 + 0.14224i
2694 0.29250 + 0.00000i -0.51518 + 0.04142i -0.51518 - 0.04142i
2695 0.94451 + 0.00000i -0.84059 + 0.00000i -0.84059 - 0.00000i
2703 Note that the eigenvectors corresponding to the eigenvalue 5.54555 +
2704 3.08545i are slightly different. This is because they differ by the
2705 multiplicative constant 0.9999984 + 0.0017674i which has magnitude 1.
2708 File: gsl-ref.info, Node: Eigenvalue and Eigenvector References, Prev: Eigenvalue and Eigenvector Examples, Up: Eigensystems
2710 14.9 References and Further Reading
2711 ===================================
2713 Further information on the algorithms described in this section can be
2714 found in the following book,
2716 G. H. Golub, C. F. Van Loan, `Matrix Computations' (3rd Ed, 1996),
2717 Johns Hopkins University Press, ISBN 0-8018-5414-8.
2719 Further information on the generalized eigensystems QZ algorithm can be
2720 found in this paper,
2722 C. Moler, G. Stewart, "An Algorithm for Generalized Matrix
2723 Eigenvalue Problems," SIAM J. Numer. Anal., Vol 10, No 2, 1973.
2725 Eigensystem routines for very large matrices can be found in the
2726 Fortran library LAPACK. The LAPACK library is described in,
2728 `LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM,
2731 `http://www.netlib.org/lapack'
2733 The LAPACK source code can be found at the website above along with an
2734 online copy of the users guide.
2737 File: gsl-ref.info, Node: Fast Fourier Transforms, Next: Numerical Integration, Prev: Eigensystems, Up: Top
2739 15 Fast Fourier Transforms (FFTs)
2740 *********************************
2742 This chapter describes functions for performing Fast Fourier Transforms
2743 (FFTs). The library includes radix-2 routines (for lengths which are a
2744 power of two) and mixed-radix routines (which work for any length). For
2745 efficiency there are separate versions of the routines for real data and
2746 for complex data. The mixed-radix routines are a reimplementation of
2747 the FFTPACK library of Paul Swarztrauber. Fortran code for FFTPACK is
2748 available on Netlib (FFTPACK also includes some routines for sine and
2749 cosine transforms but these are currently not available in GSL). For
2750 details and derivations of the underlying algorithms consult the
2751 document `GSL FFT Algorithms' (*note FFT References and Further
2756 * Mathematical Definitions::
2757 * Overview of complex data FFTs::
2758 * Radix-2 FFT routines for complex data::
2759 * Mixed-radix FFT routines for complex data::
2760 * Overview of real data FFTs::
2761 * Radix-2 FFT routines for real data::
2762 * Mixed-radix FFT routines for real data::
2763 * FFT References and Further Reading::
2766 File: gsl-ref.info, Node: Mathematical Definitions, Next: Overview of complex data FFTs, Up: Fast Fourier Transforms
2768 15.1 Mathematical Definitions
2769 =============================
2771 Fast Fourier Transforms are efficient algorithms for calculating the
2772 discrete fourier transform (DFT),
2774 x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N)
2776 The DFT usually arises as an approximation to the continuous fourier
2777 transform when functions are sampled at discrete intervals in space or
2778 time. The naive evaluation of the discrete fourier transform is a
2779 matrix-vector multiplication W\vec{z}. A general matrix-vector
2780 multiplication takes O(N^2) operations for N data-points. Fast fourier
2781 transform algorithms use a divide-and-conquer strategy to factorize the
2782 matrix W into smaller sub-matrices, corresponding to the integer
2783 factors of the length N. If N can be factorized into a product of
2784 integers f_1 f_2 ... f_n then the DFT can be computed in O(N \sum f_i)
2785 operations. For a radix-2 FFT this gives an operation count of O(N
2788 All the FFT functions offer three types of transform: forwards,
2789 inverse and backwards, based on the same mathematical definitions. The
2790 definition of the "forward fourier transform", x = FFT(z), is,
2792 x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N)
2794 and the definition of the "inverse fourier transform", x = IFFT(z), is,
2796 z_j = {1 \over N} \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N).
2798 The factor of 1/N makes this a true inverse. For example, a call to
2799 `gsl_fft_complex_forward' followed by a call to
2800 `gsl_fft_complex_inverse' should return the original data (within
2803 In general there are two possible choices for the sign of the
2804 exponential in the transform/ inverse-transform pair. GSL follows the
2805 same convention as FFTPACK, using a negative exponential for the forward
2806 transform. The advantage of this convention is that the inverse
2807 transform recreates the original function with simple fourier
2808 synthesis. Numerical Recipes uses the opposite convention, a positive
2809 exponential in the forward transform.
2811 The "backwards FFT" is simply our terminology for an unscaled
2812 version of the inverse FFT,
2814 z^{backwards}_j = \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N).
2816 When the overall scale of the result is unimportant it is often
2817 convenient to use the backwards FFT instead of the inverse to save
2818 unnecessary divisions.
2821 File: gsl-ref.info, Node: Overview of complex data FFTs, Next: Radix-2 FFT routines for complex data, Prev: Mathematical Definitions, Up: Fast Fourier Transforms
2823 15.2 Overview of complex data FFTs
2824 ==================================
2826 The inputs and outputs for the complex FFT routines are "packed arrays"
2827 of floating point numbers. In a packed array the real and imaginary
2828 parts of each complex number are placed in alternate neighboring
2829 elements. For example, the following definition of a packed array of
2833 gsl_complex_packed_array data = x;
2835 can be used to hold an array of three complex numbers, `z[3]', in the
2845 The array indices for the data have the same ordering as those in the
2846 definition of the DFT--i.e. there are no index transformations or
2847 permutations of the data.
2849 A "stride" parameter allows the user to perform transforms on the
2850 elements `z[stride*i]' instead of `z[i]'. A stride greater than 1 can
2851 be used to take an in-place FFT of the column of a matrix. A stride of
2852 1 accesses the array without any additional spacing between elements.
2854 To perform an FFT on a vector argument, such as `gsl_vector_complex
2855 * v', use the following definitions (or their equivalents) when calling
2856 the functions described in this chapter:
2858 gsl_complex_packed_array data = v->data;
2859 size_t stride = v->stride;
2862 For physical applications it is important to remember that the index
2863 appearing in the DFT does not correspond directly to a physical
2864 frequency. If the time-step of the DFT is \Delta then the
2865 frequency-domain includes both positive and negative frequencies,
2866 ranging from -1/(2\Delta) through 0 to +1/(2\Delta). The positive
2867 frequencies are stored from the beginning of the array up to the
2868 middle, and the negative frequencies are stored backwards from the end
2871 Here is a table which shows the layout of the array DATA, and the
2872 correspondence between the time-domain data z, and the frequency-domain
2878 1 z(t = 1) x(f = 1/(N Delta))
2879 2 z(t = 2) x(f = 2/(N Delta))
2880 . ........ ..................
2881 N/2 z(t = N/2) x(f = +1/(2 Delta),
2883 . ........ ..................
2884 N-3 z(t = N-3) x(f = -3/(N Delta))
2885 N-2 z(t = N-2) x(f = -2/(N Delta))
2886 N-1 z(t = N-1) x(f = -1/(N Delta))
2888 When N is even the location N/2 contains the most positive and negative
2889 frequencies (+1/(2 \Delta), -1/(2 \Delta)) which are equivalent. If N
2890 is odd then general structure of the table above still applies, but N/2
2894 File: gsl-ref.info, Node: Radix-2 FFT routines for complex data, Next: Mixed-radix FFT routines for complex data, Prev: Overview of complex data FFTs, Up: Fast Fourier Transforms
2896 15.3 Radix-2 FFT routines for complex data
2897 ==========================================
2899 The radix-2 algorithms described in this section are simple and compact,
2900 although not necessarily the most efficient. They use the Cooley-Tukey
2901 algorithm to compute in-place complex FFTs for lengths which are a power
2902 of 2--no additional storage is required. The corresponding
2903 self-sorting mixed-radix routines offer better performance at the
2904 expense of requiring additional working space.
2906 All the functions described in this section are declared in the
2907 header file `gsl_fft_complex.h'.
2909 -- Function: int gsl_fft_complex_radix2_forward
2910 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2911 -- Function: int gsl_fft_complex_radix2_transform
2912 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N,
2913 gsl_fft_direction SIGN)
2914 -- Function: int gsl_fft_complex_radix2_backward
2915 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2916 -- Function: int gsl_fft_complex_radix2_inverse
2917 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2918 These functions compute forward, backward and inverse FFTs of
2919 length N with stride STRIDE, on the packed complex array DATA
2920 using an in-place radix-2 decimation-in-time algorithm. The
2921 length of the transform N is restricted to powers of two. For the
2922 `transform' version of the function the SIGN argument can be
2923 either `forward' (-1) or `backward' (+1).
2925 The functions return a value of `GSL_SUCCESS' if no errors were
2926 detected, or `GSL_EDOM' if the length of the data N is not a power
2929 -- Function: int gsl_fft_complex_radix2_dif_forward
2930 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2931 -- Function: int gsl_fft_complex_radix2_dif_transform
2932 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N,
2933 gsl_fft_direction SIGN)
2934 -- Function: int gsl_fft_complex_radix2_dif_backward
2935 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2936 -- Function: int gsl_fft_complex_radix2_dif_inverse
2937 (gsl_complex_packed_array DATA, size_t STRIDE, size_t N)
2938 These are decimation-in-frequency versions of the radix-2 FFT
2942 Here is an example program which computes the FFT of a short pulse
2943 in a sample of length 128. To make the resulting fourier transform
2944 real the pulse is defined for equal positive and negative times (-10
2945 ... 10), where the negative times wrap around the end of the array.
2949 #include <gsl/gsl_errno.h>
2950 #include <gsl/gsl_fft_complex.h>
2952 #define REAL(z,i) ((z)[2*(i)])
2953 #define IMAG(z,i) ((z)[2*(i)+1])
2958 int i; double data[2*128];
2960 for (i = 0; i < 128; i++)
2962 REAL(data,i) = 0.0; IMAG(data,i) = 0.0;
2967 for (i = 1; i <= 10; i++)
2969 REAL(data,i) = REAL(data,128-i) = 1.0;
2972 for (i = 0; i < 128; i++)
2974 printf ("%d %e %e\n", i,
2975 REAL(data,i), IMAG(data,i));
2979 gsl_fft_complex_radix2_forward (data, 1, 128);
2981 for (i = 0; i < 128; i++)
2983 printf ("%d %e %e\n", i,
2984 REAL(data,i)/sqrt(128),
2985 IMAG(data,i)/sqrt(128));
2991 Note that we have assumed that the program is using the default error
2992 handler (which calls `abort' for any errors). If you are not using a
2993 safe error handler you would need to check the return status of
2994 `gsl_fft_complex_radix2_forward'.
2996 The transformed data is rescaled by 1/\sqrt N so that it fits on the
2997 same plot as the input. Only the real part is shown, by the choice of
2998 the input data the imaginary part is zero. Allowing for the
2999 wrap-around of negative times at t=128, and working in units of k/N,
3000 the DFT approximates the continuum fourier transform, giving a
3001 modulated sine function.
3004 File: gsl-ref.info, Node: Mixed-radix FFT routines for complex data, Next: Overview of real data FFTs, Prev: Radix-2 FFT routines for complex data, Up: Fast Fourier Transforms
3006 15.4 Mixed-radix FFT routines for complex data
3007 ==============================================
3009 This section describes mixed-radix FFT algorithms for complex data. The
3010 mixed-radix functions work for FFTs of any length. They are a
3011 reimplementation of Paul Swarztrauber's Fortran FFTPACK library. The
3012 theory is explained in the review article `Self-sorting Mixed-radix
3013 FFTs' by Clive Temperton. The routines here use the same indexing
3014 scheme and basic algorithms as FFTPACK.
3016 The mixed-radix algorithm is based on sub-transform modules--highly
3017 optimized small length FFTs which are combined to create larger FFTs.
3018 There are efficient modules for factors of 2, 3, 4, 5, 6 and 7. The
3019 modules for the composite factors of 4 and 6 are faster than combining
3020 the modules for 2*2 and 2*3.
3022 For factors which are not implemented as modules there is a
3023 fall-back to a general length-n module which uses Singleton's method for
3024 efficiently computing a DFT. This module is O(n^2), and slower than a
3025 dedicated module would be but works for any length n. Of course,
3026 lengths which use the general length-n module will still be factorized
3027 as much as possible. For example, a length of 143 will be factorized
3028 into 11*13. Large prime factors are the worst case scenario, e.g. as
3029 found in n=2*3*99991, and should be avoided because their O(n^2)
3030 scaling will dominate the run-time (consult the document `GSL FFT
3031 Algorithms' included in the GSL distribution if you encounter this
3034 The mixed-radix initialization function
3035 `gsl_fft_complex_wavetable_alloc' returns the list of factors chosen by
3036 the library for a given length N. It can be used to check how well the
3037 length has been factorized, and estimate the run-time. To a first
3038 approximation the run-time scales as N \sum f_i, where the f_i are the
3039 factors of N. For programs under user control you may wish to issue a
3040 warning that the transform will be slow when the length is poorly
3041 factorized. If you frequently encounter data lengths which cannot be
3042 factorized using the existing small-prime modules consult `GSL FFT
3043 Algorithms' for details on adding support for other factors.
3045 All the functions described in this section are declared in the
3046 header file `gsl_fft_complex.h'.
3048 -- Function: gsl_fft_complex_wavetable *
3049 gsl_fft_complex_wavetable_alloc (size_t N)
3050 This function prepares a trigonometric lookup table for a complex
3051 FFT of length N. The function returns a pointer to the newly
3052 allocated `gsl_fft_complex_wavetable' if no errors were detected,
3053 and a null pointer in the case of error. The length N is
3054 factorized into a product of subtransforms, and the factors and
3055 their trigonometric coefficients are stored in the wavetable. The
3056 trigonometric coefficients are computed using direct calls to
3057 `sin' and `cos', for accuracy. Recursion relations could be used
3058 to compute the lookup table faster, but if an application performs
3059 many FFTs of the same length then this computation is a one-off
3060 overhead which does not affect the final throughput.
3062 The wavetable structure can be used repeatedly for any transform
3063 of the same length. The table is not modified by calls to any of
3064 the other FFT functions. The same wavetable can be used for both
3065 forward and backward (or inverse) transforms of a given length.
3067 -- Function: void gsl_fft_complex_wavetable_free
3068 (gsl_fft_complex_wavetable * WAVETABLE)
3069 This function frees the memory associated with the wavetable
3070 WAVETABLE. The wavetable can be freed if no further FFTs of the
3071 same length will be needed.
3073 These functions operate on a `gsl_fft_complex_wavetable' structure
3074 which contains internal parameters for the FFT. It is not necessary to
3075 set any of the components directly but it can sometimes be useful to
3076 examine them. For example, the chosen factorization of the FFT length
3077 is given and can be used to provide an estimate of the run-time or
3078 numerical error. The wavetable structure is declared in the header file
3079 `gsl_fft_complex.h'.
3081 -- Data Type: gsl_fft_complex_wavetable
3082 This is a structure that holds the factorization and trigonometric
3083 lookup tables for the mixed radix fft algorithm. It has the
3084 following components:
3087 This is the number of complex data points
3090 This is the number of factors that the length `n' was
3094 This is the array of factors. Only the first `nf' elements
3097 `gsl_complex * trig'
3098 This is a pointer to a preallocated trigonometric lookup
3099 table of `n' complex elements.
3101 `gsl_complex * twiddle[64]'
3102 This is an array of pointers into `trig', giving the twiddle
3103 factors for each pass.
3105 The mixed radix algorithms require additional working space to hold the
3106 intermediate steps of the transform.
3108 -- Function: gsl_fft_complex_workspace *
3109 gsl_fft_complex_workspace_alloc (size_t N)
3110 This function allocates a workspace for a complex transform of
3113 -- Function: void gsl_fft_complex_workspace_free
3114 (gsl_fft_complex_workspace * WORKSPACE)
3115 This function frees the memory associated with the workspace
3116 WORKSPACE. The workspace can be freed if no further FFTs of the
3117 same length will be needed.
3119 The following functions compute the transform,
3121 -- Function: int gsl_fft_complex_forward (gsl_complex_packed_array
3122 DATA, size_t STRIDE, size_t N, const
3123 gsl_fft_complex_wavetable * WAVETABLE,
3124 gsl_fft_complex_workspace * WORK)
3125 -- Function: int gsl_fft_complex_transform (gsl_complex_packed_array
3126 DATA, size_t STRIDE, size_t N, const
3127 gsl_fft_complex_wavetable * WAVETABLE,
3128 gsl_fft_complex_workspace * WORK, gsl_fft_direction SIGN)
3129 -- Function: int gsl_fft_complex_backward (gsl_complex_packed_array
3130 DATA, size_t STRIDE, size_t N, const
3131 gsl_fft_complex_wavetable * WAVETABLE,
3132 gsl_fft_complex_workspace * WORK)
3133 -- Function: int gsl_fft_complex_inverse (gsl_complex_packed_array
3134 DATA, size_t STRIDE, size_t N, const
3135 gsl_fft_complex_wavetable * WAVETABLE,
3136 gsl_fft_complex_workspace * WORK)
3137 These functions compute forward, backward and inverse FFTs of
3138 length N with stride STRIDE, on the packed complex array DATA,
3139 using a mixed radix decimation-in-frequency algorithm. There is
3140 no restriction on the length N. Efficient modules are provided
3141 for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining
3142 factors are computed with a slow, O(n^2), general-n module. The
3143 caller must supply a WAVETABLE containing the trigonometric lookup
3144 tables and a workspace WORK. For the `transform' version of the
3145 function the SIGN argument can be either `forward' (-1) or
3148 The functions return a value of `0' if no errors were detected. The
3149 following `gsl_errno' conditions are defined for these functions:
3152 The length of the data N is not a positive integer (i.e. N is
3156 The length of the data N and the length used to compute the
3157 given WAVETABLE do not match.
3159 Here is an example program which computes the FFT of a short pulse
3160 in a sample of length 630 (=2*3*3*5*7) using the mixed-radix algorithm.
3164 #include <gsl/gsl_errno.h>
3165 #include <gsl/gsl_fft_complex.h>
3167 #define REAL(z,i) ((z)[2*(i)])
3168 #define IMAG(z,i) ((z)[2*(i)+1])
3177 gsl_fft_complex_wavetable * wavetable;
3178 gsl_fft_complex_workspace * workspace;
3180 for (i = 0; i < n; i++)
3188 for (i = 1; i <= 10; i++)
3190 REAL(data,i) = REAL(data,n-i) = 1.0;
3193 for (i = 0; i < n; i++)
3195 printf ("%d: %e %e\n", i, REAL(data,i),
3200 wavetable = gsl_fft_complex_wavetable_alloc (n);
3201 workspace = gsl_fft_complex_workspace_alloc (n);
3203 for (i = 0; i < wavetable->nf; i++)
3205 printf ("# factor %d: %d\n", i,
3206 wavetable->factor[i]);
3209 gsl_fft_complex_forward (data, 1, n,
3210 wavetable, workspace);
3212 for (i = 0; i < n; i++)
3214 printf ("%d: %e %e\n", i, REAL(data,i),
3218 gsl_fft_complex_wavetable_free (wavetable);
3219 gsl_fft_complex_workspace_free (workspace);
3223 Note that we have assumed that the program is using the default `gsl'
3224 error handler (which calls `abort' for any errors). If you are not
3225 using a safe error handler you would need to check the return status of
3226 all the `gsl' routines.
3229 File: gsl-ref.info, Node: Overview of real data FFTs, Next: Radix-2 FFT routines for real data, Prev: Mixed-radix FFT routines for complex data, Up: Fast Fourier Transforms
3231 15.5 Overview of real data FFTs
3232 ===============================
3234 The functions for real data are similar to those for complex data.
3235 However, there is an important difference between forward and inverse
3236 transforms. The fourier transform of a real sequence is not real. It
3237 is a complex sequence with a special symmetry:
3241 A sequence with this symmetry is called "conjugate-complex" or
3242 "half-complex". This different structure requires different storage
3243 layouts for the forward transform (from real to half-complex) and
3244 inverse transform (from half-complex back to real). As a consequence
3245 the routines are divided into two sets: functions in `gsl_fft_real'
3246 which operate on real sequences and functions in `gsl_fft_halfcomplex'
3247 which operate on half-complex sequences.
3249 Functions in `gsl_fft_real' compute the frequency coefficients of a
3250 real sequence. The half-complex coefficients c of a real sequence x
3251 are given by fourier analysis,
3253 c_k = \sum_{j=0}^{N-1} x_j \exp(-2 \pi i j k /N)
3255 Functions in `gsl_fft_halfcomplex' compute inverse or backwards
3256 transforms. They reconstruct real sequences by fourier synthesis from
3257 their half-complex frequency coefficients, c,
3259 x_j = {1 \over N} \sum_{k=0}^{N-1} c_k \exp(2 \pi i j k /N)
3261 The symmetry of the half-complex sequence implies that only half of the
3262 complex numbers in the output need to be stored. The remaining half can
3263 be reconstructed using the half-complex symmetry condition. This works
3264 for all lengths, even and odd--when the length is even the middle value
3265 where k=N/2 is also real. Thus only N real numbers are required to
3266 store the half-complex sequence, and the transform of a real sequence
3267 can be stored in the same size array as the original data.
3269 The precise storage arrangements depend on the algorithm, and are
3270 different for radix-2 and mixed-radix routines. The radix-2 function
3271 operates in-place, which constrains the locations where each element can
3272 be stored. The restriction forces real and imaginary parts to be stored
3273 far apart. The mixed-radix algorithm does not have this restriction,
3274 and it stores the real and imaginary parts of a given term in
3275 neighboring locations (which is desirable for better locality of memory
3279 File: gsl-ref.info, Node: Radix-2 FFT routines for real data, Next: Mixed-radix FFT routines for real data, Prev: Overview of real data FFTs, Up: Fast Fourier Transforms
3281 15.6 Radix-2 FFT routines for real data
3282 =======================================
3284 This section describes radix-2 FFT algorithms for real data. They use
3285 the Cooley-Tukey algorithm to compute in-place FFTs for lengths which
3288 The radix-2 FFT functions for real data are declared in the header
3289 files `gsl_fft_real.h'
3291 -- Function: int gsl_fft_real_radix2_transform (double DATA[], size_t
3293 This function computes an in-place radix-2 FFT of length N and
3294 stride STRIDE on the real array DATA. The output is a
3295 half-complex sequence, which is stored in-place. The arrangement
3296 of the half-complex terms uses the following scheme: for k < N/2
3297 the real part of the k-th term is stored in location k, and the
3298 corresponding imaginary part is stored in location N-k. Terms
3299 with k > N/2 can be reconstructed using the symmetry z_k =
3300 z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and
3301 count as a special case. Their real parts are stored in locations
3302 0 and N/2 respectively, while their imaginary parts which are zero
3305 The following table shows the correspondence between the output
3306 DATA and the equivalent results obtained by considering the input
3307 data as a complex sequence with zero imaginary part,
3309 complex[0].real = data[0]
3311 complex[1].real = data[1]
3312 complex[1].imag = data[N-1]
3313 ............... ................
3314 complex[k].real = data[k]
3315 complex[k].imag = data[N-k]
3316 ............... ................
3317 complex[N/2].real = data[N/2]
3318 complex[N/2].imag = 0
3319 ............... ................
3320 complex[k'].real = data[k] k' = N - k
3321 complex[k'].imag = -data[N-k]
3322 ............... ................
3323 complex[N-1].real = data[1]
3324 complex[N-1].imag = -data[N-1]
3325 Note that the output data can be converted into the full complex
3326 sequence using the function `gsl_fft_halfcomplex_unpack' described
3327 in the next section.
3329 The radix-2 FFT functions for halfcomplex data are declared in the
3330 header file `gsl_fft_halfcomplex.h'.
3332 -- Function: int gsl_fft_halfcomplex_radix2_inverse (double DATA[],
3333 size_t STRIDE, size_t N)
3334 -- Function: int gsl_fft_halfcomplex_radix2_backward (double DATA[],
3335 size_t STRIDE, size_t N)
3336 These functions compute the inverse or backwards in-place radix-2
3337 FFT of length N and stride STRIDE on the half-complex sequence
3338 DATA stored according the output scheme used by
3339 `gsl_fft_real_radix2'. The result is a real array stored in
3344 File: gsl-ref.info, Node: Mixed-radix FFT routines for real data, Next: FFT References and Further Reading, Prev: Radix-2 FFT routines for real data, Up: Fast Fourier Transforms
3346 15.7 Mixed-radix FFT routines for real data
3347 ===========================================
3349 This section describes mixed-radix FFT algorithms for real data. The
3350 mixed-radix functions work for FFTs of any length. They are a
3351 reimplementation of the real-FFT routines in the Fortran FFTPACK library
3352 by Paul Swarztrauber. The theory behind the algorithm is explained in
3353 the article `Fast Mixed-Radix Real Fourier Transforms' by Clive
3354 Temperton. The routines here use the same indexing scheme and basic
3355 algorithms as FFTPACK.
3357 The functions use the FFTPACK storage convention for half-complex
3358 sequences. In this convention the half-complex transform of a real
3359 sequence is stored with frequencies in increasing order, starting at
3360 zero, with the real and imaginary parts of each frequency in neighboring
3361 locations. When a value is known to be real the imaginary part is not
3362 stored. The imaginary part of the zero-frequency component is never
3363 stored. It is known to be zero (since the zero frequency component is
3364 simply the sum of the input data (all real)). For a sequence of even
3365 length the imaginary part of the frequency n/2 is not stored either,
3366 since the symmetry z_k = z_{N-k}^* implies that this is purely real too.
3368 The storage scheme is best shown by some examples. The table below
3369 shows the output for an odd-length sequence, n=5. The two columns give
3370 the correspondence between the 5 values in the half-complex sequence
3371 returned by `gsl_fft_real_transform', HALFCOMPLEX[] and the values
3372 COMPLEX[] that would be returned if the same real input sequence were
3373 passed to `gsl_fft_complex_backward' as a complex sequence (with
3374 imaginary parts set to `0'),
3376 complex[0].real = halfcomplex[0]
3378 complex[1].real = halfcomplex[1]
3379 complex[1].imag = halfcomplex[2]
3380 complex[2].real = halfcomplex[3]
3381 complex[2].imag = halfcomplex[4]
3382 complex[3].real = halfcomplex[3]
3383 complex[3].imag = -halfcomplex[4]
3384 complex[4].real = halfcomplex[1]
3385 complex[4].imag = -halfcomplex[2]
3387 The upper elements of the COMPLEX array, `complex[3]' and `complex[4]'
3388 are filled in using the symmetry condition. The imaginary part of the
3389 zero-frequency term `complex[0].imag' is known to be zero by the
3392 The next table shows the output for an even-length sequence, n=6 In
3393 the even case there are two values which are purely real,
3395 complex[0].real = halfcomplex[0]
3397 complex[1].real = halfcomplex[1]
3398 complex[1].imag = halfcomplex[2]
3399 complex[2].real = halfcomplex[3]
3400 complex[2].imag = halfcomplex[4]
3401 complex[3].real = halfcomplex[5]
3403 complex[4].real = halfcomplex[3]
3404 complex[4].imag = -halfcomplex[4]
3405 complex[5].real = halfcomplex[1]
3406 complex[5].imag = -halfcomplex[2]
3408 The upper elements of the COMPLEX array, `complex[4]' and `complex[5]'
3409 are filled in using the symmetry condition. Both `complex[0].imag' and
3410 `complex[3].imag' are known to be zero.
3412 All these functions are declared in the header files
3413 `gsl_fft_real.h' and `gsl_fft_halfcomplex.h'.
3415 -- Function: gsl_fft_real_wavetable * gsl_fft_real_wavetable_alloc
3417 -- Function: gsl_fft_halfcomplex_wavetable *
3418 gsl_fft_halfcomplex_wavetable_alloc (size_t N)
3419 These functions prepare trigonometric lookup tables for an FFT of
3420 size n real elements. The functions return a pointer to the newly
3421 allocated struct if no errors were detected, and a null pointer in
3422 the case of error. The length N is factorized into a product of
3423 subtransforms, and the factors and their trigonometric
3424 coefficients are stored in the wavetable. The trigonometric
3425 coefficients are computed using direct calls to `sin' and `cos',
3426 for accuracy. Recursion relations could be used to compute the
3427 lookup table faster, but if an application performs many FFTs of
3428 the same length then computing the wavetable is a one-off overhead
3429 which does not affect the final throughput.
3431 The wavetable structure can be used repeatedly for any transform
3432 of the same length. The table is not modified by calls to any of
3433 the other FFT functions. The appropriate type of wavetable must
3434 be used for forward real or inverse half-complex transforms.
3436 -- Function: void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable
3438 -- Function: void gsl_fft_halfcomplex_wavetable_free
3439 (gsl_fft_halfcomplex_wavetable * WAVETABLE)
3440 These functions free the memory associated with the wavetable
3441 WAVETABLE. The wavetable can be freed if no further FFTs of the
3442 same length will be needed.
3444 The mixed radix algorithms require additional working space to hold the
3445 intermediate steps of the transform,
3447 -- Function: gsl_fft_real_workspace * gsl_fft_real_workspace_alloc
3449 This function allocates a workspace for a real transform of length
3450 N. The same workspace can be used for both forward real and
3451 inverse halfcomplex transforms.
3453 -- Function: void gsl_fft_real_workspace_free (gsl_fft_real_workspace
3455 This function frees the memory associated with the workspace
3456 WORKSPACE. The workspace can be freed if no further FFTs of the
3457 same length will be needed.
3459 The following functions compute the transforms of real and half-complex
3462 -- Function: int gsl_fft_real_transform (double DATA[], size_t STRIDE,
3463 size_t N, const gsl_fft_real_wavetable * WAVETABLE,
3464 gsl_fft_real_workspace * WORK)
3465 -- Function: int gsl_fft_halfcomplex_transform (double DATA[], size_t
3466 STRIDE, size_t N, const gsl_fft_halfcomplex_wavetable *
3467 WAVETABLE, gsl_fft_real_workspace * WORK)
3468 These functions compute the FFT of DATA, a real or half-complex
3469 array of length N, using a mixed radix decimation-in-frequency
3470 algorithm. For `gsl_fft_real_transform' DATA is an array of
3471 time-ordered real data. For `gsl_fft_halfcomplex_transform' DATA
3472 contains fourier coefficients in the half-complex ordering
3473 described above. There is no restriction on the length N.
3474 Efficient modules are provided for subtransforms of length 2, 3, 4
3475 and 5. Any remaining factors are computed with a slow, O(n^2),
3476 general-n module. The caller must supply a WAVETABLE containing
3477 trigonometric lookup tables and a workspace WORK.
3479 -- Function: int gsl_fft_real_unpack (const double REAL_COEFFICIENT[],
3480 gsl_complex_packed_array COMPLEX_COEFFICIENT, size_t STRIDE,
3482 This function converts a single real array, REAL_COEFFICIENT into
3483 an equivalent complex array, COMPLEX_COEFFICIENT, (with imaginary
3484 part set to zero), suitable for `gsl_fft_complex' routines. The
3485 algorithm for the conversion is simply,
3487 for (i = 0; i < n; i++)
3489 complex_coefficient[i].real
3490 = real_coefficient[i];
3491 complex_coefficient[i].imag
3495 -- Function: int gsl_fft_halfcomplex_unpack (const double
3496 HALFCOMPLEX_COEFFICIENT[], gsl_complex_packed_array
3497 COMPLEX_COEFFICIENT, size_t STRIDE, size_t N)
3498 This function converts HALFCOMPLEX_COEFFICIENT, an array of
3499 half-complex coefficients as returned by `gsl_fft_real_transform',
3500 into an ordinary complex array, COMPLEX_COEFFICIENT. It fills in
3501 the complex array using the symmetry z_k = z_{N-k}^* to
3502 reconstruct the redundant elements. The algorithm for the
3505 complex_coefficient[0].real
3506 = halfcomplex_coefficient[0];
3507 complex_coefficient[0].imag
3510 for (i = 1; i < n - i; i++)
3513 = halfcomplex_coefficient[2 * i - 1];
3515 = halfcomplex_coefficient[2 * i];
3516 complex_coefficient[i].real = hc_real;
3517 complex_coefficient[i].imag = hc_imag;
3518 complex_coefficient[n - i].real = hc_real;
3519 complex_coefficient[n - i].imag = -hc_imag;
3524 complex_coefficient[i].real
3525 = halfcomplex_coefficient[n - 1];
3526 complex_coefficient[i].imag
3530 Here is an example program using `gsl_fft_real_transform' and
3531 `gsl_fft_halfcomplex_inverse'. It generates a real signal in the shape
3532 of a square pulse. The pulse is fourier transformed to frequency
3533 space, and all but the lowest ten frequency components are removed from
3534 the array of fourier coefficients returned by `gsl_fft_real_transform'.
3536 The remaining fourier coefficients are transformed back to the
3537 time-domain, to give a filtered version of the square pulse. Since
3538 fourier coefficients are stored using the half-complex symmetry both
3539 positive and negative frequencies are removed and the final filtered
3540 signal is also real.
3544 #include <gsl/gsl_errno.h>
3545 #include <gsl/gsl_fft_real.h>
3546 #include <gsl/gsl_fft_halfcomplex.h>
3554 gsl_fft_real_wavetable * real;
3555 gsl_fft_halfcomplex_wavetable * hc;
3556 gsl_fft_real_workspace * work;
3558 for (i = 0; i < n; i++)
3563 for (i = n / 3; i < 2 * n / 3; i++)
3568 for (i = 0; i < n; i++)
3570 printf ("%d: %e\n", i, data[i]);
3574 work = gsl_fft_real_workspace_alloc (n);
3575 real = gsl_fft_real_wavetable_alloc (n);
3577 gsl_fft_real_transform (data, 1, n,
3580 gsl_fft_real_wavetable_free (real);
3582 for (i = 11; i < n; i++)
3587 hc = gsl_fft_halfcomplex_wavetable_alloc (n);
3589 gsl_fft_halfcomplex_inverse (data, 1, n,
3591 gsl_fft_halfcomplex_wavetable_free (hc);
3593 for (i = 0; i < n; i++)
3595 printf ("%d: %e\n", i, data[i]);
3598 gsl_fft_real_workspace_free (work);
3603 File: gsl-ref.info, Node: FFT References and Further Reading, Prev: Mixed-radix FFT routines for real data, Up: Fast Fourier Transforms
3605 15.8 References and Further Reading
3606 ===================================
3608 A good starting point for learning more about the FFT is the review
3609 article `Fast Fourier Transforms: A Tutorial Review and A State of the
3610 Art' by Duhamel and Vetterli,
3612 P. Duhamel and M. Vetterli. Fast fourier transforms: A tutorial
3613 review and a state of the art. `Signal Processing', 19:259-299,
3616 To find out about the algorithms used in the GSL routines you may want
3617 to consult the document `GSL FFT Algorithms' (it is included in GSL, as
3618 `doc/fftalgorithms.tex'). This has general information on FFTs and
3619 explicit derivations of the implementation for each routine. There are
3620 also references to the relevant literature. For convenience some of
3621 the more important references are reproduced below.
3623 There are several introductory books on the FFT with example programs,
3624 such as `The Fast Fourier Transform' by Brigham and `DFT/FFT and
3625 Convolution Algorithms' by Burrus and Parks,
3627 E. Oran Brigham. `The Fast Fourier Transform'. Prentice Hall,
3630 C. S. Burrus and T. W. Parks. `DFT/FFT and Convolution
3631 Algorithms'. Wiley, 1984.
3633 Both these introductory books cover the radix-2 FFT in some detail.
3634 The mixed-radix algorithm at the heart of the FFTPACK routines is
3635 reviewed in Clive Temperton's paper,
3637 Clive Temperton. Self-sorting mixed-radix fast fourier transforms.
3638 `Journal of Computational Physics', 52(1):1-23, 1983.
3640 The derivation of FFTs for real-valued data is explained in the
3641 following two articles,
3643 Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C.
3644 Sidney Burrus. Real-valued fast fourier transform algorithms.
3645 `IEEE Transactions on Acoustics, Speech, and Signal Processing',
3646 ASSP-35(6):849-863, 1987.
3648 Clive Temperton. Fast mixed-radix real fourier transforms.
3649 `Journal of Computational Physics', 52:340-350, 1983.
3651 In 1979 the IEEE published a compendium of carefully-reviewed Fortran
3652 FFT programs in `Programs for Digital Signal Processing'. It is a
3653 useful reference for implementations of many different FFT algorithms,
3655 Digital Signal Processing Committee and IEEE Acoustics, Speech,
3656 and Signal Processing Committee, editors. `Programs for Digital
3657 Signal Processing'. IEEE Press, 1979.
3659 For large-scale FFT work we recommend the use of the dedicated FFTW
3660 library by Frigo and Johnson. The FFTW library is self-optimizing--it
3661 automatically tunes itself for each hardware platform in order to
3662 achieve maximum performance. It is available under the GNU GPL.
3664 FFTW Website, `http://www.fftw.org/'
3666 The source code for FFTPACK is available from Netlib,
3668 FFTPACK, `http://www.netlib.org/fftpack/'
3671 File: gsl-ref.info, Node: Numerical Integration, Next: Random Number Generation, Prev: Fast Fourier Transforms, Up: Top
3673 16 Numerical Integration
3674 ************************
3676 This chapter describes routines for performing numerical integration
3677 (quadrature) of a function in one dimension. There are routines for
3678 adaptive and non-adaptive integration of general functions, with
3679 specialised routines for specific cases. These include integration over
3680 infinite and semi-infinite ranges, singular integrals, including
3681 logarithmic singularities, computation of Cauchy principal values and
3682 oscillatory integrals. The library reimplements the algorithms used in
3683 QUADPACK, a numerical integration package written by Piessens,
3684 Doncker-Kapenga, Uberhuber and Kahaner. Fortran code for QUADPACK is
3685 available on Netlib.
3687 The functions described in this chapter are declared in the header
3688 file `gsl_integration.h'.
3692 * Numerical Integration Introduction::
3693 * QNG non-adaptive Gauss-Kronrod integration::
3694 * QAG adaptive integration::
3695 * QAGS adaptive integration with singularities::
3696 * QAGP adaptive integration with known singular points::
3697 * QAGI adaptive integration on infinite intervals::
3698 * QAWC adaptive integration for Cauchy principal values::
3699 * QAWS adaptive integration for singular functions::
3700 * QAWO adaptive integration for oscillatory functions::
3701 * QAWF adaptive integration for Fourier integrals::
3702 * Numerical integration error codes::
3703 * Numerical integration examples::
3704 * Numerical integration References and Further Reading::
3707 File: gsl-ref.info, Node: Numerical Integration Introduction, Next: QNG non-adaptive Gauss-Kronrod integration, Up: Numerical Integration
3712 Each algorithm computes an approximation to a definite integral of the
3715 I = \int_a^b f(x) w(x) dx
3717 where w(x) is a weight function (for general integrands w(x)=1). The
3718 user provides absolute and relative error bounds (epsabs, epsrel) which
3719 specify the following accuracy requirement,
3721 |RESULT - I| <= max(epsabs, epsrel |I|)
3723 where RESULT is the numerical approximation obtained by the algorithm.
3724 The algorithms attempt to estimate the absolute error ABSERR = |RESULT
3725 - I| in such a way that the following inequality holds,
3727 |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|)
3729 In short, the routines return the first approximation which has an
3730 absolute error smaller than epsabs or a relative error smaller than
3733 Note that this is an either-or constraint, not simultaneous. To
3734 compute to a specified absolute error, set epsrel to zero. To compute
3735 to a specified relative error, set epsabs to zero. The routines will
3736 fail to converge if the error bounds are too stringent, but always
3737 return the best approximation obtained up to that stage.
3739 The algorithms in QUADPACK use a naming convention based on the
3742 `Q' - quadrature routine
3744 `N' - non-adaptive integrator
3745 `A' - adaptive integrator
3747 `G' - general integrand (user-defined)
3748 `W' - weight function with integrand
3750 `S' - singularities can be more readily integrated
3751 `P' - points of special difficulty can be supplied
3752 `I' - infinite range of integration
3753 `O' - oscillatory weight function, cos or sin
3754 `F' - Fourier integral
3755 `C' - Cauchy principal value
3757 The algorithms are built on pairs of quadrature rules, a higher order
3758 rule and a lower order rule. The higher order rule is used to compute
3759 the best approximation to an integral over a small range. The
3760 difference between the results of the higher order rule and the lower
3761 order rule gives an estimate of the error in the approximation.
3765 * Integrands without weight functions::
3766 * Integrands with weight functions::
3767 * Integrands with singular weight functions::
3770 File: gsl-ref.info, Node: Integrands without weight functions, Next: Integrands with weight functions, Up: Numerical Integration Introduction
3772 16.1.1 Integrands without weight functions
3773 ------------------------------------------
3775 The algorithms for general functions (without a weight function) are
3776 based on Gauss-Kronrod rules.
3778 A Gauss-Kronrod rule begins with a classical Gaussian quadrature
3779 rule of order m. This is extended with additional points between each
3780 of the abscissae to give a higher order Kronrod rule of order 2m+1.
3781 The Kronrod rule is efficient because it reuses existing function
3782 evaluations from the Gaussian rule.
3784 The higher order Kronrod rule is used as the best approximation to
3785 the integral, and the difference between the two rules is used as an
3786 estimate of the error in the approximation.
3789 File: gsl-ref.info, Node: Integrands with weight functions, Next: Integrands with singular weight functions, Prev: Integrands without weight functions, Up: Numerical Integration Introduction
3791 16.1.2 Integrands with weight functions
3792 ---------------------------------------
3794 For integrands with weight functions the algorithms use Clenshaw-Curtis
3797 A Clenshaw-Curtis rule begins with an n-th order Chebyshev
3798 polynomial approximation to the integrand. This polynomial can be
3799 integrated exactly to give an approximation to the integral of the
3800 original function. The Chebyshev expansion can be extended to higher
3801 orders to improve the approximation and provide an estimate of the
3805 File: gsl-ref.info, Node: Integrands with singular weight functions, Prev: Integrands with weight functions, Up: Numerical Integration Introduction
3807 16.1.3 Integrands with singular weight functions
3808 ------------------------------------------------
3810 The presence of singularities (or other behavior) in the integrand can
3811 cause slow convergence in the Chebyshev approximation. The modified
3812 Clenshaw-Curtis rules used in QUADPACK separate out several common
3813 weight functions which cause slow convergence.
3815 These weight functions are integrated analytically against the
3816 Chebyshev polynomials to precompute "modified Chebyshev moments".
3817 Combining the moments with the Chebyshev approximation to the function
3818 gives the desired integral. The use of analytic integration for the
3819 singular part of the function allows exact cancellations and
3820 substantially improves the overall convergence behavior of the
3824 File: gsl-ref.info, Node: QNG non-adaptive Gauss-Kronrod integration, Next: QAG adaptive integration, Prev: Numerical Integration Introduction, Up: Numerical Integration
3826 16.2 QNG non-adaptive Gauss-Kronrod integration
3827 ===============================================
3829 The QNG algorithm is a non-adaptive procedure which uses fixed
3830 Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum
3831 of 87 points. It is provided for fast integration of smooth functions.
3833 -- Function: int gsl_integration_qng (const gsl_function * F, double
3834 A, double B, double EPSABS, double EPSREL, double * RESULT,
3835 double * ABSERR, size_t * NEVAL)
3836 This function applies the Gauss-Kronrod 10-point, 21-point,
3837 43-point and 87-point integration rules in succession until an
3838 estimate of the integral of f over (a,b) is achieved within the
3839 desired absolute and relative error limits, EPSABS and EPSREL. The
3840 function returns the final approximation, RESULT, an estimate of
3841 the absolute error, ABSERR and the number of function evaluations
3842 used, NEVAL. The Gauss-Kronrod rules are designed in such a way
3843 that each rule uses all the results of its predecessors, in order
3844 to minimize the total number of function evaluations.
3847 File: gsl-ref.info, Node: QAG adaptive integration, Next: QAGS adaptive integration with singularities, Prev: QNG non-adaptive Gauss-Kronrod integration, Up: Numerical Integration
3849 16.3 QAG adaptive integration
3850 =============================
3852 The QAG algorithm is a simple adaptive integration procedure. The
3853 integration region is divided into subintervals, and on each iteration
3854 the subinterval with the largest estimated error is bisected. This
3855 reduces the overall error rapidly, as the subintervals become
3856 concentrated around local difficulties in the integrand. These
3857 subintervals are managed by a `gsl_integration_workspace' struct, which
3858 handles the memory for the subinterval ranges, results and error
3861 -- Function: gsl_integration_workspace *
3862 gsl_integration_workspace_alloc (size_t N)
3863 This function allocates a workspace sufficient to hold N double
3864 precision intervals, their integration results and error estimates.
3866 -- Function: void gsl_integration_workspace_free
3867 (gsl_integration_workspace * W)
3868 This function frees the memory associated with the workspace W.
3870 -- Function: int gsl_integration_qag (const gsl_function * F, double
3871 A, double B, double EPSABS, double EPSREL, size_t LIMIT, int
3872 KEY, gsl_integration_workspace * WORKSPACE, double * RESULT,
3874 This function applies an integration rule adaptively until an
3875 estimate of the integral of f over (a,b) is achieved within the
3876 desired absolute and relative error limits, EPSABS and EPSREL.
3877 The function returns the final approximation, RESULT, and an
3878 estimate of the absolute error, ABSERR. The integration rule is
3879 determined by the value of KEY, which should be chosen from the
3880 following symbolic names,
3882 GSL_INTEG_GAUSS15 (key = 1)
3883 GSL_INTEG_GAUSS21 (key = 2)
3884 GSL_INTEG_GAUSS31 (key = 3)
3885 GSL_INTEG_GAUSS41 (key = 4)
3886 GSL_INTEG_GAUSS51 (key = 5)
3887 GSL_INTEG_GAUSS61 (key = 6)
3889 corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod
3890 rules. The higher-order rules give better accuracy for smooth
3891 functions, while lower-order rules save time when the function
3892 contains local difficulties, such as discontinuities.
3894 On each iteration the adaptive integration strategy bisects the
3895 interval with the largest error estimate. The subintervals and
3896 their results are stored in the memory provided by WORKSPACE. The
3897 maximum number of subintervals is given by LIMIT, which may not
3898 exceed the allocated size of the workspace.
3901 File: gsl-ref.info, Node: QAGS adaptive integration with singularities, Next: QAGP adaptive integration with known singular points, Prev: QAG adaptive integration, Up: Numerical Integration
3903 16.4 QAGS adaptive integration with singularities
3904 =================================================
3906 The presence of an integrable singularity in the integration region
3907 causes an adaptive routine to concentrate new subintervals around the
3908 singularity. As the subintervals decrease in size the successive
3909 approximations to the integral converge in a limiting fashion. This
3910 approach to the limit can be accelerated using an extrapolation
3911 procedure. The QAGS algorithm combines adaptive bisection with the Wynn
3912 epsilon-algorithm to speed up the integration of many types of
3913 integrable singularities.
3915 -- Function: int gsl_integration_qags (const gsl_function * F, double
3916 A, double B, double EPSABS, double EPSREL, size_t LIMIT,
3917 gsl_integration_workspace * WORKSPACE, double * RESULT,
3919 This function applies the Gauss-Kronrod 21-point integration rule
3920 adaptively until an estimate of the integral of f over (a,b) is
3921 achieved within the desired absolute and relative error limits,
3922 EPSABS and EPSREL. The results are extrapolated using the
3923 epsilon-algorithm, which accelerates the convergence of the
3924 integral in the presence of discontinuities and integrable
3925 singularities. The function returns the final approximation from
3926 the extrapolation, RESULT, and an estimate of the absolute error,
3927 ABSERR. The subintervals and their results are stored in the
3928 memory provided by WORKSPACE. The maximum number of subintervals
3929 is given by LIMIT, which may not exceed the allocated size of the
3934 File: gsl-ref.info, Node: QAGP adaptive integration with known singular points, Next: QAGI adaptive integration on infinite intervals, Prev: QAGS adaptive integration with singularities, Up: Numerical Integration
3936 16.5 QAGP adaptive integration with known singular points
3937 =========================================================
3939 -- Function: int gsl_integration_qagp (const gsl_function * F, double
3940 * PTS, size_t NPTS, double EPSABS, double EPSREL, size_t
3941 LIMIT, gsl_integration_workspace * WORKSPACE, double *
3942 RESULT, double * ABSERR)
3943 This function applies the adaptive integration algorithm QAGS
3944 taking account of the user-supplied locations of singular points.
3945 The array PTS of length NPTS should contain the endpoints of the
3946 integration ranges defined by the integration region and locations
3947 of the singularities. For example, to integrate over the region
3948 (a,b) with break-points at x_1, x_2, x_3 (where a < x_1 < x_2 <
3949 x_3 < b) the following PTS array should be used
3959 If you know the locations of the singular points in the integration
3960 region then this routine will be faster than `QAGS'.
3964 File: gsl-ref.info, Node: QAGI adaptive integration on infinite intervals, Next: QAWC adaptive integration for Cauchy principal values, Prev: QAGP adaptive integration with known singular points, Up: Numerical Integration
3966 16.6 QAGI adaptive integration on infinite intervals
3967 ====================================================
3969 -- Function: int gsl_integration_qagi (gsl_function * F, double
3970 EPSABS, double EPSREL, size_t LIMIT,
3971 gsl_integration_workspace * WORKSPACE, double * RESULT,
3973 This function computes the integral of the function F over the
3974 infinite interval (-\infty,+\infty). The integral is mapped onto
3975 the semi-open interval (0,1] using the transformation x = (1-t)/t,
3977 \int_{-\infty}^{+\infty} dx f(x) =
3978 \int_0^1 dt (f((1-t)/t) + f((-1+t)/t))/t^2.
3980 It is then integrated using the QAGS algorithm. The normal
3981 21-point Gauss-Kronrod rule of QAGS is replaced by a 15-point
3982 rule, because the transformation can generate an integrable
3983 singularity at the origin. In this case a lower-order rule is
3986 -- Function: int gsl_integration_qagiu (gsl_function * F, double A,
3987 double EPSABS, double EPSREL, size_t LIMIT,
3988 gsl_integration_workspace * WORKSPACE, double * RESULT,
3990 This function computes the integral of the function F over the
3991 semi-infinite interval (a,+\infty). The integral is mapped onto
3992 the semi-open interval (0,1] using the transformation x = a +
3995 \int_{a}^{+\infty} dx f(x) =
3996 \int_0^1 dt f(a + (1-t)/t)/t^2
3998 and then integrated using the QAGS algorithm.
4000 -- Function: int gsl_integration_qagil (gsl_function * F, double B,
4001 double EPSABS, double EPSREL, size_t LIMIT,
4002 gsl_integration_workspace * WORKSPACE, double * RESULT,
4004 This function computes the integral of the function F over the
4005 semi-infinite interval (-\infty,b). The integral is mapped onto
4006 the semi-open interval (0,1] using the transformation x = b -
4009 \int_{-\infty}^{b} dx f(x) =
4010 \int_0^1 dt f(b - (1-t)/t)/t^2
4012 and then integrated using the QAGS algorithm.
4015 File: gsl-ref.info, Node: QAWC adaptive integration for Cauchy principal values, Next: QAWS adaptive integration for singular functions, Prev: QAGI adaptive integration on infinite intervals, Up: Numerical Integration
4017 16.7 QAWC adaptive integration for Cauchy principal values
4018 ==========================================================
4020 -- Function: int gsl_integration_qawc (gsl_function * F, double A,
4021 double B, double C, double EPSABS, double EPSREL, size_t
4022 LIMIT, gsl_integration_workspace * WORKSPACE, double *
4023 RESULT, double * ABSERR)
4024 This function computes the Cauchy principal value of the integral
4025 of f over (a,b), with a singularity at C,
4027 I = \int_a^b dx f(x) / (x - c)
4029 The adaptive bisection algorithm of QAG is used, with
4030 modifications to ensure that subdivisions do not occur at the
4031 singular point x = c. When a subinterval contains the point x = c
4032 or is close to it then a special 25-point modified Clenshaw-Curtis
4033 rule is used to control the singularity. Further away from the
4034 singularity the algorithm uses an ordinary 15-point Gauss-Kronrod
4039 File: gsl-ref.info, Node: QAWS adaptive integration for singular functions, Next: QAWO adaptive integration for oscillatory functions, Prev: QAWC adaptive integration for Cauchy principal values, Up: Numerical Integration
4041 16.8 QAWS adaptive integration for singular functions
4042 =====================================================
4044 The QAWS algorithm is designed for integrands with algebraic-logarithmic
4045 singularities at the end-points of an integration region. In order to
4046 work efficiently the algorithm requires a precomputed table of
4049 -- Function: gsl_integration_qaws_table *
4050 gsl_integration_qaws_table_alloc (double ALPHA, double BETA, int MU,
4052 This function allocates space for a `gsl_integration_qaws_table'
4053 struct describing a singular weight function W(x) with the
4054 parameters (\alpha, \beta, \mu, \nu),
4056 W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x)
4058 where \alpha > -1, \beta > -1, and \mu = 0, 1, \nu = 0, 1. The
4059 weight function can take four different forms depending on the
4060 values of \mu and \nu,
4062 W(x) = (x-a)^alpha (b-x)^beta (mu = 0, nu = 0)
4063 W(x) = (x-a)^alpha (b-x)^beta log(x-a) (mu = 1, nu = 0)
4064 W(x) = (x-a)^alpha (b-x)^beta log(b-x) (mu = 0, nu = 1)
4065 W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1)
4067 The singular points (a,b) do not have to be specified until the
4068 integral is computed, where they are the endpoints of the
4071 The function returns a pointer to the newly allocated table
4072 `gsl_integration_qaws_table' if no errors were detected, and 0 in
4075 -- Function: int gsl_integration_qaws_table_set
4076 (gsl_integration_qaws_table * T, double ALPHA, double BETA,
4078 This function modifies the parameters (\alpha, \beta, \mu, \nu) of
4079 an existing `gsl_integration_qaws_table' struct T.
4081 -- Function: void gsl_integration_qaws_table_free
4082 (gsl_integration_qaws_table * T)
4083 This function frees all the memory associated with the
4084 `gsl_integration_qaws_table' struct T.
4086 -- Function: int gsl_integration_qaws (gsl_function * F, const double
4087 A, const double B, gsl_integration_qaws_table * T, const
4088 double EPSABS, const double EPSREL, const size_t LIMIT,
4089 gsl_integration_workspace * WORKSPACE, double * RESULT,
4091 This function computes the integral of the function f(x) over the
4092 interval (a,b) with the singular weight function (x-a)^\alpha
4093 (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x). The parameters of the
4094 weight function (\alpha, \beta, \mu, \nu) are taken from the table
4097 I = \int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x).
4099 The adaptive bisection algorithm of QAG is used. When a
4100 subinterval contains one of the endpoints then a special 25-point
4101 modified Clenshaw-Curtis rule is used to control the
4102 singularities. For subintervals which do not include the
4103 endpoints an ordinary 15-point Gauss-Kronrod integration rule is
4108 File: gsl-ref.info, Node: QAWO adaptive integration for oscillatory functions, Next: QAWF adaptive integration for Fourier integrals, Prev: QAWS adaptive integration for singular functions, Up: Numerical Integration
4110 16.9 QAWO adaptive integration for oscillatory functions
4111 ========================================================
4113 The QAWO algorithm is designed for integrands with an oscillatory
4114 factor, \sin(\omega x) or \cos(\omega x). In order to work efficiently
4115 the algorithm requires a table of Chebyshev moments which must be
4116 pre-computed with calls to the functions below.
4118 -- Function: gsl_integration_qawo_table *
4119 gsl_integration_qawo_table_alloc (double OMEGA, double L, enum
4120 gsl_integration_qawo_enum SINE, size_t N)
4121 This function allocates space for a `gsl_integration_qawo_table'
4122 struct and its associated workspace describing a sine or cosine
4123 weight function W(x) with the parameters (\omega, L),
4128 The parameter L must be the length of the interval over which the
4129 function will be integrated L = b - a. The choice of sine or
4130 cosine is made with the parameter SINE which should be chosen from
4131 one of the two following symbolic values:
4136 The `gsl_integration_qawo_table' is a table of the trigonometric
4137 coefficients required in the integration process. The parameter N
4138 determines the number of levels of coefficients that are computed.
4139 Each level corresponds to one bisection of the interval L, so that
4140 N levels are sufficient for subintervals down to the length L/2^n.
4141 The integration routine `gsl_integration_qawo' returns the error
4142 `GSL_ETABLE' if the number of levels is insufficient for the
4146 -- Function: int gsl_integration_qawo_table_set
4147 (gsl_integration_qawo_table * T, double OMEGA, double L, enum
4148 gsl_integration_qawo_enum SINE)
4149 This function changes the parameters OMEGA, L and SINE of the
4150 existing workspace T.
4152 -- Function: int gsl_integration_qawo_table_set_length
4153 (gsl_integration_qawo_table * T, double L)
4154 This function allows the length parameter L of the workspace T to
4157 -- Function: void gsl_integration_qawo_table_free
4158 (gsl_integration_qawo_table * T)
4159 This function frees all the memory associated with the workspace T.
4161 -- Function: int gsl_integration_qawo (gsl_function * F, const double
4162 A, const double EPSABS, const double EPSREL, const size_t
4163 LIMIT, gsl_integration_workspace * WORKSPACE,
4164 gsl_integration_qawo_table * WF, double * RESULT, double *
4166 This function uses an adaptive algorithm to compute the integral of
4167 f over (a,b) with the weight function \sin(\omega x) or
4168 \cos(\omega x) defined by the table WF,
4170 I = \int_a^b dx f(x) sin(omega x)
4171 I = \int_a^b dx f(x) cos(omega x)
4173 The results are extrapolated using the epsilon-algorithm to
4174 accelerate the convergence of the integral. The function returns
4175 the final approximation from the extrapolation, RESULT, and an
4176 estimate of the absolute error, ABSERR. The subintervals and
4177 their results are stored in the memory provided by WORKSPACE. The
4178 maximum number of subintervals is given by LIMIT, which may not
4179 exceed the allocated size of the workspace.
4181 Those subintervals with "large" widths d where d\omega > 4 are
4182 computed using a 25-point Clenshaw-Curtis integration rule, which
4183 handles the oscillatory behavior. Subintervals with a "small"
4184 widths where d\omega < 4 are computed using a 15-point
4185 Gauss-Kronrod integration.
4189 File: gsl-ref.info, Node: QAWF adaptive integration for Fourier integrals, Next: Numerical integration error codes, Prev: QAWO adaptive integration for oscillatory functions, Up: Numerical Integration
4191 16.10 QAWF adaptive integration for Fourier integrals
4192 =====================================================
4194 -- Function: int gsl_integration_qawf (gsl_function * F, const double
4195 A, const double EPSABS, const size_t LIMIT,
4196 gsl_integration_workspace * WORKSPACE,
4197 gsl_integration_workspace * CYCLE_WORKSPACE,
4198 gsl_integration_qawo_table * WF, double * RESULT, double *
4200 This function attempts to compute a Fourier integral of the
4201 function F over the semi-infinite interval [a,+\infty).
4203 I = \int_a^{+\infty} dx f(x) sin(omega x)
4204 I = \int_a^{+\infty} dx f(x) cos(omega x)
4206 The parameter \omega and choice of \sin or \cos is taken from the
4207 table WF (the length L can take any value, since it is overridden
4208 by this function to a value appropriate for the fourier
4209 integration). The integral is computed using the QAWO algorithm
4210 over each of the subintervals,
4213 C_2 = [a + c, a + 2 c]
4215 C_k = [a + (k-1) c, a + k c]
4217 where c = (2 floor(|\omega|) + 1) \pi/|\omega|. The width c is
4218 chosen to cover an odd number of periods so that the contributions
4219 from the intervals alternate in sign and are monotonically
4220 decreasing when F is positive and monotonically decreasing. The
4221 sum of this sequence of contributions is accelerated using the
4224 This function works to an overall absolute tolerance of ABSERR.
4225 The following strategy is used: on each interval C_k the algorithm
4226 tries to achieve the tolerance
4230 where u_k = (1 - p)p^{k-1} and p = 9/10. The sum of the geometric
4231 series of contributions from each interval gives an overall
4232 tolerance of ABSERR.
4234 If the integration of a subinterval leads to difficulties then the
4235 accuracy requirement for subsequent intervals is relaxed,
4237 TOL_k = u_k max(abserr, max_{i<k}{E_i})
4239 where E_k is the estimated error on the interval C_k.
4241 The subintervals and their results are stored in the memory
4242 provided by WORKSPACE. The maximum number of subintervals is
4243 given by LIMIT, which may not exceed the allocated size of the
4244 workspace. The integration over each subinterval uses the memory
4245 provided by CYCLE_WORKSPACE as workspace for the QAWO algorithm.
4249 File: gsl-ref.info, Node: Numerical integration error codes, Next: Numerical integration examples, Prev: QAWF adaptive integration for Fourier integrals, Up: Numerical Integration
4254 In addition to the standard error codes for invalid arguments the
4255 functions can return the following values,
4258 the maximum number of subdivisions was exceeded.
4261 cannot reach tolerance because of roundoff error, or roundoff
4262 error was detected in the extrapolation table.
4265 a non-integrable singularity or other bad integrand behavior was
4266 found in the integration interval.
4269 the integral is divergent, or too slowly convergent to be
4270 integrated numerically.
4273 File: gsl-ref.info, Node: Numerical integration examples, Next: Numerical integration References and Further Reading, Prev: Numerical integration error codes, Up: Numerical Integration
4278 The integrator `QAGS' will handle a large class of definite integrals.
4279 For example, consider the following integral, which has a
4280 algebraic-logarithmic singularity at the origin,
4282 \int_0^1 x^{-1/2} log(x) dx = -4
4284 The program below computes this integral to a relative accuracy bound of
4289 #include <gsl/gsl_integration.h>
4291 double f (double x, void * params) {
4292 double alpha = *(double *) params;
4293 double f = log(alpha*x) / sqrt(x);
4300 gsl_integration_workspace * w
4301 = gsl_integration_workspace_alloc (1000);
4303 double result, error;
4304 double expected = -4.0;
4311 gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
4312 w, &result, &error);
4314 printf ("result = % .18f\n", result);
4315 printf ("exact result = % .18f\n", expected);
4316 printf ("estimated error = % .18f\n", error);
4317 printf ("actual error = % .18f\n", result - expected);
4318 printf ("intervals = %d\n", w->size);
4320 gsl_integration_workspace_free (w);
4325 The results below show that the desired accuracy is achieved after 8
4329 result = -3.999999999999973799
4330 exact result = -4.000000000000000000
4331 estimated error = 0.000000000000246025
4332 actual error = 0.000000000000026201
4335 In fact, the extrapolation procedure used by `QAGS' produces an
4336 accuracy of almost twice as many digits. The error estimate returned by
4337 the extrapolation procedure is larger than the actual error, giving a
4338 margin of safety of one order of magnitude.
4341 File: gsl-ref.info, Node: Numerical integration References and Further Reading, Prev: Numerical integration examples, Up: Numerical Integration
4343 16.13 References and Further Reading
4344 ====================================
4346 The following book is the definitive reference for QUADPACK, and was
4347 written by the original authors. It provides descriptions of the
4348 algorithms, program listings, test programs and examples. It also
4349 includes useful advice on numerical integration and many references to
4350 the numerical integration literature used in developing QUADPACK.
4352 R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahaner.
4353 `QUADPACK A subroutine package for automatic integration' Springer
4358 File: gsl-ref.info, Node: Random Number Generation, Next: Quasi-Random Sequences, Prev: Numerical Integration, Up: Top
4360 17 Random Number Generation
4361 ***************************
4363 The library provides a large collection of random number generators
4364 which can be accessed through a uniform interface. Environment
4365 variables allow you to select different generators and seeds at runtime,
4366 so that you can easily switch between generators without needing to
4367 recompile your program. Each instance of a generator keeps track of its
4368 own state, allowing the generators to be used in multi-threaded
4369 programs. Additional functions are available for transforming uniform
4370 random numbers into samples from continuous or discrete probability
4371 distributions such as the Gaussian, log-normal or Poisson distributions.
4373 These functions are declared in the header file `gsl_rng.h'.
4377 * General comments on random numbers::
4378 * The Random Number Generator Interface::
4379 * Random number generator initialization::
4380 * Sampling from a random number generator::
4381 * Auxiliary random number generator functions::
4382 * Random number environment variables::
4383 * Copying random number generator state::
4384 * Reading and writing random number generator state::
4385 * Random number generator algorithms::
4386 * Unix random number generators::
4387 * Other random number generators::
4388 * Random Number Generator Performance::
4389 * Random Number Generator Examples::
4390 * Random Number References and Further Reading::
4391 * Random Number Acknowledgements::
4394 File: gsl-ref.info, Node: General comments on random numbers, Next: The Random Number Generator Interface, Up: Random Number Generation
4396 17.1 General comments on random numbers
4397 =======================================
4399 In 1988, Park and Miller wrote a paper entitled "Random number
4400 generators: good ones are hard to find." [Commun. ACM, 31, 1192-1201].
4401 Fortunately, some excellent random number generators are available,
4402 though poor ones are still in common use. You may be happy with the
4403 system-supplied random number generator on your computer, but you should
4404 be aware that as computers get faster, requirements on random number
4405 generators increase. Nowadays, a simulation that calls a random number
4406 generator millions of times can often finish before you can make it down
4407 the hall to the coffee machine and back.
4409 A very nice review of random number generators was written by Pierre
4410 L'Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks,
4411 ed. (Wiley, 1997). The chapter is available in postscript from
4412 L'Ecuyer's ftp site (see references). Knuth's volume on Seminumerical
4413 Algorithms (originally published in 1968) devotes 170 pages to random
4414 number generators, and has recently been updated in its 3rd edition
4415 (1997). It is brilliant, a classic. If you don't own it, you should
4416 stop reading right now, run to the nearest bookstore, and buy it.
4418 A good random number generator will satisfy both theoretical and
4419 statistical properties. Theoretical properties are often hard to obtain
4420 (they require real math!), but one prefers a random number generator
4421 with a long period, low serial correlation, and a tendency _not_ to
4422 "fall mainly on the planes." Statistical tests are performed with
4423 numerical simulations. Generally, a random number generator is used to
4424 estimate some quantity for which the theory of probability provides an
4425 exact answer. Comparison to this exact answer provides a measure of
4429 File: gsl-ref.info, Node: The Random Number Generator Interface, Next: Random number generator initialization, Prev: General comments on random numbers, Up: Random Number Generation
4431 17.2 The Random Number Generator Interface
4432 ==========================================
4434 It is important to remember that a random number generator is not a
4435 "real" function like sine or cosine. Unlike real functions, successive
4436 calls to a random number generator yield different return values. Of
4437 course that is just what you want for a random number generator, but to
4438 achieve this effect, the generator must keep track of some kind of
4439 "state" variable. Sometimes this state is just an integer (sometimes
4440 just the value of the previously generated random number), but often it
4441 is more complicated than that and may involve a whole array of numbers,
4442 possibly with some indices thrown in. To use the random number
4443 generators, you do not need to know the details of what comprises the
4444 state, and besides that varies from algorithm to algorithm.
4446 The random number generator library uses two special structs,
4447 `gsl_rng_type' which holds static information about each type of
4448 generator and `gsl_rng' which describes an instance of a generator
4449 created from a given `gsl_rng_type'.
4451 The functions described in this section are declared in the header
4455 File: gsl-ref.info, Node: Random number generator initialization, Next: Sampling from a random number generator, Prev: The Random Number Generator Interface, Up: Random Number Generation
4457 17.3 Random number generator initialization
4458 ===========================================
4460 -- Function: gsl_rng * gsl_rng_alloc (const gsl_rng_type * T)
4461 This function returns a pointer to a newly-created instance of a
4462 random number generator of type T. For example, the following
4463 code creates an instance of the Tausworthe generator,
4465 gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
4467 If there is insufficient memory to create the generator then the
4468 function returns a null pointer and the error handler is invoked
4469 with an error code of `GSL_ENOMEM'.
4471 The generator is automatically initialized with the default seed,
4472 `gsl_rng_default_seed'. This is zero by default but can be changed
4473 either directly or by using the environment variable `GSL_RNG_SEED'
4474 (*note Random number environment variables::).
4476 The details of the available generator types are described later
4479 -- Function: void gsl_rng_set (const gsl_rng * R, unsigned long int S)
4480 This function initializes (or `seeds') the random number
4481 generator. If the generator is seeded with the same value of S on
4482 two different runs, the same stream of random numbers will be
4483 generated by successive calls to the routines below. If different
4484 values of S are supplied, then the generated streams of random
4485 numbers should be completely different. If the seed S is zero
4486 then the standard seed from the original implementation is used
4487 instead. For example, the original Fortran source code for the
4488 `ranlux' generator used a seed of 314159265, and so choosing S
4489 equal to zero reproduces this when using `gsl_rng_ranlux'.
4491 -- Function: void gsl_rng_free (gsl_rng * R)
4492 This function frees all the memory associated with the generator R.
4495 File: gsl-ref.info, Node: Sampling from a random number generator, Next: Auxiliary random number generator functions, Prev: Random number generator initialization, Up: Random Number Generation
4497 17.4 Sampling from a random number generator
4498 ============================================
4500 The following functions return uniformly distributed random numbers,
4501 either as integers or double precision floating point numbers. Inline
4502 versions of these functions are used when `HAVE_INLINE' is defined. To
4503 obtain non-uniform distributions *note Random Number Distributions::.
4505 -- Function: unsigned long int gsl_rng_get (const gsl_rng * R)
4506 This function returns a random integer from the generator R. The
4507 minimum and maximum values depend on the algorithm used, but all
4508 integers in the range [MIN,MAX] are equally likely. The values of
4509 MIN and MAX can determined using the auxiliary functions
4510 `gsl_rng_max (r)' and `gsl_rng_min (r)'.
4512 -- Function: double gsl_rng_uniform (const gsl_rng * R)
4513 This function returns a double precision floating point number
4514 uniformly distributed in the range [0,1). The range includes 0.0
4515 but excludes 1.0. The value is typically obtained by dividing the
4516 result of `gsl_rng_get(r)' by `gsl_rng_max(r) + 1.0' in double
4517 precision. Some generators compute this ratio internally so that
4518 they can provide floating point numbers with more than 32 bits of
4519 randomness (the maximum number of bits that can be portably
4520 represented in a single `unsigned long int').
4522 -- Function: double gsl_rng_uniform_pos (const gsl_rng * R)
4523 This function returns a positive double precision floating point
4524 number uniformly distributed in the range (0,1), excluding both
4525 0.0 and 1.0. The number is obtained by sampling the generator
4526 with the algorithm of `gsl_rng_uniform' until a non-zero value is
4527 obtained. You can use this function if you need to avoid a
4530 -- Function: unsigned long int gsl_rng_uniform_int (const gsl_rng * R,
4531 unsigned long int N)
4532 This function returns a random integer from 0 to n-1 inclusive by
4533 scaling down and/or discarding samples from the generator R. All
4534 integers in the range [0,n-1] are produced with equal probability.
4535 For generators with a non-zero minimum value an offset is applied
4536 so that zero is returned with the correct probability.
4538 Note that this function is designed for sampling from ranges
4539 smaller than the range of the underlying generator. The parameter
4540 N must be less than or equal to the range of the generator R. If
4541 N is larger than the range of the generator then the function
4542 calls the error handler with an error code of `GSL_EINVAL' and
4545 In particular, this function is not intended for generating the
4546 full range of unsigned integer values [0,2^32-1]. Instead choose a
4547 generator with the maximal integer range and zero mimimum value,
4548 such as `gsl_rng_ranlxd1', `gsl_rng_mt19937' or `gsl_rng_taus',
4549 and sample it directly using `gsl_rng_get'. The range of each
4550 generator can be found using the auxiliary functions described in
4554 File: gsl-ref.info, Node: Auxiliary random number generator functions, Next: Random number environment variables, Prev: Sampling from a random number generator, Up: Random Number Generation
4556 17.5 Auxiliary random number generator functions
4557 ================================================
4559 The following functions provide information about an existing
4560 generator. You should use them in preference to hard-coding the
4561 generator parameters into your own code.
4563 -- Function: const char * gsl_rng_name (const gsl_rng * R)
4564 This function returns a pointer to the name of the generator. For
4567 printf ("r is a '%s' generator\n",
4570 would print something like `r is a 'taus' generator'.
4572 -- Function: unsigned long int gsl_rng_max (const gsl_rng * R)
4573 `gsl_rng_max' returns the largest value that `gsl_rng_get' can
4576 -- Function: unsigned long int gsl_rng_min (const gsl_rng * R)
4577 `gsl_rng_min' returns the smallest value that `gsl_rng_get' can
4578 return. Usually this value is zero. There are some generators
4579 with algorithms that cannot return zero, and for these generators
4580 the minimum value is 1.
4582 -- Function: void * gsl_rng_state (const gsl_rng * R)
4583 -- Function: size_t gsl_rng_size (const gsl_rng * R)
4584 These functions return a pointer to the state of generator R and
4585 its size. You can use this information to access the state
4586 directly. For example, the following code will write the state of
4587 a generator to a stream,
4589 void * state = gsl_rng_state (r);
4590 size_t n = gsl_rng_size (r);
4591 fwrite (state, n, 1, stream);
4593 -- Function: const gsl_rng_type ** gsl_rng_types_setup (void)
4594 This function returns a pointer to an array of all the available
4595 generator types, terminated by a null pointer. The function should
4596 be called once at the start of the program, if needed. The
4597 following code fragment shows how to iterate over the array of
4598 generator types to print the names of the available algorithms,
4600 const gsl_rng_type **t, **t0;
4602 t0 = gsl_rng_types_setup ();
4604 printf ("Available generators:\n");
4606 for (t = t0; *t != 0; t++)
4608 printf ("%s\n", (*t)->name);
4612 File: gsl-ref.info, Node: Random number environment variables, Next: Copying random number generator state, Prev: Auxiliary random number generator functions, Up: Random Number Generation
4614 17.6 Random number environment variables
4615 ========================================
4617 The library allows you to choose a default generator and seed from the
4618 environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED' and the
4619 function `gsl_rng_env_setup'. This makes it easy try out different
4620 generators and seeds without having to recompile your program.
4622 -- Function: const gsl_rng_type * gsl_rng_env_setup (void)
4623 This function reads the environment variables `GSL_RNG_TYPE' and
4624 `GSL_RNG_SEED' and uses their values to set the corresponding
4625 library variables `gsl_rng_default' and `gsl_rng_default_seed'.
4626 These global variables are defined as follows,
4628 extern const gsl_rng_type *gsl_rng_default
4629 extern unsigned long int gsl_rng_default_seed
4631 The environment variable `GSL_RNG_TYPE' should be the name of a
4632 generator, such as `taus' or `mt19937'. The environment variable
4633 `GSL_RNG_SEED' should contain the desired seed value. It is
4634 converted to an `unsigned long int' using the C library function
4637 If you don't specify a generator for `GSL_RNG_TYPE' then
4638 `gsl_rng_mt19937' is used as the default. The initial value of
4639 `gsl_rng_default_seed' is zero.
4642 Here is a short program which shows how to create a global generator
4643 using the environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED',
4646 #include <gsl/gsl_rng.h>
4648 gsl_rng * r; /* global generator */
4653 const gsl_rng_type * T;
4655 gsl_rng_env_setup();
4657 T = gsl_rng_default;
4658 r = gsl_rng_alloc (T);
4660 printf ("generator type: %s\n", gsl_rng_name (r));
4661 printf ("seed = %lu\n", gsl_rng_default_seed);
4662 printf ("first value = %lu\n", gsl_rng_get (r));
4668 Running the program without any environment variables uses the initial
4669 defaults, an `mt19937' generator with a seed of 0,
4672 generator type: mt19937
4674 first value = 4293858116
4676 By setting the two variables on the command line we can change the
4677 default generator and the seed,
4679 $ GSL_RNG_TYPE="taus" GSL_RNG_SEED=123 ./a.out
4682 generator type: taus
4684 first value = 2720986350
4687 File: gsl-ref.info, Node: Copying random number generator state, Next: Reading and writing random number generator state, Prev: Random number environment variables, Up: Random Number Generation
4689 17.7 Copying random number generator state
4690 ==========================================
4692 The above methods do not expose the random number `state' which changes
4693 from call to call. It is often useful to be able to save and restore
4694 the state. To permit these practices, a few somewhat more advanced
4695 functions are supplied. These include:
4697 -- Function: int gsl_rng_memcpy (gsl_rng * DEST, const gsl_rng * SRC)
4698 This function copies the random number generator SRC into the
4699 pre-existing generator DEST, making DEST into an exact copy of
4700 SRC. The two generators must be of the same type.
4702 -- Function: gsl_rng * gsl_rng_clone (const gsl_rng * R)
4703 This function returns a pointer to a newly created generator which
4704 is an exact copy of the generator R.
4707 File: gsl-ref.info, Node: Reading and writing random number generator state, Next: Random number generator algorithms, Prev: Copying random number generator state, Up: Random Number Generation
4709 17.8 Reading and writing random number generator state
4710 ======================================================
4712 The library provides functions for reading and writing the random
4713 number state to a file as binary data or formatted text.
4715 -- Function: int gsl_rng_fwrite (FILE * STREAM, const gsl_rng * R)
4716 This function writes the random number state of the random number
4717 generator R to the stream STREAM in binary format. The return
4718 value is 0 for success and `GSL_EFAILED' if there was a problem
4719 writing to the file. Since the data is written in the native
4720 binary format it may not be portable between different
4723 -- Function: int gsl_rng_fread (FILE * STREAM, gsl_rng * R)
4724 This function reads the random number state into the random number
4725 generator R from the open stream STREAM in binary format. The
4726 random number generator R must be preinitialized with the correct
4727 random number generator type since type information is not saved.
4728 The return value is 0 for success and `GSL_EFAILED' if there was a
4729 problem reading from the file. The data is assumed to have been
4730 written in the native binary format on the same architecture.
4733 File: gsl-ref.info, Node: Random number generator algorithms, Next: Unix random number generators, Prev: Reading and writing random number generator state, Up: Random Number Generation
4735 17.9 Random number generator algorithms
4736 =======================================
4738 The functions described above make no reference to the actual algorithm
4739 used. This is deliberate so that you can switch algorithms without
4740 having to change any of your application source code. The library
4741 provides a large number of generators of different types, including
4742 simulation quality generators, generators provided for compatibility
4743 with other libraries and historical generators from the past.
4745 The following generators are recommended for use in simulation. They
4746 have extremely long periods, low correlation and pass most statistical
4747 tests. For the most reliable source of uncorrelated numbers, the
4748 second-generation RANLUX generators have the strongest proof of
4751 -- Generator: gsl_rng_mt19937
4752 The MT19937 generator of Makoto Matsumoto and Takuji Nishimura is a
4753 variant of the twisted generalized feedback shift-register
4754 algorithm, and is known as the "Mersenne Twister" generator. It
4755 has a Mersenne prime period of 2^19937 - 1 (about 10^6000) and is
4756 equi-distributed in 623 dimensions. It has passed the DIEHARD
4757 statistical tests. It uses 624 words of state per generator and is
4758 comparable in speed to the other generators. The original
4759 generator used a default seed of 4357 and choosing S equal to zero
4760 in `gsl_rng_set' reproduces this. Later versions switched to 5489
4761 as the default seed, you can choose this explicitly via
4762 `gsl_rng_set' instead if you require it.
4764 For more information see,
4765 Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
4766 623-dimensionally equidistributed uniform pseudorandom number
4767 generator". `ACM Transactions on Modeling and Computer
4768 Simulation', Vol. 8, No. 1 (Jan. 1998), Pages 3-30
4770 The generator `gsl_rng_mt19937' uses the second revision of the
4771 seeding procedure published by the two authors above in 2002. The
4772 original seeding procedures could cause spurious artifacts for
4773 some seed values. They are still available through the alternative
4774 generators `gsl_rng_mt19937_1999' and `gsl_rng_mt19937_1998'.
4776 -- Generator: gsl_rng_ranlxs0
4777 -- Generator: gsl_rng_ranlxs1
4778 -- Generator: gsl_rng_ranlxs2
4779 The generator `ranlxs0' is a second-generation version of the
4780 RANLUX algorithm of Lu"scher, which produces "luxury random
4781 numbers". This generator provides single precision output (24
4782 bits) at three luxury levels `ranlxs0', `ranlxs1' and `ranlxs2',
4783 in increasing order of strength. It uses double-precision
4784 floating point arithmetic internally and can be significantly
4785 faster than the integer version of `ranlux', particularly on
4786 64-bit architectures. The period of the generator is about
4787 10^171. The algorithm has mathematically proven properties and
4788 can provide truly decorrelated numbers at a known level of
4789 randomness. The higher luxury levels provide increased
4790 decorrelation between samples as an additional safety margin.
4792 -- Generator: gsl_rng_ranlxd1
4793 -- Generator: gsl_rng_ranlxd2
4794 These generators produce double precision output (48 bits) from the
4795 RANLXS generator. The library provides two luxury levels
4796 `ranlxd1' and `ranlxd2', in increasing order of strength.
4798 -- Generator: gsl_rng_ranlux
4799 -- Generator: gsl_rng_ranlux389
4800 The `ranlux' generator is an implementation of the original
4801 algorithm developed by Lu"scher. It uses a
4802 lagged-fibonacci-with-skipping algorithm to produce "luxury random
4803 numbers". It is a 24-bit generator, originally designed for
4804 single-precision IEEE floating point numbers. This implementation
4805 is based on integer arithmetic, while the second-generation
4806 versions RANLXS and RANLXD described above provide floating-point
4807 implementations which will be faster on many platforms. The
4808 period of the generator is about 10^171. The algorithm has
4809 mathematically proven properties and it can provide truly
4810 decorrelated numbers at a known level of randomness. The default
4811 level of decorrelation recommended by Lu"scher is provided by
4812 `gsl_rng_ranlux', while `gsl_rng_ranlux389' gives the highest
4813 level of randomness, with all 24 bits decorrelated. Both types of
4814 generator use 24 words of state per generator.
4816 For more information see,
4817 M. Lu"scher, "A portable high-quality random number generator
4818 for lattice field theory calculations", `Computer Physics
4819 Communications', 79 (1994) 100-110.
4821 F. James, "RANLUX: A Fortran implementation of the
4822 high-quality pseudo-random number generator of Lu"scher",
4823 `Computer Physics Communications', 79 (1994) 111-114
4825 -- Generator: gsl_rng_cmrg
4826 This is a combined multiple recursive generator by L'Ecuyer. Its
4829 z_n = (x_n - y_n) mod m_1
4831 where the two underlying generators x_n and y_n are,
4833 x_n = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) mod m_1
4834 y_n = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) mod m_2
4836 with coefficients a_1 = 0, a_2 = 63308, a_3 = -183326, b_1 = 86098,
4837 b_2 = 0, b_3 = -539608, and moduli m_1 = 2^31 - 1 = 2147483647 and
4840 The period of this generator is lcm(m_1^3-1, m_2^3-1), which is
4841 approximately 2^185 (about 10^56). It uses 6 words of state per
4842 generator. For more information see,
4844 P. L'Ecuyer, "Combined Multiple Recursive Random Number
4845 Generators", `Operations Research', 44, 5 (1996), 816-822.
4847 -- Generator: gsl_rng_mrg
4848 This is a fifth-order multiple recursive generator by L'Ecuyer,
4849 Blouin and Coutre. Its sequence is,
4851 x_n = (a_1 x_{n-1} + a_5 x_{n-5}) mod m
4853 with a_1 = 107374182, a_2 = a_3 = a_4 = 0, a_5 = 104480 and m =
4856 The period of this generator is about 10^46. It uses 5 words of
4857 state per generator. More information can be found in the
4859 P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good
4860 multiple recursive random number generators", `ACM
4861 Transactions on Modeling and Computer Simulation' 3, 87-98
4864 -- Generator: gsl_rng_taus
4865 -- Generator: gsl_rng_taus2
4866 This is a maximally equidistributed combined Tausworthe generator
4867 by L'Ecuyer. The sequence is,
4869 x_n = (s1_n ^^ s2_n ^^ s3_n)
4873 s1_{n+1} = (((s1_n&4294967294)<<12)^^(((s1_n<<13)^^s1_n)>>19))
4874 s2_{n+1} = (((s2_n&4294967288)<< 4)^^(((s2_n<< 2)^^s2_n)>>25))
4875 s3_{n+1} = (((s3_n&4294967280)<<17)^^(((s3_n<< 3)^^s3_n)>>11))
4877 computed modulo 2^32. In the formulas above ^^ denotes
4878 "exclusive-or". Note that the algorithm relies on the properties
4879 of 32-bit unsigned integers and has been implemented using a
4880 bitmask of `0xFFFFFFFF' to make it work on 64 bit machines.
4882 The period of this generator is 2^88 (about 10^26). It uses 3
4883 words of state per generator. For more information see,
4885 P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
4886 Generators", `Mathematics of Computation', 65, 213 (1996),
4889 The generator `gsl_rng_taus2' uses the same algorithm as
4890 `gsl_rng_taus' but with an improved seeding procedure described in
4893 P. L'Ecuyer, "Tables of Maximally Equidistributed Combined
4894 LFSR Generators", `Mathematics of Computation', 68, 225
4897 The generator `gsl_rng_taus2' should now be used in preference to
4900 -- Generator: gsl_rng_gfsr4
4901 The `gfsr4' generator is like a lagged-fibonacci generator, and
4902 produces each number as an `xor''d sum of four previous values.
4904 r_n = r_{n-A} ^^ r_{n-B} ^^ r_{n-C} ^^ r_{n-D}
4906 Ziff (ref below) notes that "it is now widely known" that two-tap
4907 registers (such as R250, which is described below) have serious
4908 flaws, the most obvious one being the three-point correlation that
4909 comes from the definition of the generator. Nice mathematical
4910 properties can be derived for GFSR's, and numerics bears out the
4911 claim that 4-tap GFSR's with appropriately chosen offsets are as
4912 random as can be measured, using the author's test.
4914 This implementation uses the values suggested the example on p392
4915 of Ziff's article: A=471, B=1586, C=6988, D=9689.
4917 If the offsets are appropriately chosen (such as the one ones in
4918 this implementation), then the sequence is said to be maximal;
4919 that means that the period is 2^D - 1, where D is the longest lag.
4920 (It is one less than 2^D because it is not permitted to have all
4921 zeros in the `ra[]' array.) For this implementation with D=9689
4922 that works out to about 10^2917.
4924 Note that the implementation of this generator using a 32-bit
4925 integer amounts to 32 parallel implementations of one-bit
4926 generators. One consequence of this is that the period of this
4927 32-bit generator is the same as for the one-bit generator.
4928 Moreover, this independence means that all 32-bit patterns are
4929 equally likely, and in particular that 0 is an allowed random
4930 value. (We are grateful to Heiko Bauke for clarifying for us these
4931 properties of GFSR random number generators.)
4933 For more information see,
4934 Robert M. Ziff, "Four-tap shift-register-sequence
4935 random-number generators", `Computers in Physics', 12(4),
4936 Jul/Aug 1998, pp 385-392.
4939 File: gsl-ref.info, Node: Unix random number generators, Next: Other random number generators, Prev: Random number generator algorithms, Up: Random Number Generation
4941 17.10 Unix random number generators
4942 ===================================
4944 The standard Unix random number generators `rand', `random' and
4945 `rand48' are provided as part of GSL. Although these generators are
4946 widely available individually often they aren't all available on the
4947 same platform. This makes it difficult to write portable code using
4948 them and so we have included the complete set of Unix generators in GSL
4949 for convenience. Note that these generators don't produce high-quality
4950 randomness and aren't suitable for work requiring accurate statistics.
4951 However, if you won't be measuring statistical quantities and just want
4952 to introduce some variation into your program then these generators are
4955 -- Generator: gsl_rng_rand
4956 This is the BSD `rand' generator. Its sequence is
4958 x_{n+1} = (a x_n + c) mod m
4960 with a = 1103515245, c = 12345 and m = 2^31. The seed specifies
4961 the initial value, x_1. The period of this generator is 2^31, and
4962 it uses 1 word of storage per generator.
4964 -- Generator: gsl_rng_random_bsd
4965 -- Generator: gsl_rng_random_libc5
4966 -- Generator: gsl_rng_random_glibc2
4967 These generators implement the `random' family of functions, a set
4968 of linear feedback shift register generators originally used in BSD
4969 Unix. There are several versions of `random' in use today: the
4970 original BSD version (e.g. on SunOS4), a libc5 version (found on
4971 older GNU/Linux systems) and a glibc2 version. Each version uses a
4972 different seeding procedure, and thus produces different sequences.
4974 The original BSD routines accepted a variable length buffer for the
4975 generator state, with longer buffers providing higher-quality
4976 randomness. The `random' function implemented algorithms for
4977 buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm
4978 with the largest length that would fit into the user-supplied
4979 buffer was used. To support these algorithms additional
4980 generators are available with the following names,
4983 gsl_rng_random32_bsd
4984 gsl_rng_random64_bsd
4985 gsl_rng_random128_bsd
4986 gsl_rng_random256_bsd
4988 where the numeric suffix indicates the buffer length. The
4989 original BSD `random' function used a 128-byte default buffer and
4990 so `gsl_rng_random_bsd' has been made equivalent to
4991 `gsl_rng_random128_bsd'. Corresponding versions of the `libc5'
4992 and `glibc2' generators are also available, with the names
4993 `gsl_rng_random8_libc5', `gsl_rng_random8_glibc2', etc.
4995 -- Generator: gsl_rng_rand48
4996 This is the Unix `rand48' generator. Its sequence is
4998 x_{n+1} = (a x_n + c) mod m
5000 defined on 48-bit unsigned integers with a = 25214903917, c = 11
5001 and m = 2^48. The seed specifies the upper 32 bits of the initial
5002 value, x_1, with the lower 16 bits set to `0x330E'. The function
5003 `gsl_rng_get' returns the upper 32 bits from each term of the
5004 sequence. This does not have a direct parallel in the original
5005 `rand48' functions, but forcing the result to type `long int'
5006 reproduces the output of `mrand48'. The function
5007 `gsl_rng_uniform' uses the full 48 bits of internal state to return
5008 the double precision number x_n/m, which is equivalent to the
5009 function `drand48'. Note that some versions of the GNU C Library
5010 contained a bug in `mrand48' function which caused it to produce
5011 different results (only the lower 16-bits of the return value were
5015 File: gsl-ref.info, Node: Other random number generators, Next: Random Number Generator Performance, Prev: Unix random number generators, Up: Random Number Generation
5017 17.11 Other random number generators
5018 ====================================
5020 The generators in this section are provided for compatibility with
5021 existing libraries. If you are converting an existing program to use
5022 GSL then you can select these generators to check your new
5023 implementation against the original one, using the same random number
5024 generator. After verifying that your new program reproduces the
5025 original results you can then switch to a higher-quality generator.
5027 Note that most of the generators in this section are based on single
5028 linear congruence relations, which are the least sophisticated type of
5029 generator. In particular, linear congruences have poor properties when
5030 used with a non-prime modulus, as several of these routines do (e.g.
5031 with a power of two modulus, 2^31 or 2^32). This leads to periodicity
5032 in the least significant bits of each number, with only the higher bits
5033 having any randomness. Thus if you want to produce a random bitstream
5034 it is best to avoid using the least significant bits.
5036 -- Generator: gsl_rng_ranf
5037 This is the CRAY random number generator `RANF'. Its sequence is
5039 x_{n+1} = (a x_n) mod m
5041 defined on 48-bit unsigned integers with a = 44485709377909 and m
5042 = 2^48. The seed specifies the lower 32 bits of the initial value,
5043 x_1, with the lowest bit set to prevent the seed taking an even
5044 value. The upper 16 bits of x_1 are set to 0. A consequence of
5045 this procedure is that the pairs of seeds 2 and 3, 4 and 5, etc
5046 produce the same sequences.
5048 The generator compatible with the CRAY MATHLIB routine RANF. It
5049 produces double precision floating point numbers which should be
5050 identical to those from the original RANF.
5052 There is a subtlety in the implementation of the seeding. The
5053 initial state is reversed through one step, by multiplying by the
5054 modular inverse of a mod m. This is done for compatibility with
5055 the original CRAY implementation.
5057 Note that you can only seed the generator with integers up to
5058 2^32, while the original CRAY implementation uses non-portable
5059 wide integers which can cover all 2^48 states of the generator.
5061 The function `gsl_rng_get' returns the upper 32 bits from each term
5062 of the sequence. The function `gsl_rng_uniform' uses the full 48
5063 bits to return the double precision number x_n/m.
5065 The period of this generator is 2^46.
5067 -- Generator: gsl_rng_ranmar
5068 This is the RANMAR lagged-fibonacci generator of Marsaglia, Zaman
5069 and Tsang. It is a 24-bit generator, originally designed for
5070 single-precision IEEE floating point numbers. It was included in
5071 the CERNLIB high-energy physics library.
5073 -- Generator: gsl_rng_r250
5074 This is the shift-register generator of Kirkpatrick and Stoll. The
5075 sequence is based on the recurrence
5077 x_n = x_{n-103} ^^ x_{n-250}
5079 where ^^ denotes "exclusive-or", defined on 32-bit words. The
5080 period of this generator is about 2^250 and it uses 250 words of
5081 state per generator.
5083 For more information see,
5084 S. Kirkpatrick and E. Stoll, "A very fast shift-register
5085 sequence random number generator", `Journal of Computational
5086 Physics', 40, 517-526 (1981)
5088 -- Generator: gsl_rng_tt800
5089 This is an earlier version of the twisted generalized feedback
5090 shift-register generator, and has been superseded by the
5091 development of MT19937. However, it is still an acceptable
5092 generator in its own right. It has a period of 2^800 and uses 33
5093 words of storage per generator.
5095 For more information see,
5096 Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR
5097 Generators II", `ACM Transactions on Modelling and Computer
5098 Simulation', Vol. 4, No. 3, 1994, pages 254-266.
5100 -- Generator: gsl_rng_vax
5101 This is the VAX generator `MTH$RANDOM'. Its sequence is,
5103 x_{n+1} = (a x_n + c) mod m
5105 with a = 69069, c = 1 and m = 2^32. The seed specifies the
5106 initial value, x_1. The period of this generator is 2^32 and it
5107 uses 1 word of storage per generator.
5109 -- Generator: gsl_rng_transputer
5110 This is the random number generator from the INMOS Transputer
5111 Development system. Its sequence is,
5113 x_{n+1} = (a x_n) mod m
5115 with a = 1664525 and m = 2^32. The seed specifies the initial
5118 -- Generator: gsl_rng_randu
5119 This is the IBM `RANDU' generator. Its sequence is
5121 x_{n+1} = (a x_n) mod m
5123 with a = 65539 and m = 2^31. The seed specifies the initial value,
5124 x_1. The period of this generator was only 2^29. It has become a
5125 textbook example of a poor generator.
5127 -- Generator: gsl_rng_minstd
5128 This is Park and Miller's "minimal standard" MINSTD generator, a
5129 simple linear congruence which takes care to avoid the major
5130 pitfalls of such algorithms. Its sequence is,
5132 x_{n+1} = (a x_n) mod m
5134 with a = 16807 and m = 2^31 - 1 = 2147483647. The seed specifies
5135 the initial value, x_1. The period of this generator is about
5138 This generator is used in the IMSL Library (subroutine RNUN) and in
5139 MATLAB (the RAND function). It is also sometimes known by the
5140 acronym "GGL" (I'm not sure what that stands for).
5142 For more information see,
5143 Park and Miller, "Random Number Generators: Good ones are
5144 hard to find", `Communications of the ACM', October 1988,
5145 Volume 31, No 10, pages 1192-1201.
5147 -- Generator: gsl_rng_uni
5148 -- Generator: gsl_rng_uni32
5149 This is a reimplementation of the 16-bit SLATEC random number
5150 generator RUNIF. A generalization of the generator to 32 bits is
5151 provided by `gsl_rng_uni32'. The original source code is
5152 available from NETLIB.
5154 -- Generator: gsl_rng_slatec
5155 This is the SLATEC random number generator RAND. It is ancient.
5156 The original source code is available from NETLIB.
5158 -- Generator: gsl_rng_zuf
5159 This is the ZUFALL lagged Fibonacci series generator of Peterson.
5162 t = u_{n-273} + u_{n-607}
5165 The original source code is available from NETLIB. For more
5167 W. Petersen, "Lagged Fibonacci Random Number Generators for
5168 the NEC SX-3", `International Journal of High Speed
5171 -- Generator: gsl_rng_knuthran2
5172 This is a second-order multiple recursive generator described by
5173 Knuth in `Seminumerical Algorithms', 3rd Ed., page 108. Its
5176 x_n = (a_1 x_{n-1} + a_2 x_{n-2}) mod m
5178 with a_1 = 271828183, a_2 = 314159269, and m = 2^31 - 1.
5180 -- Generator: gsl_rng_knuthran2002
5181 -- Generator: gsl_rng_knuthran
5182 This is a second-order multiple recursive generator described by
5183 Knuth in `Seminumerical Algorithms', 3rd Ed., Section 3.6. Knuth
5184 provides its C code. The updated routine `gsl_rng_knuthran2002'
5185 is from the revised 9th printing and corrects some weaknesses in
5186 the earlier version, which is implemented as `gsl_rng_knuthran'.
5188 -- Generator: gsl_rng_borosh13
5189 -- Generator: gsl_rng_fishman18
5190 -- Generator: gsl_rng_fishman20
5191 -- Generator: gsl_rng_lecuyer21
5192 -- Generator: gsl_rng_waterman14
5193 These multiplicative generators are taken from Knuth's
5194 `Seminumerical Algorithms', 3rd Ed., pages 106-108. Their sequence
5197 x_{n+1} = (a x_n) mod m
5199 where the seed specifies the initial value, x_1. The parameters a
5200 and m are as follows, Borosh-Niederreiter: a = 1812433253, m =
5201 2^32, Fishman18: a = 62089911, m = 2^31 - 1, Fishman20: a = 48271,
5202 m = 2^31 - 1, L'Ecuyer: a = 40692, m = 2^31 - 249, Waterman: a =
5203 1566083941, m = 2^32.
5205 -- Generator: gsl_rng_fishman2x
5206 This is the L'Ecuyer-Fishman random number generator. It is taken
5207 from Knuth's `Seminumerical Algorithms', 3rd Ed., page 108. Its
5210 z_{n+1} = (x_n - y_n) mod m
5212 with m = 2^31 - 1. x_n and y_n are given by the `fishman20' and
5213 `lecuyer21' algorithms. The seed specifies the initial value, x_1.
5216 -- Generator: gsl_rng_coveyou
5217 This is the Coveyou random number generator. It is taken from
5218 Knuth's `Seminumerical Algorithms', 3rd Ed., Section 3.2.2. Its
5221 x_{n+1} = (x_n (x_n + 1)) mod m
5223 with m = 2^32. The seed specifies the initial value, x_1.
5226 File: gsl-ref.info, Node: Random Number Generator Performance, Next: Random Number Generator Examples, Prev: Other random number generators, Up: Random Number Generation
5231 The following table shows the relative performance of a selection the
5232 available random number generators. The fastest simulation quality
5233 generators are `taus', `gfsr4' and `mt19937'. The generators which
5234 offer the best mathematically-proven quality are those based on the
5237 1754 k ints/sec, 870 k doubles/sec, taus
5238 1613 k ints/sec, 855 k doubles/sec, gfsr4
5239 1370 k ints/sec, 769 k doubles/sec, mt19937
5240 565 k ints/sec, 571 k doubles/sec, ranlxs0
5241 400 k ints/sec, 405 k doubles/sec, ranlxs1
5242 490 k ints/sec, 389 k doubles/sec, mrg
5243 407 k ints/sec, 297 k doubles/sec, ranlux
5244 243 k ints/sec, 254 k doubles/sec, ranlxd1
5245 251 k ints/sec, 253 k doubles/sec, ranlxs2
5246 238 k ints/sec, 215 k doubles/sec, cmrg
5247 247 k ints/sec, 198 k doubles/sec, ranlux389
5248 141 k ints/sec, 140 k doubles/sec, ranlxd2
5250 1852 k ints/sec, 935 k doubles/sec, ran3
5251 813 k ints/sec, 575 k doubles/sec, ran0
5252 787 k ints/sec, 476 k doubles/sec, ran1
5253 379 k ints/sec, 292 k doubles/sec, ran2
5256 File: gsl-ref.info, Node: Random Number Generator Examples, Next: Random Number References and Further Reading, Prev: Random Number Generator Performance, Up: Random Number Generation
5261 The following program demonstrates the use of a random number generator
5262 to produce uniform random numbers in the range [0.0, 1.0),
5265 #include <gsl/gsl_rng.h>
5270 const gsl_rng_type * T;
5275 gsl_rng_env_setup();
5277 T = gsl_rng_default;
5278 r = gsl_rng_alloc (T);
5280 for (i = 0; i < n; i++)
5282 double u = gsl_rng_uniform (r);
5283 printf ("%.5f\n", u);
5291 Here is the output of the program,
5305 The numbers depend on the seed used by the generator. The default seed
5306 can be changed with the `GSL_RNG_SEED' environment variable to produce
5307 a different stream of numbers. The generator itself can be changed
5308 using the environment variable `GSL_RNG_TYPE'. Here is the output of
5309 the program using a seed value of 123 and the multiple-recursive
5312 $ GSL_RNG_SEED=123 GSL_RNG_TYPE=mrg ./a.out
5327 File: gsl-ref.info, Node: Random Number References and Further Reading, Next: Random Number Acknowledgements, Prev: Random Number Generator Examples, Up: Random Number Generation
5329 17.14 References and Further Reading
5330 ====================================
5332 The subject of random number generation and testing is reviewed
5333 extensively in Knuth's `Seminumerical Algorithms'.
5335 Donald E. Knuth, `The Art of Computer Programming: Seminumerical
5336 Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
5338 Further information is available in the review paper written by Pierre
5341 P. L'Ecuyer, "Random Number Generation", Chapter 4 of the Handbook
5342 on Simulation, Jerry Banks Ed., Wiley, 1998, 93-137.
5344 `http://www.iro.umontreal.ca/~lecuyer/papers.html' in the file
5347 The source code for the DIEHARD random number generator tests is also
5350 `DIEHARD source code' G. Marsaglia,
5352 `http://stat.fsu.edu/pub/diehard/'
5354 A comprehensive set of random number generator tests is available from
5357 NIST Special Publication 800-22, "A Statistical Test Suite for the
5358 Validation of Random Number Generators and Pseudo Random Number
5359 Generators for Cryptographic Applications".
5361 `http://csrc.nist.gov/rng/'
5364 File: gsl-ref.info, Node: Random Number Acknowledgements, Prev: Random Number References and Further Reading, Up: Random Number Generation
5366 17.15 Acknowledgements
5367 ======================
5369 Thanks to Makoto Matsumoto, Takuji Nishimura and Yoshiharu Kurita for
5370 making the source code to their generators (MT19937, MM&TN; TT800,
5371 MM&YK) available under the GNU General Public License. Thanks to Martin
5372 Lu"scher for providing notes and source code for the RANLXS and RANLXD
5376 File: gsl-ref.info, Node: Quasi-Random Sequences, Next: Random Number Distributions, Prev: Random Number Generation, Up: Top
5378 18 Quasi-Random Sequences
5379 *************************
5381 This chapter describes functions for generating quasi-random sequences
5382 in arbitrary dimensions. A quasi-random sequence progressively covers a
5383 d-dimensional space with a set of points that are uniformly
5384 distributed. Quasi-random sequences are also known as low-discrepancy
5385 sequences. The quasi-random sequence generators use an interface that
5386 is similar to the interface for random number generators, except that
5387 seeding is not required--each generator produces a single sequence.
5389 The functions described in this section are declared in the header
5394 * Quasi-random number generator initialization::
5395 * Sampling from a quasi-random number generator::
5396 * Auxiliary quasi-random number generator functions::
5397 * Saving and resorting quasi-random number generator state::
5398 * Quasi-random number generator algorithms::
5399 * Quasi-random number generator examples::
5400 * Quasi-random number references::
5403 File: gsl-ref.info, Node: Quasi-random number generator initialization, Next: Sampling from a quasi-random number generator, Up: Quasi-Random Sequences
5405 18.1 Quasi-random number generator initialization
5406 =================================================
5408 -- Function: gsl_qrng * gsl_qrng_alloc (const gsl_qrng_type * T,
5410 This function returns a pointer to a newly-created instance of a
5411 quasi-random sequence generator of type T and dimension D. If
5412 there is insufficient memory to create the generator then the
5413 function returns a null pointer and the error handler is invoked
5414 with an error code of `GSL_ENOMEM'.
5416 -- Function: void gsl_qrng_free (gsl_qrng * Q)
5417 This function frees all the memory associated with the generator Q.
5419 -- Function: void gsl_qrng_init (gsl_qrng * Q)
5420 This function reinitializes the generator Q to its starting point.
5421 Note that quasi-random sequences do not use a seed and always
5422 produce the same set of values.
5425 File: gsl-ref.info, Node: Sampling from a quasi-random number generator, Next: Auxiliary quasi-random number generator functions, Prev: Quasi-random number generator initialization, Up: Quasi-Random Sequences
5427 18.2 Sampling from a quasi-random number generator
5428 ==================================================
5430 -- Function: int gsl_qrng_get (const gsl_qrng * Q, double X[])
5431 This function stores the next point from the sequence generator Q
5432 in the array X. The space available for X must match the
5433 dimension of the generator. The point X will lie in the range 0 <
5434 x_i < 1 for each x_i. An inline version of this function is used
5435 when `HAVE_INLINE' is defined.
5438 File: gsl-ref.info, Node: Auxiliary quasi-random number generator functions, Next: Saving and resorting quasi-random number generator state, Prev: Sampling from a quasi-random number generator, Up: Quasi-Random Sequences
5440 18.3 Auxiliary quasi-random number generator functions
5441 ======================================================
5443 -- Function: const char * gsl_qrng_name (const gsl_qrng * Q)
5444 This function returns a pointer to the name of the generator.
5446 -- Function: size_t gsl_qrng_size (const gsl_qrng * Q)
5447 -- Function: void * gsl_qrng_state (const gsl_qrng * Q)
5448 These functions return a pointer to the state of generator R and
5449 its size. You can use this information to access the state
5450 directly. For example, the following code will write the state of
5451 a generator to a stream,
5453 void * state = gsl_qrng_state (q);
5454 size_t n = gsl_qrng_size (q);
5455 fwrite (state, n, 1, stream);
5458 File: gsl-ref.info, Node: Saving and resorting quasi-random number generator state, Next: Quasi-random number generator algorithms, Prev: Auxiliary quasi-random number generator functions, Up: Quasi-Random Sequences
5460 18.4 Saving and resorting quasi-random number generator state
5461 =============================================================
5463 -- Function: int gsl_qrng_memcpy (gsl_qrng * DEST, const gsl_qrng *
5465 This function copies the quasi-random sequence generator SRC into
5466 the pre-existing generator DEST, making DEST into an exact copy of
5467 SRC. The two generators must be of the same type.
5469 -- Function: gsl_qrng * gsl_qrng_clone (const gsl_qrng * Q)
5470 This function returns a pointer to a newly created generator which
5471 is an exact copy of the generator Q.
5474 File: gsl-ref.info, Node: Quasi-random number generator algorithms, Next: Quasi-random number generator examples, Prev: Saving and resorting quasi-random number generator state, Up: Quasi-Random Sequences
5476 18.5 Quasi-random number generator algorithms
5477 =============================================
5479 The following quasi-random sequence algorithms are available,
5481 -- Generator: gsl_qrng_niederreiter_2
5482 This generator uses the algorithm described in Bratley, Fox,
5483 Niederreiter, `ACM Trans. Model. Comp. Sim.' 2, 195 (1992). It is
5484 valid up to 12 dimensions.
5486 -- Generator: gsl_qrng_sobol
5487 This generator uses the Sobol sequence described in Antonov,
5488 Saleev, `USSR Comput. Maths. Math. Phys.' 19, 252 (1980). It is
5489 valid up to 40 dimensions.
5491 -- Generator: gsl_qrng_halton
5492 -- Generator: gsl_qrng_reverse_halton
5493 These generators use the Halton and reverse Halton sequences
5494 described in J.H. Halton, `Numerische Mathematik' 2, 84-90 (1960)
5495 and B. Vandewoestyne and R. Cools `Computational and Applied
5496 Mathematics' 189, 1&2, 341-361 (2006). They are valid up to 1229
5500 File: gsl-ref.info, Node: Quasi-random number generator examples, Next: Quasi-random number references, Prev: Quasi-random number generator algorithms, Up: Quasi-Random Sequences
5505 The following program prints the first 1024 points of the 2-dimensional
5509 #include <gsl/gsl_qrng.h>
5515 gsl_qrng * q = gsl_qrng_alloc (gsl_qrng_sobol, 2);
5517 for (i = 0; i < 1024; i++)
5520 gsl_qrng_get (q, v);
5521 printf ("%.5f %.5f\n", v[0], v[1]);
5528 Here is the output from the program,
5540 It can be seen that successive points progressively fill-in the spaces
5541 between previous points.
5544 File: gsl-ref.info, Node: Quasi-random number references, Prev: Quasi-random number generator examples, Up: Quasi-Random Sequences
5549 The implementations of the quasi-random sequence routines are based on
5550 the algorithms described in the following paper,
5552 P. Bratley and B.L. Fox and H. Niederreiter, "Algorithm 738:
5553 Programs to Generate Niederreiter's Low-discrepancy Sequences",
5554 `ACM Transactions on Mathematical Software', Vol. 20, No. 4,
5555 December, 1994, p. 494-495.
5558 File: gsl-ref.info, Node: Random Number Distributions, Next: Statistics, Prev: Quasi-Random Sequences, Up: Top
5560 19 Random Number Distributions
5561 ******************************
5563 This chapter describes functions for generating random variates and
5564 computing their probability distributions. Samples from the
5565 distributions described in this chapter can be obtained using any of the
5566 random number generators in the library as an underlying source of
5569 In the simplest cases a non-uniform distribution can be obtained
5570 analytically from the uniform distribution of a random number generator
5571 by applying an appropriate transformation. This method uses one call to
5572 the random number generator. More complicated distributions are created
5573 by the "acceptance-rejection" method, which compares the desired
5574 distribution against a distribution which is similar and known
5575 analytically. This usually requires several samples from the generator.
5577 The library also provides cumulative distribution functions and
5578 inverse cumulative distribution functions, sometimes referred to as
5579 quantile functions. The cumulative distribution functions and their
5580 inverses are computed separately for the upper and lower tails of the
5581 distribution, allowing full accuracy to be retained for small results.
5583 The functions for random variates and probability density functions
5584 described in this section are declared in `gsl_randist.h'. The
5585 corresponding cumulative distribution functions are declared in
5588 Note that the discrete random variate functions always return a
5589 value of type `unsigned int', and on most platforms this has a maximum
5590 value of 2^32-1 ~=~ 4.29e9. They should only be called with a safe
5591 range of parameters (where there is a negligible probability of a
5592 variate exceeding this limit) to prevent incorrect results due to
5597 * Random Number Distribution Introduction::
5598 * The Gaussian Distribution::
5599 * The Gaussian Tail Distribution::
5600 * The Bivariate Gaussian Distribution::
5601 * The Exponential Distribution::
5602 * The Laplace Distribution::
5603 * The Exponential Power Distribution::
5604 * The Cauchy Distribution::
5605 * The Rayleigh Distribution::
5606 * The Rayleigh Tail Distribution::
5607 * The Landau Distribution::
5608 * The Levy alpha-Stable Distributions::
5609 * The Levy skew alpha-Stable Distribution::
5610 * The Gamma Distribution::
5611 * The Flat (Uniform) Distribution::
5612 * The Lognormal Distribution::
5613 * The Chi-squared Distribution::
5614 * The F-distribution::
5615 * The t-distribution::
5616 * The Beta Distribution::
5617 * The Logistic Distribution::
5618 * The Pareto Distribution::
5619 * Spherical Vector Distributions::
5620 * The Weibull Distribution::
5621 * The Type-1 Gumbel Distribution::
5622 * The Type-2 Gumbel Distribution::
5623 * The Dirichlet Distribution::
5624 * General Discrete Distributions::
5625 * The Poisson Distribution::
5626 * The Bernoulli Distribution::
5627 * The Binomial Distribution::
5628 * The Multinomial Distribution::
5629 * The Negative Binomial Distribution::
5630 * The Pascal Distribution::
5631 * The Geometric Distribution::
5632 * The Hypergeometric Distribution::
5633 * The Logarithmic Distribution::
5634 * Shuffling and Sampling::
5635 * Random Number Distribution Examples::
5636 * Random Number Distribution References and Further Reading::
5639 File: gsl-ref.info, Node: Random Number Distribution Introduction, Next: The Gaussian Distribution, Up: Random Number Distributions
5644 Continuous random number distributions are defined by a probability
5645 density function, p(x), such that the probability of x occurring in the
5646 infinitesimal range x to x+dx is p dx.
5648 The cumulative distribution function for the lower tail P(x) is
5649 defined by the integral,
5651 P(x) = \int_{-\infty}^{x} dx' p(x')
5653 and gives the probability of a variate taking a value less than x.
5655 The cumulative distribution function for the upper tail Q(x) is
5656 defined by the integral,
5658 Q(x) = \int_{x}^{+\infty} dx' p(x')
5660 and gives the probability of a variate taking a value greater than x.
5662 The upper and lower cumulative distribution functions are related by
5663 P(x) + Q(x) = 1 and satisfy 0 <= P(x) <= 1, 0 <= Q(x) <= 1.
5665 The inverse cumulative distributions, x=P^{-1}(P) and x=Q^{-1}(Q)
5666 give the values of x which correspond to a specific value of P or Q.
5667 They can be used to find confidence limits from probability values.
5669 For discrete distributions the probability of sampling the integer
5670 value k is given by p(k), where \sum_k p(k) = 1. The cumulative
5671 distribution for the lower tail P(k) of a discrete distribution is
5674 P(k) = \sum_{i <= k} p(i)
5676 where the sum is over the allowed range of the distribution less than
5679 The cumulative distribution for the upper tail of a discrete
5680 distribution Q(k) is defined as
5682 Q(k) = \sum_{i > k} p(i)
5684 giving the sum of probabilities for all values greater than k. These
5685 two definitions satisfy the identity P(k)+Q(k)=1.
5687 If the range of the distribution is 1 to n inclusive then P(n)=1,
5688 Q(n)=0 while P(1) = p(1), Q(1)=1-p(1).
5691 File: gsl-ref.info, Node: The Gaussian Distribution, Next: The Gaussian Tail Distribution, Prev: Random Number Distribution Introduction, Up: Random Number Distributions
5693 19.2 The Gaussian Distribution
5694 ==============================
5696 -- Function: double gsl_ran_gaussian (const gsl_rng * R, double SIGMA)
5697 This function returns a Gaussian random variate, with mean zero and
5698 standard deviation SIGMA. The probability distribution for
5699 Gaussian random variates is,
5701 p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
5703 for x in the range -\infty to +\infty. Use the transformation z =
5704 \mu + x on the numbers returned by `gsl_ran_gaussian' to obtain a
5705 Gaussian distribution with mean \mu. This function uses the
5706 Box-Mueller algorithm which requires two calls to the random
5709 -- Function: double gsl_ran_gaussian_pdf (double X, double SIGMA)
5710 This function computes the probability density p(x) at X for a
5711 Gaussian distribution with standard deviation SIGMA, using the
5712 formula given above.
5715 -- Function: double gsl_ran_gaussian_ziggurat (const gsl_rng * R,
5717 -- Function: double gsl_ran_gaussian_ratio_method (const gsl_rng * R,
5719 This function computes a Gaussian random variate using the
5720 alternative Marsaglia-Tsang ziggurat and Kinderman-Monahan-Leva
5721 ratio methods. The Ziggurat algorithm is the fastest available
5722 algorithm in most cases.
5724 -- Function: double gsl_ran_ugaussian (const gsl_rng * R)
5725 -- Function: double gsl_ran_ugaussian_pdf (double X)
5726 -- Function: double gsl_ran_ugaussian_ratio_method (const gsl_rng * R)
5727 These functions compute results for the unit Gaussian
5728 distribution. They are equivalent to the functions above with a
5729 standard deviation of one, SIGMA = 1.
5731 -- Function: double gsl_cdf_gaussian_P (double X, double SIGMA)
5732 -- Function: double gsl_cdf_gaussian_Q (double X, double SIGMA)
5733 -- Function: double gsl_cdf_gaussian_Pinv (double P, double SIGMA)
5734 -- Function: double gsl_cdf_gaussian_Qinv (double Q, double SIGMA)
5735 These functions compute the cumulative distribution functions
5736 P(x), Q(x) and their inverses for the Gaussian distribution with
5737 standard deviation SIGMA.
5739 -- Function: double gsl_cdf_ugaussian_P (double X)
5740 -- Function: double gsl_cdf_ugaussian_Q (double X)
5741 -- Function: double gsl_cdf_ugaussian_Pinv (double P)
5742 -- Function: double gsl_cdf_ugaussian_Qinv (double Q)
5743 These functions compute the cumulative distribution functions
5744 P(x), Q(x) and their inverses for the unit Gaussian distribution.
5747 File: gsl-ref.info, Node: The Gaussian Tail Distribution, Next: The Bivariate Gaussian Distribution, Prev: The Gaussian Distribution, Up: Random Number Distributions
5749 19.3 The Gaussian Tail Distribution
5750 ===================================
5752 -- Function: double gsl_ran_gaussian_tail (const gsl_rng * R, double
5754 This function provides random variates from the upper tail of a
5755 Gaussian distribution with standard deviation SIGMA. The values
5756 returned are larger than the lower limit A, which must be
5757 positive. The method is based on Marsaglia's famous
5758 rectangle-wedge-tail algorithm (Ann. Math. Stat. 32, 894-899
5759 (1961)), with this aspect explained in Knuth, v2, 3rd ed, p139,586
5762 The probability distribution for Gaussian tail random variates is,
5764 p(x) dx = {1 \over N(a;\sigma) \sqrt{2 \pi \sigma^2}} \exp (- x^2/(2 \sigma^2)) dx
5766 for x > a where N(a;\sigma) is the normalization constant,
5768 N(a;\sigma) = (1/2) erfc(a / sqrt(2 sigma^2)).
5771 -- Function: double gsl_ran_gaussian_tail_pdf (double X, double A,
5773 This function computes the probability density p(x) at X for a
5774 Gaussian tail distribution with standard deviation SIGMA and lower
5775 limit A, using the formula given above.
5778 -- Function: double gsl_ran_ugaussian_tail (const gsl_rng * R, double
5780 -- Function: double gsl_ran_ugaussian_tail_pdf (double X, double A)
5781 These functions compute results for the tail of a unit Gaussian
5782 distribution. They are equivalent to the functions above with a
5783 standard deviation of one, SIGMA = 1.
5786 File: gsl-ref.info, Node: The Bivariate Gaussian Distribution, Next: The Exponential Distribution, Prev: The Gaussian Tail Distribution, Up: Random Number Distributions
5788 19.4 The Bivariate Gaussian Distribution
5789 ========================================
5791 -- Function: void gsl_ran_bivariate_gaussian (const gsl_rng * R,
5792 double SIGMA_X, double SIGMA_Y, double RHO, double * X,
5794 This function generates a pair of correlated Gaussian variates,
5795 with mean zero, correlation coefficient RHO and standard deviations
5796 SIGMA_X and SIGMA_Y in the x and y directions. The probability
5797 distribution for bivariate Gaussian random variates is,
5799 p(x,y) dx dy = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp (-(x^2/\sigma_x^2 + y^2/\sigma_y^2 - 2 \rho x y/(\sigma_x\sigma_y))/2(1-\rho^2)) dx dy
5801 for x,y in the range -\infty to +\infty. The correlation
5802 coefficient RHO should lie between 1 and -1.
5804 -- Function: double gsl_ran_bivariate_gaussian_pdf (double X, double
5805 Y, double SIGMA_X, double SIGMA_Y, double RHO)
5806 This function computes the probability density p(x,y) at (X,Y) for
5807 a bivariate Gaussian distribution with standard deviations
5808 SIGMA_X, SIGMA_Y and correlation coefficient RHO, using the
5809 formula given above.
5813 File: gsl-ref.info, Node: The Exponential Distribution, Next: The Laplace Distribution, Prev: The Bivariate Gaussian Distribution, Up: Random Number Distributions
5815 19.5 The Exponential Distribution
5816 =================================
5818 -- Function: double gsl_ran_exponential (const gsl_rng * R, double MU)
5819 This function returns a random variate from the exponential
5820 distribution with mean MU. The distribution is,
5822 p(x) dx = {1 \over \mu} \exp(-x/\mu) dx
5826 -- Function: double gsl_ran_exponential_pdf (double X, double MU)
5827 This function computes the probability density p(x) at X for an
5828 exponential distribution with mean MU, using the formula given
5832 -- Function: double gsl_cdf_exponential_P (double X, double MU)
5833 -- Function: double gsl_cdf_exponential_Q (double X, double MU)
5834 -- Function: double gsl_cdf_exponential_Pinv (double P, double MU)
5835 -- Function: double gsl_cdf_exponential_Qinv (double Q, double MU)
5836 These functions compute the cumulative distribution functions
5837 P(x), Q(x) and their inverses for the exponential distribution
5841 File: gsl-ref.info, Node: The Laplace Distribution, Next: The Exponential Power Distribution, Prev: The Exponential Distribution, Up: Random Number Distributions
5843 19.6 The Laplace Distribution
5844 =============================
5846 -- Function: double gsl_ran_laplace (const gsl_rng * R, double A)
5847 This function returns a random variate from the Laplace
5848 distribution with width A. The distribution is,
5850 p(x) dx = {1 \over 2 a} \exp(-|x/a|) dx
5852 for -\infty < x < \infty.
5854 -- Function: double gsl_ran_laplace_pdf (double X, double A)
5855 This function computes the probability density p(x) at X for a
5856 Laplace distribution with width A, using the formula given above.
5859 -- Function: double gsl_cdf_laplace_P (double X, double A)
5860 -- Function: double gsl_cdf_laplace_Q (double X, double A)
5861 -- Function: double gsl_cdf_laplace_Pinv (double P, double A)
5862 -- Function: double gsl_cdf_laplace_Qinv (double Q, double A)
5863 These functions compute the cumulative distribution functions
5864 P(x), Q(x) and their inverses for the Laplace distribution with
5868 File: gsl-ref.info, Node: The Exponential Power Distribution, Next: The Cauchy Distribution, Prev: The Laplace Distribution, Up: Random Number Distributions
5870 19.7 The Exponential Power Distribution
5871 =======================================
5873 -- Function: double gsl_ran_exppow (const gsl_rng * R, double A,
5875 This function returns a random variate from the exponential power
5876 distribution with scale parameter A and exponent B. The
5879 p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx
5881 for x >= 0. For b = 1 this reduces to the Laplace distribution.
5882 For b = 2 it has the same form as a gaussian distribution, but
5883 with a = \sqrt{2} \sigma.
5885 -- Function: double gsl_ran_exppow_pdf (double X, double A, double B)
5886 This function computes the probability density p(x) at X for an
5887 exponential power distribution with scale parameter A and exponent
5888 B, using the formula given above.
5891 -- Function: double gsl_cdf_exppow_P (double X, double A, double B)
5892 -- Function: double gsl_cdf_exppow_Q (double X, double A, double B)
5893 These functions compute the cumulative distribution functions
5894 P(x), Q(x) for the exponential power distribution with parameters
5898 File: gsl-ref.info, Node: The Cauchy Distribution, Next: The Rayleigh Distribution, Prev: The Exponential Power Distribution, Up: Random Number Distributions
5900 19.8 The Cauchy Distribution
5901 ============================
5903 -- Function: double gsl_ran_cauchy (const gsl_rng * R, double A)
5904 This function returns a random variate from the Cauchy
5905 distribution with scale parameter A. The probability distribution
5906 for Cauchy random variates is,
5908 p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx
5910 for x in the range -\infty to +\infty. The Cauchy distribution is
5911 also known as the Lorentz distribution.
5913 -- Function: double gsl_ran_cauchy_pdf (double X, double A)
5914 This function computes the probability density p(x) at X for a
5915 Cauchy distribution with scale parameter A, using the formula
5919 -- Function: double gsl_cdf_cauchy_P (double X, double A)
5920 -- Function: double gsl_cdf_cauchy_Q (double X, double A)
5921 -- Function: double gsl_cdf_cauchy_Pinv (double P, double A)
5922 -- Function: double gsl_cdf_cauchy_Qinv (double Q, double A)
5923 These functions compute the cumulative distribution functions
5924 P(x), Q(x) and their inverses for the Cauchy distribution with
5928 File: gsl-ref.info, Node: The Rayleigh Distribution, Next: The Rayleigh Tail Distribution, Prev: The Cauchy Distribution, Up: Random Number Distributions
5930 19.9 The Rayleigh Distribution
5931 ==============================
5933 -- Function: double gsl_ran_rayleigh (const gsl_rng * R, double SIGMA)
5934 This function returns a random variate from the Rayleigh
5935 distribution with scale parameter SIGMA. The distribution is,
5937 p(x) dx = {x \over \sigma^2} \exp(- x^2/(2 \sigma^2)) dx
5941 -- Function: double gsl_ran_rayleigh_pdf (double X, double SIGMA)
5942 This function computes the probability density p(x) at X for a
5943 Rayleigh distribution with scale parameter SIGMA, using the
5944 formula given above.
5947 -- Function: double gsl_cdf_rayleigh_P (double X, double SIGMA)
5948 -- Function: double gsl_cdf_rayleigh_Q (double X, double SIGMA)
5949 -- Function: double gsl_cdf_rayleigh_Pinv (double P, double SIGMA)
5950 -- Function: double gsl_cdf_rayleigh_Qinv (double Q, double SIGMA)
5951 These functions compute the cumulative distribution functions
5952 P(x), Q(x) and their inverses for the Rayleigh distribution with
5953 scale parameter SIGMA.
5956 File: gsl-ref.info, Node: The Rayleigh Tail Distribution, Next: The Landau Distribution, Prev: The Rayleigh Distribution, Up: Random Number Distributions
5958 19.10 The Rayleigh Tail Distribution
5959 ====================================
5961 -- Function: double gsl_ran_rayleigh_tail (const gsl_rng * R, double
5963 This function returns a random variate from the tail of the
5964 Rayleigh distribution with scale parameter SIGMA and a lower limit
5965 of A. The distribution is,
5967 p(x) dx = {x \over \sigma^2} \exp ((a^2 - x^2) /(2 \sigma^2)) dx
5971 -- Function: double gsl_ran_rayleigh_tail_pdf (double X, double A,
5973 This function computes the probability density p(x) at X for a
5974 Rayleigh tail distribution with scale parameter SIGMA and lower
5975 limit A, using the formula given above.
5979 File: gsl-ref.info, Node: The Landau Distribution, Next: The Levy alpha-Stable Distributions, Prev: The Rayleigh Tail Distribution, Up: Random Number Distributions
5981 19.11 The Landau Distribution
5982 =============================
5984 -- Function: double gsl_ran_landau (const gsl_rng * R)
5985 This function returns a random variate from the Landau
5986 distribution. The probability distribution for Landau random
5987 variates is defined analytically by the complex integral,
5989 p(x) = (1/(2 \pi i)) \int_{c-i\infty}^{c+i\infty} ds exp(s log(s) + x s)
5990 For numerical purposes it is more convenient to use the following
5991 equivalent form of the integral,
5993 p(x) = (1/\pi) \int_0^\infty dt \exp(-t \log(t) - x t) \sin(\pi t).
5995 -- Function: double gsl_ran_landau_pdf (double X)
5996 This function computes the probability density p(x) at X for the
5997 Landau distribution using an approximation to the formula given
6002 File: gsl-ref.info, Node: The Levy alpha-Stable Distributions, Next: The Levy skew alpha-Stable Distribution, Prev: The Landau Distribution, Up: Random Number Distributions
6004 19.12 The Levy alpha-Stable Distributions
6005 =========================================
6007 -- Function: double gsl_ran_levy (const gsl_rng * R, double C, double
6009 This function returns a random variate from the Levy symmetric
6010 stable distribution with scale C and exponent ALPHA. The symmetric
6011 stable probability distribution is defined by a fourier transform,
6013 p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^alpha)
6015 There is no explicit solution for the form of p(x) and the library
6016 does not define a corresponding `pdf' function. For \alpha = 1
6017 the distribution reduces to the Cauchy distribution. For \alpha =
6018 2 it is a Gaussian distribution with \sigma = \sqrt{2} c. For
6019 \alpha < 1 the tails of the distribution become extremely wide.
6021 The algorithm only works for 0 < alpha <= 2.
6025 File: gsl-ref.info, Node: The Levy skew alpha-Stable Distribution, Next: The Gamma Distribution, Prev: The Levy alpha-Stable Distributions, Up: Random Number Distributions
6027 19.13 The Levy skew alpha-Stable Distribution
6028 =============================================
6030 -- Function: double gsl_ran_levy_skew (const gsl_rng * R, double C,
6031 double ALPHA, double BETA)
6032 This function returns a random variate from the Levy skew stable
6033 distribution with scale C, exponent ALPHA and skewness parameter
6034 BETA. The skewness parameter must lie in the range [-1,1]. The
6035 Levy skew stable probability distribution is defined by a fourier
6038 p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^alpha (1-i beta sign(t) tan(pi alpha/2)))
6040 When \alpha = 1 the term \tan(\pi \alpha/2) is replaced by
6041 -(2/\pi)\log|t|. There is no explicit solution for the form of
6042 p(x) and the library does not define a corresponding `pdf'
6043 function. For \alpha = 2 the distribution reduces to a Gaussian
6044 distribution with \sigma = \sqrt{2} c and the skewness parameter
6045 has no effect. For \alpha < 1 the tails of the distribution
6046 become extremely wide. The symmetric distribution corresponds to
6049 The algorithm only works for 0 < alpha <= 2.
6051 The Levy alpha-stable distributions have the property that if N
6052 alpha-stable variates are drawn from the distribution p(c, \alpha,
6053 \beta) then the sum Y = X_1 + X_2 + \dots + X_N will also be
6054 distributed as an alpha-stable variate, p(N^(1/\alpha) c, \alpha,
6059 File: gsl-ref.info, Node: The Gamma Distribution, Next: The Flat (Uniform) Distribution, Prev: The Levy skew alpha-Stable Distribution, Up: Random Number Distributions
6061 19.14 The Gamma Distribution
6062 ============================
6064 -- Function: double gsl_ran_gamma (const gsl_rng * R, double A, double
6066 This function returns a random variate from the gamma
6067 distribution. The distribution function is,
6069 p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx
6073 The gamma distribution with an integer parameter A is known as the
6074 Erlang distribution.
6076 The variates are computed using the Marsaglia-Tsang fast gamma
6077 method. This function for this method was previously called
6078 `gsl_ran_gamma_mt' and can still be accessed using this name.
6080 -- Function: double gsl_ran_gamma_knuth (const gsl_rng * R, double A,
6082 This function returns a gamma variate using the algorithms from
6085 -- Function: double gsl_ran_gamma_pdf (double X, double A, double B)
6086 This function computes the probability density p(x) at X for a
6087 gamma distribution with parameters A and B, using the formula
6091 -- Function: double gsl_cdf_gamma_P (double X, double A, double B)
6092 -- Function: double gsl_cdf_gamma_Q (double X, double A, double B)
6093 -- Function: double gsl_cdf_gamma_Pinv (double P, double A, double B)
6094 -- Function: double gsl_cdf_gamma_Qinv (double Q, double A, double B)
6095 These functions compute the cumulative distribution functions
6096 P(x), Q(x) and their inverses for the gamma distribution with
6100 File: gsl-ref.info, Node: The Flat (Uniform) Distribution, Next: The Lognormal Distribution, Prev: The Gamma Distribution, Up: Random Number Distributions
6102 19.15 The Flat (Uniform) Distribution
6103 =====================================
6105 -- Function: double gsl_ran_flat (const gsl_rng * R, double A, double
6107 This function returns a random variate from the flat (uniform)
6108 distribution from A to B. The distribution is,
6110 p(x) dx = {1 \over (b-a)} dx
6112 if a <= x < b and 0 otherwise.
6114 -- Function: double gsl_ran_flat_pdf (double X, double A, double B)
6115 This function computes the probability density p(x) at X for a
6116 uniform distribution from A to B, using the formula given above.
6119 -- Function: double gsl_cdf_flat_P (double X, double A, double B)
6120 -- Function: double gsl_cdf_flat_Q (double X, double A, double B)
6121 -- Function: double gsl_cdf_flat_Pinv (double P, double A, double B)
6122 -- Function: double gsl_cdf_flat_Qinv (double Q, double A, double B)
6123 These functions compute the cumulative distribution functions
6124 P(x), Q(x) and their inverses for a uniform distribution from A to
6128 File: gsl-ref.info, Node: The Lognormal Distribution, Next: The Chi-squared Distribution, Prev: The Flat (Uniform) Distribution, Up: Random Number Distributions
6130 19.16 The Lognormal Distribution
6131 ================================
6133 -- Function: double gsl_ran_lognormal (const gsl_rng * R, double ZETA,
6135 This function returns a random variate from the lognormal
6136 distribution. The distribution function is,
6138 p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2} } \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
6142 -- Function: double gsl_ran_lognormal_pdf (double X, double ZETA,
6144 This function computes the probability density p(x) at X for a
6145 lognormal distribution with parameters ZETA and SIGMA, using the
6146 formula given above.
6149 -- Function: double gsl_cdf_lognormal_P (double X, double ZETA, double
6151 -- Function: double gsl_cdf_lognormal_Q (double X, double ZETA, double
6153 -- Function: double gsl_cdf_lognormal_Pinv (double P, double ZETA,
6155 -- Function: double gsl_cdf_lognormal_Qinv (double Q, double ZETA,
6157 These functions compute the cumulative distribution functions
6158 P(x), Q(x) and their inverses for the lognormal distribution with
6159 parameters ZETA and SIGMA.
6162 File: gsl-ref.info, Node: The Chi-squared Distribution, Next: The F-distribution, Prev: The Lognormal Distribution, Up: Random Number Distributions
6164 19.17 The Chi-squared Distribution
6165 ==================================
6167 The chi-squared distribution arises in statistics. If Y_i are n
6168 independent gaussian random variates with unit variance then the
6173 has a chi-squared distribution with n degrees of freedom.
6175 -- Function: double gsl_ran_chisq (const gsl_rng * R, double NU)
6176 This function returns a random variate from the chi-squared
6177 distribution with NU degrees of freedom. The distribution function
6180 p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx
6184 -- Function: double gsl_ran_chisq_pdf (double X, double NU)
6185 This function computes the probability density p(x) at X for a
6186 chi-squared distribution with NU degrees of freedom, using the
6187 formula given above.
6190 -- Function: double gsl_cdf_chisq_P (double X, double NU)
6191 -- Function: double gsl_cdf_chisq_Q (double X, double NU)
6192 -- Function: double gsl_cdf_chisq_Pinv (double P, double NU)
6193 -- Function: double gsl_cdf_chisq_Qinv (double Q, double NU)
6194 These functions compute the cumulative distribution functions
6195 P(x), Q(x) and their inverses for the chi-squared distribution
6196 with NU degrees of freedom.
6199 File: gsl-ref.info, Node: The F-distribution, Next: The t-distribution, Prev: The Chi-squared Distribution, Up: Random Number Distributions
6201 19.18 The F-distribution
6202 ========================
6204 The F-distribution arises in statistics. If Y_1 and Y_2 are
6205 chi-squared deviates with \nu_1 and \nu_2 degrees of freedom then the
6208 X = { (Y_1 / \nu_1) \over (Y_2 / \nu_2) }
6210 has an F-distribution F(x;\nu_1,\nu_2).
6212 -- Function: double gsl_ran_fdist (const gsl_rng * R, double NU1,
6214 This function returns a random variate from the F-distribution
6215 with degrees of freedom NU1 and NU2. The distribution function is,
6218 { \Gamma((\nu_1 + \nu_2)/2)
6219 \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) }
6220 \nu_1^{\nu_1/2} \nu_2^{\nu_2/2}
6221 x^{\nu_1/2 - 1} (\nu_2 + \nu_1 x)^{-\nu_1/2 -\nu_2/2}
6225 -- Function: double gsl_ran_fdist_pdf (double X, double NU1, double
6227 This function computes the probability density p(x) at X for an
6228 F-distribution with NU1 and NU2 degrees of freedom, using the
6229 formula given above.
6232 -- Function: double gsl_cdf_fdist_P (double X, double NU1, double NU2)
6233 -- Function: double gsl_cdf_fdist_Q (double X, double NU1, double NU2)
6234 -- Function: double gsl_cdf_fdist_Pinv (double P, double NU1, double
6236 -- Function: double gsl_cdf_fdist_Qinv (double Q, double NU1, double
6238 These functions compute the cumulative distribution functions
6239 P(x), Q(x) and their inverses for the F-distribution with NU1 and
6240 NU2 degrees of freedom.
6243 File: gsl-ref.info, Node: The t-distribution, Next: The Beta Distribution, Prev: The F-distribution, Up: Random Number Distributions
6245 19.19 The t-distribution
6246 ========================
6248 The t-distribution arises in statistics. If Y_1 has a normal
6249 distribution and Y_2 has a chi-squared distribution with \nu degrees of
6250 freedom then the ratio,
6252 X = { Y_1 \over \sqrt{Y_2 / \nu} }
6254 has a t-distribution t(x;\nu) with \nu degrees of freedom.
6256 -- Function: double gsl_ran_tdist (const gsl_rng * R, double NU)
6257 This function returns a random variate from the t-distribution.
6258 The distribution function is,
6260 p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)}
6261 (1 + x^2/\nu)^{-(\nu + 1)/2} dx
6263 for -\infty < x < +\infty.
6265 -- Function: double gsl_ran_tdist_pdf (double X, double NU)
6266 This function computes the probability density p(x) at X for a
6267 t-distribution with NU degrees of freedom, using the formula given
6271 -- Function: double gsl_cdf_tdist_P (double X, double NU)
6272 -- Function: double gsl_cdf_tdist_Q (double X, double NU)
6273 -- Function: double gsl_cdf_tdist_Pinv (double P, double NU)
6274 -- Function: double gsl_cdf_tdist_Qinv (double Q, double NU)
6275 These functions compute the cumulative distribution functions
6276 P(x), Q(x) and their inverses for the t-distribution with NU
6280 File: gsl-ref.info, Node: The Beta Distribution, Next: The Logistic Distribution, Prev: The t-distribution, Up: Random Number Distributions
6282 19.20 The Beta Distribution
6283 ===========================
6285 -- Function: double gsl_ran_beta (const gsl_rng * R, double A, double
6287 This function returns a random variate from the beta distribution.
6288 The distribution function is,
6290 p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx
6294 -- Function: double gsl_ran_beta_pdf (double X, double A, double B)
6295 This function computes the probability density p(x) at X for a
6296 beta distribution with parameters A and B, using the formula given
6300 -- Function: double gsl_cdf_beta_P (double X, double A, double B)
6301 -- Function: double gsl_cdf_beta_Q (double X, double A, double B)
6302 -- Function: double gsl_cdf_beta_Pinv (double P, double A, double B)
6303 -- Function: double gsl_cdf_beta_Qinv (double Q, double A, double B)
6304 These functions compute the cumulative distribution functions
6305 P(x), Q(x) and their inverses for the beta distribution with
6309 File: gsl-ref.info, Node: The Logistic Distribution, Next: The Pareto Distribution, Prev: The Beta Distribution, Up: Random Number Distributions
6311 19.21 The Logistic Distribution
6312 ===============================
6314 -- Function: double gsl_ran_logistic (const gsl_rng * R, double A)
6315 This function returns a random variate from the logistic
6316 distribution. The distribution function is,
6318 p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx
6320 for -\infty < x < +\infty.
6322 -- Function: double gsl_ran_logistic_pdf (double X, double A)
6323 This function computes the probability density p(x) at X for a
6324 logistic distribution with scale parameter A, using the formula
6328 -- Function: double gsl_cdf_logistic_P (double X, double A)
6329 -- Function: double gsl_cdf_logistic_Q (double X, double A)
6330 -- Function: double gsl_cdf_logistic_Pinv (double P, double A)
6331 -- Function: double gsl_cdf_logistic_Qinv (double Q, double A)
6332 These functions compute the cumulative distribution functions
6333 P(x), Q(x) and their inverses for the logistic distribution with
6337 File: gsl-ref.info, Node: The Pareto Distribution, Next: Spherical Vector Distributions, Prev: The Logistic Distribution, Up: Random Number Distributions
6339 19.22 The Pareto Distribution
6340 =============================
6342 -- Function: double gsl_ran_pareto (const gsl_rng * R, double A,
6344 This function returns a random variate from the Pareto
6345 distribution of order A. The distribution function is,
6347 p(x) dx = (a/b) / (x/b)^{a+1} dx
6351 -- Function: double gsl_ran_pareto_pdf (double X, double A, double B)
6352 This function computes the probability density p(x) at X for a
6353 Pareto distribution with exponent A and scale B, using the formula
6357 -- Function: double gsl_cdf_pareto_P (double X, double A, double B)
6358 -- Function: double gsl_cdf_pareto_Q (double X, double A, double B)
6359 -- Function: double gsl_cdf_pareto_Pinv (double P, double A, double B)
6360 -- Function: double gsl_cdf_pareto_Qinv (double Q, double A, double B)
6361 These functions compute the cumulative distribution functions
6362 P(x), Q(x) and their inverses for the Pareto distribution with
6363 exponent A and scale B.
6366 File: gsl-ref.info, Node: Spherical Vector Distributions, Next: The Weibull Distribution, Prev: The Pareto Distribution, Up: Random Number Distributions
6368 19.23 Spherical Vector Distributions
6369 ====================================
6371 The spherical distributions generate random vectors, located on a
6372 spherical surface. They can be used as random directions, for example
6373 in the steps of a random walk.
6375 -- Function: void gsl_ran_dir_2d (const gsl_rng * R, double * X,
6377 -- Function: void gsl_ran_dir_2d_trig_method (const gsl_rng * R,
6378 double * X, double * Y)
6379 This function returns a random direction vector v = (X,Y) in two
6380 dimensions. The vector is normalized such that |v|^2 = x^2 + y^2
6381 = 1. The obvious way to do this is to take a uniform random
6382 number between 0 and 2\pi and let X and Y be the sine and cosine
6383 respectively. Two trig functions would have been expensive in the
6384 old days, but with modern hardware implementations, this is
6385 sometimes the fastest way to go. This is the case for the Pentium
6386 (but not the case for the Sun Sparcstation). One can avoid the
6387 trig evaluations by choosing X and Y in the interior of a unit
6388 circle (choose them at random from the interior of the enclosing
6389 square, and then reject those that are outside the unit circle),
6390 and then dividing by \sqrt{x^2 + y^2}. A much cleverer approach,
6391 attributed to von Neumann (See Knuth, v2, 3rd ed, p140, exercise
6392 23), requires neither trig nor a square root. In this approach, U
6393 and V are chosen at random from the interior of a unit circle, and
6394 then x=(u^2-v^2)/(u^2+v^2) and y=2uv/(u^2+v^2).
6396 -- Function: void gsl_ran_dir_3d (const gsl_rng * R, double * X,
6397 double * Y, double * Z)
6398 This function returns a random direction vector v = (X,Y,Z) in
6399 three dimensions. The vector is normalized such that |v|^2 = x^2
6400 + y^2 + z^2 = 1. The method employed is due to Robert E. Knop
6401 (CACM 13, 326 (1970)), and explained in Knuth, v2, 3rd ed, p136.
6402 It uses the surprising fact that the distribution projected along
6403 any axis is actually uniform (this is only true for 3 dimensions).
6405 -- Function: void gsl_ran_dir_nd (const gsl_rng * R, size_t N, double
6407 This function returns a random direction vector v =
6408 (x_1,x_2,...,x_n) in N dimensions. The vector is normalized such
6409 that |v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1. The method uses the
6410 fact that a multivariate gaussian distribution is spherically
6411 symmetric. Each component is generated to have a gaussian
6412 distribution, and then the components are normalized. The method
6413 is described by Knuth, v2, 3rd ed, p135-136, and attributed to G.
6414 W. Brown, Modern Mathematics for the Engineer (1956).
6417 File: gsl-ref.info, Node: The Weibull Distribution, Next: The Type-1 Gumbel Distribution, Prev: Spherical Vector Distributions, Up: Random Number Distributions
6419 19.24 The Weibull Distribution
6420 ==============================
6422 -- Function: double gsl_ran_weibull (const gsl_rng * R, double A,
6424 This function returns a random variate from the Weibull
6425 distribution. The distribution function is,
6427 p(x) dx = {b \over a^b} x^{b-1} \exp(-(x/a)^b) dx
6431 -- Function: double gsl_ran_weibull_pdf (double X, double A, double B)
6432 This function computes the probability density p(x) at X for a
6433 Weibull distribution with scale A and exponent B, using the
6434 formula given above.
6437 -- Function: double gsl_cdf_weibull_P (double X, double A, double B)
6438 -- Function: double gsl_cdf_weibull_Q (double X, double A, double B)
6439 -- Function: double gsl_cdf_weibull_Pinv (double P, double A, double B)
6440 -- Function: double gsl_cdf_weibull_Qinv (double Q, double A, double B)
6441 These functions compute the cumulative distribution functions
6442 P(x), Q(x) and their inverses for the Weibull distribution with
6443 scale A and exponent B.
6446 File: gsl-ref.info, Node: The Type-1 Gumbel Distribution, Next: The Type-2 Gumbel Distribution, Prev: The Weibull Distribution, Up: Random Number Distributions
6448 19.25 The Type-1 Gumbel Distribution
6449 ====================================
6451 -- Function: double gsl_ran_gumbel1 (const gsl_rng * R, double A,
6453 This function returns a random variate from the Type-1 Gumbel
6454 distribution. The Type-1 Gumbel distribution function is,
6456 p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
6458 for -\infty < x < \infty.
6460 -- Function: double gsl_ran_gumbel1_pdf (double X, double A, double B)
6461 This function computes the probability density p(x) at X for a
6462 Type-1 Gumbel distribution with parameters A and B, using the
6463 formula given above.
6466 -- Function: double gsl_cdf_gumbel1_P (double X, double A, double B)
6467 -- Function: double gsl_cdf_gumbel1_Q (double X, double A, double B)
6468 -- Function: double gsl_cdf_gumbel1_Pinv (double P, double A, double B)
6469 -- Function: double gsl_cdf_gumbel1_Qinv (double Q, double A, double B)
6470 These functions compute the cumulative distribution functions
6471 P(x), Q(x) and their inverses for the Type-1 Gumbel distribution
6472 with parameters A and B.
6475 File: gsl-ref.info, Node: The Type-2 Gumbel Distribution, Next: The Dirichlet Distribution, Prev: The Type-1 Gumbel Distribution, Up: Random Number Distributions
6477 19.26 The Type-2 Gumbel Distribution
6478 ====================================
6480 -- Function: double gsl_ran_gumbel2 (const gsl_rng * R, double A,
6482 This function returns a random variate from the Type-2 Gumbel
6483 distribution. The Type-2 Gumbel distribution function is,
6485 p(x) dx = a b x^{-a-1} \exp(-b x^{-a}) dx
6489 -- Function: double gsl_ran_gumbel2_pdf (double X, double A, double B)
6490 This function computes the probability density p(x) at X for a
6491 Type-2 Gumbel distribution with parameters A and B, using the
6492 formula given above.
6495 -- Function: double gsl_cdf_gumbel2_P (double X, double A, double B)
6496 -- Function: double gsl_cdf_gumbel2_Q (double X, double A, double B)
6497 -- Function: double gsl_cdf_gumbel2_Pinv (double P, double A, double B)
6498 -- Function: double gsl_cdf_gumbel2_Qinv (double Q, double A, double B)
6499 These functions compute the cumulative distribution functions
6500 P(x), Q(x) and their inverses for the Type-2 Gumbel distribution
6501 with parameters A and B.
6504 File: gsl-ref.info, Node: The Dirichlet Distribution, Next: General Discrete Distributions, Prev: The Type-2 Gumbel Distribution, Up: Random Number Distributions
6506 19.27 The Dirichlet Distribution
6507 ================================
6509 -- Function: void gsl_ran_dirichlet (const gsl_rng * R, size_t K,
6510 const double ALPHA[], double THETA[])
6511 This function returns an array of K random variates from a
6512 Dirichlet distribution of order K-1. The distribution function is
6514 p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K =
6515 (1/Z) \prod_{i=1}^K \theta_i^{\alpha_i - 1} \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 ... d\theta_K
6517 for theta_i >= 0 and alpha_i >= 0. The delta function ensures
6518 that \sum \theta_i = 1. The normalization factor Z is
6520 Z = {\prod_{i=1}^K \Gamma(\alpha_i)} / {\Gamma( \sum_{i=1}^K \alpha_i)}
6522 The random variates are generated by sampling K values from gamma
6523 distributions with parameters a=alpha_i, b=1, and renormalizing.
6524 See A.M. Law, W.D. Kelton, `Simulation Modeling and Analysis'
6527 -- Function: double gsl_ran_dirichlet_pdf (size_t K, const double
6528 ALPHA[], const double THETA[])
6529 This function computes the probability density p(\theta_1, ... ,
6530 \theta_K) at THETA[K] for a Dirichlet distribution with parameters
6531 ALPHA[K], using the formula given above.
6533 -- Function: double gsl_ran_dirichlet_lnpdf (size_t K, const double
6534 ALPHA[], const double THETA[])
6535 This function computes the logarithm of the probability density
6536 p(\theta_1, ... , \theta_K) for a Dirichlet distribution with
6537 parameters ALPHA[K].
6540 File: gsl-ref.info, Node: General Discrete Distributions, Next: The Poisson Distribution, Prev: The Dirichlet Distribution, Up: Random Number Distributions
6542 19.28 General Discrete Distributions
6543 ====================================
6545 Given K discrete events with different probabilities P[k], produce a
6546 random value k consistent with its probability.
6548 The obvious way to do this is to preprocess the probability list by
6549 generating a cumulative probability array with K+1 elements:
6554 Note that this construction produces C[K]=1. Now choose a uniform
6555 deviate u between 0 and 1, and find the value of k such that C[k] <= u
6556 < C[k+1]. Although this in principle requires of order \log K steps per
6557 random number generation, they are fast steps, and if you use something
6558 like \lfloor uK \rfloor as a starting point, you can often do pretty
6561 But faster methods have been devised. Again, the idea is to
6562 preprocess the probability list, and save the result in some form of
6563 lookup table; then the individual calls for a random discrete event can
6564 go rapidly. An approach invented by G. Marsaglia (Generating discrete
6565 random numbers in a computer, Comm ACM 6, 37-38 (1963)) is very clever,
6566 and readers interested in examples of good algorithm design are
6567 directed to this short and well-written paper. Unfortunately, for
6568 large K, Marsaglia's lookup table can be quite large.
6570 A much better approach is due to Alastair J. Walker (An efficient
6571 method for generating discrete random variables with general
6572 distributions, ACM Trans on Mathematical Software 3, 253-256 (1977);
6573 see also Knuth, v2, 3rd ed, p120-121,139). This requires two lookup
6574 tables, one floating point and one integer, but both only of size K.
6575 After preprocessing, the random numbers are generated in O(1) time,
6576 even for large K. The preprocessing suggested by Walker requires
6577 O(K^2) effort, but that is not actually necessary, and the
6578 implementation provided here only takes O(K) effort. In general, more
6579 preprocessing leads to faster generation of the individual random
6580 numbers, but a diminishing return is reached pretty early. Knuth points
6581 out that the optimal preprocessing is combinatorially difficult for
6584 This method can be used to speed up some of the discrete random
6585 number generators below, such as the binomial distribution. To use it
6586 for something like the Poisson Distribution, a modification would have
6587 to be made, since it only takes a finite set of K outcomes.
6589 -- Function: gsl_ran_discrete_t * gsl_ran_discrete_preproc (size_t K,
6591 This function returns a pointer to a structure that contains the
6592 lookup table for the discrete random number generator. The array
6593 P[] contains the probabilities of the discrete events; these array
6594 elements must all be positive, but they needn't add up to one (so
6595 you can think of them more generally as "weights")--the
6596 preprocessor will normalize appropriately. This return value is
6597 used as an argument for the `gsl_ran_discrete' function below.
6599 -- Function: size_t gsl_ran_discrete (const gsl_rng * R, const
6600 gsl_ran_discrete_t * G)
6601 After the preprocessor, above, has been called, you use this
6602 function to get the discrete random numbers.
6604 -- Function: double gsl_ran_discrete_pdf (size_t K, const
6605 gsl_ran_discrete_t * G)
6606 Returns the probability P[k] of observing the variable K. Since
6607 P[k] is not stored as part of the lookup table, it must be
6608 recomputed; this computation takes O(K), so if K is large and you
6609 care about the original array P[k] used to create the lookup
6610 table, then you should just keep this original array P[k] around.
6612 -- Function: void gsl_ran_discrete_free (gsl_ran_discrete_t * G)
6613 De-allocates the lookup table pointed to by G.
6616 File: gsl-ref.info, Node: The Poisson Distribution, Next: The Bernoulli Distribution, Prev: General Discrete Distributions, Up: Random Number Distributions
6618 19.29 The Poisson Distribution
6619 ==============================
6621 -- Function: unsigned int gsl_ran_poisson (const gsl_rng * R, double
6623 This function returns a random integer from the Poisson
6624 distribution with mean MU. The probability distribution for
6625 Poisson variates is,
6627 p(k) = {\mu^k \over k!} \exp(-\mu)
6631 -- Function: double gsl_ran_poisson_pdf (unsigned int K, double MU)
6632 This function computes the probability p(k) of obtaining K from a
6633 Poisson distribution with mean MU, using the formula given above.
6636 -- Function: double gsl_cdf_poisson_P (unsigned int K, double MU)
6637 -- Function: double gsl_cdf_poisson_Q (unsigned int K, double MU)
6638 These functions compute the cumulative distribution functions
6639 P(k), Q(k) for the Poisson distribution with parameter MU.
6642 File: gsl-ref.info, Node: The Bernoulli Distribution, Next: The Binomial Distribution, Prev: The Poisson Distribution, Up: Random Number Distributions
6644 19.30 The Bernoulli Distribution
6645 ================================
6647 -- Function: unsigned int gsl_ran_bernoulli (const gsl_rng * R, double
6649 This function returns either 0 or 1, the result of a Bernoulli
6650 trial with probability P. The probability distribution for a
6657 -- Function: double gsl_ran_bernoulli_pdf (unsigned int K, double P)
6658 This function computes the probability p(k) of obtaining K from a
6659 Bernoulli distribution with probability parameter P, using the
6660 formula given above.
6664 File: gsl-ref.info, Node: The Binomial Distribution, Next: The Multinomial Distribution, Prev: The Bernoulli Distribution, Up: Random Number Distributions
6666 19.31 The Binomial Distribution
6667 ===============================
6669 -- Function: unsigned int gsl_ran_binomial (const gsl_rng * R, double
6671 This function returns a random integer from the binomial
6672 distribution, the number of successes in N independent trials with
6673 probability P. The probability distribution for binomial variates
6676 p(k) = {n! \over k! (n-k)! } p^k (1-p)^{n-k}
6680 -- Function: double gsl_ran_binomial_pdf (unsigned int K, double P,
6682 This function computes the probability p(k) of obtaining K from a
6683 binomial distribution with parameters P and N, using the formula
6687 -- Function: double gsl_cdf_binomial_P (unsigned int K, double P,
6689 -- Function: double gsl_cdf_binomial_Q (unsigned int K, double P,
6691 These functions compute the cumulative distribution functions
6692 P(k), Q(k) for the binomial distribution with parameters P and N.
6695 File: gsl-ref.info, Node: The Multinomial Distribution, Next: The Negative Binomial Distribution, Prev: The Binomial Distribution, Up: Random Number Distributions
6697 19.32 The Multinomial Distribution
6698 ==================================
6700 -- Function: void gsl_ran_multinomial (const gsl_rng * R, size_t K,
6701 unsigned int N, const double P[], unsigned int N[])
6702 This function computes a random sample N[] from the multinomial
6703 distribution formed by N trials from an underlying distribution
6704 P[K]. The distribution function for N[] is,
6706 P(n_1, n_2, ..., n_K) =
6707 (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K
6709 where (n_1, n_2, ..., n_K) are nonnegative integers with
6710 sum_{k=1}^K n_k = N, and (p_1, p_2, ..., p_K) is a probability
6711 distribution with \sum p_i = 1. If the array P[K] is not
6712 normalized then its entries will be treated as weights and
6713 normalized appropriately. The arrays N[] and P[] must both be of
6716 Random variates are generated using the conditional binomial
6717 method (see C.S. David, `The computer generation of multinomial
6718 random variates', Comp. Stat. Data Anal. 16 (1993) 205-217 for
6721 -- Function: double gsl_ran_multinomial_pdf (size_t K, const double
6722 P[], const unsigned int N[])
6723 This function computes the probability P(n_1, n_2, ..., n_K) of
6724 sampling N[K] from a multinomial distribution with parameters
6725 P[K], using the formula given above.
6727 -- Function: double gsl_ran_multinomial_lnpdf (size_t K, const double
6728 P[], const unsigned int N[])
6729 This function returns the logarithm of the probability for the
6730 multinomial distribution P(n_1, n_2, ..., n_K) with parameters
6734 File: gsl-ref.info, Node: The Negative Binomial Distribution, Next: The Pascal Distribution, Prev: The Multinomial Distribution, Up: Random Number Distributions
6736 19.33 The Negative Binomial Distribution
6737 ========================================
6739 -- Function: unsigned int gsl_ran_negative_binomial (const gsl_rng *
6740 R, double P, double N)
6741 This function returns a random integer from the negative binomial
6742 distribution, the number of failures occurring before N successes
6743 in independent trials with probability P of success. The
6744 probability distribution for negative binomial variates is,
6746 p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
6748 Note that n is not required to be an integer.
6750 -- Function: double gsl_ran_negative_binomial_pdf (unsigned int K,
6752 This function computes the probability p(k) of obtaining K from a
6753 negative binomial distribution with parameters P and N, using the
6754 formula given above.
6757 -- Function: double gsl_cdf_negative_binomial_P (unsigned int K,
6759 -- Function: double gsl_cdf_negative_binomial_Q (unsigned int K,
6761 These functions compute the cumulative distribution functions
6762 P(k), Q(k) for the negative binomial distribution with parameters
6766 File: gsl-ref.info, Node: The Pascal Distribution, Next: The Geometric Distribution, Prev: The Negative Binomial Distribution, Up: Random Number Distributions
6768 19.34 The Pascal Distribution
6769 =============================
6771 -- Function: unsigned int gsl_ran_pascal (const gsl_rng * R, double P,
6773 This function returns a random integer from the Pascal
6774 distribution. The Pascal distribution is simply a negative
6775 binomial distribution with an integer value of n.
6777 p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k
6781 -- Function: double gsl_ran_pascal_pdf (unsigned int K, double P,
6783 This function computes the probability p(k) of obtaining K from a
6784 Pascal distribution with parameters P and N, using the formula
6788 -- Function: double gsl_cdf_pascal_P (unsigned int K, double P,
6790 -- Function: double gsl_cdf_pascal_Q (unsigned int K, double P,
6792 These functions compute the cumulative distribution functions
6793 P(k), Q(k) for the Pascal distribution with parameters P and N.
6796 File: gsl-ref.info, Node: The Geometric Distribution, Next: The Hypergeometric Distribution, Prev: The Pascal Distribution, Up: Random Number Distributions
6798 19.35 The Geometric Distribution
6799 ================================
6801 -- Function: unsigned int gsl_ran_geometric (const gsl_rng * R, double
6803 This function returns a random integer from the geometric
6804 distribution, the number of independent trials with probability P
6805 until the first success. The probability distribution for
6806 geometric variates is,
6808 p(k) = p (1-p)^(k-1)
6810 for k >= 1. Note that the distribution begins with k=1 with this
6811 definition. There is another convention in which the exponent k-1
6814 -- Function: double gsl_ran_geometric_pdf (unsigned int K, double P)
6815 This function computes the probability p(k) of obtaining K from a
6816 geometric distribution with probability parameter P, using the
6817 formula given above.
6820 -- Function: double gsl_cdf_geometric_P (unsigned int K, double P)
6821 -- Function: double gsl_cdf_geometric_Q (unsigned int K, double P)
6822 These functions compute the cumulative distribution functions
6823 P(k), Q(k) for the geometric distribution with parameter P.
6826 File: gsl-ref.info, Node: The Hypergeometric Distribution, Next: The Logarithmic Distribution, Prev: The Geometric Distribution, Up: Random Number Distributions
6828 19.36 The Hypergeometric Distribution
6829 =====================================
6831 -- Function: unsigned int gsl_ran_hypergeometric (const gsl_rng * R,
6832 unsigned int N1, unsigned int N2, unsigned int T)
6833 This function returns a random integer from the hypergeometric
6834 distribution. The probability distribution for hypergeometric
6837 p(k) = C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
6839 where C(a,b) = a!/(b!(a-b)!) and t <= n_1 + n_2. The domain of k
6840 is max(0,t-n_2), ..., min(t,n_1).
6842 If a population contains n_1 elements of "type 1" and n_2 elements
6843 of "type 2" then the hypergeometric distribution gives the
6844 probability of obtaining k elements of "type 1" in t samples from
6845 the population without replacement.
6847 -- Function: double gsl_ran_hypergeometric_pdf (unsigned int K,
6848 unsigned int N1, unsigned int N2, unsigned int T)
6849 This function computes the probability p(k) of obtaining K from a
6850 hypergeometric distribution with parameters N1, N2, T, using the
6851 formula given above.
6854 -- Function: double gsl_cdf_hypergeometric_P (unsigned int K, unsigned
6855 int N1, unsigned int N2, unsigned int T)
6856 -- Function: double gsl_cdf_hypergeometric_Q (unsigned int K, unsigned
6857 int N1, unsigned int N2, unsigned int T)
6858 These functions compute the cumulative distribution functions
6859 P(k), Q(k) for the hypergeometric distribution with parameters N1,
6863 File: gsl-ref.info, Node: The Logarithmic Distribution, Next: Shuffling and Sampling, Prev: The Hypergeometric Distribution, Up: Random Number Distributions
6865 19.37 The Logarithmic Distribution
6866 ==================================
6868 -- Function: unsigned int gsl_ran_logarithmic (const gsl_rng * R,
6870 This function returns a random integer from the logarithmic
6871 distribution. The probability distribution for logarithmic random
6874 p(k) = {-1 \over \log(1-p)} {(p^k \over k)}
6878 -- Function: double gsl_ran_logarithmic_pdf (unsigned int K, double P)
6879 This function computes the probability p(k) of obtaining K from a
6880 logarithmic distribution with probability parameter P, using the
6881 formula given above.
6885 File: gsl-ref.info, Node: Shuffling and Sampling, Next: Random Number Distribution Examples, Prev: The Logarithmic Distribution, Up: Random Number Distributions
6887 19.38 Shuffling and Sampling
6888 ============================
6890 The following functions allow the shuffling and sampling of a set of
6891 objects. The algorithms rely on a random number generator as a source
6892 of randomness and a poor quality generator can lead to correlations in
6893 the output. In particular it is important to avoid generators with a
6894 short period. For more information see Knuth, v2, 3rd ed, Section
6895 3.4.2, "Random Sampling and Shuffling".
6897 -- Function: void gsl_ran_shuffle (const gsl_rng * R, void * BASE,
6898 size_t N, size_t SIZE)
6899 This function randomly shuffles the order of N objects, each of
6900 size SIZE, stored in the array BASE[0..N-1]. The output of the
6901 random number generator R is used to produce the permutation. The
6902 algorithm generates all possible n! permutations with equal
6903 probability, assuming a perfect source of random numbers.
6905 The following code shows how to shuffle the numbers from 0 to 51,
6909 for (i = 0; i < 52; i++)
6914 gsl_ran_shuffle (r, a, 52, sizeof (int));
6917 -- Function: int gsl_ran_choose (const gsl_rng * R, void * DEST,
6918 size_t K, void * SRC, size_t N, size_t SIZE)
6919 This function fills the array DEST[k] with K objects taken
6920 randomly from the N elements of the array SRC[0..N-1]. The
6921 objects are each of size SIZE. The output of the random number
6922 generator R is used to make the selection. The algorithm ensures
6923 all possible samples are equally likely, assuming a perfect source
6926 The objects are sampled _without_ replacement, thus each object can
6927 only appear once in DEST[k]. It is required that K be less than
6928 or equal to `n'. The objects in DEST will be in the same relative
6929 order as those in SRC. You will need to call `gsl_ran_shuffle(r,
6930 dest, n, size)' if you want to randomize the order.
6932 The following code shows how to select a random sample of three
6933 unique numbers from the set 0 to 99,
6935 double a[3], b[100];
6937 for (i = 0; i < 100; i++)
6942 gsl_ran_choose (r, a, 3, b, 100, sizeof (double));
6945 -- Function: void gsl_ran_sample (const gsl_rng * R, void * DEST,
6946 size_t K, void * SRC, size_t N, size_t SIZE)
6947 This function is like `gsl_ran_choose' but samples K items from
6948 the original array of N items SRC with replacement, so the same
6949 object can appear more than once in the output sequence DEST.
6950 There is no requirement that K be less than N in this case.