BLAS¶
Usage
use BLAS;
or
import BLAS;
Submodules
Support for Basic Linear Algebra Subprograms (BLAS) kernel routines.
BLAS (Basic Linear Algebra Subprograms) are the de facto standard routines for low-level matrix, vector, and scalar operations. While netlib provides the official reference version of BLAS, this documentation refers to the MKL BLAS documentation, due to interface similarities.
This module is intended to work with non-distributed dense rectangular
(DefaultRectangular
) arrays.
Compiling with BLAS¶
In order to compile a Chapel program that uses this module, the
BLAS and C_BLAS (C wrappers to BLAS) libraries must be installed on the system.
The paths to both the cblas.h
header file and BLAS library
must be passed to the -I
and -L
compiler arguments. The library name,
typically blas
, must be passed to the -l
argument as well.
The compilation command should look something like this:
chpl -I$PATH_TO_CBLAS_DIR \
-L$PATH_TO_BLAS_LIBS -lblas source.chpl
- BLAS Implementations:
There is a wide range of BLAS implementations available. This module was built and tested with netlib’s C_BLAS, but many other implementations are compatible as well. Using a version of BLAS optimized for the user’s system will yield the best performance.
There are a few known caveats with some popular implementations:
ATLAS
Compilation (linking) requires the additional flag:
-latlas
MKL
Compilation requires the additional flag in order to use the MKL header instead:
--set blasImpl=mkl
OpenBLAS
Compilation will generate warnings about incompatible pointer types, which may be ignored. These warnings are due to the header files of OpenBLAS differing from the reference C_BLAS prototypes for complex arguments by using
float*
anddouble*
pointers, instead ofvoid*
pointers.
Cray LibSci
On Cray systems with the
cray-libsci
module loaded, no compiler flags should be necessary when compiling programs that use BLAS. This is typically loaded by default, but can be manually loaded withmodule load cray-libsci
as well. Chapel programs compiled on Crays utilize thecc
wrapper as the backend compiler, which implicitly links against the libsci library. Therefore, no additional steps are required of the user.
Chapel BLAS API¶
This module provides higher-level wrappers around the BLAS functions. These provide reasonable default values for many less commonly used arguments and determine the appropriate functions to call, based on the array element types, as well as the array dimensions. The other functionality is identical to the corresponding BLAS functions.
The names of these routines are identical to the corresponding BLAS functions,
except that the type prefix is dropped. For instance, gemm
is the
wrapper for the [sdcz]gemm
routines.
The native BLAS interface can still be accessed by calling routines from the
C_BLAS
submodule.
The ldA
argument is omitted from the Chapel BLAS API. Chapel determines the
dimensions of the matrices from the arrays that are passed in, even when one is
passing in a sub-array such that the array elements are not contiguously stored
in memory.
Warning
The CHPL_LOCALE_MODEL=numa
configuration is currently not supported by
this module.
- enum BlasImpl { blas, mkl, off }¶
Available BLAS implementations for
blasImpl
- config param blasImpl = BlasImpl.blas¶
Specifies which header filename to include, based on the BLAS implementation.
Most BLAS implementations rely on
cblas.h
, which is used whenblasImpl = blas
, the default setting.blas
includescblas.h
(default)mkl
includesmkl_cblas.h
off
includes nothing
- config param blasHeader = ""¶
Manually specifies the header filename to include. This flag overrides the header determined by
blasImpl
.This flag should only be necessary if using an
BLAS
implementation with a unique header name that is not supported byblasImpl
. However, no guarantees can be made about this module working with untested implementations.
- proc isBLASType(type t) param: bool¶
Return true if type is supported by BLAS
- enum Order { Row = 101: c_int, Col }¶
Define row or column order
- enum Op { N = 111: c_int, T, H }¶
Operation of matrix : none, transpose, or adjoint
- enum Uplo { Upper = 121: c_int, Lower }¶
Storage for symmetric matrices
- enum Diag { NonUnit = 131: c_int, Unit }¶
Assume a unit or non-unit diagonal
- enum Side { Left = 141: c_int, Right }¶
Operate on the left or right side
- proc gemm(A: [?Adom], B: [?Bdom], C: [?Cdom], alpha, beta, opA: Op = Op.N, opB: Op = Op.N, order: Order = Order.Row)¶
Wrapper for the GEMM routines:
C := alpha * op(A) * op(B) + beta * C
- proc symm(A: [?Adom], B: [?Bdom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, side: Side = Side.Left, order: Order = Order.Row)¶
Wrapper for the SYMM routines:
C := alpha * A * B + beta * C
or:
C := alpha * B * A + beta * C
where
A
is a symmetric matrix.
- proc hemm(A: [?Adom], B: [?Bdom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, side: Side = Side.Left, order: Order = Order.Row)¶
Wrapper for the HEMM routines:
C := alpha * A * B + beta * C
or:
C := alpha * B * A + beta * C
where
A
is an hermitian matrix.
- proc syrk(A: [?Adom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, order: Order = Order.Row)¶
Wrapper for the SYRK routines:
C := alpha * A * A**T + beta * C
or:
C := alpha * A**T * A + beta * C
where
C
is a symmetric matrix, andA**T
is the transpose ofA
.
- proc herk(A: [?Adom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, order: Order = Order.Row)¶
Wrapper for the HERK routines:
C := alpha * A * A**H + beta * C
or:
C := alpha * A**H * A + beta * C
where
C
is an hermitian matrix, andA**H
is the conjugate transpose ofA
.
- proc syr2k(A: [?Adom], B: [?Bdom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, order: Order = Order.Row)¶
Wrapper for the SYR2K routines:
C := alpha * A * B**T + alpha * B * A**T + beta * C
or:
C := alpha * A**T * B + alpha * B**T * A + beta * C
where
C
is a symmetric matrix, andA**T
andB**T
are the transpose ofA
andB
, respectively.
- proc her2k(A: [?Adom], B: [?Bdom], C: [?Cdom], alpha, beta, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, order: Order = Order.Row)¶
Wrapper for the HER2K routines:
C := alpha * A * B**H + conjg(alpha) * B * A**H + beta * C
or:
C := alpha * A**H * B + conjg(alpha) * B**H * A + beta * C
where
C
is an hermitian matrix,B**H
is the conjugate transpose ofB
,conjg(alpha)
is the complex conjugate ofalpha
.
- proc trmm(A: [?Adom] ?eltType, B: [?Bdom] eltType, alpha, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, side: Side = Side.Left, diag: Diag = Diag.NonUnit, order: Order = Order.Row) throws¶
Wrapper for the TRMM routines:
B := alpha * op(A) * B
or:
B := alpha * B * op(A)
where
A
is a triangular matrix.- Throws
IllegalArgumentError – When B is a non-square array.
- proc trsm(A: [?Adom], B: [?Bdom], alpha, uplo: Uplo = Uplo.Upper, trans: Op = Op.N, side: Side = Side.Left, diag: Diag = Diag.NonUnit, order: Order = Order.Row) throws¶
Wrapper for the TRSM routines:
op(A) * X = alpha * B
or:
X * op(A) = alpha * B
where
A
is a triangular matrix.- Throws
IllegalArgumentError – When B is a non-square array.
- proc gemv(A: [?Adom] ?eltType, x: [?xdom] eltType, y: [?ydom] eltType, alpha, beta, opA: Op = Op.N, order: Order = Order.Row, incx: c_int = 1, incy: c_int = 1)¶
Wrapper for the GEMV routines:
y := alpha*op(A)*x + beta*y,
where
A
is anm``x``n
matrix.
- proc ger(A: [?Adom] ?eltType, X: [?Xdom] eltType, Y: [?Ydom] eltType, alpha, order: Order = Order.Row, incx: c_int = 1, incy: c_int = 1)¶
Wrapper for GER routines:
A := alpha*x*y'+ A
- proc gerc(A: [?Adom] ?eltType, X: [?Xdom] eltType, Y: [?Ydom] eltType, ref alpha: eltType, order: Order = Order.Row, incx: c_int = 1, incy: c_int = 1)¶
Wrapper for GERC routines:
A := alpha*x*conjg(y') + A
- proc geru(A: [?Adom] ?eltType, X: [?Xdom] eltType, Y: [?Ydom] eltType, ref alpha: eltType, order: Order = Order.Row, incx: c_int = 1, incy: c_int = 1)¶
Wrapper for the GERU routines:
A := alpha*x*y' + A
- proc hemv(A: [?Adom] ?eltType, X: [?vDom] eltType, Y: [vDom] eltType, ref alpha: eltType, ref beta: eltType, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1, incy: c_int = 1) throws¶
Wrapper for the HEMV routines:
y := alpha*A*x + beta*y
- Throws
IllegalArgumentError – When A is a non-square array.
- proc her(A: [?Adom] ?eltType, X: [?vDom] eltType, alpha, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1) throws¶
Wrapper for the HER routines:
A := alpha*x*conjg(x') + A
- Throws
IllegalArgumentError – When A is a non-square array.
- proc her2(A: [?Adom] ?eltType, X: [?vDom] eltType, Y: [vDom] eltType, ref alpha: eltType, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1, incy: c_int = 1) throws¶
Wrapper for HER2 routines:
A := alpha *x*conjg(y') + conjg(alpha)*y *conjg(x') + A
- Throws
IllegalArgumentError – When A is a non-square array.
- proc symv(A: [?Adom] ?eltType, X: [?vDom] eltType, Y: [vDom] eltType, alpha, beta, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1, incy: c_int = 1) throws¶
Wrapper for the SYMV routines:
y := alpha*A*x + beta*y
- Throws
IllegalArgumentError – When A is a non-square array.
- proc syr(A: [?Adom] ?eltType, X: [?vDom] eltType, alpha, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1) throws¶
Wrapper for SYR routines:
A := alpha*x*x' + A
- Throws
IllegalArgumentError – When A is a non-square array.
- proc syr2(A: [?Adom] ?eltType, X: [?vDom] eltType, Y: [vDom] eltType, alpha, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, incx: c_int = 1, incy: c_int = 1) throws¶
Wrapper for SYR2 routines:
A := alpha*x*y'+ alpha*y*x' + A
- Throws
IllegalArgumentError – When A is a non-square array.
- proc trmv(A: [?Adom] ?eltType, X: [?vDom] eltType, trans: Op = Op.N, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, diag: Diag = Diag.NonUnit, incx: c_int = 1) throws¶
Wrapper for TRMV routines:
x := op(A)*x
- Throws
IllegalArgumentError – When A is a non-square array.
- proc trsv(A: [?Adom] ?eltType, B: [?vDom] eltType, trans: Op = Op.N, order: Order = Order.Row, uplo: Uplo = Uplo.Upper, diag: Diag = Diag.NonUnit, incx: c_int = 1) throws¶
Wrapper for the TRSV routines:
A*op(x) = b
- Throws
IllegalArgumentError – When A is a non-square array.
- proc rotg(ref a: ?eltType, ref b: eltType, ref c: eltType, ref s: eltType)¶
Wrapper for the ROTG routines
Construct Givens plane rotation of point
p
defined by Cartesian coordinates(a, b)
:| c s | |a| |r| | | . | | = | | |-s c | |b| |0|
Input:
a
: x-coordinate of pointp
b
: y-coordinate of pointp
Output:
a
: stores length vector (r
) of inputs(a, b)
b
: storesz
parameter that is defined belowc
: stores value ofc
element defined aboves
: stores value ofs
element defined above
The
z
parameter (stored inb
) is defined such that:if |a| > |b| then z = s; else if c != 0 then z = 1/c; else z = 1
- proc rotmg(ref d1: ?eltType, ref d2: eltType, ref b1: eltType, b2: eltType, P: [] eltType) throws¶
Wrapper for the ROTMG routines
Generate Givens rotation of points:
|b1| |b1*sqrt(d1)| | | = H.| | |0 | |b2*sqrt(d2)|
Input:
d1
: Scaling factor forb1
(x-axis)d2
: Scaling factor forb2
(y-axis)b1
: x-coordinate of input vectorb2
: y-coordinate of input vector
Output:
d1
: Provides the first element the diagonal matrixd2
: Provides the second element the diagonal matrixb1
: Provides theb1
rotated (rotated x coordination) of the vectorP
: Parameter array of 5 elements, detailed below:|P[1] P[3]| if P[0] == -1 then H = | | |P[2] P[4]| |1.0 P[3]| if P[0] == 0 then H = | | |P[2] 1.0| |P[1] 1.0| if P[0] == 1 then H = | | |-1.0 P[4]| |1 0| if P[0] == -2 then H = | | |0 1|
- Throws
IllegalArgumentError – When P does not consist of exactly five elements.
- proc rot(X: [?D] ?eltType, Y: [D] eltType, c: eltType, s: eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for the ROT routines
Replaces the value elements of two vectors
X
andY
using the equation:X[i] = c*X[i] + s*X[i] Y[i] = c*Y[i] - s*X[i]
Input:
X
: Input vectorY
: Input vectorc
: Scalar parameters
: Scalar parameterincY
: Defines the increment for the vectorY
incX
: Defines the increment for the vectorX
Output:
X
: Vector with updated elementsY
: Vector with updated elements
- proc rotm(X: [?D] ?eltType, Y: [D] eltType, P: [] eltType, incY: c_int = 1, incX: c_int = 1) throws¶
Wrapper for the ROTM routines
Executes the modified Givens rotations with element wise substitution:
|X[i]| |X[i]| | | = H.| | |Y[i]| |Y[i]|
Input:
X
: Input vectorY
: Input vectorincY
: Defines the increment for the vectorY
incX
: Defines the increment for the vectorX
P
: Parameter array of 5 elements, detailed below:|P[1] P[3]| if P[0] == -1 then H = | | |P[2] P[4]| |1.0 P[3]| if P[0] == 0 then H = | | |P[2] 1.0| |P[1] 1.0| if P[0] == 1 then H = | | |-1.0 P[4]| |1 0| if P[0] == -2 then H = | | |0 1|
Output:
X
: Vector with updated elementsY
: Vector with updated elements
- Throws
IllegalArgumentError – When P does not consist of exactly five elements.
- proc swap(X: [?D] ?eltType, Y: [D] eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for the SWAP routines
Exchanges elements of two vectors.
Input:
X
: Input vectorY
: Input vectorincY
: Defines the increment for the vectorY
incX
: Defines the increment for the vectorX
Output:
X
: Contains inputY
elementsY
: Contains inputX
elements
- proc scal(X: [?D] ?eltType, ref alpha: eltType, incX: c_int = 1)¶
Wrapper for the SCAL routines
Calculates the product of a vector
X
with scalar alpha:X' = alpha*X
Input:
X
: Input vectoralpha
: Scalar parameter
Output:
X
: Vector updated by the equation:X[i] = alpha*X[i]
- proc copy(X: [?D] ?eltType, Y: [D] eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for the COPY routines
Copies one vector (
X
the source) to another (Y
the destination):Y = X
Input:
X
: Input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
Output:
Y
: Contains the values copied fromX
vector
- proc axpy(X: [?D] ?eltType, Y: [D] eltType, ref alpha: eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for the AXPY routines
Computes the vector-scalar product of apha and
X
and adds the result toY
:Y = alpha*X + Y
Input:
X
: Input vectoralpha
: Scalar parameterincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
Output:
Y
: Vector updated by the equation:Y[i] = alpha*X[i]+Y[i]
- proc dot(X: [?xD] ?eltType, Y: [?yD] eltType, incY: c_int = 1, incX: c_int = 1): eltType¶
Wrapper for DOT routines
Returns the dot product of two vectors:
X*Y
Input:
X
: Input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
- Returns
Scalar value of dot product
- proc dotu(X: [?D] ?eltType, Y: [D] eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for DOTU routines (
DOTU_SUB
)Obtains the dot product of two complex vectors:
X*Y
Input:
X
: Input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
- Returns
The complex dot product
- proc dotc(X: [?D] ?eltType, Y: [D] eltType, incY: c_int = 1, incX: c_int = 1)¶
Wrapper for DOTC routines (
DOTC_SUB
)Obtains the dot product of conjugated X vector with Y vector:
conj(X)*Y
Input:
X
: Conjugated input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
- Returns
The complex dot product
- proc dsdot(X: [?D] real(32), Y: [D] real(32), incY: c_int = 1, incX: c_int = 1): real(64)¶
Wrapper for SDOT routines (
DSDOT
variant)Returns the dot product of two
real(32)
vectors as areal(64)
, usingreal(64)
precision internally:X*Y
Input:
X
: Input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
- Returns
Scalar value of dot product
- proc sdsdot(X: [?D] real(32), Y: [D] real(32), incY: c_int = 1, incX: c_int = 1): real(32)¶
Wrapper for SDOT routines (
SDSDOT
variant)Returns the dot product of two
real(32)
vectors as areal(32)
, usingreal(64)
precision internally:X*Y
Input:
X
: Input vectorY
: Input vectorincX
: Defines the increment for the vectorX
incY
: Defines the increment for the vectorY
- Returns
Scalar value of dot product
- proc nrm2(X: [?D] ?eltType, incX: c_int = 1)¶
Wrapper for NRM2 routines
Returns the Euclidean norm of vector
X
:||X||
Input:
X
: Input vectorincX
: Defines the increment for the vectorX
- Returns
The 2-norm of
X
vector
- proc asum(X: [?D] ?eltType, incX: c_int = 1)¶
Wrapper for the ASUM routines
Returns the sum of the magnitude values of X elements:
|Re X[1]| + |Im X[1]| + |Re X[2]| + |Im X[2]|+ ... + |Re X[N]| + |Im X[N]|.
Input:
X
: Input vectorincX
: Defines the increment for the vectorX
- Returns
The 1-norm of
X
vector