htswanalysis/MACS/lib/gsl/gsl-1.11/doc/blas.texi

   1 @cindex linear algebra, BLAS
   2 @cindex matrix, operations
   3 @cindex vector, operations
   4 @cindex BLAS
   5 @cindex CBLAS
   6 @cindex Basic Linear Algebra Subroutines (BLAS)
   7
   8 The Basic Linear Algebra Subprograms (@sc{blas}) define a set of fundamental
   9 operations on vectors and matrices which can be used to create optimized
  10 higher-level linear algebra functionality.
  11
  12 The library provides a low-level layer which corresponds directly to the
  13 C-language @sc{blas} standard, referred to here as ``@sc{cblas}'', and a
  14 higher-level interface for operations on GSL vectors and matrices.
  15 Users who are interested in simple operations on GSL vector and matrix
  16 objects should use the high-level layer, which is declared in the file
  17 @code{gsl_blas.h}.  This should satisfy the needs of most users.  Note
  18 that GSL matrices are implemented using dense-storage so the interface
  19 only includes the corresponding dense-storage @sc{blas} functions.  The full
  20 @sc{blas} functionality for band-format and packed-format matrices is
  21 available through the low-level @sc{cblas} interface.
  22
  23 The interface for the @code{gsl_cblas} layer is specified in the file
  24 @code{gsl_cblas.h}.  This interface corresponds to the @sc{blas} Technical
  25 Forum's draft standard for the C interface to legacy @sc{blas}
  26 implementations. Users who have access to other conforming @sc{cblas}
  27 implementations can use these in place of the version provided by the
  28 library.  Note that users who have only a Fortran @sc{blas} library can
  29 use a @sc{cblas} conformant wrapper to convert it into a @sc{cblas}
  30 library.  A reference @sc{cblas} wrapper for legacy Fortran
  31 implementations exists as part of the draft @sc{cblas} standard and can
  32 be obtained from Netlib.  The complete set of @sc{cblas} functions is
  33 listed in an appendix (@pxref{GSL CBLAS Library}).
  34
  35 There are three levels of @sc{blas} operations,
  36
  37 @table @b
  38 @item Level 1
  39 Vector operations, e.g. @math{y = \alpha x + y}
  40 @item Level 2
  41 Matrix-vector operations, e.g. @math{y = \alpha A x + \beta y}
  42 @item Level 3
  43 Matrix-matrix operations, e.g. @math{C = \alpha A B + C}
  44 @end table
  45
  46 @noindent
  47 Each routine has a name which specifies the operation, the type of
  48 matrices involved and their precisions.  Some of the most common
  49 operations and their names are given below,
  50
  51 @table @b
  52 @item DOT
  53 scalar product, @math{x^T y}
  54 @item AXPY
  55 vector sum, @math{\alpha x + y}
  56 @item MV
  57 matrix-vector product, @math{A x}
  58 @item SV
  59 matrix-vector solve, @math{inv(A) x}
  60 @item MM
  61 matrix-matrix product, @math{A B}
  62 @item SM
  63 matrix-matrix solve, @math{inv(A) B}
  64 @end table
  65
  66 @noindent
  67 The types of matrices are,
  68
  69 @table @b
  70 @item GE
  71 general
  72 @item GB
  73 general band
  74 @item SY
  75 symmetric
  76 @item SB
  77 symmetric band
  78 @item SP
  79 symmetric packed
  80 @item HE
  81 hermitian
  82 @item HB
  83 hermitian band
  84 @item HP
  85 hermitian packed
  86 @item TR
  87 triangular
  88 @item TB
  89 triangular band
  90 @item TP
  91 triangular packed
  92 @end table
  93
  94 @noindent
  95 Each operation is defined for four precisions,
  96
  97 @table @b
  98 @item S
  99 single real
 100 @item D
 101 double real
 102 @item C
 103 single complex
 104 @item Z
 105 double complex
 106 @end table
 107
 108 @noindent
 109 Thus, for example, the name @sc{sgemm} stands for ``single-precision
 110 general matrix-matrix multiply'' and @sc{zgemm} stands for
 111 ``double-precision complex matrix-matrix multiply''.
 112
 113 @menu
 114 * GSL BLAS Interface::
 115 * BLAS Examples::
 116 * BLAS References and Further Reading::
 117 @end menu
 118
 119 @node GSL BLAS Interface
 120 @section GSL BLAS Interface
 121
 122 GSL provides dense vector and matrix objects, based on the relevant
 123 built-in types.  The library provides an interface to the @sc{blas}
 124 operations which apply to these objects.  The interface to this
 125 functionality is given in the file @code{gsl_blas.h}.
 126
 127 @comment CblasNoTrans, CblasTrans, CblasConjTrans
 128 @comment CblasUpper, CblasLower
 129 @comment CblasNonUnit, CblasUnit
 130 @comment CblasLeft, CblasRight
 131
 132 @menu
 133 * Level 1 GSL BLAS Interface::
 134 * Level 2 GSL BLAS Interface::
 135 * Level 3 GSL BLAS Interface::
 136 @end menu
 137
 138 @node Level 1 GSL BLAS Interface
 139 @subsection Level 1
 140
 141 @deftypefun int gsl_blas_sdsdot (float @var{alpha}, const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, float * @var{result})
 142 @cindex DOT, Level-1 BLAS
 143 This function computes the sum @math{\alpha + x^T y} for the vectors
 144 @var{x} and @var{y}, returning the result in @var{result}.
 145 @end deftypefun
 146
 147 @deftypefun int gsl_blas_sdot (const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, float * @var{result})
 148 @deftypefunx int gsl_blas_dsdot (const gsl_vector_float * @var{x}, const gsl_vector_float * @var{y}, double * @var{result})
 149 @deftypefunx int gsl_blas_ddot (const gsl_vector * @var{x}, const gsl_vector * @var{y}, double * @var{result})
 150 These functions compute the scalar product @math{x^T y} for the vectors
 151 @var{x} and @var{y}, returning the result in @var{result}.
 152 @end deftypefun
 153
 154 @deftypefun int gsl_blas_cdotu (const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_complex_float * @var{dotu})
 155 @deftypefunx int gsl_blas_zdotu (const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_complex * @var{dotu})
 156 These functions compute the complex scalar product @math{x^T y} for the
 157 vectors @var{x} and @var{y}, returning the result in @var{result}
 158 @end deftypefun
 159
 160 @deftypefun int gsl_blas_cdotc (const gsl_vector_complex_float * @var{x}, const gsl_vector_complex_float * @var{y}, gsl_complex_float * @var{dotc})
 161 @deftypefunx int gsl_blas_zdotc (const gsl_vector_complex * @var{x}, const gsl_vector_complex * @var{y}, gsl_complex * @var{dotc})
 162 These functions compute the complex conjugate scalar product @math{x^H
 163 y} for the vectors @var{x} and @var{y}, returning the result in
 164 @var{result}
 165 @end deftypefun
 166
 167 @deftypefun float gsl_blas_snrm2 (const gsl_vector_float * @var{x})
 168 @deftypefunx double gsl_blas_dnrm2 (const gsl_vector * @var{x})
 169 @cindex NRM2, Level-1 BLAS
 170 These functions compute the Euclidean norm
 171 @c{$||x||_2 = \sqrt{\sum x_i^2}$}
 172 @math{||x||_2 = \sqrt @{\sum x_i^2@}} of the vector @var{x}.
 173 @end deftypefun
 174
 175 @deftypefun float gsl_blas_scnrm2 (const gsl_vector_complex_float * @var{x})
 176 @deftypefunx double gsl_blas_dznrm2 (const gsl_vector_complex * @var{x})
 177 These functions compute the Euclidean norm of the complex vector @var{x},
 178 @tex
 179 \beforedisplay
 180 $$
 181 ||x||_2 = \sqrt{\sum (\Re(x_i)^2 + \Im(x_i)^2)}.
 182 $$
 183 \afterdisplay
 184 @end tex
 185 @ifinfo
 186
 187 @example
 188 ||x||_2 = \sqrt @{\sum (\Re(x_i)^2 + \Im(x_i)^2)@}.
 189 @end example
 190 @end ifinfo
 191 @end deftypefun
 192
 193 @deftypefun float gsl_blas_sasum (const gsl_vector_float * @var{x})
 194 @deftypefunx double gsl_blas_dasum (const gsl_vector * @var{x})
 195 @cindex ASUM, Level-1 BLAS
 196 These functions compute the absolute sum @math{\sum |x_i|} of the
 197 elements of the vector @var{x}.
 198 @end deftypefun
 199
 200 @deftypefun float gsl_blas_scasum (const gsl_vector_complex_float * @var{x})
 201 @deftypefunx double gsl_blas_dzasum (const gsl_vector_complex * @var{x})
 202 These functions compute the sum of the magnitudes of the real and
 203 imaginary parts of the complex vector @var{x},
 204 @c{$\sum \left( |\Re(x_i)| + |\Im(x_i)| \right)$}
 205 @math{\sum |\Re(x_i)| + |\Im(x_i)|}.
 206 @end deftypefun
 207
 208 @deftypefun CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * @var{x})
 209 @deftypefunx CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * @var{x})
 210 @deftypefunx CBLAS_INDEX_t gsl_blas_icamax (const gsl_vector_complex_float * @var{x})
 211 @deftypefunx CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * @var{x})
 212 @cindex AMAX, Level-1 BLAS
 213 These functions return the index of the largest element of the vector
 214 @var{x}. The largest element is determined by its absolute magnitude for
 215 real vectors and by the sum of the magnitudes of the real and imaginary
 216 parts @math{|\Re(x_i)| + |\Im(x_i)|} for complex vectors.  If the
 217 largest value occurs several times then the index of the first
 218 occurrence is returned.
 219 @end deftypefun
 220
 221 @deftypefun int gsl_blas_sswap (gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
 222 @deftypefunx int gsl_blas_dswap (gsl_vector * @var{x}, gsl_vector * @var{y})
 223 @deftypefunx int gsl_blas_cswap (gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
 224 @deftypefunx int gsl_blas_zswap (gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
 225 @cindex SWAP, Level-1 BLAS
 226 These functions exchange the elements of the vectors @var{x} and @var{y}.
 227 @end deftypefun
 228
 229 @deftypefun int gsl_blas_scopy (const gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
 230 @deftypefunx int gsl_blas_dcopy (const gsl_vector * @var{x}, gsl_vector * @var{y})
 231 @deftypefunx int gsl_blas_ccopy (const gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
 232 @deftypefunx int gsl_blas_zcopy (const gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
 233 @cindex COPY, Level-1 BLAS
 234 These functions copy the elements of the vector @var{x} into the vector
 235 @var{y}.
 236 @end deftypefun
 237
 238
 239 @deftypefun int gsl_blas_saxpy (float @var{alpha}, const gsl_vector_float * @var{x}, gsl_vector_float * @var{y})
 240 @deftypefunx int gsl_blas_daxpy (double @var{alpha}, const gsl_vector * @var{x}, gsl_vector * @var{y})
 241 @deftypefunx int gsl_blas_caxpy (const gsl_complex_float @var{alpha}, const gsl_vector_complex_float * @var{x}, gsl_vector_complex_float * @var{y})
 242 @deftypefunx int gsl_blas_zaxpy (const gsl_complex @var{alpha}, const gsl_vector_complex * @var{x}, gsl_vector_complex * @var{y})
 243 @cindex AXPY, Level-1 BLAS
 244 @cindex DAXPY, Level-1 BLAS
 245 @cindex SAXPY, Level-1 BLAS
 246 These functions compute the sum @math{y = \alpha x + y} for the vectors
 247 @var{x} and @var{y}.
 248 @end deftypefun
 249
 250 @deftypefun void gsl_blas_sscal (float @var{alpha}, gsl_vector_float * @var{x})
 251 @deftypefunx void gsl_blas_dscal (double @var{alpha}, gsl_vector * @var{x})
 252 @deftypefunx void gsl_blas_cscal (const gsl_complex_float @var{alpha}, gsl_vector_complex_float * @var{x})
 253 @deftypefunx void gsl_blas_zscal (const gsl_complex @var{alpha}, gsl_vector_complex * @var{x})
 254 @deftypefunx void gsl_blas_csscal (float @var{alpha}, gsl_vector_complex_float * @var{x})
 255 @deftypefunx void gsl_blas_zdscal (double @var{alpha}, gsl_vector_complex * @var{x})
 256 @cindex SCAL, Level-1 BLAS
 257 These functions rescale the vector @var{x} by the multiplicative factor
 258 @var{alpha}.
 259 @end deftypefun
 260
 261 @deftypefun int gsl_blas_srotg (float @var{a}[], float @var{b}[], float @var{c}[], float @var{s}[])
 262 @deftypefunx int gsl_blas_drotg (double @var{a}[], double @var{b}[], double @var{c}[], double @var{s}[])
 263 @cindex ROTG, Level-1 BLAS
 264 @cindex Givens Rotation, BLAS
 265 These functions compute a Givens rotation @math{(c,s)} which zeroes the
 266 vector @math{(a,b)},
 267 @tex
 268 \beforedisplay
 269 $$
 270 \left(
 271 \matrix{c&s\cr
 272 -s&c\cr}
 273 \right)
 274 \left(
 275 \matrix{a\cr
 276 b\cr}
 277 \right)
 278 =
 279 \left(
 280 \matrix{r'\cr
 281 0\cr}
 282 \right)
 283 $$
 284 \afterdisplay
 285 @end tex
 286 @ifinfo
 287
 288 @example
 289 [  c  s ] [ a ] = [ r ]
 290 [ -s  c ] [ b ]   [ 0 ]
 291 @end example
 292
 293 @end ifinfo
 294 @noindent
 295 The variables @var{a} and @var{b} are overwritten by the routine.
 296 @end deftypefun
 297
 298 @deftypefun int gsl_blas_srot (gsl_vector_float * @var{x}, gsl_vector_float * @var{y}, float @var{c}, float @var{s})
 299 @deftypefunx int gsl_blas_drot (gsl_vector * @var{x}, gsl_vector * @var{y}, const double @var{c}, const double @var{s})
 300 These functions apply a Givens rotation @math{(x', y') = (c x + s y, -s
 301 x + c y)} to the vectors @var{x}, @var{y}.
 302 @end deftypefun
 303
 304 @deftypefun int gsl_blas_srotmg (float @var{d1}[], float @var{d2}[], float @var{b1}[], float @var{b2}, float @var{P}[])
 305 @deftypefunx int gsl_blas_drotmg (double @var{d1}[], double @var{d2}[], double @var{b1}[], double @var{b2}, double @var{P}[])
 306 @cindex Modified Givens Rotation, BLAS
 307 @cindex Givens Rotation, Modified, BLAS
 308 These functions compute a modified Givens transformation.  The modified
 309 Givens transformation is defined in the original Level-1 @sc{blas}
 310 specification, given in the references.
 311 @end deftypefun
 312
 313 @deftypefun int gsl_blas_srotm (gsl_vector_float * @var{x}, gsl_vector_float * @var{y}, const float @var{P}[])
 314 @deftypefunx int gsl_blas_drotm (gsl_vector * @var{x}, gsl_vector * @var{y}, const double @var{P}[])
 315 These functions apply a modified Givens transformation.
 316 @end deftypefun
 317
 318 @node Level 2 GSL BLAS Interface
 319 @subsection Level 2
 320
 321 @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})
 322 @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})
 323 @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})
 324 @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})
 325 @cindex GEMV, Level-2 BLAS
 326 These functions compute the matrix-vector product and sum @math{y =
 327 \alpha op(A) x + \beta y}, where @math{op(A) = A},
 328 @math{A^T}, @math{A^H} for @var{TransA} = @code{CblasNoTrans},
 329 @code{CblasTrans}, @code{CblasConjTrans}.
 330 @end deftypefun
 331
 332
 333 @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})
 334 @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})
 335 @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})
 336 @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})
 337 @cindex TRMV, Level-2 BLAS
 338 These functions compute the matrix-vector product
 339 @math{x = op(A) x} for the triangular matrix @var{A}, where
 340 @math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
 341 @code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}.  When
 342 @var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
 343 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 344 of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
 345 diagonal of the matrix is used, but if @var{Diag} is @code{CblasUnit}
 346 then the diagonal elements of the matrix @var{A} are taken as unity and
 347 are not referenced.
 348 @end deftypefun
 349
 350
 351 @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})
 352 @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})
 353 @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})
 354 @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})
 355 @cindex TRSV, Level-2 BLAS
 356 These functions compute @math{inv(op(A)) x} for @var{x}, where
 357 @math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
 358 @code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}.  When
 359 @var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
 360 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 361 of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
 362 diagonal of the matrix is used, but if @var{Diag} is @code{CblasUnit}
 363 then the diagonal elements of the matrix @var{A} are taken as unity and
 364 are not referenced.
 365 @end deftypefun
 366
 367
 368 @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})
 369 @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})
 370 @cindex SYMV, Level-2 BLAS
 371 These functions compute the matrix-vector product and sum @math{y =
 372 \alpha A x + \beta y} for the symmetric matrix @var{A}.  Since the
 373 matrix @var{A} is symmetric only its upper half or lower half need to be
 374 stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
 375 and diagonal of @var{A} are used, and when @var{Uplo} is
 376 @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 377 used.
 378 @end deftypefun
 379
 380 @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})
 381 @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})
 382 @cindex HEMV, Level-2 BLAS
 383 These functions compute the matrix-vector product and sum @math{y =
 384 \alpha A x + \beta y} for the hermitian matrix @var{A}.  Since the
 385 matrix @var{A} is hermitian only its upper half or lower half need to be
 386 stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
 387 and diagonal of @var{A} are used, and when @var{Uplo} is
 388 @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 389 used.  The imaginary elements of the diagonal are automatically assumed
 390 to be zero and are not referenced.
 391 @end deftypefun
 392
 393 @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})
 394 @deftypefunx int gsl_blas_dger (double @var{alpha}, const gsl_vector * @var{x}, const gsl_vector * @var{y}, gsl_matrix * @var{A})
 395 @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})
 396 @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})
 397 @cindex GER, Level-2 BLAS
 398 @cindex GERU, Level-2 BLAS
 399 These functions compute the rank-1 update @math{A = \alpha x y^T + A} of
 400 the matrix @var{A}.
 401 @end deftypefun
 402
 403 @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})
 404 @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})
 405 @cindex GERC, Level-2 BLAS
 406 These functions compute the conjugate rank-1 update @math{A = \alpha x
 407 y^H + A} of the matrix @var{A}.
 408 @end deftypefun
 409
 410 @deftypefun int gsl_blas_ssyr (CBLAS_UPLO_t @var{Uplo}, float @var{alpha}, const gsl_vector_float * @var{x}, gsl_matrix_float * @var{A})
 411 @deftypefunx int gsl_blas_dsyr (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_vector * @var{x}, gsl_matrix * @var{A})
 412 @cindex SYR, Level-2 BLAS
 413 These functions compute the symmetric rank-1 update @math{A = \alpha x
 414 x^T + A} of the symmetric matrix @var{A}.  Since the matrix @var{A} is
 415 symmetric only its upper half or lower half need to be stored.  When
 416 @var{Uplo} is @code{CblasUpper} then the upper triangle and diagonal of
 417 @var{A} are used, and when @var{Uplo} is @code{CblasLower} then the
 418 lower triangle and diagonal of @var{A} are used.
 419 @end deftypefun
 420
 421 @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})
 422 @deftypefunx int gsl_blas_zher (CBLAS_UPLO_t @var{Uplo}, double @var{alpha}, const gsl_vector_complex * @var{x}, gsl_matrix_complex * @var{A})
 423 @cindex HER, Level-2 BLAS
 424 These functions compute the hermitian rank-1 update @math{A = \alpha x
 425 x^H + A} of the hermitian matrix @var{A}.  Since the matrix @var{A} is
 426 hermitian only its upper half or lower half need to be stored.  When
 427 @var{Uplo} is @code{CblasUpper} then the upper triangle and diagonal of
 428 @var{A} are used, and when @var{Uplo} is @code{CblasLower} then the
 429 lower triangle and diagonal of @var{A} are used.  The imaginary elements
 430 of the diagonal are automatically set to zero.
 431 @end deftypefun
 432
 433 @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})
 434 @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})
 435 @cindex SYR2, Level-2 BLAS
 436 These functions compute the symmetric rank-2 update @math{A = \alpha x
 437 y^T + \alpha y x^T + A} of the symmetric matrix @var{A}.  Since the
 438 matrix @var{A} is symmetric only its upper half or lower half need to be
 439 stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
 440 and diagonal of @var{A} are used, and when @var{Uplo} is
 441 @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 442 used.
 443 @end deftypefun
 444
 445 @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})
 446 @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})
 447 @cindex HER2, Level-2 BLAS
 448 These functions compute the hermitian rank-2 update @math{A = \alpha x
 449 y^H + \alpha^* y x^H A} of the hermitian matrix @var{A}.  Since the
 450 matrix @var{A} is hermitian only its upper half or lower half need to be
 451 stored.  When @var{Uplo} is @code{CblasUpper} then the upper triangle
 452 and diagonal of @var{A} are used, and when @var{Uplo} is
 453 @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 454 used.  The imaginary elements of the diagonal are automatically set to zero.
 455 @end deftypefun
 456
 457 @node Level 3 GSL BLAS Interface
 458 @subsection Level 3
 459
 460
 461 @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})
 462 @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})
 463 @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})
 464 @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})
 465 @cindex GEMM, Level-3 BLAS
 466 These functions compute the matrix-matrix product and sum @math{C =
 467 \alpha op(A) op(B) + \beta C} where @math{op(A) = A}, @math{A^T},
 468 @math{A^H} for @var{TransA} = @code{CblasNoTrans}, @code{CblasTrans},
 469 @code{CblasConjTrans} and similarly for the parameter @var{TransB}.
 470 @end deftypefun
 471
 472
 473 @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})
 474 @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})
 475 @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})
 476 @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})
 477 @cindex SYMM, Level-3 BLAS
 478 These functions compute the matrix-matrix product and sum @math{C =
 479 \alpha A B + \beta C} for @var{Side} is @code{CblasLeft} and @math{C =
 480 \alpha B A + \beta C} for @var{Side} is @code{CblasRight}, where the
 481 matrix @var{A} is symmetric.  When @var{Uplo} is @code{CblasUpper} then
 482 the upper triangle and diagonal of @var{A} are used, and when @var{Uplo}
 483 is @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 484 used.
 485 @end deftypefun
 486
 487 @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})
 488 @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})
 489 @cindex HEMM, Level-3 BLAS
 490 These functions compute the matrix-matrix product and sum @math{C =
 491 \alpha A B + \beta C} for @var{Side} is @code{CblasLeft} and @math{C =
 492 \alpha B A + \beta C} for @var{Side} is @code{CblasRight}, where the
 493 matrix @var{A} is hermitian.  When @var{Uplo} is @code{CblasUpper} then
 494 the upper triangle and diagonal of @var{A} are used, and when @var{Uplo}
 495 is @code{CblasLower} then the lower triangle and diagonal of @var{A} are
 496 used.  The imaginary elements of the diagonal are automatically set to
 497 zero.
 498 @end deftypefun
 499
 500 @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})
 501 @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})
 502 @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})
 503 @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})
 504 @cindex TRMM, Level-3 BLAS
 505 These functions compute the matrix-matrix product @math{B = \alpha op(A)
 506 B} for @var{Side} is @code{CblasLeft} and @math{B = \alpha B op(A)} for
 507 @var{Side} is @code{CblasRight}.  The matrix @var{A} is triangular and
 508 @math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
 509 @code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}. When
 510 @var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
 511 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 512 of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
 513 diagonal of @var{A} is used, but if @var{Diag} is @code{CblasUnit} then
 514 the diagonal elements of the matrix @var{A} are taken as unity and are
 515 not referenced.
 516 @end deftypefun
 517
 518
 519 @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})
 520 @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})
 521 @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})
 522 @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})
 523 @cindex TRSM, Level-3 BLAS
 524 These functions compute the inverse-matrix matrix product
 525 @math{B = \alpha op(inv(A))B} for @var{Side} is
 526 @code{CblasLeft} and @math{B = \alpha B op(inv(A))} for
 527 @var{Side} is @code{CblasRight}.  The matrix @var{A} is triangular and
 528 @math{op(A) = A}, @math{A^T}, @math{A^H} for @var{TransA} =
 529 @code{CblasNoTrans}, @code{CblasTrans}, @code{CblasConjTrans}. When
 530 @var{Uplo} is @code{CblasUpper} then the upper triangle of @var{A} is
 531 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 532 of @var{A} is used.  If @var{Diag} is @code{CblasNonUnit} then the
 533 diagonal of @var{A} is used, but if @var{Diag} is @code{CblasUnit} then
 534 the diagonal elements of the matrix @var{A} are taken as unity and are
 535 not referenced.
 536 @end deftypefun
 537
 538 @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})
 539 @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})
 540 @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})
 541 @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})
 542 @cindex SYRK, Level-3 BLAS
 543 These functions compute a rank-k update of the symmetric matrix @var{C},
 544 @math{C = \alpha A A^T + \beta C} when @var{Trans} is
 545 @code{CblasNoTrans} and @math{C = \alpha A^T A + \beta C} when
 546 @var{Trans} is @code{CblasTrans}.  Since the matrix @var{C} is symmetric
 547 only its upper half or lower half need to be stored.  When @var{Uplo} is
 548 @code{CblasUpper} then the upper triangle and diagonal of @var{C} are
 549 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 550 and diagonal of @var{C} are used.
 551 @end deftypefun
 552
 553 @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})
 554 @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})
 555 @cindex HERK, Level-3 BLAS
 556 These functions compute a rank-k update of the hermitian matrix @var{C},
 557 @math{C = \alpha A A^H + \beta C} when @var{Trans} is
 558 @code{CblasNoTrans} and @math{C = \alpha A^H A + \beta C} when
 559 @var{Trans} is @code{CblasTrans}.  Since the matrix @var{C} is hermitian
 560 only its upper half or lower half need to be stored.  When @var{Uplo} is
 561 @code{CblasUpper} then the upper triangle and diagonal of @var{C} are
 562 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 563 and diagonal of @var{C} are used.  The imaginary elements of the
 564 diagonal are automatically set to zero.
 565 @end deftypefun
 566
 567 @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})
 568 @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})
 569 @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})
 570 @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})
 571 @cindex SYR2K, Level-3 BLAS
 572 These functions compute a rank-2k update of the symmetric matrix @var{C},
 573 @math{C = \alpha A B^T + \alpha B A^T + \beta C} when @var{Trans} is
 574 @code{CblasNoTrans} and @math{C = \alpha A^T B + \alpha B^T A + \beta C} when
 575 @var{Trans} is @code{CblasTrans}.  Since the matrix @var{C} is symmetric
 576 only its upper half or lower half need to be stored.  When @var{Uplo} is
 577 @code{CblasUpper} then the upper triangle and diagonal of @var{C} are
 578 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 579 and diagonal of @var{C} are used.
 580 @end deftypefun
 581
 582 @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})
 583 @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})
 584 @cindex HER2K, Level-3 BLAS
 585 These functions compute a rank-2k update of the hermitian matrix @var{C},
 586 @math{C = \alpha A B^H + \alpha^* B A^H + \beta C} when @var{Trans} is
 587 @code{CblasNoTrans} and @math{C = \alpha A^H B + \alpha^* B^H A + \beta C} when
 588 @var{Trans} is @code{CblasConjTrans}.  Since the matrix @var{C} is hermitian
 589 only its upper half or lower half need to be stored.  When @var{Uplo} is
 590 @code{CblasUpper} then the upper triangle and diagonal of @var{C} are
 591 used, and when @var{Uplo} is @code{CblasLower} then the lower triangle
 592 and diagonal of @var{C} are used.  The imaginary elements of the
 593 diagonal are automatically set to zero.
 594 @end deftypefun
 595
 596 @node BLAS Examples
 597 @section Examples
 598
 599 The following program computes the product of two matrices using the
 600 Level-3 @sc{blas} function @sc{dgemm},
 601 @tex
 602 \beforedisplay
 603 $$
 604 \left(
 605 \matrix{0.11&0.12&0.13\cr
 606 0.21&0.22&0.23\cr}
 607 \right)
 608 \left(
 609 \matrix{1011&1012\cr
 610 1021&1022\cr
 611 1031&1031\cr}
 612 \right)
 613 =
 614 \left(
 615 \matrix{367.76&368.12\cr
 616 674.06&674.72\cr}
 617 \right)
 618 $$
 619 \afterdisplay
 620 @end tex
 621 @ifinfo
 622
 623 @example
 624 [ 0.11 0.12 0.13 ]  [ 1011 1012 ]     [ 367.76 368.12 ]
 625 [ 0.21 0.22 0.23 ]  [ 1021 1022 ]  =  [ 674.06 674.72 ]
 626                     [ 1031 1032 ]
 627 @end example
 628
 629 @end ifinfo
 630 @noindent
 631 The matrices are stored in row major order, according to the C convention
 632 for arrays.
 633
 634 @example
 635 @verbatiminclude examples/blas.c
 636 @end example
 637
 638 @noindent
 639 Here is the output from the program,
 640
 641 @example
 642 $ ./a.out
 643 @verbatiminclude examples/blas.out
 644 @end example
 645
 646 @node BLAS References and Further Reading
 647 @section References and Further Reading
 648
 649 Information on the @sc{blas} standards, including both the legacy and
 650 draft interface standards, is available online from the @sc{blas}
 651 Homepage and @sc{blas} Technical Forum web-site.
 652
 653 @itemize @asis
 654 @item
 655 @cite{BLAS Homepage} @*
 656 @uref{http://www.netlib.org/blas/}
 657 @item
 658 @cite{BLAS Technical Forum} @*
 659 @uref{http://www.netlib.org/cgi-bin/checkout/blast/blast.pl}
 660 @end itemize
 661
 662 @noindent
 663 The following papers contain the specifications for Level 1, Level 2 and
 664 Level 3 @sc{blas}.
 665
 666 @itemize @asis
 667 @item
 668 C. Lawson, R. Hanson, D. Kincaid, F. Krogh, ``Basic Linear Algebra
 669 Subprograms for Fortran Usage'', @cite{ACM Transactions on Mathematical
 670 Software}, Vol.@: 5 (1979), Pages 308--325.
 671
 672 @item
 673 J.J. Dongarra, J. DuCroz, S. Hammarling, R. Hanson, ``An Extended Set of
 674 Fortran Basic Linear Algebra Subprograms'', @cite{ACM Transactions on
 675 Mathematical Software}, Vol.@: 14, No.@: 1 (1988), Pages 1--32.
 676
 677 @item
 678 J.J. Dongarra, I. Duff, J. DuCroz, S. Hammarling, ``A Set of
 679 Level 3 Basic Linear Algebra Subprograms'', @cite{ACM Transactions on
 680 Mathematical Software}, Vol.@: 16 (1990), Pages 1--28.
 681 @end itemize
 682
 683 @noindent
 684 Postscript versions of the latter two papers are available from
 685 @uref{http://www.netlib.org/blas/}. A @sc{cblas} wrapper for Fortran @sc{blas}
 686 libraries is available from the same location.