## SugarBLAS.jl

Syntactic sugar for BLAS polynomials
Popularity
29 Stars
Updated Last
2 Years Ago
Started In
August 2016

# SugarBLAS

`BLAS` functions are unaesthetic and annoying without good knowledge of the positional arguments. This package provides macros for `BLAS` functions representing polynomials. The main macro of the package is `@blas!` for most of the use cases: `copy!`, `scale!` and `axpy!`. Non mutable versions of this operator are already very easy to write so they are not included.

The macros will output a function from `BASE` module, this allows defining new behavior for custom types. Note that the output won't necessarily belong to the julia `BLAS` API, e.g. `copy!` is used instead of `BASE.LinAlg.BLAS.blascopy!` for better performance.

For now the package supports the most common BLAS functions from the internal API. The access for these functions is private since the official API is private aswell and may change in the future.

This documentation offers great examples but is by no means super extensive, for more examples check the test folder of the repository.

## Installing

To install the package, use the following command inside Julia's REPL:

`Pkg.add("SugarBLAS")`

## Usage

`@blas!` matches the expression and decides which function to call. As long as it is correctly parenthesized putting more variables won't be an issue.

```julia> macroexpand(SugarBLAS, :(@blas! Y = (a*b +c)*(X*Z) + Y))
:(Base.LinAlg.axpy!(a * b + c,X * Z,Y))

julia> macroexpand(SugarBLAS, :(@blas! X = (a+c)*X))
:(scale!(a + c,X))```

When doing this just imagine the BLAS expression.

```Y = a*X + Y
->
a := (a*b +c); X := (X*Z)
->
Y = (a*b +c)*(X*Z) + Y```

### Updating operators

Both `*=` and `+=` are supported. `*=` can only be used for scaling given that is pretty unambigous.

```julia> macroexpand(SugarBLAS, :(@blas! Y += X)) == macroexpand(SugarBLAS, :(@blas! Y = Y + X))
true```

### Commutativity

`+` is assumed as the only commutative operator, it is important to note here that `*` is not treated as commutative and therefore some expressions will lead to errors.

```julia> a = 2.3;

julia> X = rand(10,10);

julia> Y = rand(10,10);

julia> @blas! Y += X*a
ERROR: MethodError: `axpy!` has no method matching axpy!(::Array{Float64,2}, ::Float64, ::Array{Float64,2})```

The package assumes types by its position in the multiplication, this doesn't happen with addition and that's why it conserves its property.

```julia> macroexpand(SugarBLAS, :(@blas! Y = X + Y)) == macroexpand(SugarBLAS, :(@blas! Y = Y + X))
true```

## blas!

Macro for most of the use cases: `copy!`, `scale!` and `axpy!`

### scale!

Scale an array `X` by a scalar `a` overwriting `X` in-place.

Polynomials

• `X *= a`
• `X = a*X`

Example

```julia> macroexpand(SugarBLAS, :(@blas! X *= a))
:(scale!(a,X))```

### axpy!

Overwrite `Y` with `a*X + Y`. Return `Y`.

Polynomials

• `Y += X`
• `Y += a*X`

Example

```julia> macroexpand(SugarBLAS, :(@blas! Y += X))
:(Base.LinAlg.axpy!(1.0,X,Y))

julia> macroexpand(SugarBLAS, :(@blas! Y += a*X))
:(Base.LinAlg.axpy!(a,X,Y))```

### copy!

Copy all elements from collection `Y` to array `X`. Return `X`.

Polynomials

• `X = Y`

Example

```julia> macroexpand(SugarBLAS, :(@blas! X = Y))
:(copy!(X,Y))```

## Internal API

Macro for most of the functions available in the JuliaLang internal BLAS API.

### scale!

Scale an array `X` by a scalar `a` overwriting `X` in-place.

Polynomials

• `X *= a`
• `X = a*X`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@scale! X *= a))
:(scale!(a,X))```

### axpy!

Overwrite `Y` with `a*X + Y`. Return `Y`.

Polynomials

• `Y ±= X`
• `Y ±= a*X`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@axpy! Y += X))
:(Base.LinAlg.axpy!(1.0,X,Y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@axpy! Y += a*X))
:(Base.LinAlg.axpy!(a,X,Y))```

### copy!

Copy all elements from collection `Y` to array `X`. Return `X`.

Polynomials

• `X = Y`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@copy! X = Y))
:(copy!(X,Y))```

### ger!

Rank-1 update of the matrix `A` with vectors `x` and `y` as `alpha*x*y' + A`.

Polynomials

• `A ±= alpha*x*y'`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@ger! A -= alpha*x*y'))
:(Base.LinAlg.BLAS.ger!(-alpha,x,y,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@ger! A += alpha*x*y'))
:(Base.LinAlg.BLAS.ger!(alpha,x,y,A))```

### syr!

Rank-1 update of the symmetric matrix `A` with vector `x` as `alpha*x*(x)ᵀ + A`. When left side has `A['U']` the upper triangle of `A` is updated (`'L'` for lower triangle). Return `A`.

Polynomials

• `A[uplo] ±= alpha*x*(x)ᵀ`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@syr! A['U'] -= alpha*x*(x)ᵀ))
:(Base.LinAlg.BLAS.syr!('U',-alpha,x,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@syr! A['L'] += alpha*x*(x)ᵀ))
:(Base.LinAlg.BLAS.syr!('L',alpha,x,A))```

### syrk

Return either the upper triangle or the lower triangle, depending on (`'U'` or `'L'`), of `alpha*A*(A)ᵀ` or `alpha*(A)ᵀ*A`.

Polynomials

• `alpha*A*(A)ᵀ uplo=ul`
• `alpha*(A)ᵀ*A uplo=ul`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk alpha*A*(A)ᵀ uplo='U'))
:(Base.LinAlg.BLAS.syrk('U','N',alpha,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk alpha*(A)ᵀ*A uplo='L'))
:(Base.LinAlg.BLAS.syrk('L','T',alpha,A))```

### syrk!

Rank-k update of the symmetric matrix `C` as `alpha*A*(A)ᵀ + beta*C` or `alpha*(A)ᵀ*A + beta*C`. When the left hand side is`C['U']` the upper triangle of `C` is updated (`'L'` for lower triangle). Return `C`.

Polynomials

• `C[uplo] ±= alpha*A*(A)ᵀ`
• `C[uplo] = beta*C ± alpha*(A)ᵀ*A`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk! C['U'] -= alpha*A*(A)ᵀ))
:(Base.LinAlg.BLAS.syrk!('U','N',-alpha,A,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk! C['L'] = beta*C - alpha*(A)ᵀ*A))
:(Base.LinAlg.BLAS.syrk!('L','T',-alpha,A,beta,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk! C['U'] += alpha*A*(A)ᵀ))
:(Base.LinAlg.BLAS.syrk!('U','N',alpha,A,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@syrk! C['L'] = alpha*(A)ᵀ*A + beta*C))
:(Base.LinAlg.BLAS.syrk!('L','T',alpha,A,beta,C))```

### her!

Methods for complex arrays only. Rank-1 update of the Hermitian matrix `A` with vector `x` as `alpha*x*x' + A`. Whenthe left hand side is `A['U']` the upper triangle of `A` is updated (`'L'` for lower triangle). Return `A`.

Polynomials

• `A[uplo] ±= alpha*x*x'`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@her! A['U'] -= alpha*x*x'))
:(Base.LinAlg.BLAS.her!('U',-alpha,x,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@her! A['L'] = A - alpha*x*x'))
:(Base.LinAlg.BLAS.her!('L',-alpha,x,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@her! A['U'] += alpha*x*x'))
:(Base.LinAlg.BLAS.her!('U',alpha,x,A))```

### herk

Methods for complex arrays only. Returns either the upper triangle or the lower triangle, according to uplo ('U' or 'L'), of alphaAA' or alpha*A'*A, according to trans ('N' or 'T').

Polynomials

• `alpha*A*A' uplo=ul`
• `alpha*A'*A uplo=ul`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk alpha*A*A' uplo='U'))
:(Base.LinAlg.BLAS.herk('U','N',alpha,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk alpha*A'*A uplo='U'))
:(Base.LinAlg.BLAS.herk('U','T',alpha,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk alpha*A*A' uplo='L'))
:(Base.LinAlg.BLAS.herk('L','N',alpha,A))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk alpha*A'*A uplo='L'))
:(Base.LinAlg.BLAS.herk('L','T',alpha,A))```

### herk!

Methods for complex arrays only. Rank-k update of the Hermitian matrix `C` as `alpha*A*A' + beta*C` or `alpha*A'*A + beta*C`. When the left hand side is `C['U']` the upper triangle of `C` is updated (`'L'` for lower triangle). Return `C`.

Polynomials

• `C[uplo] ±= alpha*A*A'`
• `C[uplo] = beta*C ± alpha*A'*A`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk! C['L'] -= alpha*A'*A))
:(Base.LinAlg.BLAS.herk!('L','T',-alpha,A,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk! C['U'] = C - alpha*A*A'))
:(Base.LinAlg.BLAS.herk!('U','N',-alpha,A,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk! C['L'] = beta*C - alpha*A'*A))
:(Base.LinAlg.BLAS.herk!('L','T',-alpha,A,beta,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk! C['U'] += alpha*A*A'))
:(Base.LinAlg.BLAS.herk!('U','N',alpha,A,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@herk! C['L'] = alpha*A'*A + beta*C))
:(Base.LinAlg.BLAS.herk!('L','T',alpha,A,beta,C))```

### gbmv

Return `alpha*A*x` or `alpha*A'*x`. The matrix `A` is a general band matrix of dimension `m` by `size(A,2)` with `kl` sub-diagonals and `ku` super-diagonals.

Polynomials

• `alpha*A[kl:ku,h=m]*x`
• `alpha*A[h=m,kl:ku]'*x`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@gbmv alpha*A[0:ku,h=2]*x))
:(Base.LinAlg.BLAS.gbmv('N',2,0,ku,alpha,A,x))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gbmv alpha*A[h=m,-kl:ku]*x))
:(Base.LinAlg.BLAS.gbmv('N',m,kl,ku,alpha,A,x))```

### gbmv!

Update vector `y` as `alpha*A*x + beta*y` or `alpha*A'*x + beta*y`. The matrix `A` is a general band matrix of dimension `m` by `size(A,2)` with `kl` sub-diagonals and `ku` super-diagonals. Return the updated `y`.

Polynomials

• `y ±= alpha*A[kl:ku,h=m]*x`
• `y = beta*y ± alpha*A[h=m,kl:ku]'*x`

Example

```@test macroexpand(SugarBLAS, :(SugarBLAS.@gbmv! y -= alpha*A[h=m,-kl:ku]*x))
:(Base.LinAlg.BLAS.gbmv!('N',m,kl,ku,-alpha,A,x,1.0,y))

@test macroexpand(SugarBLAS, :(SugarBLAS.@gbmv! y = beta*y - alpha*A[h=2, 0:ku]'*x))
:(Base.LinAlg.BLAS.gbmv!('T',2,0,ku,-alpha,A,x,beta,y))

@test macroexpand(SugarBLAS, :(SugarBLAS.@gbmv! y = alpha*A[0:ku,h=2]*x + y))
:(Base.LinAlg.BLAS.gbmv!('N',2,0,ku,alpha,A,x,1.0,y))

@test macroexpand(SugarBLAS, :(SugarBLAS.@gbmv! y += alpha*A[h=m,-kl:ku]*x))
:(Base.LinAlg.BLAS.gbmv!('N',m,kl,ku,alpha,A,x,1.0,y))

@test macroexpand(SugarBLAS, :(SugarBLAS.@gbmv! y = alpha*A[kl:ku, h=m]'*x + beta*y))
:(Base.LinAlg.BLAS.gbmv!('T',m,-kl,ku,alpha,A,x,beta,y))```

### sbmv

Return `alpha*A*x` where `A` is a symmetric band matrix of order `size(A,2)` with `k` super-diagonals stored in the argument `A`.

Polynomials

• `A[0:k,uplo]*xv`
• `alpha*A[0:k,uplo]*x`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv A['U',0:k]*x))
:(Base.LinAlg.BLAS.sbmv('U',k,A,x))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv alpha*A[0:k,'L']*x))
:(Base.LinAlg.BLAS.sbmv('L',k,alpha,A,x))```

### sbmv!

Update vector `y` as `alpha*A*x + beta*y` where `A` is a a symmetric band matrix of order `size(A,2)` with `k` super-diagonals stored in the argument `A`. If `A[...,'U']` is used multiplication is done with `A`'s upper triangle, `L` is for the lower triangle. Return updated `y`.

Polynomials

• `y ±= alpha*A[0:k,uplo]*x`
• `y = beta*y ± alpha*A[0:k,uplo]*x`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv! y -= alpha*A['U',0:k]*x))
:(Base.LinAlg.BLAS.sbmv!('U',k,-alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv! y = beta*y - alpha*A[0:k,'U']*x))
:(Base.LinAlg.BLAS.sbmv!('U',k,-alpha,A,x,beta,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv! y = beta*y - alpha*A[0:k,'L']*x))
:(Base.LinAlg.BLAS.sbmv!('L',k,-alpha,A,x,beta,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv! y += alpha*A[0:k,'L']*x))
:(Base.LinAlg.BLAS.sbmv!('L',k,alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@sbmv! y = alpha*A['L',0:k]*x + beta*y))
:(Base.LinAlg.BLAS.sbmv!('L',k,alpha,A,x,beta,y))```

### gemm

Return `alpha*A*B`, `alpha*A'*B`, `alpha*A*B'` or `alpha*A'*B'`.

Polynomials

• `A*B`
• `A'*B`
• `A*B'`
• `A'*B'`
• `alpha*A*B`
• `alpha*A'*B`
• `alpha*A*B'`
• `alpha*A'*B'`

Example

```@test macroexpand(SugarBLAS, :(SugarBLAS.@gemm alpha*A*B))
:(Base.LinAlg.BLAS.gemm('N','N',alpha,A,B))

@test macroexpand(SugarBLAS, :(SugarBLAS.@gemm A*B'))
:(Base.LinAlg.BLAS.gemm('N','T',A,B))```

### gemm!

Update `C` as `alpha*A*B + beta*C` or the other three variants according to the combination of transposes of `A` and `B`. Return updated C.

Polynomials

• `C ±= alpha*A*B`
• `C ±= alpha*A'*B`
• `C ±= alpha*A*B'`
• `C ±= alpha*A'*B'`
• `C = beta*C ± alpha*A*B`
• `C = beta*C ± alpha*A'*B`
• `C = beta*C ± alpha*A*B'`
• `C = beta*C ± alpha*A'*B'`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemm! C -= alpha*A*B))
:(Base.LinAlg.BLAS.gemm!('N','N',-alpha,A,B,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemm! C = beta*C - alpha*A*B))
:(Base.LinAlg.BLAS.gemm!('N','N',-alpha,A,B,beta,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemm! C += alpha*A*B))
:(Base.LinAlg.BLAS.gemm!('N','N',alpha,A,B,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemm! C = 3.4*C - alpha*A'*B'))
:(Base.LinAlg.BLAS.gemm!('T','T',-alpha,A,B,3.4,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemm! C = alpha*A'*B + beta*C))
:(Base.LinAlg.BLAS.gemm!('T','N',alpha,A,B,beta,C))```

### gemv

Return `alpha*A*x` or `alpha*A'*x`.

Polynomials

• `A*x`
• `A'*x`
• `alpha*A*x`
• `alpha*A'*x`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv A'*x))
:(Base.LinAlg.BLAS.gemv('T',A,x))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv alpha*A*x))
:(Base.LinAlg.BLAS.gemv('N',alpha,A,x))```

### gemv!

Update the vector `y` as `alpha*A*x + beta*y` or `alpha*A'*x + beta*y`. Return updated `y`.

Polynomials

• `y ±= alpha*A*x`
• `y ±= alpha*A'*x`
• `y = beta*y ± alpha*A*x`
• `y = beta*y ± alpha*A'*x`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv! y -= alpha*A*x))
:(Base.LinAlg.BLAS.gemv!('N',-alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv! y = beta*y - alpha*A*x))
:(Base.LinAlg.BLAS.gemv!('N',-alpha,A,x,beta,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv! y = beta*y - 1.5*A'*x))
:(Base.LinAlg.BLAS.gemv!('T',-1.5,A,x,beta,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv! y += alpha*A*x))
:(Base.LinAlg.BLAS.gemv!('N',alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@gemv! y = alpha*A*x + beta*y))
:(Base.LinAlg.BLAS.gemv!('N',alpha,A,x,beta,y))```

### symm

Return `alpha*A*B` or `alpha*B*A` according to `"symm"`. `A` is assumed to be symmetric. Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper).

Polynomials

• `A["symm", uplo]*B`
• `A*B["symm", uplo]`
• `alpha*A["symm", uplo]*B `
• `alpha*A*B["symm", uplo]`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@symm alpha*A["symm", 'L']*B))
:(Base.LinAlg.BLAS.symm('L','L',alpha,A,B))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symm A*B["symm", 'U']))
:(Base.LinAlg.BLAS.symm('R','U',A,B))```

### symm!

Update `C` as `alpha*A*B + beta*C` or `alpha*B*A + beta*C` according to `"symm"`. `A` is assumed to be symmetric. Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper). Return updated `C`.

Polynomials

• `C = alpha*A["symm",uplo]*B`
• `C = alpha*A*B["symm",uplo]`
• `C = beta*C ± alpha*A["symm",uplo]*B`
• `C = beta*C ± alpha*A*B["symm",uplo]`

Example

```julia> macroexpand(SugarBLAS, :(SugarBLAS.@symm! C -= alpha*A["symm", 'L']*B))
:(Base.LinAlg.BLAS.symm!('L','L',-alpha,A,B,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symm! C = C - alpha*A["symm", 'U']*B))
:(Base.LinAlg.BLAS.symm!('L','U',-alpha,A,B,1.0,C))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symm! C = beta*C - alpha*A["symm", 'L']*B))
:(Base.LinAlg.BLAS.symm!('L','L',-alpha,A,B,beta,C))```

### symv

Return `alpha*A*x`. `A` is assumed to be symmetric. Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper).

Polynomials

• `A[uplo]*x`
• `alpha*A[uplo]*x`
```julia> macroexpand(SugarBLAS, :(SugarBLAS.@symv alpha*A['U']*x))
:(Base.LinAlg.BLAS.symv('U',alpha,A,x))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symv A['L']*x))
:(Base.LinAlg.BLAS.symv('L',A,x))```

### symv!

Update the vector `y` as `alpha*A*x + beta*y`. `A` is assumed to be symmetric. Only the `uplo` triangle of `A` is used (`'L'` for lower and `'U'` for upper). Return updated y.

Polynomials

• `y ±= alpha*A[uplo]*x`
• `y = beta*y ± alpha*A[uplo]*x`
```julia> macroexpand(SugarBLAS, :(SugarBLAS.@symv! y -= alpha*A['U']*x))
:(Base.LinAlg.BLAS.symv!('U',-alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symv! y = y - alpha*A['L']*x))
(Base.LinAlg.BLAS.symv!('L',-alpha,A,x,1.0,y))

julia> macroexpand(SugarBLAS, :(SugarBLAS.@symv! y = beta*y + alpha*A['U']*x))
:(Base.LinAlg.BLAS.symv!('U',-alpha,A,x,beta,y))```