@cindex linear algebra, BLAS
@cindex matrix, operations
@cindex vector, operations
@cindex BLAS
@cindex CBLAS
@cindex Basic Linear Algebra Subroutines (BLAS)

The Basic Linear Algebra Subprograms (@sc{blas}) define a set of fundamental
operations on vectors and matrices which can be used to create optimized
higher-level linear algebra functionality.

The library provides a low-level layer which corresponds directly to the
C-language @sc{blas} standard, referred to here as ``@sc{cblas}'', and a
higher-level interface for operations on GSL vectors and matrices.
Users who are interested in simple operations on GSL vector and matrix
objects should use the high-level layer described
in this chapter.  The functions are declared in the file
@file{gsl_blas.h} and should satisfy the needs of most users.  

Note that GSL matrices are implemented using dense-storage so the
interface only includes the corresponding dense-storage @sc{blas}
functions.  The full @sc{blas} functionality for band-format and
packed-format matrices is available through the low-level @sc{cblas}
interface.  Similarly, GSL vectors are restricted to positive strides,
whereas the low-level @sc{cblas} interface supports negative
strides as specified in the @sc{blas} standard.@footnote{In the low-level
@sc{cblas} interface, a negative stride accesses the vector elements
in reverse order, i.e. the @math{i}-th element is given by
@math{(N-i)*|incx|} for @math{incx < 0}.}
 
The interface for the @code{gsl_cblas} layer is specified in the file
@file{gsl_cblas.h}.  This interface corresponds to the @sc{blas} Technical
Forum's standard for the C interface to legacy @sc{blas}
implementations. Users who have access to other conforming @sc{cblas}
implementations can use these in place of the version provided by the
library.  Note that users who have only a Fortran @sc{blas} library can
use a @sc{cblas} conformant wrapper to convert it into a @sc{cblas}
library.  A reference @sc{cblas} wrapper for legacy Fortran
implementations exists as part of the @sc{cblas} standard and can
be obtained from Netlib.  The complete set of @sc{cblas} functions is
listed in an appendix (@pxref{GSL CBLAS Library}).

There are three levels of @sc{blas} operations,

@table @b
@item Level 1
Vector operations, e.g. @math{y = \alpha x + y}
@item Level 2
Matrix-vector operations, e.g. @math{y = \alpha A x + \beta y}
@item Level 3
Matrix-matrix operations, e.g. @math{C = \alpha A B + C}
@end table

@noindent
Each routine has a name which specifies the operation, the type of
matrices involved and their precisions.  Some of the most common
operations and their names are given below,

@table @b
@item DOT
scalar product, @math{x^T y}
@item AXPY
vector sum, @math{\alpha x + y}
@item MV
matrix-vector product, @math{A x}
@item SV
matrix-vector solve, @math{inv(A) x}
@item MM
matrix-matrix product, @math{A B}
@item SM
matrix-matrix solve, @math{inv(A) B}
@end table

@noindent
The types of matrices are,

@table @b
@item GE
general
@item GB
general band
@item SY
symmetric
@item SB
symmetric band
@item SP
symmetric packed
@item HE
hermitian
@item HB
hermitian band
@item HP
hermitian packed
@item TR
triangular 
@item TB
triangular band
@item TP
triangular packed
@end table

@noindent
Each operation is defined for four precisions,

@table @b
@item S
single real
@item D
double real
@item C
single complex
@item Z
double complex
@end table

@noindent
Thus, for example, the name @sc{sgemm} stands for ``single-precision
general matrix-matrix multiply'' and @sc{zgemm} stands for
``double-precision complex matrix-matrix multiply''.

Note that the vector and matrix arguments to BLAS functions must not
be aliased, as the results are undefined when the underlying arrays
overlap (@pxref{Aliasing of arrays}). 

@menu
* GSL BLAS Interface::          
* BLAS Examples::               
* BLAS References and Further Reading::  
@end menu

@node GSL BLAS Interface
@section GSL BLAS Interface

GSL provides dense vector and matrix objects, based on the relevant
built-in types.  The library provides an interface to the @sc{blas}
operations which apply to these objects.  The interface to this
functionality is given in the file @file{gsl_blas.h}.

@comment CblasNoTrans, CblasTrans, CblasConjTrans
@comment CblasUpper, CblasLower
@comment CblasNonUnit, CblasUnit
@comment CblasLeft, CblasRight

@menu
* Level 1 GSL BLAS Interface::  
* Level 2 GSL BLAS Interface::  
* Level 3 GSL BLAS Interface::  
@end menu

@node Level 1 GSL BLAS Interface
@subsection Level 1 

@deftypefun int gsl_blas_sdsdot (float @var{alpha}, const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, float * @var{result})
@cindex DOT, Level-1 BLAS
This function computes the sum @math{\alpha + x^T y} for the vectors
@var{x} and @var{y}, returning the result in @var{result}.
@end deftypefun

@deftypefun int gsl_blas_sdot (const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, float * @var{result})
@deftypefunx int gsl_blas_dsdot (const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, double * @var{result})
@deftypefunx int gsl_blas_ddot (const gsl_vector * @var{x}, const gsl_vector * @var{y}, double * @var{result})
These functions compute the scalar product @math{x^T y} for the vectors
@var{x} and @var{y}, returning the result in @var{result}.
@end deftypefun

@deftypefun int gsl_blas_cdotu (const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_complex_float * @var{dotu})
@deftypefunx int gsl_blas_zdotu (const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_complex * @var{dotu})
These functions compute the complex scalar product @math{x^T y} for the
vectors @var{x} and @var{y}, returning the result in @var{dotu}
@end deftypefun

@deftypefun int gsl_blas_cdotc (const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_complex_float * @var{dotc})
@deftypefunx int gsl_blas_zdotc (const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_complex * @var{dotc})
These functions compute the complex conjugate scalar product @math{x^H
y} for the vectors @var{x} and @var{y}, returning the result in
@var{dotc}
@end deftypefun

@deftypefun float gsl_blas_snrm2 (const gsl_vector_float * @var{x})
@deftypefunx double gsl_blas_dnrm2 (const gsl_vector * @var{x})
@cindex NRM2, Level-1 BLAS
These functions compute the Euclidean norm 
@c{$||x||_2 = \sqrt{\sum x_i^2}$}
@math{||x||_2 = \sqrt @{\sum x_i^2@}} of the vector @var{x}.
@end deftypefun

@deftypefun float gsl_blas_scnrm2 (const gsl_vector_complex_float * @var{x})
@deftypefunx double gsl_blas_dznrm2 (const gsl_vector_complex * @var{x})
These functions compute the Euclidean norm of the complex vector @var{x},
@tex
\beforedisplay
$$
||x||_2 = \sqrt{\sum (\Re(x_i)^2 + \Im(x_i)^2)}.
$$
\afterdisplay
@end tex
@ifinfo

@example
||x||_2 = \sqrt @{\sum (\Re(x_i)^2 + \Im(x_i)^2)@}.
@end example
@end ifinfo
@end deftypefun

@deftypefun float gsl_blas_sasum (const gsl_vector_float * @var{x})
@deftypefunx double gsl_blas_dasum (const gsl_vector * @var{x})
@cindex ASUM, Level-1 BLAS
These functions compute the absolute sum @math{\sum |x_i|} of the
elements of the vector @var{x}.
@end deftypefun

@deftypefun float gsl_blas_scasum (const gsl_vector_complex_float * @var{x})
@deftypefunx double gsl_blas_dzasum (const gsl_vector_complex * @var{x})
These functions compute the sum of the magnitudes of the real and
imaginary parts of the complex vector @var{x}, 
@c{$\sum \left( |\Re(x_i)| + |\Im(x_i)| \right)$}
@math{\sum |\Re(x_i)| + |\Im(x_i)|}.
@end deftypefun

@deftypefun CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * @var{x})
@deftypefunx CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * @var{x})
@deftypefunx CBLAS_INDEX_t gsl_blas_icamax (const gsl_vector_complex_float * @var{x})
@deftypefunx CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * @var{x})
@cindex AMAX, Level-1 BLAS
These functions return the index of the largest element of the vector
@var{x}. The largest element is determined by its absolute magnitude for
real vectors and by the sum of the magnitudes of the real and imaginary
parts @math{|\Re(x_i)| + |\Im(x_i)|} for complex vectors.  If the
largest value occurs several times then the index of the first
occurrence is returned.
@end deftypefun

@deftypefun int gsl_blas_sswap (gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
@deftypefunx int gsl_blas_dswap (gsl_vector * @var{x}, gsl_vector * @var{y})
@deftypefunx int gsl_blas_cswap (gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
@deftypefunx int gsl_blas_zswap (gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
@cindex SWAP, Level-1 BLAS
These functions exchange the elements of the vectors @var{x} and @var{y}.
@end deftypefun

@deftypefun int gsl_blas_scopy (const gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
@deftypefunx int gsl_blas_dcopy (const gsl_vector * @var{x}, gsl_vector * @var{y})
@deftypefunx int gsl_blas_ccopy (const gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
@deftypefunx int gsl_blas_zcopy (const gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
@cindex COPY, Level-1 BLAS
These functions copy the elements of the vector @var{x} into the vector
@var{y}.
@end deftypefun


@deftypefun int gsl_blas_saxpy (float @var{alpha}, const gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
@deftypefunx int gsl_blas_daxpy (double @var{alpha}, const gsl_vector * @var{x}, gsl_vector * @var{y})
@deftypefunx int gsl_blas_caxpy (const gsl_complex_float @var{alpha}, const gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
@deftypefunx int gsl_blas_zaxpy (const gsl_complex @var{alpha}, const gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
@cindex AXPY, Level-1 BLAS
@cindex DAXPY, Level-1 BLAS
@cindex SAXPY, Level-1 BLAS
These functions compute the sum @math{y = \alpha x + y} for the vectors
@var{x} and @var{y}.
@end deftypefun

@deftypefun void gsl_blas_sscal (float @var{alpha}, gsl_vector_float * @var{x})
@deftypefunx void gsl_blas_dscal (double @var{alpha}, gsl_vector * @var{x})
@deftypefunx void gsl_blas_cscal (const gsl_complex_float @var{alpha}, gsl_vector_complex_float * @var{x})
@deftypefunx void gsl_blas_zscal (const gsl_complex @var{alpha}, gsl_vector_complex * @var{x})
@deftypefunx void gsl_blas_csscal (float @var{alpha}, gsl_vector_complex_float * @var{x})
@deftypefunx void gsl_blas_zdscal (double @var{alpha}, gsl_vector_complex * @var{x})
@cindex SCAL, Level-1 BLAS
These functions rescale the vector @var{x} by the multiplicative factor
@var{alpha}.
@end deftypefun

@deftypefun int gsl_blas_srotg (float @var{a}[], float @var{b}[], float @var{c}[], float @var{s}[])
@deftypefunx int gsl_blas_drotg (double @var{a}[], double @var{b}[], double @var{c}[], double @var{s}[])
@cindex ROTG, Level-1 BLAS
@cindex Givens Rotation, BLAS
These functions compute a Givens rotation @math{(c,s)} which zeroes the
vector @math{(a,b)},
@tex
\beforedisplay
$$
\left(
\matrix{c&s\cr
-s&c\cr}
\right)
\left(
\matrix{a\cr
b\cr}
\right)
=
\left(
\matrix{r'\cr
0\cr}
\right)
$$
\afterdisplay
@end tex
@ifinfo

@example
[  c  s ] [ a ] = [ r ]
[ -s  c ] [ b ]   [ 0 ]
@end example

@end ifinfo
@noindent
The variables @var{a} and @var{b} are overwritten by the routine.
@end deftypefun

@deftypefun int gsl_blas_srot (gsl_vector_float * @var{x}, gsl_vector_float * @var{y}, float @var{c}, float @var{s})
@deftypefunx int gsl_blas_drot (gsl_vector * @var{x}, gsl_vector * @var{y}, const double @var{c}, const double @var{s})
These functions apply a Givens rotation @math{(x', y') = (c x + s y, -s
x + c y)} to the vectors @var{x}, @var{y}.
@end deftypefun

@deftypefun int gsl_blas_srotmg (float @var{d1}[], float @var{d2}[], float @var{b1}[], float @var{b2}, float @var{P}[])
@deftypefunx int gsl_blas_drotmg (double @var{d1}[], double @var{d2}[], double @var{b1}[], double @var{b2}, double @var{P}[])
@cindex Modified Givens Rotation, BLAS
@cindex Givens Rotation, Modified, BLAS
These functions compute a modified Givens transformation.  The modified
Givens transformation is defined in the original Level-1 @sc{blas}
specification, given in the references.
@end deftypefun

@deftypefun int gsl_blas_srotm (gsl_vector_float * @var{x}, gsl_vector_float * @var{y}, const float @var{P}[])
@deftypefunx int gsl_blas_drotm (gsl_vector * @var{x}, gsl_vector * @var{y}, const double @var{P}[])
These functions apply a modified Givens transformation.  
@end deftypefun

@node Level 2 GSL BLAS Interface
@subsection Level 2 

@deftypefun int gsl_blas_sgemv (CBLAS_TRANSPOSE_t @var{TransA}, float @var{alpha}, const gsl_matrix_float * @var{A}, const gsl_vector_float * @var{x}, float @var{beta}, gsl_vector_float * @var{y})
@deftypefunx int gsl_blas_dgemv (CBLAS_TRANSPOSE_t @var{TransA}, double @var{alpha}, const gsl_matrix * @var{A}, const gsl_vector * @var{x}, double @var{beta}, gsl_vector * @var{y})
@deftypefunx int gsl_blas_cgemv (CBLAS_TRANSPOSE_t @var{TransA}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_vector_complex_float * @var{x}, const gsl_complex_float @var{beta}, gsl_vector_complex_float * @var{y})
@deftypefunx int gsl_blas_zgemv (CBLAS_TRANSPOSE_t @var{TransA}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_vector_complex * @var{x}, const gsl_complex @var{beta}, gsl_vector_complex * @var{y})
@cindex GEMV, Level-2 BLAS
These functions compute the matrix-vector product and sum @math{y =
\alpha op(A) x + \beta y}, where @math{op(A) = A},
@math{A^T}, @math{A^H} for @var{TransA} = @code{CblasNoTrans},
@code{CblasTrans}, @code{CblasConjTrans}.
@end deftypefun


@deftypefun int gsl_blas_strmv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_float * @var{A}, gsl_vector_float * @var{x})
@deftypefunx int gsl_blas_dtrmv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix * @var{A}, gsl_vector * @var{x})
@deftypefunx int gsl_blas_ctrmv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_complex_float * @var{A}, gsl_vector_complex_float * @var{x})
@deftypefunx int gsl_blas_ztrmv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_complex * @var{A}, gsl_vector_complex * @var{x})
@cindex TRMV, Level-2 BLAS
These functions compute the matrix-vector product 
@math{x = op(A) x} for the triangular matrix @var{A}, where
@math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
@code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}.  When
@var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
diagonal of the matrix is used, but if @var{Diag} is @code{CblasUnit}
then the diagonal elements of the matrix @var{A} are taken as unity and
are not referenced.
@end deftypefun


@deftypefun int gsl_blas_strsv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_float * @var{A}, gsl_vector_float * @var{x})
@deftypefunx int gsl_blas_dtrsv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix * @var{A}, gsl_vector * @var{x})
@deftypefunx int gsl_blas_ctrsv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_complex_float * @var{A}, gsl_vector_complex_float * @var{x})
@deftypefunx int gsl_blas_ztrsv (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_matrix_complex * @var{A}, gsl_vector_complex * @var{x})
@cindex TRSV, Level-2 BLAS
These functions compute @math{inv(op(A)) x} for @var{x}, where
@math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
@code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}.  When
@var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
diagonal of the matrix is used, but if @var{Diag} is @code{CblasUnit}
then the diagonal elements of the matrix @var{A} are taken as unity and
are not referenced.
@end deftypefun


@deftypefun int gsl_blas_ssymv (CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_matrix_float * @var{A}, const gsl_vector_float * @var{x}, float @var{beta}, gsl_vector_float * @var{y})
@deftypefunx int gsl_blas_dsymv (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_matrix * @var{A}, const gsl_vector * @var{x}, double @var{beta}, gsl_vector * @var{y})
@cindex SYMV, Level-2 BLAS
These functions compute the matrix-vector product and sum @math{y =
\alpha A x + \beta y} for the symmetric matrix @var{A}.  Since the
matrix @var{A} is symmetric only its upper half or lower half need to be
stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
and diagonal of @var{A} are used, and when @var{Uplo} is
@code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.
@end deftypefun

@deftypefun int gsl_blas_chemv (CBLAS_UPLO_t @var{Uplo}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_vector_complex_float * @var{x}, const gsl_complex_float @var{beta}, gsl_vector_complex_float * @var{y})
@deftypefunx int gsl_blas_zhemv (CBLAS_UPLO_t @var{Uplo}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_vector_complex * @var{x}, const gsl_complex @var{beta}, gsl_vector_complex * @var{y})
@cindex HEMV, Level-2 BLAS
These functions compute the matrix-vector product and sum @math{y =
\alpha A x + \beta y} for the hermitian matrix @var{A}.  Since the
matrix @var{A} is hermitian only its upper half or lower half need to be
stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
and diagonal of @var{A} are used, and when @var{Uplo} is
@code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.  The imaginary elements of the diagonal are automatically assumed
to be zero and are not referenced.
@end deftypefun

@deftypefun int gsl_blas_sger (float @var{alpha}, const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, gsl_matrix_float * @var{A})
@deftypefunx int gsl_blas_dger (double @var{alpha}, const gsl_vector * @var{x}, const gsl_vector * @var{y}, gsl_matrix * @var{A})
@deftypefunx int gsl_blas_cgeru (const gsl_complex_float @var{alpha}, const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_matrix_complex_float * @var{A})
@deftypefunx int gsl_blas_zgeru (const gsl_complex @var{alpha}, const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_matrix_complex * @var{A})
@cindex GER, Level-2 BLAS
@cindex GERU, Level-2 BLAS
These functions compute the rank-1 update @math{A = \alpha x y^T + A} of
the matrix @var{A}.
@end deftypefun

@deftypefun int gsl_blas_cgerc (const gsl_complex_float @var{alpha}, const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_matrix_complex_float * @var{A})
@deftypefunx int gsl_blas_zgerc (const gsl_complex @var{alpha}, const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_matrix_complex * @var{A})
@cindex GERC, Level-2 BLAS
These functions compute the conjugate rank-1 update @math{A = \alpha x
y^H + A} of the matrix @var{A}.
@end deftypefun

@deftypefun int gsl_blas_ssyr (CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_vector_float * @var{x}, gsl_matrix_float * @var{A})
@deftypefunx int gsl_blas_dsyr (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_vector * @var{x}, gsl_matrix * @var{A})
@cindex SYR, Level-2 BLAS
These functions compute the symmetric rank-1 update @math{A = \alpha x
x^T + A} of the symmetric matrix @var{A}.  Since the matrix @var{A} is
symmetric only its upper half or lower half need to be stored.  When
@var{Uplo} is @code{CblasUpper} then the upper triangle and diagonal of
@var{A} are used, and when @var{Uplo} is @code{CblasLower} then the
lower triangle and diagonal of @var{A} are used.
@end deftypefun

@deftypefun int gsl_blas_cher (CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_vector_complex_float * @var{x}, gsl_matrix_complex_float * @var{A})
@deftypefunx int gsl_blas_zher (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_vector_complex * @var{x}, gsl_matrix_complex * @var{A})
@cindex HER, Level-2 BLAS
These functions compute the hermitian rank-1 update @math{A = \alpha x
x^H + A} of the hermitian matrix @var{A}.  Since the matrix @var{A} is
hermitian only its upper half or lower half need to be stored.  When
@var{Uplo} is @code{CblasUpper} then the upper triangle and diagonal of
@var{A} are used, and when @var{Uplo} is @code{CblasLower} then the
lower triangle and diagonal of @var{A} are used.  The imaginary elements
of the diagonal are automatically set to zero.
@end deftypefun

@deftypefun int gsl_blas_ssyr2 (CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, gsl_matrix_float * @var{A})
@deftypefunx int gsl_blas_dsyr2 (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_vector * @var{x}, const gsl_vector * @var{y}, gsl_matrix * @var{A})
@cindex SYR2, Level-2 BLAS
These functions compute the symmetric rank-2 update @math{A = \alpha x
y^T + \alpha y x^T + A} of the symmetric matrix @var{A}.  Since the
matrix @var{A} is symmetric only its upper half or lower half need to be
stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
and diagonal of @var{A} are used, and when @var{Uplo} is
@code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.
@end deftypefun

@deftypefun int gsl_blas_cher2 (CBLAS_UPLO_t @var{Uplo}, const gsl_complex_float @var{alpha}, const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_matrix_complex_float * @var{A})
@deftypefunx int gsl_blas_zher2 (CBLAS_UPLO_t @var{Uplo}, const gsl_complex @var{alpha}, const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_matrix_complex * @var{A})
@cindex HER2, Level-2 BLAS
These functions compute the hermitian rank-2 update @math{A = \alpha x
y^H + \alpha^* y x^H + A} of the hermitian matrix @var{A}.  Since the
matrix @var{A} is hermitian only its upper half or lower half need to be
stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
and diagonal of @var{A} are used, and when @var{Uplo} is
@code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.  The imaginary elements of the diagonal are automatically set to zero.
@end deftypefun

@node Level 3 GSL BLAS Interface
@subsection Level 3


@deftypefun int gsl_blas_sgemm (CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_TRANSPOSE_t @var{TransB}, float @var{alpha}, const gsl_matrix_float * @var{A}, const gsl_matrix_float * @var{B}, float @var{beta}, gsl_matrix_float * @var{C})
@deftypefunx int gsl_blas_dgemm (CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_TRANSPOSE_t @var{TransB}, double @var{alpha}, const gsl_matrix * @var{A}, const gsl_matrix * @var{B}, double @var{beta}, gsl_matrix * @var{C})
@deftypefunx int gsl_blas_cgemm (CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_TRANSPOSE_t @var{TransB}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_matrix_complex_float * @var{B}, const gsl_complex_float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zgemm (CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_TRANSPOSE_t @var{TransB}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{B}, const gsl_complex @var{beta}, gsl_matrix_complex * @var{C})
@cindex GEMM, Level-3 BLAS
These functions compute the matrix-matrix product and sum @math{C =
\alpha op(A) op(B) + \beta C} where @math{op(A) = A}, @math{A^T},
@math{A^H} for @var{TransA} = @code{CblasNoTrans}, @code{CblasTrans},
@code{CblasConjTrans} and similarly for the parameter @var{TransB}.
@end deftypefun


@deftypefun int gsl_blas_ssymm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_matrix_float * @var{A}, const gsl_matrix_float * @var{B}, float @var{beta}, gsl_matrix_float * @var{C})
@deftypefunx int gsl_blas_dsymm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_matrix * @var{A}, const gsl_matrix * @var{B}, double @var{beta}, gsl_matrix * @var{C})
@deftypefunx int gsl_blas_csymm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_matrix_complex_float * @var{B}, const gsl_complex_float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zsymm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{B}, const gsl_complex @var{beta}, gsl_matrix_complex * @var{C})
@cindex SYMM, Level-3 BLAS
These functions compute the matrix-matrix product and sum @math{C =
\alpha A B + \beta C} for @var{Side} is @code{CblasLeft} and @math{C =
\alpha B A + \beta C} for @var{Side} is @code{CblasRight}, where the
matrix @var{A} is symmetric.  When @var{Uplo} is @code{CblasUpper} then
the upper triangle and diagonal of @var{A} are used, and when @var{Uplo}
is @code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.
@end deftypefun

@deftypefun int gsl_blas_chemm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_matrix_complex_float * @var{B}, const gsl_complex_float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zhemm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{B}, const gsl_complex @var{beta}, gsl_matrix_complex * @var{C})
@cindex HEMM, Level-3 BLAS
These functions compute the matrix-matrix product and sum @math{C =
\alpha A B + \beta C} for @var{Side} is @code{CblasLeft} and @math{C =
\alpha B A + \beta C} for @var{Side} is @code{CblasRight}, where the
matrix @var{A} is hermitian.  When @var{Uplo} is @code{CblasUpper} then
the upper triangle and diagonal of @var{A} are used, and when @var{Uplo}
is @code{CblasLower} then the lower triangle and diagonal of @var{A} are
used.  The imaginary elements of the diagonal are automatically set to
zero.
@end deftypefun

@deftypefun int gsl_blas_strmm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, float @var{alpha}, const gsl_matrix_float * @var{A}, gsl_matrix_float * @var{B})
@deftypefunx int gsl_blas_dtrmm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, double @var{alpha}, const gsl_matrix * @var{A}, gsl_matrix * @var{B})
@deftypefunx int gsl_blas_ctrmm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, gsl_matrix_complex_float * @var{B})
@deftypefunx int gsl_blas_ztrmm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, gsl_matrix_complex * @var{B})
@cindex TRMM, Level-3 BLAS
These functions compute the matrix-matrix product @math{B = \alpha op(A)
B} for @var{Side} is @code{CblasLeft} and @math{B = \alpha B op(A)} for
@var{Side} is @code{CblasRight}.  The matrix @var{A} is triangular and
@math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
@code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}. When
@var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
diagonal of @var{A} is used, but if @var{Diag} is @code{CblasUnit} then
the diagonal elements of the matrix @var{A} are taken as unity and are
not referenced.
@end deftypefun


@deftypefun int gsl_blas_strsm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, float @var{alpha}, const gsl_matrix_float * @var{A}, gsl_matrix_float * @var{B})
@deftypefunx int gsl_blas_dtrsm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, double @var{alpha}, const gsl_matrix * @var{A}, gsl_matrix * @var{B})
@deftypefunx int gsl_blas_ctrsm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, gsl_matrix_complex_float * @var{B})
@deftypefunx int gsl_blas_ztrsm (CBLAS_SIDE_t @var{Side}, CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{TransA}, CBLAS_DIAG_t @var{Diag}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, gsl_matrix_complex * @var{B})
@cindex TRSM, Level-3 BLAS
These functions compute the inverse-matrix matrix product 
@math{B = \alpha op(inv(A))B} for @var{Side} is 
@code{CblasLeft} and @math{B = \alpha B op(inv(A))} for
@var{Side} is @code{CblasRight}.  The matrix @var{A} is triangular and
@math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
@code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}. When
@var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
diagonal of @var{A} is used, but if @var{Diag} is @code{CblasUnit} then
the diagonal elements of the matrix @var{A} are taken as unity and are
not referenced.
@end deftypefun

@deftypefun int gsl_blas_ssyrk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, float @var{alpha}, const gsl_matrix_float * @var{A}, float @var{beta}, gsl_matrix_float * @var{C})
@deftypefunx int gsl_blas_dsyrk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, double @var{alpha}, const gsl_matrix * @var{A}, double @var{beta}, gsl_matrix * @var{C})
@deftypefunx int gsl_blas_csyrk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_complex_float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zsyrk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_complex @var{beta}, gsl_matrix_complex * @var{C})
@cindex SYRK, Level-3 BLAS
These functions compute a rank-k update of the symmetric matrix @var{C},
@math{C = \alpha A A^T + \beta C} when @var{Trans} is
@code{CblasNoTrans} and @math{C = \alpha A^T A + \beta C} when
@var{Trans} is @code{CblasTrans}.  Since the matrix @var{C} is symmetric
only its upper half or lower half need to be stored.  When @var{Uplo} is
@code{CblasUpper} then the upper triangle and diagonal of @var{C} are
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
and diagonal of @var{C} are used.
@end deftypefun

@deftypefun int gsl_blas_cherk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, float @var{alpha}, const gsl_matrix_complex_float * @var{A}, float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zherk (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, double @var{alpha}, const gsl_matrix_complex * @var{A}, double @var{beta}, gsl_matrix_complex * @var{C})
@cindex HERK, Level-3 BLAS
These functions compute a rank-k update of the hermitian matrix @var{C},
@math{C = \alpha A A^H + \beta C} when @var{Trans} is
@code{CblasNoTrans} and @math{C = \alpha A^H A + \beta C} when
@var{Trans} is @code{CblasConjTrans}.  Since the matrix @var{C} is hermitian
only its upper half or lower half need to be stored.  When @var{Uplo} is
@code{CblasUpper} then the upper triangle and diagonal of @var{C} are
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
and diagonal of @var{C} are used.  The imaginary elements of the
diagonal are automatically set to zero.
@end deftypefun

@deftypefun int gsl_blas_ssyr2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, float @var{alpha}, const gsl_matrix_float * @var{A}, const gsl_matrix_float * @var{B}, float @var{beta}, gsl_matrix_float * @var{C})
@deftypefunx int gsl_blas_dsyr2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, double @var{alpha}, const gsl_matrix * @var{A}, const gsl_matrix * @var{B}, double @var{beta}, gsl_matrix * @var{C})
@deftypefunx int gsl_blas_csyr2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_matrix_complex_float * @var{B}, const gsl_complex_float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zsyr2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{B}, const gsl_complex @var{beta}, gsl_matrix_complex * @var{C})
@cindex SYR2K, Level-3 BLAS
These functions compute a rank-2k update of the symmetric matrix @var{C},
@math{C = \alpha A B^T + \alpha B A^T + \beta C} when @var{Trans} is
@code{CblasNoTrans} and @math{C = \alpha A^T B + \alpha B^T A + \beta C} when
@var{Trans} is @code{CblasTrans}.  Since the matrix @var{C} is symmetric
only its upper half or lower half need to be stored.  When @var{Uplo} is
@code{CblasUpper} then the upper triangle and diagonal of @var{C} are
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
and diagonal of @var{C} are used.
@end deftypefun

@deftypefun int gsl_blas_cher2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex_float @var{alpha}, const gsl_matrix_complex_float * @var{A}, const gsl_matrix_complex_float * @var{B}, float @var{beta}, gsl_matrix_complex_float * @var{C})
@deftypefunx int gsl_blas_zher2k (CBLAS_UPLO_t @var{Uplo}, CBLAS_TRANSPOSE_t @var{Trans}, const gsl_complex @var{alpha}, const gsl_matrix_complex * @var{A}, const gsl_matrix_complex * @var{B}, double @var{beta}, gsl_matrix_complex * @var{C})
@cindex HER2K, Level-3 BLAS
These functions compute a rank-2k update of the hermitian matrix @var{C},
@math{C = \alpha A B^H + \alpha^* B A^H + \beta C} when @var{Trans} is
@code{CblasNoTrans} and @math{C = \alpha A^H B + \alpha^* B^H A + \beta C} when
@var{Trans} is @code{CblasConjTrans}.  Since the matrix @var{C} is hermitian
only its upper half or lower half need to be stored.  When @var{Uplo} is
@code{CblasUpper} then the upper triangle and diagonal of @var{C} are
used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
and diagonal of @var{C} are used.  The imaginary elements of the
diagonal are automatically set to zero.
@end deftypefun

@node BLAS Examples
@section Examples

The following program computes the product of two matrices using the
Level-3 @sc{blas} function @sc{dgemm},
@tex
\beforedisplay
$$
\left(
\matrix{0.11&0.12&0.13\cr
0.21&0.22&0.23\cr}
\right)
\left(
\matrix{1011&1012\cr
1021&1022\cr
1031&1031\cr}
\right)
=
\left(
\matrix{367.76&368.12\cr
674.06&674.72\cr}
\right)
$$
\afterdisplay
@end tex
@ifinfo

@example
[ 0.11 0.12 0.13 ]  [ 1011 1012 ]     [ 367.76 368.12 ]
[ 0.21 0.22 0.23 ]  [ 1021 1022 ]  =  [ 674.06 674.72 ]
                    [ 1031 1032 ]
@end example

@end ifinfo
@noindent
The matrices are stored in row major order, according to the C convention 
for arrays.

@example
@verbatiminclude examples/blas.c
@end example

@noindent
Here is the output from the program,

@example
$ ./a.out
@verbatiminclude examples/blas.txt
@end example

@node BLAS References and Further Reading
@section References and Further Reading

Information on the @sc{blas} standards, including both the legacy and
updated interface standards, is available online from the @sc{blas}
Homepage and @sc{blas} Technical Forum web-site.

@itemize @w{}
@item
@cite{BLAS Homepage} @*
@uref{http://www.netlib.org/blas/}
@item
@cite{BLAS Technical Forum} @*
@uref{http://www.netlib.org/blas/blast-forum/}
@end itemize

@noindent
The following papers contain the specifications for Level 1, Level 2 and
Level 3 @sc{blas}.

@itemize @w{}
@item
C. Lawson, R. Hanson, D. Kincaid, F. Krogh, ``Basic Linear Algebra
Subprograms for Fortran Usage'', @cite{ACM Transactions on Mathematical
Software}, Vol.@: 5 (1979), Pages 308--325.

@item
J.J. Dongarra, J. DuCroz, S. Hammarling, R. Hanson, ``An Extended Set of
Fortran Basic Linear Algebra Subprograms'', @cite{ACM Transactions on
Mathematical Software}, Vol.@: 14, No.@: 1 (1988), Pages 1--32.

@item
J.J. Dongarra, I. Duff, J. DuCroz, S. Hammarling, ``A Set of
Level 3 Basic Linear Algebra Subprograms'', @cite{ACM Transactions on
Mathematical Software}, Vol.@: 16 (1990), Pages 1--28.
@end itemize

@noindent
Postscript versions of the latter two papers are available from
@uref{http://www.netlib.org/blas/}. A @sc{cblas} wrapper for Fortran @sc{blas}
libraries is available from the same location.