The functions described in this chapter provide a simple vector and matrix interface to ordinary C arrays. By writing your functions in terms of vectors and matrices you can pass a single structure containing both data and dimensions without needing additional function arguments. The structures are compatible with the vector and matrix formats used by BLAS routines.

Data types
==========

All the functions are available in versions 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 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.

The block struct
================

For consistency all memory is allocated through a gsl_block structure. The structure contains two members, 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.

The functions for allocating and deallocating blocks are defined in gsl_block.h

Block allocation
================

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.

 -- Function: gsl_block * gsl_block_alloc (size_t N)
     This function allocates memory for a block of 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 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.

 -- Function: gsl_block * gsl_block_calloc (size_t N)
     This function allocates memory for a block and initializes all the elements of the block to zero.

 -- Function: void gsl_block_free (gsl_block * B)
     This function frees the memory used by a block B previously allocated with gsl_block_alloc or gsl_block_calloc.

Reading and writing blocks
==========================

The library provides functions for reading and writing blocks to a file as binary data or formatted text.

 -- Function: int gsl_block_fwrite (FILE * STREAM, const gsl_block * B)
     This function writes the elements of the block B to the stream STREAM in binary format. The return value is 0 for success and 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.

 -- Function: int gsl_block_fread (FILE * STREAM, gsl_block * B)
     This function reads into the block B from the open stream STREAM in binary format. The block B must be preallocated with the correct length since the function uses the size of B to determine how many bytes to read. The return value is 0 for success and 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.

 -- Function: int gsl_block_fprintf (FILE * STREAM, const gsl_block * B, const char * FORMAT)
     This function writes the elements of the block B line-by-line to the stream STREAM using the format specifier FORMAT, which should be one of the %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.

 -- Function: int gsl_block_fscanf (FILE * STREAM, gsl_block * B)
     This function reads formatted data from the stream STREAM into the block B. The block B must be preallocated with the correct length since the function uses the size of B to determine how many numbers to read.  The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.

Example programs for blocks
===========================

The following program shows how to allocate a block,

     #include <stdio.h>
     #include <gsl_block.h>

     int main ()
     {
       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);
     }

Here is the output from the program,

     length of block = 100
     block data address = 0x804b0d8

The vector struct
=================

Vectors are defined by a gsl_vector structure which describes a slice of a block. A vector represents a view of an area of memory, and different vectors can be created which point to the same block. For example, it is possible to make two vectors which point to the even and odd elements of a third vector, or to make vectors which correspond to the rows and columns of a matrix.

The gsl_vector structure contains four members, the size, the stride, a pointer to the memory where the elements are stored, DATA, and a pointer to the block owned by the vector, BLOCK, if any. The structure is very simple and looks like this,

     typedef struct
     {
       size_t size;
       size_t stride;
       double * data;
       gsl_block * block;
     } 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.  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 owned by the vector (if any). This block will be deallocated when the vector is freed.  If the vector is only a view of another object then the pointer BLOCK is null.

The functions for allocating and accessing vectors are defined in gsl_vector.h

Vector allocation
=================

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.

 -- Function: gsl_vector * gsl_vector_alloc (size_t N)
     This function creates a vector of length 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 BLOCK member of the vector struct. The block is ``owned'' by the vector, and will be deallocated when the vector is deallocated.

 -- Function: gsl_vector * gsl_vector_calloc (size_t N)
     This function allocates memory for a vector of length N and initializes all the elements of the vector to zero.

 -- Function: gsl_vector * gsl_vector_alloc_from_block (gsl_block * B, size_t OFFSET, size_t N, size_t STRIDE)
     This function creates a vector as a slice of an existing block B, returning a pointer to a newly initialized vector struct. The start of the vector is offset by OFFSET elements from the start of the block. The vector has N elements, with a step-size of STRIDE from one element to the next.  Mathematically, the I-th element of the vector is given by,

          v(i) = b->data[offset + i*stride]

     where the index I runs from 0 to n-1.

     A null pointer is returned if the combined parameters (OFFSET,N,STRIDE) overrun the end of the block or if insufficient memory is available to store the vector.

     The vector is only a view of the block B, and the block is not owned by the vector. When the vector is deallocated the block B will continue to exist. This memory can only be deallocated by freeing the block itself. Of course, this block should not be deallocated while the vector is still in use.

 -- Function: gsl_vector * gsl_vector_alloc_from_vector (gsl_vector * V, size_t OFFSET, size_t N, size_t STRIDE)
     This function creates a vector as a slice of another vector V, returning a pointer to a newly initialized vector struct. The start of the new vector is offset by OFFSET elements from the start of the original vector. 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.

     A null pointer is returned if the combined parameters (OFFSET,N,STRIDE) overrun the end of the original vector or if insufficient memory is available store the new vector.

     The new vector is only a view of the block underlying the original vector, V. The block is not owned by the new vector. When the new vector is deallocated 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 new vector is still in use.

 -- Function: void gsl_vector_free (gsl_vector * V)
     This function frees a previously allocated vector V. If the vector was created using gsl_vector_alloc then the block underlying the vector will also be deallocated.  If the vector has been created from another object (using gsl_vector_alloc_from_block or gsl_vector_alloc_from_vector) then the memory is still owned by that object and will not be deallocated.

Accessing vector elements
=========================

Unlike FORTRAN compilers, C compilers do not usually provide support for range checking of vectors and matrices. However, the functions gsl_vector_get and gsl_vector_set can perform 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. 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 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.

 -- Function: double gsl_vector_get (const gsl_vector * V, size_t I)
     This function returns the I-th element of a vector V. If I lies outside the allowed range of 0 to N-1 then the error handler is invoked and 0 is returned.

 -- Function: void gsl_vector_set (gsl_vector * V, size_t I, double X)
     This function sets the value of the I-th element of a vector V to X. If I lies outside the allowed range of 0 to N-1 then the error handler is invoked.

Reading and writing vectors
===========================

The library provides functions for reading and writing vectors to a file as binary data or formatted text.

 -- Function: int gsl_vector_fwrite (FILE * STREAM, const gsl_vector * V)
     This function writes the elements of the vector V to the stream STREAM in binary format. The return value is 0 for success and 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.

 -- Function: int gsl_vector_fread (FILE * STREAM, gsl_vector * V)
     This function reads into the vector V from the open stream STREAM in binary format. The vector V must be preallocated with the correct length since the function uses the size of V to determine how many bytes to read. The return value is 0 for success and 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.

 -- Function: int gsl_vector_fprintf (FILE * STREAM, const gsl_vector * V, const char * FORMAT)
     This function writes the elements of the vector V line-by-line to the stream STREAM using the format specifier FORMAT, which should be one of the %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.

 -- Function: int gsl_vector_fscanf (FILE * STREAM, gsl_vector * V)
     This function reads formatted data from the stream STREAM into the vector V. The vector V must be preallocated with the correct length since the function uses the size of V to determine how many numbers to read.  The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.

Copying vectors
===============

Common operations on vectors such as addition and multiplication are available in the BLAS part of the library (see BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following function falls into this category.

 -- Function: int gsl_vector_copy (gsl_vector * DEST, const gsl_vector * SRC)
     This function copies the elements of the vector SRC into the vector DEST. The two vectors must be the same length.

Exchanging elements
===================

The following function can be used to exchange, or permute, the elements of a vector.

 -- Function: int gsl_vector_swap (gsl_vector * V, size_t i, size_t j)
     This function exchanges the I-th and J-th elements of the vector V in-place.

Vector properties
=================

 -- Function: int gsl_vector_isnull (gsl_vector * V)
     This function returns 1 if all the elements of the vector V are zero, and 0 otherwise.

Example programs for 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_vector.h>

     int main ()
     {
       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));
         }
     }

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

     int main ()
     {
       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);
       }
     }

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

The matrix struct
=================

Matrices are defined by a gsl_matrix structure which describes a generalized slice of a block. Like a vector it represents a view of an area of memory, but uses two indices instead of one.

The gsl_matrix structure contains five members, the two dimensions of the matrix, a physical dimension, a pointer to the memory where the elements of the matrix are stored, DATA, and a pointer to the block owned by the matrix BLOCK, if any.  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 dim2;
       double * data;
       gsl_block * block;
     } gsl_matrix ;

Matrices are stored in row-major order, meaning that each row of elements forms a contiguous block in memory.  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 DIM2 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 DIM2 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 .. .. .. ..
     10 11 12 13 .. .. .. ..
     20 21 22 23 .. .. .. ..

Each unused memory location is represented by ``..''.  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 owned by the matrix (if any). This block will be deallocated when the matrix is freed.  If the matrix is only a view of another object then the pointer BLOCK is null.

The functions for allocating and accessing matrices are defined in gsl_matrix.h

Matrix allocation
=================

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.

 -- Function: gsl_matrix * gsl_matrix_alloc (size_t N1, size_t N2)
     This function creates a matrix of size N1 rows by 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 BLOCK member of the matrix struct. The block is ``owned'' by the matrix, and will be deallocated when the matrix is deallocated.

 -- Function: gsl_matrix * gsl_matrix_calloc (size_t N1, size_t N2)
     This function allocates memory for a matrix of size N1 rows by N2 columns and initializes all the elements of the matrix to zero.

 -- Function: gsl_matrix * gsl_matrix_alloc_from_block (gsl_block * B, size_t OFFSET, size_t N1, size_t N2, size_t DIM2)
     This function creates a matrix as a slice of the block B, returning a pointer to a newly initialized matrix struct. The start of the matrix is offset by OFFSET elements from the start of the block. The matrix has N1 rows and N2 columns, with the physical number of columns in memory given by DIM2. Mathematically, the (I,J)-th element of the matrix is given by,

          m(i,j) = b->data[offset + i*dim2 + j]

     where the index I runs from 0 to n1-1 and the index J runs from 0 to n2-1.

     A null pointer is returned if the combined parameters (OFFSET,N1,N2,DIM2) overrun the end of the block or if insufficient memory is available to store the matrix.

     The matrix is only a view of the block B, and the block is not owned by the matrix. When the matrix is deallocated the block B will continue to exist. This memory can only be deallocated by freeing the block itself. Of course, this block should not be deallocated while the matrix is still in use.

 -- Function: gsl_matrix * gsl_matrix_alloc_from_matrix (gsl_matrix * M, size_t K1, size_t K2, size_t N1, size_t N2)

     This function creates a matrix as a submatrix of the matrix M, returning a pointer to a newly initialized matrix struct. The upper-left element of the submatrix is the element (K1,K2) of the original matrix. The submatrix has N1 rows and N2 columns. The physical number of columns in memory given by DIM2 is unchanged.  Mathematically, the (I,J)-th element of the new matrix is given by,

          m'(i,j) = m->data[(k1*m->dim2 + k2) + i*m->dim2 + j]

     where the index I runs from 0 to n1-1 and the index J runs from 0 to n2-1.

     A null pointer is returned if the combined parameters (K1,K2,N1,N2,DIM2) overrun the end of the original matrix or if insufficient memory is available to store the matrix.

     The new matrix is only a view of the block underlying the existing matrix, M. The block is not owned by the new matrix. When the new matrix is deallocated 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 new matrix is still in use.

 -- Function: void gsl_matrix_free (gsl_matrix * M)
     This function frees a previously allocated matrix M. If the matrix was created using gsl_matrix_alloc then the block underlying the matrix will also be deallocated.  If the matrix has been created from another object (using gsl_matrix_alloc_from_block or gsl_matrix_alloc_from_matrix) then the memory is still owned by that object and will not be deallocated.

Accessing matrix elements
=========================

The functions for accessing the elements of a matrix use the same range checking system as vectors. You 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->dim2 + j]

where DIM2 is the physical row-length of the matrix.

 -- Function: double gsl_matrix_get (const gsl_matrix * M, size_t I, size_t J)
     These functions return the (I,J)th element of a matrix M. If I or J lie outside the allowed range of 0 to N1-1 and 0 to N2-1 then the error handler is invoked and 0 is returned.

 -- Function: void gsl_matrix_set (gsl_matrix * M, size_t I, size_t J, double X)
     These functions set the value of the (I,J)th element of a matrix M to X. If I or J lies outside the allowed range of 0 to N1-1 and 0 to N2-1 then the error handler is invoked.

Reading and writing matrices
============================

The library provides functions for reading and writing matrices to a file as binary data or formatted text.

 -- Function: int gsl_matrix_fwrite (FILE * STREAM, const gsl_matrix * M)
     This function writes the elements of the matrix M to the stream STREAM in binary format. The return value is 0 for success and 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.

 -- Function: int gsl_matrix_fread (FILE * STREAM, gsl_matrix * M)
     This function reads into the matrix M from the open stream STREAM in binary format. The matrix M must be preallocated with the correct length since the function uses the size of M to determine how many bytes to read. The return value is 0 for success and 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.

 -- Function: int gsl_matrix_fprintf (FILE * STREAM, const gsl_matrix * M, const char * FORMAT)
     This function writes the elements of the matrix M line-by-line to the stream STREAM using the format specifier FORMAT, which should be one of the %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.

 -- Function: int gsl_matrix_fscanf (FILE * STREAM, gsl_matrix * M)
     This function reads formatted data from the stream STREAM into the matrix M. The matrix M must be preallocated with the correct length since the function uses the size of M to determine how many numbers to read.  The function returns 0 for success and GSL_EFAILED if there was a problem reading from the file.

Slicing rows and columns
========================

In general there are two ways to access an object, by reference or by copying. The functions described in this section create vectors which allow access to a row or column of a matrix by reference. Modifying elements of the vector is equivalent to modifying the matrix, since both the vector and the matrix point to the same memory block.

 -- Function: gsl_vector * gsl_vector_alloc_row_from_matrix (gsl_matrix * M, size_t I)
     This function allocates a new gsl_vector struct which points to the I-th row of the matrix M.

 -- Function: gsl_vector * gsl_vector_alloc_col_from_matrix (gsl_matrix * M, size_t J)
     This function allocates a new gsl_vector struct which points to the J-th column of the matrix M.

Copying matrices
================

Common operations on matrices such as addition and multiplication are available in the BLAS part of the library (see BLAS Support). However, it is useful to have a small number of utility functions which do not require the full BLAS code. The following function falls into this category.

 -- Function: int gsl_matrix_copy (gsl_matrix * DEST, const gsl_matrix * SRC)
     This function copies the elements of the matrix SRC into the matrix DEST. The two matrix must have the same size.

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.

 -- Function: int gsl_matrix_copy_row (gsl_vector * V, const gsl_matrix * M, size_t I)
     This function copies the elements of the I-th row of the matrix M into the vector V. The length of the vector must be the same as the length of the row.

 -- Function: int gsl_matrix_copy_col (gsl_vector * V, const gsl_matrix * M, size_t J)
     This function copies the elements of the I-th column of the matrix M into the vector V. The length of the vector must be the same as the length of the column.

 -- Function: int gsl_matrix_set_row (gsl_matrix * M, size_t I, const gsl_vector * V)
     This function copies the elements of the vector V into the I-th row of the matrix M. The length of the vector must be the same as the length of the row.

 -- Function: int gsl_matrix_set_col (gsl_matrix * M, size_t J, const gsl_vector * V)
     This function copies the elements of the vector V into the I-th column of the matrix M. The length of the vector must be the same as the length of the column.

Exchanging rows and columns
===========================

The following functions can be used to exchange the rows and columns of a matrix.

 -- Function: int gsl_matrix_swap_rows (gsl_matrix * M, size_t I, size_t J)
     This function exchanges the I-th and J-th rows of the matrix M in-place.

 -- Function: int gsl_matrix_swap_cols (gsl_matrix * M, size_t I, size_t J)
     This function exchanges the I-th and J-th columns of the matrix M in-place.

 -- Function: int gsl_matrix_swap_rowcol (gsl_matrix * M, size_t I, size_t J)
     This function exchanges the I-th row and J-th column of the matrix M in-place. The matrix must be square for this operation to be possible.

Example programs for matrices
=============================

This program 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_matrix.h>

     int main ()
     {
       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 < 100; i++)
         for (j = 0; j < 3; j++)
           printf("m(%d,%d) = %g\n", i, j, gsl_matrix_get (m, i, j));
     }

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

     int main ()
     {
       int i, j, differences = 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", "w") ;
          gsl_matrix_fwrite (f, m);
          fclose (f);
       }

       {
          FILE * f = fopen("test.dat", "r") ;
          gsl_matrix_fread (f, a);
          fclose (f);
       }

       for (i = 0; i < 100; i++)
         for (j = 0 ; j < 100; j++)
             if (gsl_matrix_get(m, i, j) != gsl_matrix_get(a, i, j))
                differences ++ ;

       printf("differences = %d (should be zero)\n", differences) ;

     }

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.