## blas

An interface to level 1, 2 and 3 BLAS linear algebra routines.

## TOC »

- blas
- Usage
- Documentation
- Naming conventions for routines
- Vector copy routines
- BLAS level 1 routines
- Conventions
- Apply plane rotation
- Scale vector
- Swap the elements of two vectors
- Real vector dot product
- Complex vector dot product
- Hermitian vector dot product
- Vector multiply-add
- Vector multiply-add with optional offset
- Euclidean norm of a vector
- Sum of absolute values of the elements in a vector
- Sum of absolute values of the elements in a vector

- BLAS level 2 routines
- Conventions
- General matrix-vector multiply-add
- Banded matrix-vector multiply-add
- Hermitian matrix-vector multiply-add
- Hermitian banded matrix-vector multiply-add
- Symmetric matrix-vector multiply-add
- Banded symmetric matrix-vector multiply-add
- Triangular matrix-vector multiply-add
- Banded triangular matrix-vector multiply-add
- Triangular matrix equation solve
- Banded triangular matrix equation solve
- Rank 1 operation
- Rank 1 operation with optional offset
- Rank 1 operation on complex matrices and vectors
- Rank 1 operation on complex matrices and vectors
- Hermitian rank 1 operation
- Hermitian rank 2 operation
- Symmetric rank 1 operation
- Symmetric rank 2 operation

- BLAS level 3 routines

- Examples
- About this egg

## Usage

(require-extension blas)

## Documentation

### Naming conventions for routines

Every routine in the BLAS library comes in four flavors, each prefixed by the letters S, D, C, and Z, respectively. Each letter indicates the format of input data:

- S stands for single-precision (32-bit IEEE floating point numbers),
- D stands for double-precision (64-bit IEEE floating point numbers),
- C stands for complex numbers (represented by pairs of 32-bit IEEE floating point numbers),
- Z stands for double complex numbers (represented by pairs of 64-bit IEEE floating point numbers)

In addition, each BLAS routine in this egg comes in three flavors:

- Safe, pure

*Safe* routines check the sizes of their input arguments. For example, if a routine is supplied arguments that indicate that an input matrix is of dimensions *M*-by-*N*, then the argument corresponding to that matrix is checked that it is of size *M * N*.

*Pure* routines do not alter their arguments in any way. A new matrix or vector is allocated for the return value of the routine.

- Safe, destructive (suffix: !)

*Safe* routines check the sizes of their input arguments. For example, if a routine is supplied arguments that indicate that an input matrix is of dimensions *M*-by-*N*, then the argument corresponding to that matrix is checked that it is of size *M * N*.

*Destructive* routines can modify some or all of their arguments. They are given names ending in exclamation mark. Please consult the BLAS documentation to determine which functions modify their input arguments.

- Unsafe, destructive (prefix:
**unsafe-**, suffix: !)

*Unsafe* routines do not check the sizes of their input arguments. They invoke the corresponding BLAS routines directly. Unsafe routines do not have pure variants.

For example, function *xGEMM* (matrix-matrix multiplication) comes in the following variants:

BLAS name | Safe, pure | Safe, destructive | Unsafe, destructive |
---|---|---|---|

SGEMM | sgemm | sgemm! | unsafe-sgemm! |

DGEMM | dgemm | dgemm! | unsafe-dgemm! |

CGEMM | cgemm | cgemm! | unsafe-cgemm! |

ZGEMM | zgemm | zgemm! | unsafe-zgemm! |

### Vector copy routines

`scopy:`procedure

`dcopy:`procedure

`ccopy:`procedure

`zcopy:`procedureThese procedures return a copy of given input SRFI-4 vector. The returned vector is allocated with the corresponding SRFI-4 constructor, and the input vector is copied to it by the corresponding BLAS copy procedure.

### BLAS level 1 routines

#### Conventions

The BLAS level 1 procedures in this egg differ from the actual routines they invoke by the position of the vector increment arguments (`INCX` and `INCY`). In this egg, these arguments are optional; the default value of `INCX` and `INCY` is 1.

In the procedure signatures below, these optional arguments are indicated by [ and ] (square brackets).

#### Apply plane rotation

`srot:`procedure

`drot:`procedure`xROT`applies a plane rotation matrix to a sequence of ordered pairs:`(x_i , y_i)`, for`i = 1, 2, ..., n`.`X`and`Y`are vector of dimensions`(N-1) * abs(incx) + 1`and`(N-1) * abs(incy) + 1`, respectively.`C`and`S`are respectively the cosine and sine of the plane of rotation.

#### Scale vector

`sscal:`procedure

`dscal:`procedure

`cscal:`procedure

`zscal:`procedure`xSCAL`scales a vector with a scalar:`x := alpha * x`.

#### Swap the elements of two vectors

`sswap:`procedure

`dswap:`procedure

`cswap:`procedure

`zswap:`procedure`xSWAP`interchanges the elements of two vectors:`x <-> y`.

#### Real vector dot product

`sdot:`procedure

`ddot:`procedure`xDOT`computes the dot product of two vectors of real values:`dot := x'*y = \Sum_{i=1}^{n} (x_i * y_i)`.

#### Complex vector dot product

`cdotu:`procedure

`zdotu:`procedure`xDOTU`computes the dot product of two vectors of complex values:`dotu := x'*y = \Sum_{i=1}^{n} (x_i * y_i)`.

#### Hermitian vector dot product

`cdotc:`procedure

`zdotc:`procedure`xDOTC`computes the dot product of the conjugates of two complex vectors:`dotu := conjg(x')*y = \Sum_{i=1}^{n} (conjg(x_i) * y_i)`, for`i = 1, 2, ..., n`.

#### Vector multiply-add

`saxpy:`procedure

`daxpy:`procedure

`caxpy:`procedure

`zaxpy:`procedure`xAXPY`adds a scalar multiple of a vector to another vector:`y := alpha * x + y`.

#### Vector multiply-add with optional offset

`siaxpy:`procedure

`diaxpy:`procedure

`ciaxpy:`procedure

`ziaxpy:`procedure`xIAXPY`adds a scalar multiple of a vector to another vector, where the beginning of each vector argument can be offset:`y[yofs:n] := alpha * x[xofs:n] + y[yofs:n]`.

#### Euclidean norm of a vector

`snrm2:`procedure

`dnrm2:`procedure

`cnrm2:`procedure

`znrm2:`procedure`xNRM2`computes the Euclidean (L2) norm of a vector.

#### Sum of absolute values of the elements in a vector

`sasum:`procedure

`dasum:`procedure

`casum:`procedure

`zasum:`procedure`xASUM`sums the absolute values of the elements in a vector.

#### Sum of absolute values of the elements in a vector

`samax:`procedure

`damax:`procedure

`camax:`procedure

`zamax:`procedure`xAMAX`searches a vector for the first occurrence of its maximum absolute value, and returns the index of that element.

### BLAS level 2 routines

#### Conventions

The BLAS level 2 procedures in this egg differ from the actual routines they invoke by the position of the leading dimension argument (`LDA`) and the vector increment arguments (`INCX` and `INCY`). In this egg, these arguments are optional; the default value of `LDA`is the largest matrix dimension, depending on the semantics of the respective operation, and the default value of `INCX` and `INCY` is 1.

In the procedure signatures below, these optional arguments are indicated by [ and ] (square brackets).

Argument `ORDER` is one of `RowMajor` or `ColMajor` to indicate that the input and output matrices are in row-major or column-major form, respectively.

Where present, argument `TRANS` can be one of `NoTrans` or `Trans` to indicate whether the input matrix is to be transposed or not.

Where present, argument `UPLO` can be one of `Upper` or `Lower` to indicate whether the upper or lower triangular part of an input symmetric matrix is to referenced,or to specify the type of an input triangular matrix.

Where present, argument `DIAG` can be one of `NonUnit` or `Unit` to indicate whether an input triangular matrix is unit triangular or not.

#### General matrix-vector multiply-add

`sgemv:`procedure

`dgemv:`procedure

`cgemv:`procedure

`zgemv:`procedure`xGEMV`performs the matrix-vector multiply-add operation of the form`y := alpha*op( A )*x + beta*y`, where`op( X )`is one of`op( A ) = A`or`op( A ) = A'`.`ALPHA`and`BETA`are scalars, and`A`is an`M x N`matrix.`X`is a vector of size`(1 + ( N - 1 ) * abs(INCX))`when argument`TRANS`is`NoTrans`, and`(1 + ( M - 1 ) * abs(INCX))`otherwise.`Y`is a vector of size`(1 + ( M - 1 ) * abs(INCY))`when argument`TRANS`is`NoTrans`, and`(1 + ( N - 1 ) * abs(INCY))`otherwise.

#### Banded matrix-vector multiply-add

`sgbmv:`procedure

`dgbmv:`procedure

`cgbmv:`procedure

`zgbmv:`procedure`xGBMV`performs the matrix-vector multiply-add operation of the form`y := alpha*op( A )*x + beta*y`, where`op( X )`is one of`op( A ) = A`or`op( A ) = A'`.`ALPHA`and`BETA`are scalars, and`A`is an`M x N`banded matrix, with`KL`sub-diagonals and`KU`super-diagonals.`X`is a vector of size`(1 + ( N - 1 ) * abs(INCX))`when argument`TRANS`is`NoTrans`, and`(1 + ( M - 1 ) * abs(INCX))`otherwise.`Y`is a vector of size`(1 + ( M - 1 ) * abs(INCY))`when argument`TRANS`is`NoTrans`, and`(1 + ( N - 1 ) * abs(INCY))`otherwise.

#### Hermitian matrix-vector multiply-add

`chemv:`procedure

`zhemv:`procedure`xHEMV`performs the matrix-vector multiply-add operation of the form`y := alpha*op( A )*x + beta*y`, where`op( X )`is one of`op( A ) = A`or`op( A ) = A'`.`ALPHA`and`BETA`are scalars, and`A`is an`N x N`Hermitian matrix.`X`and`Y`are`N`element vectors.

#### Hermitian banded matrix-vector multiply-add

`chbmv:`procedure

`zhbmv:`procedure`xHBMV`performs the matrix-vector multiply-add operation of the form`y := alpha*op( A )*x + beta*y`, where`op( X )`is one of`op( A ) = A`or`op( A ) = A'`.`ALPHA`and`BETA`are scalars, and`A`is an`N x N`Hermitian banded matrix, with`K`super-diagonals.`X`and`Y`are`N`element vectors.

#### Symmetric matrix-vector multiply-add

`ssymv:`procedure

`dsymv:`procedure`xSYMV`performs matrix-vector multiply-add operation of the form`y := alpha*A*x + beta*y`.`ALPHA`and`BETA`are scalars, and`A`is an`N x N`symmetric matrix.`X`and`Y`are`N`element vectors.

#### Banded symmetric matrix-vector multiply-add

`ssbmv:`procedure

`dsbmv:`procedure`xSBMV`performs matrix-vector multiply-add operation of the form`y := alpha*A*B + beta*y`.`ALPHA`and`BETA`are scalars, and`A`is an`N x N`symmetric banded matrix, with`K`super-diagonals.`X`and`Y`are`N`element vectors.

#### Triangular matrix-vector multiply-add

`strmv:`procedure

`dtrmv:`procedure

`ctrmv:`procedure

`ztrmv:`procedure`xTRMV`performs matrix-vector multiply-add operation of the form`y := alpha*op( A )*x`, where`op ( A ) = A`or`op ( A ) = A'``ALPHA`and`BETA`are scalars, and`A`is an`N x N`upper or lower triangular matrix.`X`is a vector of length`(1 + (n - 1) * abs(INCX))`.

#### Banded triangular matrix-vector multiply-add

`stbmv:`procedure

`dtbmv:`procedure

`ctbmv:`procedure

`ztbmv:`procedure`xTBMV`performs matrix-vector multiply-add operation of the form`y := alpha*A*B + beta*y`, where`op ( A ) = A`or`op ( A ) = A'``ALPHA`and`BETA`are scalars, and`A`is an`N x N`upper or lower triangular banded matrix, with`K+1`diagonals.`X`is a vector of length`(1 + (n - 1) * abs(INCX))`.

#### Triangular matrix equation solve

`strsv:`procedure

`dtrsv:`procedure

`ctrsv:`procedure

`ztrsv:`procedure`xTRSV`solves one of the systems of equations`A*x = b`or`A'*x = b`.`ALPHA`and`BETA`are scalars,`A`is a upper or lower triangular matrix, and`B`is a`N`element vector.

#### Banded triangular matrix equation solve

`stbsv:`procedure

`dtbsv:`procedure

`ctbsv:`procedure

`ztbsv:`procedure`xTBSV`solves one of the systems of equations`A*x = b`or`A'*x = b`.`ALPHA`and`BETA`are scalars,`A`is a upper or lower banded triangular matrix with`K+1`diagonals, and`B`is a`N`element vector.

#### Rank 1 operation

`sger:`procedure

`dger:`procedure`xGER`performs the rank 1 operation`A := alpha*x*y' + A`.`ALPHA`is a scalar,`X`is an`M`element vector,`Y`is an`N`element vector, and`A`is an`M x N`matrix.

#### Rank 1 operation with optional offset

`siger:`procedure

`diger:`procedure`xIGER`performs the rank 1 operation`A := alpha*x[xofs:M]*y'[yofs:N] + A`.`ALPHA`is a scalar,`X`is an`M`element vector,`Y`is an`N`element vector, and`A`is an`M x N`matrix.

#### Rank 1 operation on complex matrices and vectors

`cgeru:`procedure

`zgeru:`procedure`xGERU`performs the rank 1 operation`A := alpha*x*y' + A`.`ALPHA`is a scalar,`X`is an`M`element vector,`Y`is an`N`element vector, and`A`is an`M x N`matrix.

#### Rank 1 operation on complex matrices and vectors

`cgerc:`procedure

`zgerc:`procedure`xGERC`performs the rank 1 operation`A := alpha*x*conjg(y') + A`.`ALPHA`is a scalar,`X`is an`M`element vector,`Y`is an`N`element vector, and`A`is an`M x N`matrix.

#### Hermitian rank 1 operation

`cher:`procedure

`zher:`procedure`xHER`performs the Hermitian rank 1 operation`A := alpha*x*conjg(x') + A`.`ALPHA`is a scalar,`X`is an`N`element vector, and`A`is an`N x N`Hermitian matrix.

#### Hermitian rank 2 operation

`cher2:`procedure

`zher2:`procedure`xHER2`performs the Hermitian rank 2 operation`A := alpha*x*conjg(y') + conjg(alpha)*y*conjg(x') + A`.`ALPHA`is a scalar,`X`and`Y`are`N`element vectors, and`A`is an`N x N`Hermitian matrix.

#### Symmetric rank 1 operation

`ssyr:`procedure

`dsyr:`procedure`xSYR`performs the symmetric rank 1 operation`A := alpha*x*x' + A`.`ALPHA`is a scalar,`X`is an`N`element vector, and`A`is an`N x N`symmetric matrix.

#### Symmetric rank 2 operation

`ssyr2:`procedure

`dsyr2:`procedure`xSYR2`performs the symmetric rank 2 operation`A := alpha*x*y' + alpha*y*x' + A`.`ALPHA`is a scalar,`X`and`Y`are`N`element vectors, and`A`is an`N x N`symmetric matrix.

### BLAS level 3 routines

#### Conventions

The BLAS level 3 procedures in this egg differ from the actual routines they invoke by the position of the leading dimension arguments (`LDA`, `LDB`, and `LDC`). In this egg, these arguments are optional, and their default values are set to the largest matrix dimension, depending on the semantics of the respective operation.

In the procedure signatures below, these optional arguments are indicated by [ and ] (square brackets).

Argument `ORDER` is one of `RowMajor` or `ColMajor` to indicate that the input and output matrices are in row-major or column-major form, respectively.

Where present, arguments `TRANS`, `TRANSA`, `TRANSB` can be one of `NoTrans` or `Trans` to indicate whether the respective input matrices are to be transposed or not.

Where present, argument `SIDE` can be one of `Left` or `Right` to indicate whether an input symmetric matrix appears on the left or right in the respective operation.

Where present, argument `UPLO` can be one of `Upper` or `Lower` to indicate whether the upper or lower triangular part of an input symmetric matrix is to referenced,or to specify the type of an input triangular matrix.

Where present, argument `DIAG` can be one of `NonUnit` or `Unit` to indicate whether an input triangular matrix is unit triangular or not.

#### General matrix multiply-add

`sgemm:`procedure

`dgemm:`procedure

`cgemm:`procedure

`zgemm:`procedure`xGEMM`performs matrix-matrix multiply-add operation of the form`C := alpha*op( A )*op( B ) + beta*C`, where`op( X )`is one of`op( X ) = X`or`op( X ) = X'`.`ALPHA`and`BETA`are scalars, and`A`,`B`and`C`are matrices, with`op( A )`an`M x K`matrix,`op( B )`a`K x N`matrix and`C`an`M x N`matrix.

#### Symmetric matrix multiply-add

`ssymm:`procedure

`dsymm:`procedure

`csymm:`procedure

`zsymm:`procedure`xSYMM`performs matrix-matrix multiply-add operation of the form`C := alpha*A*B + beta*C`or`C := alpha*B*A + beta*C`.`ALPHA`and`BETA`are scalars,`A`is a symmetric matrix, and`B`and`C`are`M x N`matrices.

#### Symmetric rank k operation

`ssyrk:`procedure

`dsyrk:`procedure

`csyrk:`procedure

`zsyrk:`procedure`xSYRK`performs one of the symmetric rank k operations{{C := alpha*A*A' + beta*C}} or {{C := alpha*A'*A + beta*C}}.

`ALPHA`and`BETA`are scalars,`A`is an`N x K`or`K x N`matrix, and`C`is an`N x N`symmetric matrix.

#### Hermitian rank k operation

`cherk:`procedure

`zherk:`procedure`xHERK`performs one of the hermitian rank k operations`C := alpha*A*conjg(A') + beta*C`or`C := alpha*conjg(A')*A + beta*C`.`ALPHA`and`BETA`are scalars,`A`is an`N x K`or`K x N`matrix, and`C`is an`N x N`hermitian matrix.

#### Symmetric rank 2k operation

`ssyr2k:`procedure

`dsyr2k:`procedure

`csyr2k:`procedure

`zsyr2k:`procedure`xSYR2K`performs one of the symmetric rank 2k operations`C := alpha*A*B' + beta*C`or`C := alpha*B'*A + beta*C`.`ALPHA`and`BETA`are scalars,`A`and`B`are`N x K`or`K x N`matrices, and`C`is an`N x N`symmetric matrix.

#### Hermitian rank 2k operation

`cher2k:`procedure

`zher2k:`procedure`xHER2K`performs one of the hermitian rank 2k operations`C := alpha*A*conjg(B') + beta*C`or`C := alpha*conjg(B')*A + beta*C`.`ALPHA`and`BETA`are scalars,`A`and`B`are`N x K`or`K x N`matrices, and`C`is an`N x N`hermitian matrix.

#### Triangular matrix multiply

`strmm:`procedure

`dtrmm:`procedure

`ctrmm:`procedure

`ztrmm:`procedure`xTRMM`performs matrix-matrix multiply operation of the form`B := alpha*op( A )*B`or`B := alpha*B*op( A )`.`ALPHA`is a scalar,`A`is an upper or lower triangular matrix, and`B`is an`M x N`matrix.

#### Triangular matrix equation solve

`strsm:`procedure

`dtrsm:`procedure

`ctrsm:`procedure

`ztrsm:`procedure`xTRSM`solves one of the matrix equations{{op( A )*X = alpha*B}} or {{X*op( A ) = alpha*B}}.

`op( A )`is one of`op( A ) = A`or`op( A ) = A'`.`ALPHA`and`BETA`are scalars,`A`is a upper or lower triangular matrix, and`B`is a`M x N`matrix.

## Examples

(use srfi-4 blas) (define order ColMajor) (define transa NoTrans) (define m 4) (define n 4) (define alpha 1) (define beta 0) (define a ; column-major order! (f64vector 1 2 3 4 1 1 1 1 3 4 5 6 5 6 7 8) ) (define x (f64vector 1 2 1 1)) (define y (f64vector 0 0 0 0)) (dgemv! order transa m n alpha a x beta y) (print y)

## About this egg

### Author

Felix Winkelmann and Ivan Raikov

### Version history

- 4.1
- Fixed use of memq [thanks to Moritz Heidkamp]
- 4.0
- Using bind instead of easyffi
- 3.1
- Updated test script to return proper exit code
- 3.0
- Eliminated reduntant blas: prefix from names of exported symbols
- 2.8
- Eliminated use of noop
- 2.7
- Switched to wiki documentation
- 2.6
- Ported to Chicken 4
- 2.5
- Build script updated for better cross-platform compatibility
- 2.4
- Added iger procedures; fixed a bug in the default arguments of level 2 routines
- 2.3
- Added iaxpy procedures; fixed a bug in the default arguments of level 1 routines
- 2.2
- Fixed a bug in the renaming of C routines
- 2.1
- Added eggdoc property to meta file
- 2.0
- An overhaul of the library to introduce safe, unsafe, and pure variants of each routine
- 1.8
- Added icopy procedures [by Ivan Raikov]
- 1.7
- Support for Mac OS X added [by Will Farr]
- 1.6
- Fixed bug in blas library test code
- 1.5
- Added support for CLAPACK [by Ivan Raikov]
- 1.4
- Added support for atlas CBLAS library [by Stu Glaser]
- 1.3
- Tries to find a proper CBLAS library (currently only GSL)
- 1.2
- Adapted to new FFI macro-names
- 1.1
- Adapted to new setup scheme

### License

Copyright (c) 2003-2006, Felix L. Winkelmann Copyright (c) 2007-2014 Ivan Raikov All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. Neither the name of the author nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.