5 The functions described in this chapter provide a simple vector and
6 matrix interface to ordinary C arrays. The memory management of these
7 arrays is implemented using a single underlying type, known as a
8 block. By writing your functions in terms of vectors and matrices you
9 can pass a single structure containing both data and dimensions as an
10 argument without needing additional function parameters. The structures
11 are compatible with the vector and matrix formats used by @sc{blas}
19 * Vector and Matrix References and Further Reading::
25 All the functions are available for each of the standard data-types.
26 The versions for @code{double} have the prefix @code{gsl_block},
27 @code{gsl_vector} and @code{gsl_matrix}. Similarly the versions for
28 single-precision @code{float} arrays have the prefix
29 @code{gsl_block_float}, @code{gsl_vector_float} and
30 @code{gsl_matrix_float}. The full list of available types is given
36 gsl_block_long_double long double
38 gsl_block_uint unsigned int
40 gsl_block_ulong unsigned long
42 gsl_block_ushort unsigned short
44 gsl_block_uchar unsigned char
45 gsl_block_complex complex double
46 gsl_block_complex_float complex float
47 gsl_block_complex_long_double complex long double
51 Corresponding types exist for the @code{gsl_vector} and
52 @code{gsl_matrix} functions.
59 For consistency all memory is allocated through a @code{gsl_block}
60 structure. The structure contains two components, the size of an area of
61 memory and a pointer to the memory. The @code{gsl_block} structure looks
74 Vectors and matrices are made by @dfn{slicing} an underlying block. A
75 slice is a set of elements formed from an initial offset and a
76 combination of indices and step-sizes. In the case of a matrix the
77 step-size for the column index represents the row-length. The step-size
78 for a vector is known as the @dfn{stride}.
80 The functions for allocating and deallocating blocks are defined in
85 * Reading and writing blocks::
86 * Example programs for blocks::
89 @node Block allocation
90 @subsection Block allocation
92 The functions for allocating memory to a block follow the style of
93 @code{malloc} and @code{free}. In addition they also perform their own
94 error checking. If there is insufficient memory available to allocate a
95 block then the functions call the GSL error handler (with an error
96 number of @code{GSL_ENOMEM}) in addition to returning a null
97 pointer. Thus if you use the library error handler to abort your program
98 then it isn't necessary to check every @code{alloc}.
100 @deftypefun {gsl_block *} gsl_block_alloc (size_t @var{n})
101 This function allocates memory for a block of @var{n} double-precision
102 elements, returning a pointer to the block struct. The block is not
103 initialized and so the values of its elements are undefined. Use the
104 function @code{gsl_block_calloc} if you want to ensure that all the
105 elements are initialized to zero.
107 A null pointer is returned if insufficient memory is available to create
111 @deftypefun {gsl_block *} gsl_block_calloc (size_t @var{n})
112 This function allocates memory for a block and initializes all the
113 elements of the block to zero.
116 @deftypefun void gsl_block_free (gsl_block * @var{b})
117 This function frees the memory used by a block @var{b} previously
118 allocated with @code{gsl_block_alloc} or @code{gsl_block_calloc}. The
119 block @var{b} must be a valid block object (a null pointer is not
123 @node Reading and writing blocks
124 @subsection Reading and writing blocks
126 The library provides functions for reading and writing blocks to a file
127 as binary data or formatted text.
129 @deftypefun int gsl_block_fwrite (FILE * @var{stream}, const gsl_block * @var{b})
130 This function writes the elements of the block @var{b} to the stream
131 @var{stream} in binary format. The return value is 0 for success and
132 @code{GSL_EFAILED} if there was a problem writing to the file. Since the
133 data is written in the native binary format it may not be portable
134 between different architectures.
137 @deftypefun int gsl_block_fread (FILE * @var{stream}, gsl_block * @var{b})
138 This function reads into the block @var{b} from the open stream
139 @var{stream} in binary format. The block @var{b} must be preallocated
140 with the correct length since the function uses the size of @var{b} to
141 determine how many bytes to read. The return value is 0 for success and
142 @code{GSL_EFAILED} if there was a problem reading from the file. The
143 data is assumed to have been written in the native binary format on the
147 @deftypefun int gsl_block_fprintf (FILE * @var{stream}, const gsl_block * @var{b}, const char * @var{format})
148 This function writes the elements of the block @var{b} line-by-line to
149 the stream @var{stream} using the format specifier @var{format}, which
150 should be one of the @code{%g}, @code{%e} or @code{%f} formats for
151 floating point numbers and @code{%d} for integers. The function returns
152 0 for success and @code{GSL_EFAILED} if there was a problem writing to
156 @deftypefun int gsl_block_fscanf (FILE * @var{stream}, gsl_block * @var{b})
157 This function reads formatted data from the stream @var{stream} into the
158 block @var{b}. The block @var{b} must be preallocated with the correct
159 length since the function uses the size of @var{b} to determine how many
160 numbers to read. The function returns 0 for success and
161 @code{GSL_EFAILED} if there was a problem reading from the file.
165 @node Example programs for blocks
166 @subsection Example programs for blocks
168 The following program shows how to allocate a block,
171 @verbatiminclude examples/block.c
176 Here is the output from the program,
179 @verbatiminclude examples/block.out
186 @cindex stride, of vector index
188 Vectors are defined by a @code{gsl_vector} structure which describes a
189 slice of a block. Different vectors can be created which point to the
190 same block. A vector slice is a set of equally-spaced elements of an
193 The @code{gsl_vector} structure contains five components, the
194 @dfn{size}, the @dfn{stride}, a pointer to the memory where the elements
195 are stored, @var{data}, a pointer to the block owned by the vector,
196 @var{block}, if any, and an ownership flag, @var{owner}. The structure
197 is very simple and looks like this,
212 The @var{size} is simply the number of vector elements. The range of
213 valid indices runs from 0 to @code{size-1}. The @var{stride} is the
214 step-size from one element to the next in physical memory, measured in
215 units of the appropriate datatype. The pointer @var{data} gives the
216 location of the first element of the vector in memory. The pointer
217 @var{block} stores the location of the memory block in which the vector
218 elements are located (if any). If the vector owns this block then the
219 @var{owner} field is set to one and the block will be deallocated when the
220 vector is freed. If the vector points to a block owned by another
221 object then the @var{owner} field is zero and any underlying block will not be
222 deallocated with the vector.
224 The functions for allocating and accessing vectors are defined in
228 * Vector allocation::
229 * Accessing vector elements::
230 * Initializing vector elements::
231 * Reading and writing vectors::
234 * Exchanging elements::
235 * Vector operations::
236 * Finding maximum and minimum elements of vectors::
237 * Vector properties::
238 * Example programs for vectors::
241 @node Vector allocation
242 @subsection Vector allocation
244 The functions for allocating memory to a vector follow the style of
245 @code{malloc} and @code{free}. In addition they also perform their own
246 error checking. If there is insufficient memory available to allocate a
247 vector then the functions call the GSL error handler (with an error
248 number of @code{GSL_ENOMEM}) in addition to returning a null
249 pointer. Thus if you use the library error handler to abort your program
250 then it isn't necessary to check every @code{alloc}.
252 @deftypefun {gsl_vector *} gsl_vector_alloc (size_t @var{n})
253 This function creates a vector of length @var{n}, returning a pointer to
254 a newly initialized vector struct. A new block is allocated for the
255 elements of the vector, and stored in the @var{block} component of the
256 vector struct. The block is ``owned'' by the vector, and will be
257 deallocated when the vector is deallocated.
260 @deftypefun {gsl_vector *} gsl_vector_calloc (size_t @var{n})
261 This function allocates memory for a vector of length @var{n} and
262 initializes all the elements of the vector to zero.
265 @deftypefun void gsl_vector_free (gsl_vector * @var{v})
266 This function frees a previously allocated vector @var{v}. If the
267 vector was created using @code{gsl_vector_alloc} then the block
268 underlying the vector will also be deallocated. If the vector has
269 been created from another object then the memory is still owned by
270 that object and will not be deallocated. The vector @var{v} must be a
271 valid vector object (a null pointer is not allowed).
274 @node Accessing vector elements
275 @subsection Accessing vector elements
276 @cindex vectors, range-checking
277 @cindex range-checking for vectors
278 @cindex bounds checking, extension to GCC
279 @cindex gcc extensions, range-checking
280 @cindex Fortran range checking, equivalent in gcc
282 Unlike @sc{fortran} compilers, C compilers do not usually provide
283 support for range checking of vectors and matrices. Range checking is
284 available in the GNU C Compiler bounds-checking extension, but it is not
285 part of the default installation of GCC. The functions
286 @code{gsl_vector_get} and @code{gsl_vector_set} can perform portable
287 range checking for you and report an error if you attempt to access
288 elements outside the allowed range.
290 The functions for accessing the elements of a vector or matrix are
291 defined in @file{gsl_vector.h} and declared @code{extern inline} to
292 eliminate function-call overhead. You must compile your program with
293 the macro @code{HAVE_INLINE} defined to use these functions.
295 If necessary you can turn off range checking completely without
296 modifying any source files by recompiling your program with the
297 preprocessor definition @code{GSL_RANGE_CHECK_OFF}. Provided your
298 compiler supports inline functions the effect of turning off range
299 checking is to replace calls to @code{gsl_vector_get(v,i)} by
300 @code{v->data[i*v->stride]} and calls to @code{gsl_vector_set(v,i,x)} by
301 @code{v->data[i*v->stride]=x}. Thus there should be no performance
302 penalty for using the range checking functions when range checking is
305 @deftypefun double gsl_vector_get (const gsl_vector * @var{v}, size_t @var{i})
306 This function returns the @var{i}-th element of a vector @var{v}. If
307 @var{i} lies outside the allowed range of 0 to @math{@var{n}-1} then the error
308 handler is invoked and 0 is returned. @inlinefn{}
311 @deftypefun void gsl_vector_set (gsl_vector * @var{v}, size_t @var{i}, double @var{x})
312 This function sets the value of the @var{i}-th element of a vector
313 @var{v} to @var{x}. If @var{i} lies outside the allowed range of 0 to
314 @math{@var{n}-1} then the error handler is invoked. @inlinefn{}
317 @deftypefun {double *} gsl_vector_ptr (gsl_vector * @var{v}, size_t @var{i})
318 @deftypefunx {const double *} gsl_vector_const_ptr (const gsl_vector * @var{v}, size_t @var{i})
319 These functions return a pointer to the @var{i}-th element of a vector
320 @var{v}. If @var{i} lies outside the allowed range of 0 to @math{@var{n}-1}
321 then the error handler is invoked and a null pointer is returned. @inlinefns{}
324 @node Initializing vector elements
325 @subsection Initializing vector elements
326 @cindex vectors, initializing
327 @cindex initializing vectors
329 @deftypefun void gsl_vector_set_all (gsl_vector * @var{v}, double @var{x})
330 This function sets all the elements of the vector @var{v} to the value
334 @deftypefun void gsl_vector_set_zero (gsl_vector * @var{v})
335 This function sets all the elements of the vector @var{v} to zero.
338 @deftypefun int gsl_vector_set_basis (gsl_vector * @var{v}, size_t @var{i})
339 This function makes a basis vector by setting all the elements of the
340 vector @var{v} to zero except for the @var{i}-th element which is set to
344 @node Reading and writing vectors
345 @subsection Reading and writing vectors
347 The library provides functions for reading and writing vectors to a file
348 as binary data or formatted text.
350 @deftypefun int gsl_vector_fwrite (FILE * @var{stream}, const gsl_vector * @var{v})
351 This function writes the elements of the vector @var{v} to the stream
352 @var{stream} in binary format. The return value is 0 for success and
353 @code{GSL_EFAILED} if there was a problem writing to the file. Since the
354 data is written in the native binary format it may not be portable
355 between different architectures.
358 @deftypefun int gsl_vector_fread (FILE * @var{stream}, gsl_vector * @var{v})
359 This function reads into the vector @var{v} from the open stream
360 @var{stream} in binary format. The vector @var{v} must be preallocated
361 with the correct length since the function uses the size of @var{v} to
362 determine how many bytes to read. The return value is 0 for success and
363 @code{GSL_EFAILED} if there was a problem reading from the file. The
364 data is assumed to have been written in the native binary format on the
368 @deftypefun int gsl_vector_fprintf (FILE * @var{stream}, const gsl_vector * @var{v}, const char * @var{format})
369 This function writes the elements of the vector @var{v} line-by-line to
370 the stream @var{stream} using the format specifier @var{format}, which
371 should be one of the @code{%g}, @code{%e} or @code{%f} formats for
372 floating point numbers and @code{%d} for integers. The function returns
373 0 for success and @code{GSL_EFAILED} if there was a problem writing to
377 @deftypefun int gsl_vector_fscanf (FILE * @var{stream}, gsl_vector * @var{v})
378 This function reads formatted data from the stream @var{stream} into the
379 vector @var{v}. The vector @var{v} must be preallocated with the correct
380 length since the function uses the size of @var{v} to determine how many
381 numbers to read. The function returns 0 for success and
382 @code{GSL_EFAILED} if there was a problem reading from the file.
386 @subsection Vector views
388 In addition to creating vectors from slices of blocks it is also
389 possible to slice vectors and create vector views. For example, a
390 subvector of another vector can be described with a view, or two views
391 can be made which provide access to the even and odd elements of a
394 A vector view is a temporary object, stored on the stack, which can be
395 used to operate on a subset of vector elements. Vector views can be
396 defined for both constant and non-constant vectors, using separate types
397 that preserve constness. A vector view has the type
398 @code{gsl_vector_view} and a constant vector view has the type
399 @code{gsl_vector_const_view}. In both cases the elements of the view
400 can be accessed as a @code{gsl_vector} using the @code{vector} component
401 of the view object. A pointer to a vector of type @code{gsl_vector *}
402 or @code{const gsl_vector *} can be obtained by taking the address of
403 this component with the @code{&} operator.
405 When using this pointer it is important to ensure that the view itself
406 remains in scope---the simplest way to do so is by always writing the
407 pointer as @code{&}@var{view}@code{.vector}, and never storing this value
410 @deftypefun gsl_vector_view gsl_vector_subvector (gsl_vector * @var{v}, size_t @var{offset}, size_t @var{n})
411 @deftypefunx {gsl_vector_const_view} gsl_vector_const_subvector (const gsl_vector * @var{v}, size_t @var{offset}, size_t @var{n})
412 These functions return a vector view of a subvector of another vector
413 @var{v}. The start of the new vector is offset by @var{offset} elements
414 from the start of the original vector. The new vector has @var{n}
415 elements. Mathematically, the @var{i}-th element of the new vector
416 @var{v'} is given by,
419 v'(i) = v->data[(offset + i)*v->stride]
423 where the index @var{i} runs from 0 to @code{n-1}.
425 The @code{data} pointer of the returned vector struct is set to null if
426 the combined parameters (@var{offset},@var{n}) overrun the end of the
429 The new vector is only a view of the block underlying the original
430 vector, @var{v}. The block containing the elements of @var{v} is not
431 owned by the new vector. When the view goes out of scope the original
432 vector @var{v} and its block will continue to exist. The original
433 memory can only be deallocated by freeing the original vector. Of
434 course, the original vector should not be deallocated while the view is
437 The function @code{gsl_vector_const_subvector} is equivalent to
438 @code{gsl_vector_subvector} but can be used for vectors which are
439 declared @code{const}.
443 @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})
444 @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})
445 These functions return a vector view of a subvector of another vector
446 @var{v} with an additional stride argument. The subvector is formed in
447 the same way as for @code{gsl_vector_subvector} but the new vector has
448 @var{n} elements with a step-size of @var{stride} from one element to
449 the next in the original vector. Mathematically, the @var{i}-th element
450 of the new vector @var{v'} is given by,
453 v'(i) = v->data[(offset + i*stride)*v->stride]
457 where the index @var{i} runs from 0 to @code{n-1}.
459 Note that subvector views give direct access to the underlying elements
460 of the original vector. For example, the following code will zero the
461 even elements of the vector @code{v} of length @code{n}, while leaving the
462 odd elements untouched,
465 gsl_vector_view v_even
466 = gsl_vector_subvector_with_stride (v, 0, 2, n/2);
467 gsl_vector_set_zero (&v_even.vector);
471 A vector view can be passed to any subroutine which takes a vector
472 argument just as a directly allocated vector would be, using
473 @code{&}@var{view}@code{.vector}. For example, the following code
474 computes the norm of the odd elements of @code{v} using the @sc{blas}
478 gsl_vector_view v_odd
479 = gsl_vector_subvector_with_stride (v, 1, 2, n/2);
480 double r = gsl_blas_dnrm2 (&v_odd.vector);
483 The function @code{gsl_vector_const_subvector_with_stride} is equivalent
484 to @code{gsl_vector_subvector_with_stride} but can be used for vectors
485 which are declared @code{const}.
488 @deftypefun gsl_vector_view gsl_vector_complex_real (gsl_vector_complex * @var{v})
489 @deftypefunx {gsl_vector_const_view} gsl_vector_complex_const_real (const gsl_vector_complex * @var{v})
490 These functions return a vector view of the real parts of the complex
493 The function @code{gsl_vector_complex_const_real} is equivalent to
494 @code{gsl_vector_complex_real} but can be used for vectors which are
495 declared @code{const}.
498 @deftypefun gsl_vector_view gsl_vector_complex_imag (gsl_vector_complex * @var{v})
499 @deftypefunx {gsl_vector_const_view} gsl_vector_complex_const_imag (const gsl_vector_complex * @var{v})
500 These functions return a vector view of the imaginary parts of the
501 complex vector @var{v}.
503 The function @code{gsl_vector_complex_const_imag} is equivalent to
504 @code{gsl_vector_complex_imag} but can be used for vectors which are
505 declared @code{const}.
508 @deftypefun gsl_vector_view gsl_vector_view_array (double * @var{base}, size_t @var{n})
509 @deftypefunx {gsl_vector_const_view} gsl_vector_const_view_array (const double * @var{base}, size_t @var{n})
510 These functions return a vector view of an array. The start of the new
511 vector is given by @var{base} and has @var{n} elements. Mathematically,
512 the @var{i}-th element of the new vector @var{v'} is given by,
519 where the index @var{i} runs from 0 to @code{n-1}.
521 The array containing the elements of @var{v} is not owned by the new
522 vector view. When the view goes out of scope the original array will
523 continue to exist. The original memory can only be deallocated by
524 freeing the original pointer @var{base}. Of course, the original array
525 should not be deallocated while the view is still in use.
527 The function @code{gsl_vector_const_view_array} is equivalent to
528 @code{gsl_vector_view_array} but can be used for arrays which are
529 declared @code{const}.
532 @deftypefun gsl_vector_view gsl_vector_view_array_with_stride (double * @var{base}, size_t @var{stride}, size_t @var{n})
533 @deftypefunx gsl_vector_const_view gsl_vector_const_view_array_with_stride (const double * @var{base}, size_t @var{stride}, size_t @var{n})
534 These functions return a vector view of an array @var{base} with an
535 additional stride argument. The subvector is formed in the same way as
536 for @code{gsl_vector_view_array} but the new vector has @var{n} elements
537 with a step-size of @var{stride} from one element to the next in the
538 original array. Mathematically, the @var{i}-th element of the new
539 vector @var{v'} is given by,
542 v'(i) = base[i*stride]
546 where the index @var{i} runs from 0 to @code{n-1}.
548 Note that the view gives direct access to the underlying elements of the
549 original array. A vector view can be passed to any subroutine which
550 takes a vector argument just as a directly allocated vector would be,
551 using @code{&}@var{view}@code{.vector}.
553 The function @code{gsl_vector_const_view_array_with_stride} is
554 equivalent to @code{gsl_vector_view_array_with_stride} but can be used
555 for arrays which are declared @code{const}.
559 @comment @node Modifying subvector views
560 @comment @subsection Modifying subvector views
562 @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})
563 @comment This function modifies and existing vector view @var{v} to form a new
564 @comment view of a vector @var{base}, starting from element @var{offset}. The
565 @comment vector has @var{n} elements separated by stride @var{stride}. Any
566 @comment existing view in @var{v} will be lost as a result of this function.
567 @comment @end deftypefun
569 @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})
570 @comment This function modifies and existing vector view @var{v} to form a new
571 @comment view of an array @var{base}, starting from element @var{offset}. The
572 @comment vector has @var{n} elements separated by stride @var{stride}. Any
573 @comment existing view in @var{v} will be lost as a result of this function.
574 @comment @end deftypefun
576 @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})
577 @comment This function creates a vector as a slice of an existing block @var{b},
578 @comment returning a pointer to a newly initialized vector struct. The start of
579 @comment the vector is offset by @var{offset} elements from the start of the
580 @comment block. The vector has @var{n} elements, with a step-size of @var{stride}
581 @comment from one element to the next. Mathematically, the @var{i}-th element of
582 @comment the vector is given by,
585 @comment v(i) = b->data[offset + i*stride]
586 @comment @end example
588 @comment where the index @var{i} runs from 0 to @code{n-1}.
590 @comment A null pointer is returned if the combined parameters
591 @comment (@var{offset},@var{n},@var{stride}) overrun the end of the block or if
592 @comment insufficient memory is available to store the vector.
594 @comment The vector is only a view of the block @var{b}, and the block is not
595 @comment owned by the vector. When the vector is deallocated the block @var{b}
596 @comment will continue to exist. This memory can only be deallocated by freeing
597 @comment the block itself. Of course, this block should not be deallocated while
598 @comment the vector is still in use.
599 @comment @end deftypefun
601 @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})
602 @comment This function creates a vector as a slice of another vector @var{v},
603 @comment returning a pointer to a newly initialized vector struct. The start of
604 @comment the new vector is offset by @var{offset} elements from the start of the
605 @comment original vector. The new vector has @var{n} elements, with a step-size
606 @comment of @var{stride} from one element to the next in the original vector.
607 @comment Mathematically, the @var{i}-th element of the new vector @var{v'} is
611 @comment v'(i) = v->data[(offset + i*stride)*v->stride]
612 @comment @end example
614 @comment where the index @var{i} runs from 0 to @code{n-1}.
616 @comment A null pointer is returned if the combined parameters
617 @comment (@var{offset},@var{n},@var{stride}) overrun the end of the original
618 @comment vector or if insufficient memory is available store the new vector.
620 @comment The new vector is only a view of the block underlying the original
621 @comment vector, @var{v}. The block is not owned by the new vector. When the new
622 @comment vector is deallocated the original vector @var{v} and its block will
623 @comment continue to exist. The original memory can only be deallocated by
624 @comment freeing the original vector. Of course, the original vector should not
625 @comment be deallocated while the new vector is still in use.
626 @comment @end deftypefun
628 @node Copying vectors
629 @subsection Copying vectors
631 Common operations on vectors such as addition and multiplication are
632 available in the @sc{blas} part of the library (@pxref{BLAS
633 Support}). However, it is useful to have a small number of utility
634 functions which do not require the full @sc{blas} code. The following
635 functions fall into this category.
637 @deftypefun int gsl_vector_memcpy (gsl_vector * @var{dest}, const gsl_vector * @var{src})
638 This function copies the elements of the vector @var{src} into the
639 vector @var{dest}. The two vectors must have the same length.
642 @deftypefun int gsl_vector_swap (gsl_vector * @var{v}, gsl_vector * @var{w})
643 This function exchanges the elements of the vectors @var{v} and @var{w}
644 by copying. The two vectors must have the same length.
648 @node Exchanging elements
649 @subsection Exchanging elements
651 The following function can be used to exchange, or permute, the elements
654 @deftypefun int gsl_vector_swap_elements (gsl_vector * @var{v}, size_t @var{i}, size_t @var{j})
655 This function exchanges the @var{i}-th and @var{j}-th elements of the
656 vector @var{v} in-place.
659 @deftypefun int gsl_vector_reverse (gsl_vector * @var{v})
660 This function reverses the order of the elements of the vector @var{v}.
664 @node Vector operations
665 @subsection Vector operations
667 The following operations are only defined for real vectors.
669 @deftypefun int gsl_vector_add (gsl_vector * @var{a}, const gsl_vector * @var{b})
670 This function adds the elements of vector @var{b} to the elements of
671 vector @var{a}, @math{a'_i = a_i + b_i}. The two vectors must have the
675 @deftypefun int gsl_vector_sub (gsl_vector * @var{a}, const gsl_vector * @var{b})
676 This function subtracts the elements of vector @var{b} from the elements of
677 vector @var{a}, @math{a'_i = a_i - b_i}. The two vectors must have the
681 @deftypefun int gsl_vector_mul (gsl_vector * @var{a}, const gsl_vector * @var{b})
682 This function multiplies the elements of vector @var{a} by the elements of
683 vector @var{b}, @math{a'_i = a_i * b_i}. The two vectors must have the
687 @deftypefun int gsl_vector_div (gsl_vector * @var{a}, const gsl_vector * @var{b})
688 This function divides the elements of vector @var{a} by the elements of
689 vector @var{b}, @math{a'_i = a_i / b_i}. The two vectors must have the
693 @deftypefun int gsl_vector_scale (gsl_vector * @var{a}, const double @var{x})
694 This function multiplies the elements of vector @var{a} by the constant
695 factor @var{x}, @math{a'_i = x a_i}.
698 @deftypefun int gsl_vector_add_constant (gsl_vector * @var{a}, const double @var{x})
699 This function adds the constant value @var{x} to the elements of the
700 vector @var{a}, @math{a'_i = a_i + x}.
703 @node Finding maximum and minimum elements of vectors
704 @subsection Finding maximum and minimum elements of vectors
706 @deftypefun double gsl_vector_max (const gsl_vector * @var{v})
707 This function returns the maximum value in the vector @var{v}.
710 @deftypefun double gsl_vector_min (const gsl_vector * @var{v})
711 This function returns the minimum value in the vector @var{v}.
714 @deftypefun void gsl_vector_minmax (const gsl_vector * @var{v}, double * @var{min_out}, double * @var{max_out})
715 This function returns the minimum and maximum values in the vector
716 @var{v}, storing them in @var{min_out} and @var{max_out}.
719 @deftypefun size_t gsl_vector_max_index (const gsl_vector * @var{v})
720 This function returns the index of the maximum value in the vector @var{v}.
721 When there are several equal maximum elements then the lowest index is
725 @deftypefun size_t gsl_vector_min_index (const gsl_vector * @var{v})
726 This function returns the index of the minimum value in the vector @var{v}.
727 When there are several equal minimum elements then the lowest index is
731 @deftypefun void gsl_vector_minmax_index (const gsl_vector * @var{v}, size_t * @var{imin}, size_t * @var{imax})
732 This function returns the indices of the minimum and maximum values in
733 the vector @var{v}, storing them in @var{imin} and @var{imax}. When
734 there are several equal minimum or maximum elements then the lowest
735 indices are returned.
738 @node Vector properties
739 @subsection Vector properties
741 @deftypefun int gsl_vector_isnull (const gsl_vector * @var{v})
742 @deftypefunx int gsl_vector_ispos (const gsl_vector * @var{v})
743 @deftypefunx int gsl_vector_isneg (const gsl_vector * @var{v})
744 @deftypefunx int gsl_vector_isnonneg (const gsl_vector * @var{v})
745 These functions return 1 if all the elements of the vector @var{v} are
746 zero, strictly positive, strictly negative, or non-negative
747 respectively, and 0 otherwise.
750 @node Example programs for vectors
751 @subsection Example programs for vectors
753 This program shows how to allocate, initialize and read from a vector
754 using the functions @code{gsl_vector_alloc}, @code{gsl_vector_set} and
755 @code{gsl_vector_get}.
758 @verbatiminclude examples/vector.c
763 Here is the output from the program. The final loop attempts to read
764 outside the range of the vector @code{v}, and the error is trapped by
765 the range-checking code in @code{gsl_vector_get}.
772 gsl: vector_source.c:12: ERROR: index out of range
773 Default GSL error handler invoked.
774 Aborted (core dumped)
779 The next program shows how to write a vector to a file.
782 @verbatiminclude examples/vectorw.c
787 After running this program the file @file{test.dat} should contain the
788 elements of @code{v}, written using the format specifier
789 @code{%.5g}. The vector could then be read back in using the function
790 @code{gsl_vector_fscanf (f, v)} as follows:
793 @verbatiminclude examples/vectorr.c
800 @cindex physical dimension, matrices
801 @cindex trailing dimension, matrices
802 @cindex leading dimension, matrices
804 Matrices are defined by a @code{gsl_matrix} structure which describes a
805 generalized slice of a block. Like a vector it represents a set of
806 elements in an area of memory, but uses two indices instead of one.
808 The @code{gsl_matrix} structure contains six components, the two
809 dimensions of the matrix, a physical dimension, a pointer to the memory
810 where the elements of the matrix are stored, @var{data}, a pointer to
811 the block owned by the matrix @var{block}, if any, and an ownership
812 flag, @var{owner}. The physical dimension determines the memory layout
813 and can differ from the matrix dimension to allow the use of
814 submatrices. The @code{gsl_matrix} structure is very simple and looks
831 Matrices are stored in row-major order, meaning that each row of
832 elements forms a contiguous block in memory. This is the standard
833 ``C-language ordering'' of two-dimensional arrays. Note that @sc{fortran}
834 stores arrays in column-major order. The number of rows is @var{size1}.
835 The range of valid row indices runs from 0 to @code{size1-1}. Similarly
836 @var{size2} is the number of columns. The range of valid column indices
837 runs from 0 to @code{size2-1}. The physical row dimension @var{tda}, or
838 @dfn{trailing dimension}, specifies the size of a row of the matrix as
841 For example, in the following matrix @var{size1} is 3, @var{size2} is 4,
842 and @var{tda} is 8. The physical memory layout of the matrix begins in
843 the top left hand-corner and proceeds from left to right along each row
848 00 01 02 03 XX XX XX XX
849 10 11 12 13 XX XX XX XX
850 20 21 22 23 XX XX XX XX
855 Each unused memory location is represented by ``@code{XX}''. The
856 pointer @var{data} gives the location of the first element of the matrix
857 in memory. The pointer @var{block} stores the location of the memory
858 block in which the elements of the matrix are located (if any). If the
859 matrix owns this block then the @var{owner} field is set to one and the
860 block will be deallocated when the matrix is freed. If the matrix is
861 only a slice of a block owned by another object then the @var{owner} field is
862 zero and any underlying block will not be freed.
864 The functions for allocating and accessing matrices are defined in
868 * Matrix allocation::
869 * Accessing matrix elements::
870 * Initializing matrix elements::
871 * Reading and writing matrices::
873 * Creating row and column views::
875 * Copying rows and columns::
876 * Exchanging rows and columns::
877 * Matrix operations::
878 * Finding maximum and minimum elements of matrices::
879 * Matrix properties::
880 * Example programs for matrices::
883 @node Matrix allocation
884 @subsection Matrix allocation
886 The functions for allocating memory to a matrix follow the style of
887 @code{malloc} and @code{free}. They also perform their own error
888 checking. If there is insufficient memory available to allocate a matrix
889 then the functions call the GSL error handler (with an error number of
890 @code{GSL_ENOMEM}) in addition to returning a null pointer. Thus if you
891 use the library error handler to abort your program then it isn't
892 necessary to check every @code{alloc}.
894 @deftypefun {gsl_matrix *} gsl_matrix_alloc (size_t @var{n1}, size_t @var{n2})
895 This function creates a matrix of size @var{n1} rows by @var{n2}
896 columns, returning a pointer to a newly initialized matrix struct. A new
897 block is allocated for the elements of the matrix, and stored in the
898 @var{block} component of the matrix struct. The block is ``owned'' by the
899 matrix, and will be deallocated when the matrix is deallocated.
902 @deftypefun {gsl_matrix *} gsl_matrix_calloc (size_t @var{n1}, size_t @var{n2})
903 This function allocates memory for a matrix of size @var{n1} rows by
904 @var{n2} columns and initializes all the elements of the matrix to zero.
907 @deftypefun void gsl_matrix_free (gsl_matrix * @var{m})
908 This function frees a previously allocated matrix @var{m}. If the
909 matrix was created using @code{gsl_matrix_alloc} then the block
910 underlying the matrix will also be deallocated. If the matrix has
911 been created from another object then the memory is still owned by
912 that object and will not be deallocated. The matrix @var{m} must be a
913 valid matrix object (a null pointer is not allowed).
916 @node Accessing matrix elements
917 @subsection Accessing matrix elements
918 @cindex matrices, range-checking
919 @cindex range-checking for matrices
921 The functions for accessing the elements of a matrix use the same range
922 checking system as vectors. You can turn off range checking by recompiling
923 your program with the preprocessor definition
924 @code{GSL_RANGE_CHECK_OFF}.
926 The elements of the matrix are stored in ``C-order'', where the second
927 index moves continuously through memory. More precisely, the element
928 accessed by the function @code{gsl_matrix_get(m,i,j)} and
929 @code{gsl_matrix_set(m,i,j,x)} is
932 m->data[i * m->tda + j]
937 where @var{tda} is the physical row-length of the matrix.
939 @deftypefun double gsl_matrix_get (const gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
940 This function returns the @math{(i,j)}-th element of a matrix
941 @var{m}. If @var{i} or @var{j} lie outside the allowed range of 0 to
942 @math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error handler is invoked and 0
943 is returned. @inlinefn{}
946 @deftypefun void gsl_matrix_set (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j}, double @var{x})
947 This function sets the value of the @math{(i,j)}-th element of a
948 matrix @var{m} to @var{x}. If @var{i} or @var{j} lies outside the
949 allowed range of 0 to @math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error
950 handler is invoked. @inlinefn{}
953 @deftypefun {double *} gsl_matrix_ptr (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
954 @deftypefunx {const double *} gsl_matrix_const_ptr (const gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
955 These functions return a pointer to the @math{(i,j)}-th element of a
956 matrix @var{m}. If @var{i} or @var{j} lie outside the allowed range of
957 0 to @math{@var{n1}-1} and 0 to @math{@var{n2}-1} then the error handler is invoked
958 and a null pointer is returned. @inlinefns{}
961 @node Initializing matrix elements
962 @subsection Initializing matrix elements
963 @cindex matrices, initializing
964 @cindex initializing matrices
965 @cindex identity matrix
966 @cindex matrix, identity
969 @cindex constant matrix
970 @cindex matrix, constant
972 @deftypefun void gsl_matrix_set_all (gsl_matrix * @var{m}, double @var{x})
973 This function sets all the elements of the matrix @var{m} to the value
977 @deftypefun void gsl_matrix_set_zero (gsl_matrix * @var{m})
978 This function sets all the elements of the matrix @var{m} to zero.
981 @deftypefun void gsl_matrix_set_identity (gsl_matrix * @var{m})
982 This function sets the elements of the matrix @var{m} to the
983 corresponding elements of the identity matrix, @math{m(i,j) =
984 \delta(i,j)}, i.e. a unit diagonal with all off-diagonal elements zero.
985 This applies to both square and rectangular matrices.
988 @node Reading and writing matrices
989 @subsection Reading and writing matrices
991 The library provides functions for reading and writing matrices to a file
992 as binary data or formatted text.
994 @deftypefun int gsl_matrix_fwrite (FILE * @var{stream}, const gsl_matrix * @var{m})
995 This function writes the elements of the matrix @var{m} to the stream
996 @var{stream} in binary format. The return value is 0 for success and
997 @code{GSL_EFAILED} if there was a problem writing to the file. Since the
998 data is written in the native binary format it may not be portable
999 between different architectures.
1002 @deftypefun int gsl_matrix_fread (FILE * @var{stream}, gsl_matrix * @var{m})
1003 This function reads into the matrix @var{m} from the open stream
1004 @var{stream} in binary format. The matrix @var{m} must be preallocated
1005 with the correct dimensions since the function uses the size of @var{m} to
1006 determine how many bytes to read. The return value is 0 for success and
1007 @code{GSL_EFAILED} if there was a problem reading from the file. The
1008 data is assumed to have been written in the native binary format on the
1012 @deftypefun int gsl_matrix_fprintf (FILE * @var{stream}, const gsl_matrix * @var{m}, const char * @var{format})
1013 This function writes the elements of the matrix @var{m} line-by-line to
1014 the stream @var{stream} using the format specifier @var{format}, which
1015 should be one of the @code{%g}, @code{%e} or @code{%f} formats for
1016 floating point numbers and @code{%d} for integers. The function returns
1017 0 for success and @code{GSL_EFAILED} if there was a problem writing to
1021 @deftypefun int gsl_matrix_fscanf (FILE * @var{stream}, gsl_matrix * @var{m})
1022 This function reads formatted data from the stream @var{stream} into the
1023 matrix @var{m}. The matrix @var{m} must be preallocated with the correct
1024 dimensions since the function uses the size of @var{m} to determine how many
1025 numbers to read. The function returns 0 for success and
1026 @code{GSL_EFAILED} if there was a problem reading from the file.
1030 @subsection Matrix views
1032 A matrix view is a temporary object, stored on the stack, which can be
1033 used to operate on a subset of matrix elements. Matrix views can be
1034 defined for both constant and non-constant matrices using separate types
1035 that preserve constness. A matrix view has the type
1036 @code{gsl_matrix_view} and a constant matrix view has the type
1037 @code{gsl_matrix_const_view}. In both cases the elements of the view
1038 can by accessed using the @code{matrix} component of the view object. A
1039 pointer @code{gsl_matrix *} or @code{const gsl_matrix *} can be obtained
1040 by taking the address of the @code{matrix} component with the @code{&}
1041 operator. In addition to matrix views it is also possible to create
1042 vector views of a matrix, such as row or column views.
1044 @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})
1045 @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})
1046 These functions return a matrix view of a submatrix of the matrix
1047 @var{m}. The upper-left element of the submatrix is the element
1048 (@var{k1},@var{k2}) of the original matrix. The submatrix has @var{n1}
1049 rows and @var{n2} columns. The physical number of columns in memory
1050 given by @var{tda} is unchanged. Mathematically, the
1051 @math{(i,j)}-th element of the new matrix is given by,
1054 m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j]
1058 where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1059 runs from 0 to @code{n2-1}.
1061 The @code{data} pointer of the returned matrix struct is set to null if
1062 the combined parameters (@var{i},@var{j},@var{n1},@var{n2},@var{tda})
1063 overrun the ends of the original matrix.
1065 The new matrix view is only a view of the block underlying the existing
1066 matrix, @var{m}. The block containing the elements of @var{m} is not
1067 owned by the new matrix view. When the view goes out of scope the
1068 original matrix @var{m} and its block will continue to exist. The
1069 original memory can only be deallocated by freeing the original matrix.
1070 Of course, the original matrix should not be deallocated while the view
1073 The function @code{gsl_matrix_const_submatrix} is equivalent to
1074 @code{gsl_matrix_submatrix} but can be used for matrices which are
1075 declared @code{const}.
1079 @deftypefun gsl_matrix_view gsl_matrix_view_array (double * @var{base}, size_t @var{n1}, size_t @var{n2})
1080 @deftypefunx gsl_matrix_const_view gsl_matrix_const_view_array (const double * @var{base}, size_t @var{n1}, size_t @var{n2})
1081 These functions return a matrix view of the array @var{base}. The
1082 matrix has @var{n1} rows and @var{n2} columns. The physical number of
1083 columns in memory is also given by @var{n2}. Mathematically, the
1084 @math{(i,j)}-th element of the new matrix is given by,
1087 m'(i,j) = base[i*n2 + j]
1091 where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1092 runs from 0 to @code{n2-1}.
1094 The new matrix is only a view of the array @var{base}. When the view
1095 goes out of scope the original array @var{base} will continue to exist.
1096 The original memory can only be deallocated by freeing the original
1097 array. Of course, the original array should not be deallocated while
1098 the view is still in use.
1100 The function @code{gsl_matrix_const_view_array} is equivalent to
1101 @code{gsl_matrix_view_array} but can be used for matrices which are
1102 declared @code{const}.
1106 @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})
1107 @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})
1108 These functions return a matrix view of the array @var{base} with a
1109 physical number of columns @var{tda} which may differ from the corresponding
1110 dimension of the matrix. The matrix has @var{n1} rows and @var{n2}
1111 columns, and the physical number of columns in memory is given by
1112 @var{tda}. Mathematically, the @math{(i,j)}-th element of the new
1116 m'(i,j) = base[i*tda + j]
1120 where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1121 runs from 0 to @code{n2-1}.
1123 The new matrix is only a view of the array @var{base}. When the view
1124 goes out of scope the original array @var{base} will continue to exist.
1125 The original memory can only be deallocated by freeing the original
1126 array. Of course, the original array should not be deallocated while
1127 the view is still in use.
1129 The function @code{gsl_matrix_const_view_array_with_tda} is equivalent
1130 to @code{gsl_matrix_view_array_with_tda} but can be used for matrices
1131 which are declared @code{const}.
1134 @deftypefun gsl_matrix_view gsl_matrix_view_vector (gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2})
1135 @deftypefunx gsl_matrix_const_view gsl_matrix_const_view_vector (const gsl_vector * @var{v}, size_t @var{n1}, size_t @var{n2})
1136 These functions return a matrix view of the vector @var{v}. The matrix
1137 has @var{n1} rows and @var{n2} columns. The vector must have unit
1138 stride. The physical number of columns in memory is also given by
1139 @var{n2}. Mathematically, the @math{(i,j)}-th element of the new
1143 m'(i,j) = v->data[i*n2 + j]
1147 where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1148 runs from 0 to @code{n2-1}.
1150 The new matrix is only a view of the vector @var{v}. When the view
1151 goes out of scope the original vector @var{v} will continue to exist.
1152 The original memory can only be deallocated by freeing the original
1153 vector. Of course, the original vector should not be deallocated while
1154 the view is still in use.
1156 The function @code{gsl_matrix_const_view_vector} is equivalent to
1157 @code{gsl_matrix_view_vector} but can be used for matrices which are
1158 declared @code{const}.
1162 @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})
1163 @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})
1164 These functions return a matrix view of the vector @var{v} with a
1165 physical number of columns @var{tda} which may differ from the
1166 corresponding matrix dimension. The vector must have unit stride. The
1167 matrix has @var{n1} rows and @var{n2} columns, and the physical number
1168 of columns in memory is given by @var{tda}. Mathematically, the
1169 @math{(i,j)}-th element of the new matrix is given by,
1172 m'(i,j) = v->data[i*tda + j]
1176 where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1177 runs from 0 to @code{n2-1}.
1179 The new matrix is only a view of the vector @var{v}. When the view
1180 goes out of scope the original vector @var{v} will continue to exist.
1181 The original memory can only be deallocated by freeing the original
1182 vector. Of course, the original vector should not be deallocated while
1183 the view is still in use.
1185 The function @code{gsl_matrix_const_view_vector_with_tda} is equivalent
1186 to @code{gsl_matrix_view_vector_with_tda} but can be used for matrices
1187 which are declared @code{const}.
1191 @comment @node Modifying matrix views
1192 @comment @subsection Modifying matrix views
1194 @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})
1195 @comment This function modifies and existing matrix view @var{m} to form a new
1196 @comment view of a matrix @var{mm}, starting from element (@var{k1},@var{k2}).
1197 @comment The matrix view has @var{n1} rows and @var{n2} columns. Any existing
1198 @comment view in @var{m} will be lost as a result of this function.
1199 @comment @end deftypefun
1201 @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})
1202 @comment This function modifies and existing matrix view @var{m} to form a new
1203 @comment view of a vector @var{v}, starting from element @var{offset}. The
1204 @comment vector has @var{n1} rows and @var{n2} columns. Any
1205 @comment existing view in @var{m} will be lost as a result of this function.
1206 @comment @end deftypefun
1208 @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})
1209 @comment This function modifies and existing matrix view @var{m} to form a new
1210 @comment view of an array @var{base}, starting from element @var{offset}. The
1211 @comment matrix has @var{n1} rows and @var{n2} columns. Any
1212 @comment existing view in @var{m} will be lost as a result of this function.
1213 @comment @end deftypefun
1215 @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})
1216 @comment This function creates a matrix as a slice of the block @var{b},
1217 @comment returning a pointer to a newly initialized matrix struct. The start of
1218 @comment the matrix is offset by @var{offset} elements from the start of the
1219 @comment block. The matrix has @var{n1} rows and @var{n2} columns, with the
1220 @comment physical number of columns in memory given by @var{tda}.
1221 @comment Mathematically, the (@var{i},@var{j})-th element of the matrix is given by,
1224 @comment m(i,j) = b->data[offset + i*tda + j]
1225 @comment @end example
1227 @comment where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1228 @comment runs from 0 to @code{n2-1}.
1230 @comment A null pointer is returned if the combined parameters
1231 @comment (@var{offset},@var{n1},@var{n2},@var{tda}) overrun the end of the block
1232 @comment or if insufficient memory is available to store the matrix.
1234 @comment The matrix is only a view of the block @var{b}, and the block is not
1235 @comment owned by the matrix. When the matrix is deallocated the block @var{b}
1236 @comment will continue to exist. This memory can only be deallocated by freeing
1237 @comment the block itself. Of course, this block should not be deallocated while
1238 @comment the matrix is still in use.
1239 @comment @end deftypefun
1241 @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})
1243 @comment This function creates a matrix as a submatrix of the matrix @var{m},
1244 @comment returning a pointer to a newly initialized matrix struct. The upper-left
1245 @comment element of the submatrix is the element (@var{k1},@var{k2}) of the
1246 @comment original matrix. The submatrix has @var{n1} rows and @var{n2} columns.
1247 @comment The physical number of columns in memory given by @var{tda} is
1248 @comment unchanged. Mathematically, the (@var{i},@var{j})-th element of the
1249 @comment new matrix is given by,
1252 @comment m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j]
1253 @comment @end example
1255 @comment where the index @var{i} runs from 0 to @code{n1-1} and the index @var{j}
1256 @comment runs from 0 to @code{n2-1}.
1258 @comment A null pointer is returned if the combined parameters
1259 @comment (@var{k1},@var{k2},@var{n1},@var{n2},@var{tda}) overrun the end of the
1260 @comment original matrix or if insufficient memory is available to store the matrix.
1262 @comment The new matrix is only a view of the block underlying the existing
1263 @comment matrix, @var{m}. The block is not owned by the new matrix. When the new
1264 @comment matrix is deallocated the original matrix @var{m} and its block will
1265 @comment continue to exist. The original memory can only be deallocated by
1266 @comment freeing the original matrix. Of course, the original matrix should not
1267 @comment be deallocated while the new matrix is still in use.
1268 @comment @end deftypefun
1270 @node Creating row and column views
1271 @subsection Creating row and column views
1273 In general there are two ways to access an object, by reference or by
1274 copying. The functions described in this section create vector views
1275 which allow access to a row or column of a matrix by reference.
1276 Modifying elements of the view is equivalent to modifying the matrix,
1277 since both the vector view and the matrix point to the same memory
1280 @deftypefun gsl_vector_view gsl_matrix_row (gsl_matrix * @var{m}, size_t @var{i})
1281 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_row (const gsl_matrix * @var{m}, size_t @var{i})
1282 These functions return a vector view of the @var{i}-th row of the matrix
1283 @var{m}. The @code{data} pointer of the new vector is set to null if
1284 @var{i} is out of range.
1286 The function @code{gsl_vector_const_row} is equivalent to
1287 @code{gsl_matrix_row} but can be used for matrices which are declared
1291 @deftypefun gsl_vector_view gsl_matrix_column (gsl_matrix * @var{m}, size_t @var{j})
1292 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_column (const gsl_matrix * @var{m}, size_t @var{j})
1293 These functions return a vector view of the @var{j}-th column of the
1294 matrix @var{m}. The @code{data} pointer of the new vector is set to
1295 null if @var{j} is out of range.
1297 The function @code{gsl_vector_const_column} is equivalent to
1298 @code{gsl_matrix_column} but can be used for matrices which are declared
1302 @deftypefun gsl_vector_view gsl_matrix_subrow (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{offset}, size_t @var{n})
1303 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_subrow (const gsl_matrix * @var{m}, size_t @var{i}, size_t @var{offset}, size_t @var{n})
1304 These functions return a vector view of the @var{i}-th row of the matrix
1305 @var{m} beginning at @var{offset} elements past the first column and
1306 containing @var{n} elements. The @code{data} pointer of the new vector
1307 is set to null if @var{i}, @var{offset}, or @var{n} are out of range.
1309 The function @code{gsl_vector_const_subrow} is equivalent to
1310 @code{gsl_matrix_subrow} but can be used for matrices which are declared
1314 @deftypefun gsl_vector_view gsl_matrix_subcolumn (gsl_matrix * @var{m}, size_t @var{j}, size_t @var{offset}, size_t @var{n})
1315 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_subcolumn (const gsl_matrix * @var{m}, size_t @var{j}, size_t @var{offset}, size_t @var{n})
1316 These functions return a vector view of the @var{j}-th column of the matrix
1317 @var{m} beginning at @var{offset} elements past the first row and
1318 containing @var{n} elements. The @code{data} pointer of the new vector
1319 is set to null if @var{j}, @var{offset}, or @var{n} are out of range.
1321 The function @code{gsl_vector_const_subcolumn} is equivalent to
1322 @code{gsl_matrix_subcolumn} but can be used for matrices which are declared
1326 @cindex matrix diagonal
1327 @cindex diagonal, of a matrix
1328 @deftypefun gsl_vector_view gsl_matrix_diagonal (gsl_matrix * @var{m})
1329 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_diagonal (const gsl_matrix * @var{m})
1330 These functions returns a vector view of the diagonal of the matrix
1331 @var{m}. The matrix @var{m} is not required to be square. For a
1332 rectangular matrix the length of the diagonal is the same as the smaller
1333 dimension of the matrix.
1335 The function @code{gsl_matrix_const_diagonal} is equivalent to
1336 @code{gsl_matrix_diagonal} but can be used for matrices which are
1337 declared @code{const}.
1340 @cindex matrix subdiagonal
1341 @cindex subdiagonal, of a matrix
1342 @deftypefun gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * @var{m}, size_t @var{k})
1343 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_subdiagonal (const gsl_matrix * @var{m}, size_t @var{k})
1344 These functions return a vector view of the @var{k}-th subdiagonal of
1345 the matrix @var{m}. The matrix @var{m} is not required to be square.
1346 The diagonal of the matrix corresponds to @math{k = 0}.
1348 The function @code{gsl_matrix_const_subdiagonal} is equivalent to
1349 @code{gsl_matrix_subdiagonal} but can be used for matrices which are
1350 declared @code{const}.
1353 @cindex matrix superdiagonal
1354 @cindex superdiagonal, matrix
1355 @deftypefun gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * @var{m}, size_t @var{k})
1356 @deftypefunx {gsl_vector_const_view} gsl_matrix_const_superdiagonal (const gsl_matrix * @var{m}, size_t @var{k})
1357 These functions return a vector view of the @var{k}-th superdiagonal of
1358 the matrix @var{m}. The matrix @var{m} is not required to be square. The
1359 diagonal of the matrix corresponds to @math{k = 0}.
1361 The function @code{gsl_matrix_const_superdiagonal} is equivalent to
1362 @code{gsl_matrix_superdiagonal} but can be used for matrices which are
1363 declared @code{const}.
1366 @comment @deftypefun {gsl_vector *} gsl_vector_alloc_row_from_matrix (gsl_matrix * @var{m}, size_t @var{i})
1367 @comment This function allocates a new @code{gsl_vector} struct which points to
1368 @comment the @var{i}-th row of the matrix @var{m}.
1369 @comment @end deftypefun
1371 @comment @deftypefun {gsl_vector *} gsl_vector_alloc_col_from_matrix (gsl_matrix * @var{m}, size_t @var{j})
1372 @comment This function allocates a new @code{gsl_vector} struct which points to
1373 @comment the @var{j}-th column of the matrix @var{m}.
1374 @comment @end deftypefun
1376 @node Copying matrices
1377 @subsection Copying matrices
1379 @deftypefun int gsl_matrix_memcpy (gsl_matrix * @var{dest}, const gsl_matrix * @var{src})
1380 This function copies the elements of the matrix @var{src} into the
1381 matrix @var{dest}. The two matrices must have the same size.
1384 @deftypefun int gsl_matrix_swap (gsl_matrix * @var{m1}, gsl_matrix * @var{m2})
1385 This function exchanges the elements of the matrices @var{m1} and
1386 @var{m2} by copying. The two matrices must have the same size.
1389 @node Copying rows and columns
1390 @subsection Copying rows and columns
1392 The functions described in this section copy a row or column of a matrix
1393 into a vector. This allows the elements of the vector and the matrix to
1394 be modified independently. Note that if the matrix and the vector point
1395 to overlapping regions of memory then the result will be undefined. The
1396 same effect can be achieved with more generality using
1397 @code{gsl_vector_memcpy} with vector views of rows and columns.
1399 @deftypefun int gsl_matrix_get_row (gsl_vector * @var{v}, const gsl_matrix * @var{m}, size_t @var{i})
1400 This function copies the elements of the @var{i}-th row of the matrix
1401 @var{m} into the vector @var{v}. The length of the vector must be the
1402 same as the length of the row.
1405 @deftypefun int gsl_matrix_get_col (gsl_vector * @var{v}, const gsl_matrix * @var{m}, size_t @var{j})
1406 This function copies the elements of the @var{j}-th column of the matrix
1407 @var{m} into the vector @var{v}. The length of the vector must be the
1408 same as the length of the column.
1411 @deftypefun int gsl_matrix_set_row (gsl_matrix * @var{m}, size_t @var{i}, const gsl_vector * @var{v})
1412 This function copies the elements of the vector @var{v} into the
1413 @var{i}-th row of the matrix @var{m}. The length of the vector must be
1414 the same as the length of the row.
1417 @deftypefun int gsl_matrix_set_col (gsl_matrix * @var{m}, size_t @var{j}, const gsl_vector * @var{v})
1418 This function copies the elements of the vector @var{v} into the
1419 @var{j}-th column of the matrix @var{m}. The length of the vector must be
1420 the same as the length of the column.
1423 @node Exchanging rows and columns
1424 @subsection Exchanging rows and columns
1426 The following functions can be used to exchange the rows and columns of
1429 @deftypefun int gsl_matrix_swap_rows (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
1430 This function exchanges the @var{i}-th and @var{j}-th rows of the matrix
1434 @deftypefun int gsl_matrix_swap_columns (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
1435 This function exchanges the @var{i}-th and @var{j}-th columns of the
1436 matrix @var{m} in-place.
1439 @deftypefun int gsl_matrix_swap_rowcol (gsl_matrix * @var{m}, size_t @var{i}, size_t @var{j})
1440 This function exchanges the @var{i}-th row and @var{j}-th column of the
1441 matrix @var{m} in-place. The matrix must be square for this operation to
1445 @deftypefun int gsl_matrix_transpose_memcpy (gsl_matrix * @var{dest}, const gsl_matrix * @var{src})
1446 This function makes the matrix @var{dest} the transpose of the matrix
1447 @var{src} by copying the elements of @var{src} into @var{dest}. This
1448 function works for all matrices provided that the dimensions of the matrix
1449 @var{dest} match the transposed dimensions of the matrix @var{src}.
1452 @deftypefun int gsl_matrix_transpose (gsl_matrix * @var{m})
1453 This function replaces the matrix @var{m} by its transpose by copying
1454 the elements of the matrix in-place. The matrix must be square for this
1455 operation to be possible.
1458 @node Matrix operations
1459 @subsection Matrix operations
1461 The following operations are defined for real and complex matrices.
1463 @deftypefun int gsl_matrix_add (gsl_matrix * @var{a}, const gsl_matrix * @var{b})
1464 This function adds the elements of matrix @var{b} to the elements of
1465 matrix @var{a}, @math{a'(i,j) = a(i,j) + b(i,j)}. The two matrices must have the
1469 @deftypefun int gsl_matrix_sub (gsl_matrix * @var{a}, const gsl_matrix * @var{b})
1470 This function subtracts the elements of matrix @var{b} from the elements of
1471 matrix @var{a}, @math{a'(i,j) = a(i,j) - b(i,j)}. The two matrices must have the
1475 @deftypefun int gsl_matrix_mul_elements (gsl_matrix * @var{a}, const gsl_matrix * @var{b})
1476 This function multiplies the elements of matrix @var{a} by the elements of
1477 matrix @var{b}, @math{a'(i,j) = a(i,j) * b(i,j)}. The two matrices must have the
1481 @deftypefun int gsl_matrix_div_elements (gsl_matrix * @var{a}, const gsl_matrix * @var{b})
1482 This function divides the elements of matrix @var{a} by the elements of
1483 matrix @var{b}, @math{a'(i,j) = a(i,j) / b(i,j)}. The two matrices must have the
1487 @deftypefun int gsl_matrix_scale (gsl_matrix * @var{a}, const double @var{x})
1488 This function multiplies the elements of matrix @var{a} by the constant
1489 factor @var{x}, @math{a'(i,j) = x a(i,j)}.
1492 @deftypefun int gsl_matrix_add_constant (gsl_matrix * @var{a}, const double @var{x})
1493 This function adds the constant value @var{x} to the elements of the
1494 matrix @var{a}, @math{a'(i,j) = a(i,j) + x}.
1497 @node Finding maximum and minimum elements of matrices
1498 @subsection Finding maximum and minimum elements of matrices
1500 The following operations are only defined for real matrices.
1502 @deftypefun double gsl_matrix_max (const gsl_matrix * @var{m})
1503 This function returns the maximum value in the matrix @var{m}.
1506 @deftypefun double gsl_matrix_min (const gsl_matrix * @var{m})
1507 This function returns the minimum value in the matrix @var{m}.
1510 @deftypefun void gsl_matrix_minmax (const gsl_matrix * @var{m}, double * @var{min_out}, double * @var{max_out})
1511 This function returns the minimum and maximum values in the matrix
1512 @var{m}, storing them in @var{min_out} and @var{max_out}.
1515 @deftypefun void gsl_matrix_max_index (const gsl_matrix * @var{m}, size_t * @var{imax}, size_t * @var{jmax})
1516 This function returns the indices of the maximum value in the matrix
1517 @var{m}, storing them in @var{imax} and @var{jmax}. When there are
1518 several equal maximum elements then the first element found is returned,
1519 searching in row-major order.
1522 @deftypefun void gsl_matrix_min_index (const gsl_matrix * @var{m}, size_t * @var{imin}, size_t * @var{jmin})
1523 This function returns the indices of the minimum value in the matrix
1524 @var{m}, storing them in @var{imin} and @var{jmin}. When there are
1525 several equal minimum elements then the first element found is returned,
1526 searching in row-major order.
1529 @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})
1530 This function returns the indices of the minimum and maximum values in
1531 the matrix @var{m}, storing them in (@var{imin},@var{jmin}) and
1532 (@var{imax},@var{jmax}). When there are several equal minimum or maximum
1533 elements then the first elements found are returned, searching in
1537 @node Matrix properties
1538 @subsection Matrix properties
1540 @deftypefun int gsl_matrix_isnull (const gsl_matrix * @var{m})
1541 @deftypefunx int gsl_matrix_ispos (const gsl_matrix * @var{m})
1542 @deftypefunx int gsl_matrix_isneg (const gsl_matrix * @var{m})
1543 @deftypefunx int gsl_matrix_isnonneg (const gsl_matrix * @var{m})
1544 These functions return 1 if all the elements of the matrix @var{m} are
1545 zero, strictly positive, strictly negative, or non-negative
1546 respectively, and 0 otherwise. To test whether a matrix is
1547 positive-definite, use the Cholesky decomposition (@pxref{Cholesky
1551 @node Example programs for matrices
1552 @subsection Example programs for matrices
1554 The program below shows how to allocate, initialize and read from a matrix
1555 using the functions @code{gsl_matrix_alloc}, @code{gsl_matrix_set} and
1556 @code{gsl_matrix_get}.
1559 @verbatiminclude examples/matrix.c
1564 Here is the output from the program. The final loop attempts to read
1565 outside the range of the matrix @code{m}, and the error is trapped by
1566 the range-checking code in @code{gsl_matrix_get}.
1578 gsl: matrix_source.c:13: ERROR: first index out of range
1579 Default GSL error handler invoked.
1580 Aborted (core dumped)
1585 The next program shows how to write a matrix to a file.
1588 @verbatiminclude examples/matrixw.c
1593 After running this program the file @file{test.dat} should contain the
1594 elements of @code{m}, written in binary format. The matrix which is read
1595 back in using the function @code{gsl_matrix_fread} should be exactly
1596 equal to the original matrix.
1598 The following program demonstrates the use of vector views. The program
1599 computes the column norms of a matrix.
1602 @verbatiminclude examples/vectorview.c
1606 Here is the output of the program,
1610 @verbatiminclude examples/vectorview.out
1614 The results can be confirmed using @sc{gnu octave},
1618 GNU Octave, version 2.0.16.92
1619 octave> m = sin(0:9)' * ones(1,10)
1620 + ones(10,1) * cos(0:9);
1621 octave> sqrt(sum(m.^2))
1623 4.3146 3.1205 2.1932 3.2611 2.5342 2.5728
1624 4.2047 3.6520 2.0852 3.0731
1628 @node Vector and Matrix References and Further Reading
1629 @section References and Further Reading
1631 The block, vector and matrix objects in GSL follow the @code{valarray}
1632 model of C++. A description of this model can be found in the following
1638 @cite{The C++ Programming Language} (3rd Ed),
1639 Section 22.4 Vector Arithmetic.
1640 Addison-Wesley 1997, ISBN 0-201-88954-4.