Go to the first, previous, next, last section, table of contents.
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 BLAS routines.
All the functions are available for each of the standard data-types.
The versions for double have the prefix gsl_block,
gsl_vector and gsl_matrix. Similarly the versions for
single-precision float arrays have the prefix
gsl_block_float, gsl_vector_float and
gsl_matrix_float. The full list of available types is given
below,
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
Corresponding types exist for the gsl_vector and
gsl_matrix functions.
For consistency all memory is allocated through a gsl_block
structure. The structure contains two components, the size of an area of
memory and a pointer to the memory. The gsl_block structure looks
like this,
typedef struct
{
size_t size;
double * data;
} gsl_block;
Vectors and matrices are made by 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 stride.
The functions for allocating and deallocating blocks are defined in 'gsl_block.h'
The functions for allocating memory to a block follow the style of
malloc and 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 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 alloc.
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.
gsl_block_alloc or gsl_block_calloc.
The library provides functions for reading and writing blocks to a file as binary data or formatted text.
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.
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.
%g, %e or %f formats for
floating point numbers and %d for integers. The function returns
0 for success and GSL_EFAILED if there was a problem writing to
the file.
GSL_EFAILED if there was a problem reading from the file.
The following program shows how to allocate a block,
#include <stdio.h>
#include <gsl/gsl_block.h>
int
main (void)
{
gsl_block * b = gsl_block_alloc (100);
printf ("length of block = %u\n", b->size);
printf ("block data address = %#x\n", b->data);
gsl_block_free (b);
return 0;
}
Here is the output from the program,
length of block = 100 block data address = 0x804b0d8
Vectors are defined by a 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 gsl_vector structure contains five components, the
size, the stride, a pointer to the memory where the elements
are stored, data, a pointer to the block owned by the vector,
block, if any, and an ownership flag, owner. The structure
is very simple and looks like this,
typedef struct
{
size_t size;
size_t stride;
double * data;
gsl_block * block;
int owner;
} gsl_vector;
The size is simply the number of vector elements. The range of
valid indices runs from 0 to size-1. The stride is the
step-size from one element to the next in physical memory, measured in
units of the appropriate datatype. The pointer data gives the
location of the first element of the vector in memory. The pointer
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
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 owner field is zero and any underlying block will not be
deallocated.
The functions for allocating and accessing vectors are defined in 'gsl_vector.h'
The functions for allocating memory to a vector follow the style of
malloc and 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 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 alloc.
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.
Unlike 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 extension checkergcc but it is
not available on every platform. The functions gsl_vector_get
and 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 'gsl_vector.h' and declared extern inline to
eliminate function-call overhead. You must compile your program with
the macro 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 GSL_RANGE_CHECK_OFF. Provided your compiler supports
inline functions the effect of turning off range checking is to replace
calls to gsl_vector_get(v,i) by v->data[i*v->stride] and
calls to gsl_vector_set(v,i,x) by 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.
The library provides functions for reading and writing vectors to a file as binary data or formatted text.
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.
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.
%g, %e or %f formats for
floating point numbers and %d for integers. The function returns
0 for success and GSL_EFAILED if there was a problem writing to
the file.
GSL_EFAILED if there was a problem reading from the file.
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
gsl_vector_view and a constant vector view has the type
gsl_vector_const_view. In both cases the elements of the view
can be accessed as a gsl_vector using the vector component
of the view object. A pointer to a vector of type gsl_vector *
or const gsl_vector * can be obtained by taking the address of
this component with the & operator.
v'(i) = v->data[(offset + i)*v->stride]
where the index i runs from 0 to n-1.
The data pointer of the returned vector struct is set to null if
the combined parameters (offset,n) overrun the end of the
original vector.
The new vector is only a view of the block underlying the original vector, v. The block containing the elements of v is not owned by the new vector. When the view goes out of scope the original vector 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 gsl_vector_const_subvector is equivalent to
gsl_vector_subvector but can be used for vectors which are
declared const.
gsl_vector_subvector but the new vector has
n elements with a step-size of stride from one element to
the next in the original vector. Mathematically, the i-th element
of the new vector v' is given by,
v'(i) = v->data[(offset + i*stride)*v->stride]
where the index i runs from 0 to 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 v of length n, while leaving the
odd elements untouched,
gsl_vector_view v_even = gsl_vector_subvector_with_stride (v, 0, 2, n/2); gsl_vector_set_zero (&v_even.vector);
A vector view can be passed to any subroutine which takes a vector
argument just as a directly allocated vector would be, using
&view.vector. For example, the following code
computes the norm of odd elements of v using the BLAS
routine DNRM2,
gsl_vector_view v_odd = gsl_vector_subvector_with_stride (v, 1, 2, n/2); double r = gsl_blas_dnrm2 (&v_odd.vector);
The function gsl_vector_const_subvector_with_stride is equivalent
to gsl_vector_subvector_with_stride but can be used for vectors
which are declared const.
The function gsl_vector_complex_const_real is equivalent to
gsl_vector_complex_real but can be used for vectors which are
declared const.
The function gsl_vector_complex_const_imag is equivalent to
gsl_vector_complex_imag but can be used for vectors which are
declared const.
v'(i) = base[i]
where the index i runs from 0 to n-1.
The array containing the elements of 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 base. Of course, the original array should not be deallocated while the view is still in use.
The function gsl_vector_const_view_array is equivalent to
gsl_vector_view_array but can be used for arrays which are
declared const.
gsl_vector_view_array but the new vector has n elements
with a step-size of stride from one element to the next in the
original array. Mathematically, the i-th element of the new
vector v' is given by,
v'(i) = base[i*stride]
where the index i runs from 0 to 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 &view.vector.
The function gsl_vector_const_view_array_with_stride is
equivalent to gsl_vector_view_array_with_stride but can be used
for arrays which are declared const.
Common operations on vectors such as addition and multiplication are available in the BLAS part of the library (see section BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following functions fall into this category.
The following function can be used to exchange, or permute, the elements of a vector.
The following operations are only defined for real vectors.
This program shows how to allocate, initialize and read from a vector
using the functions gsl_vector_alloc, gsl_vector_set and
gsl_vector_get.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (3);
for (i = 0; i < 3; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
for (i = 0; i < 100; i++)
{
printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
}
return 0;
}
Here is the output from the program. The final loop attempts to read
outside the range of the vector v, and the error is trapped by
the range-checking code in gsl_vector_get.
v_0 = 1.23 v_1 = 2.23 v_2 = 3.23 gsl: vector_source.c:12: ERROR: index out of range IOT trap/Abort (core dumped)
The next program shows how to write a vector to a file.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (100);
for (i = 0; i < 100; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, v, "%.5g");
fclose (f);
}
return 0;
}
After running this program the file 'test.dat' should contain the
elements of v, written using the format specifier
%.5g. The vector could then be read back in using the function
gsl_vector_fscanf (f, v) as follows:
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (10);
{
FILE * f = fopen ("test.dat", "r");
gsl_vector_fscanf (f, v);
fclose (f);
}
for (i = 0; i < 10; i++)
{
printf ("%g\n", gsl_vector_get(v, i));
}
return 0;
}
Matrices are defined by a 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 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, data, a pointer to
the block owned by the matrix block, if any, and an ownership
flag, owner. The physical dimension determines the memory layout
and can differ from the matrix dimension to allow the use of
submatrices. The gsl_matrix structure is very simple and looks
like this,
typedef struct
{
size_t size1;
size_t size2;
size_t tda;
double * data;
gsl_block * block;
int owner;
} gsl_matrix;
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 FORTRAN
stores arrays in column-major order. The number of rows is size1.
The range of valid row indices runs from 0 to size1-1. Similarly
size2 is the number of columns. The range of valid column indices
runs from 0 to size2-1. The physical row dimension tda, or
trailing dimension, specifies the size of a row of the matrix as
laid out in memory.
For example, in the following matrix size1 is 3, size2 is 4, and 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.
00 01 02 03 XX XX XX XX 10 11 12 13 XX XX XX XX 20 21 22 23 XX XX XX XX
Each unused memory location is represented by "XX". The
pointer data gives the location of the first element of the matrix
in memory. The pointer 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 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 owner field is
zero and any underlying block will not be freed.
The functions for allocating and accessing matrices are defined in 'gsl_matrix.h'
The functions for allocating memory to a matrix follow the style of
malloc and 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
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 alloc.
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 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
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 gsl_matrix_get(m,i,j) and
gsl_matrix_set(m,i,j,x) is
m->data[i * m->tda + j]
where tda is the physical row-length of the matrix.
The library provides functions for reading and writing matrices to a file as binary data or formatted text.
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.
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.
%g, %e or %f formats for
floating point numbers and %d for integers. The function returns
0 for success and GSL_EFAILED if there was a problem writing to
the file.
GSL_EFAILED if there was a problem reading from the file.
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
gsl_matrix_view and a constant matrix view has the type
gsl_matrix_const_view. In both cases the elements of the view
can by accessed using the matrix component of the view object. A
pointer gsl_matrix * or const gsl_matrix * can be obtained
by taking the address of the matrix component with the &
operator. In addition to matrix views it is also possible to create
vector views of a matrix, such as row or column views.
m'(i,j) = m->data[(k1*m->tda + k1) + i*m->tda + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The data pointer of the returned matrix struct is set to null if
the combined parameters (i,j,n1,n2,tda)
overrun the ends of the original matrix.
The new matrix view is only a view of the block underlying the existing matrix, m. The block containing the elements of m is not owned by the new matrix view. When the view goes out of scope the original matrix 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 gsl_matrix_const_submatrix is equivalent to
gsl_matrix_submatrix but can be used for matrices which are
declared const.
m'(i,j) = base[i*n2 + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the array base. When the view goes out of scope the original array 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 gsl_matrix_const_view_array is equivalent to
gsl_matrix_view_array but can be used for matrices which are
declared const.
m'(i,j) = base[i*tda + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the array base. When the view goes out of scope the original array 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 gsl_matrix_const_view_array_with_tda is equivalent
to gsl_matrix_view_array_with_tda but can be used for matrices
which are declared const.
m'(i,j) = v->data[i*n2 + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the vector v. When the view goes out of scope the original vector 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 gsl_matrix_const_view_vector is equivalent to
gsl_matrix_view_vector but can be used for matrices which are
declared const.
m'(i,j) = v->data[i*tda + j]
where the index i runs from 0 to n1-1 and the index j
runs from 0 to n2-1.
The new matrix is only a view of the vector v. When the view goes out of scope the original vector 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 gsl_matrix_const_view_vector_with_tda is equivalent
to gsl_matrix_view_vector_with_tda but can be used for matrices
which are declared const.
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.
data pointer of the new vector is set to null if
i is out of range.
The function gsl_vector_const_row is equivalent to
gsl_matrix_row but can be used for matrices which are declared
const.
data pointer of the new vector is set to
null if j is out of range.
The function gsl_vector_const_column equivalent to
gsl_matrix_column but can be used for matrices which are declared
const.
The function gsl_matrix_const_diagonal is equivalent to
gsl_matrix_diagonal but can be used for matrices which are
declared const.
The function gsl_matrix_const_subdiagonal is equivalent to
gsl_matrix_subdiagonal but can be used for matrices which are
declared const.
The function gsl_matrix_const_superdiagonal is equivalent to
gsl_matrix_superdiagonal but can be used for matrices which are
declared const.
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
gsl_vector_memcpy with vector views of rows and columns.
The following functions can be used to exchange the rows and columns of a matrix.
The following operations are defined for real and complex matrices.
The following operations are only defined for real matrices.
The program below shows how to allocate, initialize and read from a matrix
using the functions gsl_matrix_alloc, gsl_matrix_set and
gsl_matrix_get.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int
main (void)
{
int i, j;
gsl_matrix * m = gsl_matrix_alloc (10, 3);
for (i = 0; i < 10; i++)
for (j = 0; j < 3; j++)
gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
for (i = 0; i < 10; i++)
for (j = 0; j < 3; j++)
printf ("m(%d,%d) = %g\n", i, j,
gsl_matrix_get (m, i, j));
return 0;
}
Here is the output from the program. The final loop attempts to read
outside the range of the matrix m, and the error is trapped by
the range-checking code in gsl_matrix_get.
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 IOT trap/Abort (core dumped)
The next program shows how to write a matrix to a file.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int
main (void)
{
int i, j, k = 0;
gsl_matrix * m = gsl_matrix_alloc (100, 100);
gsl_matrix * a = gsl_matrix_alloc (100, 100);
for (i = 0; i < 100; i++)
for (j = 0; j < 100; j++)
gsl_matrix_set (m, i, j, 0.23 + i + j);
{
FILE * f = fopen ("test.dat", "wb");
gsl_matrix_fwrite (f, m);
fclose (f);
}
{
FILE * f = fopen ("test.dat", "rb");
gsl_matrix_fread (f, a);
fclose (f);
}
for (i = 0; i < 100; i++)
for (j = 0; j < 100; j++)
{
double mij = gsl_matrix_get (m, i, j);
double aij = gsl_matrix_get (a, i, j);
if (mij != aij) k++;
}
printf ("differences = %d (should be zero)\n", k);
return (k > 0);
}
After running this program the file 'test.dat' should contain the
elements of m, written in binary format. The matrix which is read
back in using the function 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.
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
int
main (void)
{
size_t i,j;
gsl_matrix *m = gsl_matrix_alloc (10, 10);
for (i = 0; i < 10; i++)
for (j = 0; j < 10; j++)
gsl_matrix_set (m, i, j, sin (i) + cos (j));
for (j = 0; j < 10; j++)
{
gsl_vector_view column = gsl_matrix_column (m, j);
double d;
d = gsl_blas_dnrm2 (&column.vector);
printf ("matrix column %d, norm = %g\n", j, d);
}
gsl_matrix_free (m);
return 0;
}
Here is the output of the program, which can be confirmed using GNU OCTAVE,
$ ./a.out matrix column 0, norm = 4.31461 matrix column 1, norm = 3.1205 matrix column 2, norm = 2.19316 matrix column 3, norm = 3.26114 matrix column 4, norm = 2.53416 matrix column 5, norm = 2.57281 matrix column 6, norm = 4.20469 matrix column 7, norm = 3.65202 matrix column 8, norm = 2.08524 matrix column 9, norm = 3.07313
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
The block, vector and matrix objects in GSL follow the valarray
model of C++. A description of this model can be found in the following
reference,
Go to the first, previous, next, last section, table of contents.