diff options
Diffstat (limited to 'gsl-1.9/doc/vectors.texi')
| -rw-r--r-- | gsl-1.9/doc/vectors.texi | 1616 |
1 files changed, 1616 insertions, 0 deletions
diff --git a/gsl-1.9/doc/vectors.texi b/gsl-1.9/doc/vectors.texi new file mode 100644 index 0000000..e17a1af --- /dev/null +++ b/gsl-1.9/doc/vectors.texi @@ -0,0 +1,1616 @@ +@cindex blocks +@cindex vectors +@cindex matrices + +The functions described in this chapter provide a simple vector and +matrix interface to ordinary C arrays. The memory management of these +arrays is implemented using a single underlying type, known as a +block. By writing your functions in terms of vectors and matrices you +can pass a single structure containing both data and dimensions as an +argument without needing additional function parameters. The structures +are compatible with the vector and matrix formats used by @sc{blas} +routines. + +@menu +* Data types:: +* Blocks:: +* Vectors:: +* Matrices:: +* Vector and Matrix References and Further Reading:: +@end menu + +@node Data types +@section Data types + +All the functions are available for each of the standard data-types. +The versions for @code{double} have the prefix @code{gsl_block}, +@code{gsl_vector} and @code{gsl_matrix}. Similarly the versions for +single-precision @code{float} arrays have the prefix +@code{gsl_block_float}, @code{gsl_vector_float} and +@code{gsl_matrix_float}. The full list of available types is given +below, + +@example +gsl_block double +gsl_block_float float +gsl_block_long_double long double +gsl_block_int int +gsl_block_uint unsigned int +gsl_block_long long +gsl_block_ulong unsigned long +gsl_block_short short +gsl_block_ushort unsigned short +gsl_block_char char +gsl_block_uchar unsigned char +gsl_block_complex complex double +gsl_block_complex_float complex float +gsl_block_complex_long_double complex long double +@end example + +@noindent +Corresponding types exist for the @code{gsl_vector} and +@code{gsl_matrix} functions. + + + +@node Blocks +@section Blocks + +For consistency all memory is allocated through a @code{gsl_block} +structure. The structure contains two components, the size of an area of +memory and a pointer to the memory. The @code{gsl_block} structure looks +like this, + +@example +typedef struct +@{ + size_t size; + double * data; +@} gsl_block; +@end example +@comment + +@noindent +Vectors and matrices are made by @dfn{slicing} an underlying block. A +slice is a set of elements formed from an initial offset and a +combination of indices and step-sizes. In the case of a matrix the +step-size for the column index represents the row-length. The step-size +for a vector is known as the @dfn{stride}. + +The functions for allocating and deallocating blocks are defined in +@file{gsl_block.h} + +@menu +* Block allocation:: +* Reading and writing blocks:: +* Example programs for blocks:: +@end menu + +@node Block allocation +@subsection Block allocation + +The functions for allocating memory to a block follow the style of +@code{malloc} and @code{free}. In addition they also perform their own +error checking. If there is insufficient memory available to allocate a +block then the functions call the GSL error handler (with an error +number of @code{GSL_ENOMEM}) in addition to returning a null +pointer. Thus if you use the library error handler to abort your program +then it isn't necessary to check every @code{alloc}. + +@deftypefun {gsl_block *} gsl_block_alloc (size_t @var{n}) +This function allocates memory for a block of @var{n} double-precision +elements, returning a pointer to the block struct. The block is not +initialized and so the values of its elements are undefined. Use the +function @code{gsl_block_calloc} if you want to ensure that all the +elements are initialized to zero. + +A null pointer is returned if insufficient memory is available to create +the block. +@end deftypefun + +@deftypefun {gsl_block *} gsl_block_calloc (size_t @var{n}) +This function allocates memory for a block and initializes all the +elements of the block to zero. +@end deftypefun + +@deftypefun void gsl_block_free (gsl_block * @var{b}) +This function frees the memory used by a block @var{b} previously +allocated with @code{gsl_block_alloc} or @code{gsl_block_calloc}. The +block @var{b} must be a valid block object (a null pointer is not +allowed). +@end deftypefun + +@node Reading and writing blocks +@subsection Reading and writing blocks + +The library provides functions for reading and writing blocks to a file +as binary data or formatted text. + +@deftypefun int gsl_block_fwrite (FILE * @var{stream}, const gsl_block * @var{b}) +This function writes the elements of the block @var{b} to the stream +@var{stream} in binary format. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem writing to the file. Since the +data is written in the native binary format it may not be portable +between different architectures. +@end deftypefun + +@deftypefun int gsl_block_fread (FILE * @var{stream}, gsl_block * @var{b}) +This function reads into the block @var{b} from the open stream +@var{stream} in binary format. The block @var{b} must be preallocated +with the correct length since the function uses the size of @var{b} to +determine how many bytes to read. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. The +data is assumed to have been written in the native binary format on the +same architecture. +@end deftypefun + +@deftypefun int gsl_block_fprintf (FILE * @var{stream}, const gsl_block * @var{b}, const char * @var{format}) +This function writes the elements of the block @var{b} line-by-line to +the stream @var{stream} using the format specifier @var{format}, which +should be one of the @code{%g}, @code{%e} or @code{%f} formats for +floating point numbers and @code{%d} for integers. The function returns +0 for success and @code{GSL_EFAILED} if there was a problem writing to +the file. +@end deftypefun + +@deftypefun int gsl_block_fscanf (FILE * @var{stream}, gsl_block * @var{b}) +This function reads formatted data from the stream @var{stream} into the +block @var{b}. The block @var{b} must be preallocated with the correct +length since the function uses the size of @var{b} to determine how many +numbers to read. The function returns 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. +@end deftypefun +@comment + +@node Example programs for blocks +@subsection Example programs for blocks + +The following program shows how to allocate a block, + +@example +@verbatiminclude examples/block.c +@end example +@comment + +@noindent +Here is the output from the program, + +@example +@verbatiminclude examples/block.out +@end example +@comment + +@node Vectors +@section Vectors +@cindex vectors +@cindex stride, of vector index + +Vectors are defined by a @code{gsl_vector} structure which describes a +slice of a block. Different vectors can be created which point to the +same block. A vector slice is a set of equally-spaced elements of an +area of memory. + +The @code{gsl_vector} structure contains five components, the +@dfn{size}, the @dfn{stride}, a pointer to the memory where the elements +are stored, @var{data}, a pointer to the block owned by the vector, +@var{block}, if any, and an ownership flag, @var{owner}. The structure +is very simple and looks like this, + +@example +typedef struct +@{ + size_t size; + size_t stride; + double * data; + gsl_block * block; + int owner; +@} gsl_vector; +@end example +@comment + +@noindent +The @var{size} is simply the number of vector elements. The range of +valid indices runs from 0 to @code{size-1}. The @var{stride} is the +step-size from one element to the next in physical memory, measured in +units of the appropriate datatype. The pointer @var{data} gives the +location of the first element of the vector in memory. The pointer +@var{block} stores the location of the memory block in which the vector +elements are located (if any). If the vector owns this block then the +@var{owner} field is set to one and the block will be deallocated when the +vector is freed. If the vector points to a block owned by another +object then the @var{owner} field is zero and any underlying block will not be +deallocated with the vector. + +The functions for allocating and accessing vectors are defined in +@file{gsl_vector.h} + +@menu +* Vector allocation:: +* Accessing vector elements:: +* Initializing vector elements:: +* Reading and writing vectors:: +* Vector views:: +* Copying vectors:: +* Exchanging elements:: +* Vector operations:: +* Finding maximum and minimum elements of vectors:: +* Vector properties:: +* Example programs for vectors:: +@end menu + +@node Vector allocation +@subsection Vector allocation + +The functions for allocating memory to a vector follow the style of +@code{malloc} and @code{free}. In addition they also perform their own +error checking. If there is insufficient memory available to allocate a +vector then the functions call the GSL error handler (with an error +number of @code{GSL_ENOMEM}) in addition to returning a null +pointer. Thus if you use the library error handler to abort your program +then it isn't necessary to check every @code{alloc}. + +@deftypefun {gsl_vector *} gsl_vector_alloc (size_t @var{n}) +This function creates a vector of length @var{n}, returning a pointer to +a newly initialized vector struct. A new block is allocated for the +elements of the vector, and stored in the @var{block} component of the +vector struct. The block is ``owned'' by the vector, and will be +deallocated when the vector is deallocated. +@end deftypefun + +@deftypefun {gsl_vector *} gsl_vector_calloc (size_t @var{n}) +This function allocates memory for a vector of length @var{n} and +initializes all the elements of the vector to zero. +@end deftypefun + +@deftypefun void gsl_vector_free (gsl_vector * @var{v}) +This function frees a previously allocated vector @var{v}. If the +vector was created using @code{gsl_vector_alloc} then the block +underlying the vector will also be deallocated. If the vector has +been created from another object then the memory is still owned by +that object and will not be deallocated. The vector @var{v} must be a +valid vector object (a null pointer is not allowed). +@end deftypefun + +@node Accessing vector elements +@subsection Accessing vector elements +@cindex vectors, range-checking +@cindex range-checking for vectors +@cindex bounds checking, extension to GCC +@cindex gcc extensions, range-checking +@cindex Fortran range checking, equivalent in gcc + +Unlike @sc{fortran} compilers, C compilers do not usually provide +support for range checking of vectors and matrices. Range checking is +available in the GNU C Compiler bounds-checking extension, but it is not +part of the default installation of GCC. The functions +@code{gsl_vector_get} and @code{gsl_vector_set} can perform portable +range checking for you and report an error if you attempt to access +elements outside the allowed range. + +The functions for accessing the elements of a vector or matrix are +defined in @file{gsl_vector.h} and declared @code{extern inline} to +eliminate function-call overhead. You must compile your program with +the macro @code{HAVE_INLINE} defined to use these functions. + +If necessary you can turn off range checking completely without +modifying any source files by recompiling your program with the +preprocessor definition @code{GSL_RANGE_CHECK_OFF}. Provided your +compiler supports inline functions the effect of turning off range +checking is to replace calls to @code{gsl_vector_get(v,i)} by +@code{v->data[i*v->stride]} and calls to @code{gsl_vector_set(v,i,x)} by +@code{v->data[i*v->stride]=x}. Thus there should be no performance +penalty for using the range checking functions when range checking is +turned off. + +@deftypefun double gsl_vector_get (const gsl_vector * @var{v}, size_t @var{i}) +This function returns the @var{i}-th element of a vector @var{v}. If +@var{i} lies outside the allowed range of 0 to @math{@var{n}-1} then the error +handler is invoked and 0 is returned. +@end deftypefun + +@deftypefun void gsl_vector_set (gsl_vector * @var{v}, size_t @var{i}, double @var{x}) +This function sets the value of the @var{i}-th element of a vector +@var{v} to @var{x}. If @var{i} lies outside the allowed range of 0 to +@math{@var{n}-1} then the error handler is invoked. +@end deftypefun + +@deftypefun {double *} gsl_vector_ptr (gsl_vector * @var{v}, size_t @var{i}) +@deftypefunx {const double *} gsl_vector_const_ptr (const gsl_vector * @var{v}, size_t @var{i}) +These functions return a pointer to the @var{i}-th element of a vector +@var{v}. If @var{i} lies outside the allowed range of 0 to @math{@var{n}-1} +then the error handler is invoked and a null pointer is returned. +@end deftypefun + +@node Initializing vector elements +@subsection Initializing vector elements +@cindex vectors, initializing +@cindex initializing vectors + +@deftypefun void gsl_vector_set_all (gsl_vector * @var{v}, double @var{x}) +This function sets all the elements of the vector @var{v} to the value +@var{x}. +@end deftypefun + +@deftypefun void gsl_vector_set_zero (gsl_vector * @var{v}) +This function sets all the elements of the vector @var{v} to zero. +@end deftypefun + +@deftypefun int gsl_vector_set_basis (gsl_vector * @var{v}, size_t @var{i}) +This function makes a basis vector by setting all the elements of the +vector @var{v} to zero except for the @var{i}-th element which is set to +one. +@end deftypefun + +@node Reading and writing vectors +@subsection Reading and writing vectors + +The library provides functions for reading and writing vectors to a file +as binary data or formatted text. + +@deftypefun int gsl_vector_fwrite (FILE * @var{stream}, const gsl_vector * @var{v}) +This function writes the elements of the vector @var{v} to the stream +@var{stream} in binary format. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem writing to the file. Since the +data is written in the native binary format it may not be portable +between different architectures. +@end deftypefun + +@deftypefun int gsl_vector_fread (FILE * @var{stream}, gsl_vector * @var{v}) +This function reads into the vector @var{v} from the open stream +@var{stream} in binary format. The vector @var{v} must be preallocated +with the correct length since the function uses the size of @var{v} to +determine how many bytes to read. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. The +data is assumed to have been written in the native binary format on the +same architecture. +@end deftypefun + +@deftypefun int gsl_vector_fprintf (FILE * @var{stream}, const gsl_vector * @var{v}, const char * @var{format}) +This function writes the elements of the vector @var{v} line-by-line to +the stream @var{stream} using the format specifier @var{format}, which +should be one of the @code{%g}, @code{%e} or @code{%f} formats for +floating point numbers and @code{%d} for integers. The function returns +0 for success and @code{GSL_EFAILED} if there was a problem writing to +the file. +@end deftypefun + +@deftypefun int gsl_vector_fscanf (FILE * @var{stream}, gsl_vector * @var{v}) +This function reads formatted data from the stream @var{stream} into the +vector @var{v}. The vector @var{v} must be preallocated with the correct +length since the function uses the size of @var{v} to determine how many +numbers to read. The function returns 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. +@end deftypefun + +@node Vector views +@subsection Vector views + +In addition to creating vectors from slices of blocks it is also +possible to slice vectors and create vector views. For example, a +subvector of another vector can be described with a view, or two views +can be made which provide access to the even and odd elements of a +vector. + +A vector view is a temporary object, stored on the stack, which can be +used to operate on a subset of vector elements. Vector views can be +defined for both constant and non-constant vectors, using separate types +that preserve constness. A vector view has the type +@code{gsl_vector_view} and a constant vector view has the type +@code{gsl_vector_const_view}. In both cases the elements of the view +can be accessed as a @code{gsl_vector} using the @code{vector} component +of the view object. A pointer to a vector of type @code{gsl_vector *} +or @code{const gsl_vector *} can be obtained by taking the address of +this component with the @code{&} operator. + +When using this pointer it is important to ensure that the view itself +remains in scope---the simplest way to do so is by always writing the +pointer as @code{&}@var{view}@code{.vector}, and never storing this value +in another variable. + +@deftypefun gsl_vector_view gsl_vector_subvector (gsl_vector * @var{v}, size_t @var{offset}, size_t @var{n}) +@deftypefunx {gsl_vector_const_view} gsl_vector_const_subvector (const gsl_vector * @var{v}, size_t @var{offset}, size_t @var{n}) +These functions return a vector view of a subvector of another vector +@var{v}. The start of the new vector is offset by @var{offset} elements +from the start of the original vector. The new vector has @var{n} +elements. Mathematically, the @var{i}-th element of the new vector +@var{v'} is given by, + +@example +v'(i) = v->data[(offset + i)*v->stride] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n-1}. + +The @code{data} pointer of the returned vector struct is set to null if +the combined parameters (@var{offset},@var{n}) overrun the end of the +original vector. + +The new vector is only a view of the block underlying the original +vector, @var{v}. The block containing the elements of @var{v} is not +owned by the new vector. When the view goes out of scope the original +vector @var{v} and its block will continue to exist. The original +memory can only be deallocated by freeing the original vector. Of +course, the original vector should not be deallocated while the view is +still in use. + +The function @code{gsl_vector_const_subvector} is equivalent to +@code{gsl_vector_subvector} but can be used for vectors which are +declared @code{const}. +@end deftypefun + + +@deftypefun gsl_vector_view gsl_vector_subvector_with_stride (gsl_vector * @var{v}, size_t @var{offset}, size_t @var{stride}, size_t @var{n}) +@deftypefunx {gsl_vector_const_view} gsl_vector_const_subvector_with_stride (const gsl_vector * @var{v}, size_t @var{offset}, size_t @var{stride}, size_t @var{n}) +These functions return a vector view of a subvector of another vector +@var{v} with an additional stride argument. The subvector is formed in +the same way as for @code{gsl_vector_subvector} but the new vector has +@var{n} elements with a step-size of @var{stride} from one element to +the next in the original vector. Mathematically, the @var{i}-th element +of the new vector @var{v'} is given by, + +@example +v'(i) = v->data[(offset + i*stride)*v->stride] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n-1}. + +Note that subvector views give direct access to the underlying elements +of the original vector. For example, the following code will zero the +even elements of the vector @code{v} of length @code{n}, while leaving the +odd elements untouched, + +@example +gsl_vector_view v_even + = gsl_vector_subvector_with_stride (v, 0, 2, n/2); +gsl_vector_set_zero (&v_even.vector); +@end example + +@noindent +A vector view can be passed to any subroutine which takes a vector +argument just as a directly allocated vector would be, using +@code{&}@var{view}@code{.vector}. For example, the following code +computes the norm of the odd elements of @code{v} using the @sc{blas} +routine @sc{dnrm2}, + +@example +gsl_vector_view v_odd + = gsl_vector_subvector_with_stride (v, 1, 2, n/2); +double r = gsl_blas_dnrm2 (&v_odd.vector); +@end example + +The function @code{gsl_vector_const_subvector_with_stride} is equivalent +to @code{gsl_vector_subvector_with_stride} but can be used for vectors +which are declared @code{const}. +@end deftypefun + +@deftypefun gsl_vector_view gsl_vector_complex_real (gsl_vector_complex * @var{v}) +@deftypefunx {gsl_vector_const_view} gsl_vector_complex_const_real (const gsl_vector_complex * @var{v}) +These functions return a vector view of the real parts of the complex +vector @var{v}. + +The function @code{gsl_vector_complex_const_real} is equivalent to +@code{gsl_vector_complex_real} but can be used for vectors which are +declared @code{const}. +@end deftypefun + +@deftypefun gsl_vector_view gsl_vector_complex_imag (gsl_vector_complex * @var{v}) +@deftypefunx {gsl_vector_const_view} gsl_vector_complex_const_imag (const gsl_vector_complex * @var{v}) +These functions return a vector view of the imaginary parts of the +complex vector @var{v}. + +The function @code{gsl_vector_complex_const_imag} is equivalent to +@code{gsl_vector_complex_imag} but can be used for vectors which are +declared @code{const}. +@end deftypefun + +@deftypefun gsl_vector_view gsl_vector_view_array (double * @var{base}, size_t @var{n}) +@deftypefunx {gsl_vector_const_view} gsl_vector_const_view_array (const double * @var{base}, size_t @var{n}) +These functions return a vector view of an array. The start of the new +vector is given by @var{base} and has @var{n} elements. Mathematically, +the @var{i}-th element of the new vector @var{v'} is given by, + +@example +v'(i) = base[i] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n-1}. + +The array containing the elements of @var{v} is not owned by the new +vector view. When the view goes out of scope the original array will +continue to exist. The original memory can only be deallocated by +freeing the original pointer @var{base}. Of course, the original array +should not be deallocated while the view is still in use. + +The function @code{gsl_vector_const_view_array} is equivalent to +@code{gsl_vector_view_array} but can be used for arrays which are +declared @code{const}. +@end deftypefun + +@deftypefun gsl_vector_view gsl_vector_view_array_with_stride (double * @var{base}, size_t @var{stride}, size_t @var{n}) +@deftypefunx gsl_vector_const_view gsl_vector_const_view_array_with_stride (const double * @var{base}, size_t @var{stride}, size_t @var{n}) +These functions return a vector view of an array @var{base} with an +additional stride argument. The subvector is formed in the same way as +for @code{gsl_vector_view_array} but the new vector has @var{n} elements +with a step-size of @var{stride} from one element to the next in the +original array. Mathematically, the @var{i}-th element of the new +vector @var{v'} is given by, + +@example +v'(i) = base[i*stride] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n-1}. + +Note that the view gives direct access to the underlying elements of the +original array. A vector view can be passed to any subroutine which +takes a vector argument just as a directly allocated vector would be, +using @code{&}@var{view}@code{.vector}. + +The function @code{gsl_vector_const_view_array_with_stride} is +equivalent to @code{gsl_vector_view_array_with_stride} but can be used +for arrays which are declared @code{const}. +@end deftypefun + + +@comment @node Modifying subvector views +@comment @subsection Modifying subvector views +@comment +@comment @deftypefun int gsl_vector_view_from_vector (gsl_vector * @var{v}, gsl_vector * @var{base}, size_t @var{offset}, size_t @var{n}, size_t @var{stride}) +@comment This function modifies and existing vector view @var{v} to form a new +@comment view of a vector @var{base}, starting from element @var{offset}. The +@comment vector has @var{n} elements separated by stride @var{stride}. Any +@comment existing view in @var{v} will be lost as a result of this function. +@comment @end deftypefun +@comment +@comment @deftypefun int gsl_vector_view_from_array (gsl_vector * @var{v}, double * @var{base}, size_t @var{offset}, size_t @var{n}, size_t @var{stride}) +@comment This function modifies and existing vector view @var{v} to form a new +@comment view of an array @var{base}, starting from element @var{offset}. The +@comment vector has @var{n} elements separated by stride @var{stride}. Any +@comment existing view in @var{v} will be lost as a result of this function. +@comment @end deftypefun + +@comment @deftypefun {gsl_vector *} gsl_vector_alloc_from_block (gsl_block * @var{b}, size_t @var{offset}, size_t @var{n}, size_t @var{stride}) +@comment This function creates a vector as a slice of an existing block @var{b}, +@comment returning a pointer to a newly initialized vector struct. The start of +@comment the vector is offset by @var{offset} elements from the start of the +@comment block. The vector has @var{n} elements, with a step-size of @var{stride} +@comment from one element to the next. Mathematically, the @var{i}-th element of +@comment the vector is given by, +@comment +@comment @example +@comment v(i) = b->data[offset + i*stride] +@comment @end example +@comment @noindent +@comment where the index @var{i} runs from 0 to @code{n-1}. +@comment +@comment A null pointer is returned if the combined parameters +@comment (@var{offset},@var{n},@var{stride}) overrun the end of the block or if +@comment insufficient memory is available to store the vector. +@comment +@comment The vector is only a view of the block @var{b}, and the block is not +@comment owned by the vector. When the vector is deallocated the block @var{b} +@comment will continue to exist. This memory can only be deallocated by freeing +@comment the block itself. Of course, this block should not be deallocated while +@comment the vector is still in use. +@comment @end deftypefun +@comment +@comment @deftypefun {gsl_vector *} gsl_vector_alloc_from_vector (gsl_vector * @var{v}, size_t @var{offset}, size_t @var{n}, size_t @var{stride}) +@comment This function creates a vector as a slice of another vector @var{v}, +@comment returning a pointer to a newly initialized vector struct. The start of +@comment the new vector is offset by @var{offset} elements from the start of the +@comment original vector. The new vector has @var{n} elements, with a step-size +@comment of @var{stride} from one element to the next in the original vector. +@comment Mathematically, the @var{i}-th element of the new vector @var{v'} is +@comment given by, +@comment +@comment @example +@comment v'(i) = v->data[(offset + i*stride)*v->stride] +@comment @end example +@comment @noindent +@comment where the index @var{i} runs from 0 to @code{n-1}. +@comment +@comment A null pointer is returned if the combined parameters +@comment (@var{offset},@var{n},@var{stride}) overrun the end of the original +@comment vector or if insufficient memory is available store the new vector. +@comment +@comment The new vector is only a view of the block underlying the original +@comment vector, @var{v}. The block is not owned by the new vector. When the new +@comment vector is deallocated the original vector @var{v} and its block will +@comment continue to exist. The original memory can only be deallocated by +@comment freeing the original vector. Of course, the original vector should not +@comment be deallocated while the new vector is still in use. +@comment @end deftypefun + +@node Copying vectors +@subsection Copying vectors + +Common operations on vectors such as addition and multiplication are +available in the @sc{blas} part of the library (@pxref{BLAS +Support}). However, it is useful to have a small number of utility +functions which do not require the full @sc{blas} code. The following +functions fall into this category. + +@deftypefun int gsl_vector_memcpy (gsl_vector * @var{dest}, const gsl_vector * @var{src}) +This function copies the elements of the vector @var{src} into the +vector @var{dest}. The two vectors must have the same length. +@end deftypefun + +@deftypefun int gsl_vector_swap (gsl_vector * @var{v}, gsl_vector * @var{w}) +This function exchanges the elements of the vectors @var{v} and @var{w} +by copying. The two vectors must have the same length. +@end deftypefun + + +@node Exchanging elements +@subsection Exchanging elements + +The following function can be used to exchange, or permute, the elements +of a vector. + +@deftypefun int gsl_vector_swap_elements (gsl_vector * @var{v}, size_t @var{i}, size_t @var{j}) +This function exchanges the @var{i}-th and @var{j}-th elements of the +vector @var{v} in-place. +@end deftypefun + +@deftypefun int gsl_vector_reverse (gsl_vector * @var{v}) +This function reverses the order of the elements of the vector @var{v}. +@end deftypefun + + +@node Vector operations +@subsection Vector operations + +The following operations are only defined for real vectors. + +@deftypefun int gsl_vector_add (gsl_vector * @var{a}, const gsl_vector * @var{b}) +This function adds the elements of vector @var{b} to the elements of +vector @var{a}, @math{a'_i = a_i + b_i}. The two vectors must have the +same length. +@end deftypefun + +@deftypefun int gsl_vector_sub (gsl_vector * @var{a}, const gsl_vector * @var{b}) +This function subtracts the elements of vector @var{b} from the elements of +vector @var{a}, @math{a'_i = a_i - b_i}. The two vectors must have the +same length. +@end deftypefun + +@deftypefun int gsl_vector_mul (gsl_vector * @var{a}, const gsl_vector * @var{b}) +This function multiplies the elements of vector @var{a} by the elements of +vector @var{b}, @math{a'_i = a_i * b_i}. The two vectors must have the +same length. +@end deftypefun + +@deftypefun int gsl_vector_div (gsl_vector * @var{a}, const gsl_vector * @var{b}) +This function divides the elements of vector @var{a} by the elements of +vector @var{b}, @math{a'_i = a_i / b_i}. The two vectors must have the +same length. +@end deftypefun + +@deftypefun int gsl_vector_scale (gsl_vector * @var{a}, const double @var{x}) +This function multiplies the elements of vector @var{a} by the constant +factor @var{x}, @math{a'_i = x a_i}. +@end deftypefun + +@deftypefun int gsl_vector_add_constant (gsl_vector * @var{a}, const double @var{x}) +This function adds the constant value @var{x} to the elements of the +vector @var{a}, @math{a'_i = a_i + x}. +@end deftypefun + +@node Finding maximum and minimum elements of vectors +@subsection Finding maximum and minimum elements of vectors + +@deftypefun double gsl_vector_max (const gsl_vector * @var{v}) +This function returns the maximum value in the vector @var{v}. +@end deftypefun + +@deftypefun double gsl_vector_min (const gsl_vector * @var{v}) +This function returns the minimum value in the vector @var{v}. +@end deftypefun + +@deftypefun void gsl_vector_minmax (const gsl_vector * @var{v}, double * @var{min_out}, double * @var{max_out}) +This function returns the minimum and maximum values in the vector +@var{v}, storing them in @var{min_out} and @var{max_out}. +@end deftypefun + +@deftypefun size_t gsl_vector_max_index (const gsl_vector * @var{v}) +This function returns the index of the maximum value in the vector @var{v}. +When there are several equal maximum elements then the lowest index is +returned. +@end deftypefun + +@deftypefun size_t gsl_vector_min_index (const gsl_vector * @var{v}) +This function returns the index of the minimum value in the vector @var{v}. +When there are several equal minimum elements then the lowest index is +returned. +@end deftypefun + +@deftypefun void gsl_vector_minmax_index (const gsl_vector * @var{v}, size_t * @var{imin}, size_t * @var{imax}) +This function returns the indices of the minimum and maximum values in +the vector @var{v}, storing them in @var{imin} and @var{imax}. When +there are several equal minimum or maximum elements then the lowest +indices are returned. +@end deftypefun + +@node Vector properties +@subsection Vector properties + +@deftypefun int gsl_vector_isnull (const gsl_vector * @var{v}) +@deftypefunx int gsl_vector_ispos (const gsl_vector * @var{v}) +@deftypefunx int gsl_vector_isneg (const gsl_vector * @var{v}) +These functions return 1 if all the elements of the vector @var{v} +are zero, strictly positive, or strictly negative respectively, and 0 +otherwise. To test for a non-negative vector, use the expression +@code{!gsl_vector_isneg(v)}. +@end deftypefun + +@node Example programs for vectors +@subsection Example programs for vectors + +This program shows how to allocate, initialize and read from a vector +using the functions @code{gsl_vector_alloc}, @code{gsl_vector_set} and +@code{gsl_vector_get}. + +@example +@verbatiminclude examples/vector.c +@end example +@comment + +@noindent +Here is the output from the program. The final loop attempts to read +outside the range of the vector @code{v}, and the error is trapped by +the range-checking code in @code{gsl_vector_get}. + +@example +$ ./a.out +v_0 = 1.23 +v_1 = 2.23 +v_2 = 3.23 +gsl: vector_source.c:12: ERROR: index out of range +Default GSL error handler invoked. +Aborted (core dumped) +@end example +@comment + +@noindent +The next program shows how to write a vector to a file. + +@example +@verbatiminclude examples/vectorw.c +@end example +@comment + +@noindent +After running this program the file @file{test.dat} should contain the +elements of @code{v}, written using the format specifier +@code{%.5g}. The vector could then be read back in using the function +@code{gsl_vector_fscanf (f, v)} as follows: + +@example +@verbatiminclude examples/vectorr.c +@end example + + +@node Matrices +@section Matrices +@cindex matrices +@cindex physical dimension, matrices +@cindex trailing dimension, matrices +@cindex leading dimension, matrices + +Matrices are defined by a @code{gsl_matrix} structure which describes a +generalized slice of a block. Like a vector it represents a set of +elements in an area of memory, but uses two indices instead of one. + +The @code{gsl_matrix} structure contains six components, the two +dimensions of the matrix, a physical dimension, a pointer to the memory +where the elements of the matrix are stored, @var{data}, a pointer to +the block owned by the matrix @var{block}, if any, and an ownership +flag, @var{owner}. The physical dimension determines the memory layout +and can differ from the matrix dimension to allow the use of +submatrices. The @code{gsl_matrix} structure is very simple and looks +like this, + +@example +typedef struct +@{ + size_t size1; + size_t size2; + size_t tda; + double * data; + gsl_block * block; + int owner; +@} gsl_matrix; +@end example +@comment + +@noindent +Matrices are stored in row-major order, meaning that each row of +elements forms a contiguous block in memory. This is the standard +``C-language ordering'' of two-dimensional arrays. Note that @sc{fortran} +stores arrays in column-major order. The number of rows is @var{size1}. +The range of valid row indices runs from 0 to @code{size1-1}. Similarly +@var{size2} is the number of columns. The range of valid column indices +runs from 0 to @code{size2-1}. The physical row dimension @var{tda}, or +@dfn{trailing dimension}, specifies the size of a row of the matrix as +laid out in memory. + +For example, in the following matrix @var{size1} is 3, @var{size2} is 4, +and @var{tda} is 8. The physical memory layout of the matrix begins in +the top left hand-corner and proceeds from left to right along each row +in turn. + +@example +@group +00 01 02 03 XX XX XX XX +10 11 12 13 XX XX XX XX +20 21 22 23 XX XX XX XX +@end group +@end example + +@noindent +Each unused memory location is represented by ``@code{XX}''. The +pointer @var{data} gives the location of the first element of the matrix +in memory. The pointer @var{block} stores the location of the memory +block in which the elements of the matrix are located (if any). If the +matrix owns this block then the @var{owner} field is set to one and the +block will be deallocated when the matrix is freed. If the matrix is +only a slice of a block owned by another object then the @var{owner} field is +zero and any underlying block will not be freed. + +The functions for allocating and accessing matrices are defined in +@file{gsl_matrix.h} + +@menu +* Matrix allocation:: +* Accessing matrix elements:: +* Initializing matrix elements:: +* Reading and writing matrices:: +* Matrix views:: +* Creating row and column views:: +* Copying matrices:: +* Copying rows and columns:: +* Exchanging rows and columns:: +* Matrix operations:: +* Finding maximum and minimum elements of matrices:: +* Matrix properties:: +* Example programs for matrices:: +@end menu + +@node Matrix allocation +@subsection Matrix allocation + +The functions for allocating memory to a matrix follow the style of +@code{malloc} and @code{free}. They also perform their own error +checking. If there is insufficient memory available to allocate a vector +then the functions call the GSL error handler (with an error number of +@code{GSL_ENOMEM}) in addition to returning a null pointer. Thus if you +use the library error handler to abort your program then it isn't +necessary to check every @code{alloc}. + +@deftypefun {gsl_matrix *} gsl_matrix_alloc (size_t @var{n1}, size_t @var{n2}) +This function creates a matrix of size @var{n1} rows by @var{n2} +columns, returning a pointer to a newly initialized matrix struct. A new +block is allocated for the elements of the matrix, and stored in the +@var{block} component of the matrix struct. The block is ``owned'' by the +matrix, and will be deallocated when the matrix is deallocated. +@end deftypefun + +@deftypefun {gsl_matrix *} gsl_matrix_calloc (size_t @var{n1}, size_t @var{n2}) +This function allocates memory for a matrix of size @var{n1} rows by +@var{n2} columns and initializes all the elements of the matrix to zero. +@end deftypefun + +@deftypefun void gsl_matrix_free (gsl_matrix * @var{m}) +This function frees a previously allocated matrix @var{m}. If the +matrix was created using @code{gsl_matrix_alloc} then the block +underlying the matrix will also be deallocated. If the matrix has +been created from another object then the memory is still owned by +that object and will not be deallocated. The matrix @var{m} must be a +valid matrix object (a null pointer is not allowed). +@end deftypefun + +@node Accessing matrix elements +@subsection Accessing matrix elements +@cindex matrices, range-checking +@cindex range-checking for matrices + +The functions for accessing the elements of a matrix use the same range +checking system as vectors. You can turn off range checking by recompiling +your program with the preprocessor definition +@code{GSL_RANGE_CHECK_OFF}. + +The elements of the matrix are stored in ``C-order'', where the second +index moves continuously through memory. More precisely, the element +accessed by the function @code{gsl_matrix_get(m,i,j)} and +@code{gsl_matrix_set(m,i,j,x)} is + +@example +m->data[i * m->tda + j] +@end example +@comment + +@noindent +where @var{tda} is the physical row-length of the matrix. + +@deftypefun double gsl_matrix_get (const gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +This function returns the @math{(i,j)}-th element of a matrix +@var{m}. If @var{i} or @var{j} lie outside the allowed range of 0 to +@math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error handler is invoked and 0 +is returned. +@end deftypefun + +@deftypefun void gsl_matrix_set (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}, double @var{x}) +This function sets the value of the @math{(i,j)}-th element of a +matrix @var{m} to @var{x}. If @var{i} or @var{j} lies outside the +allowed range of 0 to @math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error +handler is invoked. +@end deftypefun + +@deftypefun {double *} gsl_matrix_ptr (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +@deftypefunx {const double *} gsl_matrix_const_ptr (const gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +These functions return a pointer to the @math{(i,j)}-th element of a +matrix @var{m}. If @var{i} or @var{j} lie outside the allowed range of +0 to @math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error handler is invoked +and a null pointer is returned. +@end deftypefun + +@node Initializing matrix elements +@subsection Initializing matrix elements +@cindex matrices, initializing +@cindex initializing matrices +@cindex identity matrix +@cindex matrix, identity +@cindex zero matrix +@cindex matrix, zero +@cindex constant matrix +@cindex matrix, constant + +@deftypefun void gsl_matrix_set_all (gsl_matrix * @var{m}, double @var{x}) +This function sets all the elements of the matrix @var{m} to the value +@var{x}. +@end deftypefun + +@deftypefun void gsl_matrix_set_zero (gsl_matrix * @var{m}) +This function sets all the elements of the matrix @var{m} to zero. +@end deftypefun + +@deftypefun void gsl_matrix_set_identity (gsl_matrix * @var{m}) +This function sets the elements of the matrix @var{m} to the +corresponding elements of the identity matrix, @math{m(i,j) = +\delta(i,j)}, i.e. a unit diagonal with all off-diagonal elements zero. +This applies to both square and rectangular matrices. +@end deftypefun + +@node Reading and writing matrices +@subsection Reading and writing matrices + +The library provides functions for reading and writing matrices to a file +as binary data or formatted text. + +@deftypefun int gsl_matrix_fwrite (FILE * @var{stream}, const gsl_matrix * @var{m}) +This function writes the elements of the matrix @var{m} to the stream +@var{stream} in binary format. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem writing to the file. Since the +data is written in the native binary format it may not be portable +between different architectures. +@end deftypefun + +@deftypefun int gsl_matrix_fread (FILE * @var{stream}, gsl_matrix * @var{m}) +This function reads into the matrix @var{m} from the open stream +@var{stream} in binary format. The matrix @var{m} must be preallocated +with the correct dimensions since the function uses the size of @var{m} to +determine how many bytes to read. The return value is 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. The +data is assumed to have been written in the native binary format on the +same architecture. +@end deftypefun + +@deftypefun int gsl_matrix_fprintf (FILE * @var{stream}, const gsl_matrix * @var{m}, const char * @var{format}) +This function writes the elements of the matrix @var{m} line-by-line to +the stream @var{stream} using the format specifier @var{format}, which +should be one of the @code{%g}, @code{%e} or @code{%f} formats for +floating point numbers and @code{%d} for integers. The function returns +0 for success and @code{GSL_EFAILED} if there was a problem writing to +the file. +@end deftypefun + +@deftypefun int gsl_matrix_fscanf (FILE * @var{stream}, gsl_matrix * @var{m}) +This function reads formatted data from the stream @var{stream} into the +matrix @var{m}. The matrix @var{m} must be preallocated with the correct +dimensions since the function uses the size of @var{m} to determine how many +numbers to read. The function returns 0 for success and +@code{GSL_EFAILED} if there was a problem reading from the file. +@end deftypefun + +@node Matrix views +@subsection Matrix views + +A matrix view is a temporary object, stored on the stack, which can be +used to operate on a subset of matrix elements. Matrix views can be +defined for both constant and non-constant matrices using separate types +that preserve constness. A matrix view has the type +@code{gsl_matrix_view} and a constant matrix view has the type +@code{gsl_matrix_const_view}. In both cases the elements of the view +can by accessed using the @code{matrix} component of the view object. A +pointer @code{gsl_matrix *} or @code{const gsl_matrix *} can be obtained +by taking the address of the @code{matrix} component with the @code{&} +operator. In addition to matrix views it is also possible to create +vector views of a matrix, such as row or column views. + +@deftypefun gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * @var{m}, size_t @var{k1}, size_t @var{k2}, size_t @var{n1}, size_t @var{n2}) +@deftypefunx gsl_matrix_const_view gsl_matrix_const_submatrix (const gsl_matrix * @var{m}, size_t @var{k1}, size_t @var{k2}, size_t @var{n1}, size_t @var{n2}) +These functions return a matrix view of a submatrix of the matrix +@var{m}. The upper-left element of the submatrix is the element +(@var{k1},@var{k2}) of the original matrix. The submatrix has @var{n1} +rows and @var{n2} columns. The physical number of columns in memory +given by @var{tda} is unchanged. Mathematically, the +@math{(i,j)}-th element of the new matrix is given by, + +@example +m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +runs from 0 to @code{n2-1}. + +The @code{data} pointer of the returned matrix struct is set to null if +the combined parameters (@var{i},@var{j},@var{n1},@var{n2},@var{tda}) +overrun the ends of the original matrix. + +The new matrix view is only a view of the block underlying the existing +matrix, @var{m}. The block containing the elements of @var{m} is not +owned by the new matrix view. When the view goes out of scope the +original matrix @var{m} and its block will continue to exist. The +original memory can only be deallocated by freeing the original matrix. +Of course, the original matrix should not be deallocated while the view +is still in use. + +The function @code{gsl_matrix_const_submatrix} is equivalent to +@code{gsl_matrix_submatrix} but can be used for matrices which are +declared @code{const}. +@end deftypefun + + +@deftypefun gsl_matrix_view gsl_matrix_view_array (double * @var{base}, size_t @var{n1}, size_t @var{n2}) +@deftypefunx gsl_matrix_const_view gsl_matrix_const_view_array (const double * @var{base}, size_t @var{n1}, size_t @var{n2}) +These functions return a matrix view of the array @var{base}. The +matrix has @var{n1} rows and @var{n2} columns. The physical number of +columns in memory is also given by @var{n2}. Mathematically, the +@math{(i,j)}-th element of the new matrix is given by, + +@example +m'(i,j) = base[i*n2 + j] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +runs from 0 to @code{n2-1}. + +The new matrix is only a view of the array @var{base}. When the view +goes out of scope the original array @var{base} will continue to exist. +The original memory can only be deallocated by freeing the original +array. Of course, the original array should not be deallocated while +the view is still in use. + +The function @code{gsl_matrix_const_view_array} is equivalent to +@code{gsl_matrix_view_array} but can be used for matrices which are +declared @code{const}. +@end deftypefun + + +@deftypefun gsl_matrix_view gsl_matrix_view_array_with_tda (double * @var{base}, size_t @var{n1}, size_t @var{n2}, size_t @var{tda}) +@deftypefunx gsl_matrix_const_view gsl_matrix_const_view_array_with_tda (const double * @var{base}, size_t @var{n1}, size_t @var{n2}, size_t @var{tda}) +These functions return a matrix view of the array @var{base} with a +physical number of columns @var{tda} which may differ from the corresponding +dimension of the matrix. The matrix has @var{n1} rows and @var{n2} +columns, and the physical number of columns in memory is given by +@var{tda}. Mathematically, the @math{(i,j)}-th element of the new +matrix is given by, + +@example +m'(i,j) = base[i*tda + j] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +runs from 0 to @code{n2-1}. + +The new matrix is only a view of the array @var{base}. When the view +goes out of scope the original array @var{base} will continue to exist. +The original memory can only be deallocated by freeing the original +array. Of course, the original array should not be deallocated while +the view is still in use. + +The function @code{gsl_matrix_const_view_array_with_tda} is equivalent +to @code{gsl_matrix_view_array_with_tda} but can be used for matrices +which are declared @code{const}. +@end deftypefun + +@deftypefun gsl_matrix_view gsl_matrix_view_vector (gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2}) +@deftypefunx gsl_matrix_const_view gsl_matrix_const_view_vector (const gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2}) +These functions return a matrix view of the vector @var{v}. The matrix +has @var{n1} rows and @var{n2} columns. The vector must have unit +stride. The physical number of columns in memory is also given by +@var{n2}. Mathematically, the @math{(i,j)}-th element of the new +matrix is given by, + +@example +m'(i,j) = v->data[i*n2 + j] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +runs from 0 to @code{n2-1}. + +The new matrix is only a view of the vector @var{v}. When the view +goes out of scope the original vector @var{v} will continue to exist. +The original memory can only be deallocated by freeing the original +vector. Of course, the original vector should not be deallocated while +the view is still in use. + +The function @code{gsl_matrix_const_view_vector} is equivalent to +@code{gsl_matrix_view_vector} but can be used for matrices which are +declared @code{const}. +@end deftypefun + + +@deftypefun gsl_matrix_view gsl_matrix_view_vector_with_tda (gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2}, size_t @var{tda}) +@deftypefunx gsl_matrix_const_view gsl_matrix_const_view_vector_with_tda (const gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2}, size_t @var{tda}) +These functions return a matrix view of the vector @var{v} with a +physical number of columns @var{tda} which may differ from the +corresponding matrix dimension. The vector must have unit stride. The +matrix has @var{n1} rows and @var{n2} columns, and the physical number +of columns in memory is given by @var{tda}. Mathematically, the +@math{(i,j)}-th element of the new matrix is given by, + +@example +m'(i,j) = v->data[i*tda + j] +@end example + +@noindent +where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +runs from 0 to @code{n2-1}. + +The new matrix is only a view of the vector @var{v}. When the view +goes out of scope the original vector @var{v} will continue to exist. +The original memory can only be deallocated by freeing the original +vector. Of course, the original vector should not be deallocated while +the view is still in use. + +The function @code{gsl_matrix_const_view_vector_with_tda} is equivalent +to @code{gsl_matrix_view_vector_with_tda} but can be used for matrices +which are declared @code{const}. +@end deftypefun + + +@comment @node Modifying matrix views +@comment @subsection Modifying matrix views +@comment +@comment @deftypefun int gsl_matrix_view_from_matrix (gsl_matrix * @var{m}, gsl_matrix * @var{mm}, const size_t @var{k1}, const size_t @var{k2}, const size_t @var{n1}, const size_t @var{n2}) +@comment This function modifies and existing matrix view @var{m} to form a new +@comment view of a matrix @var{mm}, starting from element (@var{k1},@var{k2}). +@comment The matrix view has @var{n1} rows and @var{n2} columns. Any existing +@comment view in @var{m} will be lost as a result of this function. +@comment @end deftypefun +@comment +@comment @deftypefun int gsl_matrix_view_from_vector (gsl_matrix * @var{m}, gsl_vector * @var{v}, const size_t @var{offset}, const size_t @var{n1}, const size_t @var{n2}) +@comment This function modifies and existing matrix view @var{m} to form a new +@comment view of a vector @var{v}, starting from element @var{offset}. The +@comment vector has @var{n1} rows and @var{n2} columns. Any +@comment existing view in @var{m} will be lost as a result of this function. +@comment @end deftypefun +@comment +@comment @deftypefun int gsl_matrix_view_from_array (gsl_matrix * @var{m}, double * @var{base}, const size_t @var{offset}, const size_t @var{n1}, const size_t @var{n2}) +@comment This function modifies and existing matrix view @var{m} to form a new +@comment view of an array @var{base}, starting from element @var{offset}. The +@comment matrix has @var{n1} rows and @var{n2} columns. Any +@comment existing view in @var{m} will be lost as a result of this function. +@comment @end deftypefun +@comment +@comment @deftypefun {gsl_matrix *} gsl_matrix_alloc_from_block (gsl_block * @var{b}, size_t @var{offset}, size_t @var{n1}, size_t @var{n2}, size_t @var{tda}) +@comment This function creates a matrix as a slice of the block @var{b}, +@comment returning a pointer to a newly initialized matrix struct. The start of +@comment the matrix is offset by @var{offset} elements from the start of the +@comment block. The matrix has @var{n1} rows and @var{n2} columns, with the +@comment physical number of columns in memory given by @var{tda}. +@comment Mathematically, the (@var{i},@var{j})-th element of the matrix is given by, +@comment +@comment @example +@comment m(i,j) = b->data[offset + i*tda + j] +@comment @end example +@comment @noindent +@comment where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +@comment runs from 0 to @code{n2-1}. +@comment +@comment A null pointer is returned if the combined parameters +@comment (@var{offset},@var{n1},@var{n2},@var{tda}) overrun the end of the block +@comment or if insufficient memory is available to store the matrix. +@comment +@comment The matrix is only a view of the block @var{b}, and the block is not +@comment owned by the matrix. When the matrix is deallocated the block @var{b} +@comment will continue to exist. This memory can only be deallocated by freeing +@comment the block itself. Of course, this block should not be deallocated while +@comment the matrix is still in use. +@comment @end deftypefun +@comment +@comment @deftypefun {gsl_matrix *} gsl_matrix_alloc_from_matrix (gsl_matrix * @var{m}, size_t @var{k1}, size_t @var{k2}, size_t @var{n1}, size_t @var{n2}) +@comment +@comment This function creates a matrix as a submatrix of the matrix @var{m}, +@comment returning a pointer to a newly initialized matrix struct. The upper-left +@comment element of the submatrix is the element (@var{k1},@var{k2}) of the +@comment original matrix. The submatrix has @var{n1} rows and @var{n2} columns. +@comment The physical number of columns in memory given by @var{tda} is +@comment unchanged. Mathematically, the (@var{i},@var{j})-th element of the +@comment new matrix is given by, +@comment +@comment @example +@comment m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j] +@comment @end example +@comment @noindent +@comment where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j} +@comment runs from 0 to @code{n2-1}. +@comment +@comment A null pointer is returned if the combined parameters +@comment (@var{k1},@var{k2},@var{n1},@var{n2},@var{tda}) overrun the end of the +@comment original matrix or if insufficient memory is available to store the matrix. +@comment +@comment The new matrix is only a view of the block underlying the existing +@comment matrix, @var{m}. The block is not owned by the new matrix. When the new +@comment matrix is deallocated the original matrix @var{m} and its block will +@comment continue to exist. The original memory can only be deallocated by +@comment freeing the original matrix. Of course, the original matrix should not +@comment be deallocated while the new matrix is still in use. +@comment @end deftypefun + +@node Creating row and column views +@subsection Creating row and column views + +In general there are two ways to access an object, by reference or by +copying. The functions described in this section create vector views +which allow access to a row or column of a matrix by reference. +Modifying elements of the view is equivalent to modifying the matrix, +since both the vector view and the matrix point to the same memory +block. + +@deftypefun gsl_vector_view gsl_matrix_row (gsl_matrix * @var{m}, size_t @var{i}) +@deftypefunx {gsl_vector_const_view} gsl_matrix_const_row (const gsl_matrix * @var{m}, size_t @var{i}) +These functions return a vector view of the @var{i}-th row of the matrix +@var{m}. The @code{data} pointer of the new vector is set to null if +@var{i} is out of range. + +The function @code{gsl_vector_const_row} is equivalent to +@code{gsl_matrix_row} but can be used for matrices which are declared +@code{const}. +@end deftypefun + +@deftypefun gsl_vector_view gsl_matrix_column (gsl_matrix * @var{m}, size_t @var{j}) +@deftypefunx {gsl_vector_const_view} gsl_matrix_const_column (const gsl_matrix * @var{m}, size_t @var{j}) +These functions return a vector view of the @var{j}-th column of the +matrix @var{m}. The @code{data} pointer of the new vector is set to +null if @var{j} is out of range. + +The function @code{gsl_vector_const_column} is equivalent to +@code{gsl_matrix_column} but can be used for matrices which are declared +@code{const}. +@end deftypefun + +@cindex matrix diagonal +@cindex diagonal, of a matrix +@deftypefun gsl_vector_view gsl_matrix_diagonal (gsl_matrix * @var{m}) +@deftypefunx {gsl_vector_const_view} gsl_matrix_const_diagonal (const gsl_matrix * @var{m}) +These functions returns a vector view of the diagonal of the matrix +@var{m}. The matrix @var{m} is not required to be square. For a +rectangular matrix the length of the diagonal is the same as the smaller +dimension of the matrix. + +The function @code{gsl_matrix_const_diagonal} is equivalent to +@code{gsl_matrix_diagonal} but can be used for matrices which are +declared @code{const}. +@end deftypefun + +@cindex matrix subdiagonal +@cindex subdiagonal, of a matrix +@deftypefun gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * @var{m}, size_t @var{k}) +@deftypefunx {gsl_vector_const_view} gsl_matrix_const_subdiagonal (const gsl_matrix * @var{m}, size_t @var{k}) +These functions return a vector view of the @var{k}-th subdiagonal of +the matrix @var{m}. The matrix @var{m} is not required to be square. +The diagonal of the matrix corresponds to @math{k = 0}. + +The function @code{gsl_matrix_const_subdiagonal} is equivalent to +@code{gsl_matrix_subdiagonal} but can be used for matrices which are +declared @code{const}. +@end deftypefun + +@cindex matrix superdiagonal +@cindex superdiagonal, matrix +@deftypefun gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * @var{m}, size_t @var{k}) +@deftypefunx {gsl_vector_const_view} gsl_matrix_const_superdiagonal (const gsl_matrix * @var{m}, size_t @var{k}) +These functions return a vector view of the @var{k}-th superdiagonal of +the matrix @var{m}. The matrix @var{m} is not required to be square. The +diagonal of the matrix corresponds to @math{k = 0}. + +The function @code{gsl_matrix_const_superdiagonal} is equivalent to +@code{gsl_matrix_superdiagonal} but can be used for matrices which are +declared @code{const}. +@end deftypefun + +@comment @deftypefun {gsl_vector *} gsl_vector_alloc_row_from_matrix (gsl_matrix * @var{m}, size_t @var{i}) +@comment This function allocates a new @code{gsl_vector} struct which points to +@comment the @var{i}-th row of the matrix @var{m}. +@comment @end deftypefun +@comment +@comment @deftypefun {gsl_vector *} gsl_vector_alloc_col_from_matrix (gsl_matrix * @var{m}, size_t @var{j}) +@comment This function allocates a new @code{gsl_vector} struct which points to +@comment the @var{j}-th column of the matrix @var{m}. +@comment @end deftypefun + +@node Copying matrices +@subsection Copying matrices + +@deftypefun int gsl_matrix_memcpy (gsl_matrix * @var{dest}, const gsl_matrix * @var{src}) +This function copies the elements of the matrix @var{src} into the +matrix @var{dest}. The two matrices must have the same size. +@end deftypefun + +@deftypefun int gsl_matrix_swap (gsl_matrix * @var{m1}, gsl_matrix * @var{m2}) +This function exchanges the elements of the matrices @var{m1} and +@var{m2} by copying. The two matrices must have the same size. +@end deftypefun + +@node Copying rows and columns +@subsection Copying rows and columns + +The functions described in this section copy a row or column of a matrix +into a vector. This allows the elements of the vector and the matrix to +be modified independently. Note that if the matrix and the vector point +to overlapping regions of memory then the result will be undefined. The +same effect can be achieved with more generality using +@code{gsl_vector_memcpy} with vector views of rows and columns. + +@deftypefun int gsl_matrix_get_row (gsl_vector * @var{v}, const gsl_matrix * @var{m}, size_t @var{i}) +This function copies the elements of the @var{i}-th row of the matrix +@var{m} into the vector @var{v}. The length of the vector must be the +same as the length of the row. +@end deftypefun + +@deftypefun int gsl_matrix_get_col (gsl_vector * @var{v}, const gsl_matrix * @var{m}, size_t @var{j}) +This function copies the elements of the @var{j}-th column of the matrix +@var{m} into the vector @var{v}. The length of the vector must be the +same as the length of the column. +@end deftypefun + +@deftypefun int gsl_matrix_set_row (gsl_matrix * @var{m}, size_t @var{i}, const gsl_vector * @var{v}) +This function copies the elements of the vector @var{v} into the +@var{i}-th row of the matrix @var{m}. The length of the vector must be +the same as the length of the row. +@end deftypefun + +@deftypefun int gsl_matrix_set_col (gsl_matrix * @var{m}, size_t @var{j}, const gsl_vector * @var{v}) +This function copies the elements of the vector @var{v} into the +@var{j}-th column of the matrix @var{m}. The length of the vector must be +the same as the length of the column. +@end deftypefun + +@node Exchanging rows and columns +@subsection Exchanging rows and columns + +The following functions can be used to exchange the rows and columns of +a matrix. + +@deftypefun int gsl_matrix_swap_rows (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +This function exchanges the @var{i}-th and @var{j}-th rows of the matrix +@var{m} in-place. +@end deftypefun + +@deftypefun int gsl_matrix_swap_columns (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +This function exchanges the @var{i}-th and @var{j}-th columns of the +matrix @var{m} in-place. +@end deftypefun + +@deftypefun int gsl_matrix_swap_rowcol (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}) +This function exchanges the @var{i}-th row and @var{j}-th column of the +matrix @var{m} in-place. The matrix must be square for this operation to +be possible. +@end deftypefun + +@deftypefun int gsl_matrix_transpose_memcpy (gsl_matrix * @var{dest}, const gsl_matrix * @var{src}) +This function makes the matrix @var{dest} the transpose of the matrix +@var{src} by copying the elements of @var{src} into @var{dest}. This +function works for all matrices provided that the dimensions of the matrix +@var{dest} match the transposed dimensions of the matrix @var{src}. +@end deftypefun + +@deftypefun int gsl_matrix_transpose (gsl_matrix * @var{m}) +This function replaces the matrix @var{m} by its transpose by copying +the elements of the matrix in-place. The matrix must be square for this +operation to be possible. +@end deftypefun + +@node Matrix operations +@subsection Matrix operations + +The following operations are defined for real and complex matrices. + +@deftypefun int gsl_matrix_add (gsl_matrix * @var{a}, const gsl_matrix * @var{b}) +This function adds the elements of matrix @var{b} to the elements of +matrix @var{a}, @math{a'(i,j) = a(i,j) + b(i,j)}. The two matrices must have the +same dimensions. +@end deftypefun + +@deftypefun int gsl_matrix_sub (gsl_matrix * @var{a}, const gsl_matrix * @var{b}) +This function subtracts the elements of matrix @var{b} from the elements of +matrix @var{a}, @math{a'(i,j) = a(i,j) - b(i,j)}. The two matrices must have the +same dimensions. +@end deftypefun + +@deftypefun int gsl_matrix_mul_elements (gsl_matrix * @var{a}, const gsl_matrix * @var{b}) +This function multiplies the elements of matrix @var{a} by the elements of +matrix @var{b}, @math{a'(i,j) = a(i,j) * b(i,j)}. The two matrices must have the +same dimensions. +@end deftypefun + +@deftypefun int gsl_matrix_div_elements (gsl_matrix * @var{a}, const gsl_matrix * @var{b}) +This function divides the elements of matrix @var{a} by the elements of +matrix @var{b}, @math{a'(i,j) = a(i,j) / b(i,j)}. The two matrices must have the +same dimensions. +@end deftypefun + +@deftypefun int gsl_matrix_scale (gsl_matrix * @var{a}, const double @var{x}) +This function multiplies the elements of matrix @var{a} by the constant +factor @var{x}, @math{a'(i,j) = x a(i,j)}. +@end deftypefun + +@deftypefun int gsl_matrix_add_constant (gsl_matrix * @var{a}, const double @var{x}) +This function adds the constant value @var{x} to the elements of the +matrix @var{a}, @math{a'(i,j) = a(i,j) + x}. +@end deftypefun + +@node Finding maximum and minimum elements of matrices +@subsection Finding maximum and minimum elements of matrices + +The following operations are only defined for real matrices. + +@deftypefun double gsl_matrix_max (const gsl_matrix * @var{m}) +This function returns the maximum value in the matrix @var{m}. +@end deftypefun + +@deftypefun double gsl_matrix_min (const gsl_matrix * @var{m}) +This function returns the minimum value in the matrix @var{m}. +@end deftypefun + +@deftypefun void gsl_matrix_minmax (const gsl_matrix * @var{m}, double * @var{min_out}, double * @var{max_out}) +This function returns the minimum and maximum values in the matrix +@var{m}, storing them in @var{min_out} and @var{max_out}. +@end deftypefun + +@deftypefun void gsl_matrix_max_index (const gsl_matrix * @var{m}, size_t * @var{imax}, size_t * @var{jmax}) +This function returns the indices of the maximum value in the matrix +@var{m}, storing them in @var{imax} and @var{jmax}. When there are +several equal maximum elements then the first element found is returned, +searching in row-major order. +@end deftypefun + +@deftypefun void gsl_matrix_min_index (const gsl_matrix * @var{m}, size_t * @var{imin}, size_t * @var{jmin}) +This function returns the indices of the minimum value in the matrix +@var{m}, storing them in @var{imin} and @var{jmin}. When there are +several equal minimum elements then the first element found is returned, +searching in row-major order. +@end deftypefun + +@deftypefun void gsl_matrix_minmax_index (const gsl_matrix * @var{m}, size_t * @var{imin}, size_t * @var{jmin}, size_t * @var{imax}, size_t * @var{jmax}) +This function returns the indices of the minimum and maximum values in +the matrix @var{m}, storing them in (@var{imin},@var{jmin}) and +(@var{imax},@var{jmax}). When there are several equal minimum or maximum +elements then the first elements found are returned, searching in +row-major order. +@end deftypefun + +@node Matrix properties +@subsection Matrix properties + +@deftypefun int gsl_matrix_isnull (const gsl_matrix * @var{m}) +@deftypefunx int gsl_matrix_ispos (const gsl_matrix * @var{m}) +@deftypefunx int gsl_matrix_isneg (const gsl_matrix * @var{m}) +These functions return 1 if all the elements of the matrix @var{m} are +zero, strictly positive, or strictly negative respectively, and 0 +otherwise. To test for a non-negative matrix, use the expression +@code{!gsl_matrix_isneg(m)}. To test whether a matrix is +positive-definite, use the Cholesky decomposition (@pxref{Cholesky Decomposition}). +@end deftypefun + +@node Example programs for matrices +@subsection Example programs for matrices + +The program below shows how to allocate, initialize and read from a matrix +using the functions @code{gsl_matrix_alloc}, @code{gsl_matrix_set} and +@code{gsl_matrix_get}. + +@example +@verbatiminclude examples/matrix.c +@end example +@comment + +@noindent +Here is the output from the program. The final loop attempts to read +outside the range of the matrix @code{m}, and the error is trapped by +the range-checking code in @code{gsl_matrix_get}. + +@example +$ ./a.out +m(0,0) = 0.23 +m(0,1) = 1.23 +m(0,2) = 2.23 +m(1,0) = 100.23 +m(1,1) = 101.23 +m(1,2) = 102.23 +... +m(9,2) = 902.23 +gsl: matrix_source.c:13: ERROR: first index out of range +Default GSL error handler invoked. +Aborted (core dumped) +@end example +@comment + +@noindent +The next program shows how to write a matrix to a file. + +@example +@verbatiminclude examples/matrixw.c +@end example +@comment + +@noindent +After running this program the file @file{test.dat} should contain the +elements of @code{m}, written in binary format. The matrix which is read +back in using the function @code{gsl_matrix_fread} should be exactly +equal to the original matrix. + +The following program demonstrates the use of vector views. The program +computes the column norms of a matrix. + +@example +@verbatiminclude examples/vectorview.c +@end example + +@noindent +Here is the output of the program, + +@example +$ ./a.out +@verbatiminclude examples/vectorview.out +@end example + +@noindent +The results can be confirmed using @sc{gnu octave}, + +@example +$ octave +GNU Octave, version 2.0.16.92 +octave> m = sin(0:9)' * ones(1,10) + + ones(10,1) * cos(0:9); +octave> sqrt(sum(m.^2)) +ans = + 4.3146 3.1205 2.1932 3.2611 2.5342 2.5728 + 4.2047 3.6520 2.0852 3.0731 +@end example + + +@node Vector and Matrix References and Further Reading +@section References and Further Reading + +The block, vector and matrix objects in GSL follow the @code{valarray} +model of C++. A description of this model can be found in the following +reference, + +@itemize @asis +@item +B. Stroustrup, +@cite{The C++ Programming Language} (3rd Ed), +Section 22.4 Vector Arithmetic. +Addison-Wesley 1997, ISBN 0-201-88954-4. +@end itemize |