BLAS and LAPACK reference
StandardPacked.jl makes several BLAS and LAPACK functions related to packed matrices available that are not provided by Base.
The general rule of thumb is: Whenever a function is provided in Base for general matrices, this package makes its packed counterpart available. Note that there are a lot of additional BLAS and LAPACK functions with respect to packed matrices for which no interface is provided because the corresponding general function is not made available in Base.
If you want to have one of these functions included, just open an issue or pull request.
It is sometimes rather mysterious why a function has the s (for symmetric) or h (for Hermitian) prefix: for example,
spr!is the symmetric rank-1 update, and since it is also implemented (though not made available in Base) for symmetric complex-valued matrices, there is a separate instancehpr!that does the Hermitian equivalent.spev!is the symmetric or Hermitian eigensolver (depending on the type), and it does not change its name tohpev!even though LAPACK does. There is no symmetric complex eigensolver.hptrd!is the Hermitian or symmetric tridiagonalized (depending on the type), and it does not change its name tosptrd!even though LAPACK does.
This mystery comes from following Base's convention of names.
Note that for every function, this package provides an undocumented "raw" variant that does almost no checking apart from whether the uplo parameter is correct plus a check of the return value. This raw version is usually not made available directly by the Julia Base.
Additionally, StandardPacked.jl provide a convenience wrapper that works with the AbstractArray interface, automatically infers dimensions and checks for the compatibility of all the sizes - this is what is commonly provided by Julia. And finally, the vector representing the packed matrix can also be replaced by a SPMatrixUnscaled - then, the uplo argument is not present, as it is determined from the type.
Even the convenience wrappers give more control than you might be used to. Whenever some allocations might occur due to the need for temporaries or additional output arrays, these wrappers will allow to pass preallocated arrays if they are sized appropriately (where output arguments must match exactly in size, but temporaries may also be larger than necessary).
Whenever a parameter is documented with a missing default, this is such a parameter that might be preallocated if necessary. The parameters are all passed by position, not by keyword, which allows Julia to remove unnecessary checks by compiling the correct version.
BLAS Level 2
Some of these functions are already present in Base.BLAS and are re-exported in StandardPacked. They are listed here for completeness only.
LinearAlgebra.BLAS.spmv! — Functionspmv!(α, AP::SPMatrixUnscaled, x, β)LinearAlgebra.BLAS.hpmv! — Functionhpmv!(α, AP::SPMatrixUnscaled, x, β, y)LinearAlgebra.BLAS.spr! — Functionspr!(α, x, AP::SPMatrix)The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.hpr! — Functionhpr!(uplo, α, x, AP::AbstractVector)
hpr!(α, x, AP::SPMatrix)Update matrix $A$ as $A + \alpha x x'$, where $A$ is a Hermitian matrix provided in packed format AP and x is a vector.
With uplo = 'U', the vector AP must contain the upper triangular part of the Hermitian matrix packed sequentially, column by column, so that AP[1] contains A[1, 1], AP[2] and AP[3] contain A[1, 2] and A[2, 2] respectively, and so on.
With uplo = 'L', the vector AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so that AP[1] contains A[1, 1], AP[2] and AP[3] contain A[2, 1] and A[3, 1] respectively, and so on.
The scalar input α must be real.
The array inputs x and AP must all be of ComplexF32 or ComplexF64 type. Return the updated AP.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
BLAS Level 3
StandardPacked.gemmt! — Functiongemmt!(uplo, transA, transB, alpha, A, B, beta, C)gemmt! computes a matrix-matrix product with general matrices but updates only the upper or lower triangular part of the result matrix.
This function is a recent BLAS extension; for OpenBLAS, it requires at least version 0.3.22 (which is not yet shipped with Julia). If the currently available BLAS does not offer gemmt, the function falls back to gemm.
LAPACK
Conversion
StandardPacked.trttp! — Functiontrttp!(uplo, A, AP::AbstractVector)
trttp!(A, AP::SPMatrixUnscaled)trttp! copies a triangular matrix from the standard full format (TR) to the standard packed format (TP).
This function is also implemented in plain Julia and therefore works with arbitrary element types.
StandardPacked.tpttr! — Functiontpttr!(uplo, AP::AbstractVector, A)
tpttr!(AP::SPMatrixUnscaled, A)tpttr! copies a triangular matrix from the standard packed format (TP) to the standard full format (TR).
This function is also implemented in plain Julia and therefore works with arbitrary element types.
Linear systems - Cholesky
StandardPacked.pptrf! — Functionpptrf!(uplo, AP::AbstractVector) -> (AP, info)
pptrf!(AP::SPMatrix) -> (AP, info)pptrf! computes the Cholesky factorization of a real symmetric or complex Hermitian positive definite matrix $A$ stored in packed format AP.
The factorization has the form
- $A = U' U$, if
uplo == 'U', or - $A = L L'$, if
uplo == 'L',
where $U$ is an upper triangular matrix and $L$ is lower triangular.
If info > 0, the leading minor of order info is not positive definite, and the factorization could not be completed.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.pptrs! — Functionpptrs!(uplo, AP::AbstractVector, B)
pptrs!(AP::SPMatrixUnscaled, B)pptrs! solves a system of linear equations $A X = B$ with a symmetric or Hermitian positive definite matrix $A$ in packed storage AP using the Cholesky factorization $A = U' U$ or $A = L L'$ computed by pptrf!.
StandardPacked.pptri! — Functionpptri!(uplo, AP::AbstractVector)
pptri!(AP::SPMatrixUnscaled)pptri! computes the inverse of a real symmetric or complex Hermitian positive definite matrix A using the Cholesky factorization $A = U' U$ or $A = L L'$ computed by pptrf!.
Linear systems - symmetric indefinite
StandardPacked.spsv! — Functionspsv!(uplo, AP::AbstractVector, ipiv=missing, B) -> (B, AP, ipiv)
spsv!(AP::SPMatrix, ipiv=missing, B) -> (B, AP, ipiv)spsv! computes the solution to a real or complex system of linear equations $A X = B$, where $A$ is an $N$-by-$N$ symmetric matrix stored in packed format AP and $X$ and B are $N$-by-$N_{\mathrm{rhs}}$ matrices.
The diagonal pivoting method is used to factor $A$ as
- $A = U D U^\top$, if
UPLO = 'U', or - $A = L D L^\top$, if
UPLO = 'L',
where U (or L) is a product of permutation and unit upper (lower) triangular matrices, D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of $A$ is then used to solve the system of equations $A X = B$.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.hpsv! — Functionhpsv!(uplo, AP::AbstractVector, ipiv=missing, B) -> (B, AP, ipiv)
hpsv!(AP::SPMatrix, ipiv=missing, B) -> (B, AP, ipiv)hpsv! computes the solution to a complex system of linear equations $A X = B$, where $A$ is an $N$-by-$N$ Hermitian matrix stored in packed format AP and $X$ and B are $N$-by-$N_{\mathrm{rhs}}$ matrices.
The diagonal pivoting method is used to factor $A$ as
- $A = U D U'$, if
UPLO = 'U', or - $A = L D L'$, if
UPLO = 'L',
where U (or L) is a product of permutation and unit upper (lower) triangular matrices, D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of $A$ is then used to solve the system of equations $A X = B$.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.sptrf! — Functionsptrf!(uplo, AP::AbstractVector, ipiv=missing) -> (AP, ipiv, info)
sptrf!(AP::SPMatrix, ipiv=missing) -> (AP, ipiv, info)sptrf! computes the factorization of a real or complex symmetric matrix $A$ stored in packed format AP using the Bunch-Kaufman diagonal pivoting method:
$A = U D U^\top$ or $A = L D L^\top$
where $U$ (or $L$) is a product of permutation and unit upper (lower) triangular matrices, and $D$ is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.hptrf! — Functionhptrf!(uplo, AP::AbstractVector, ipiv=missing) -> (AP, ipiv, info)
hptrf!(AP::SPMatrix, ipiv=missing) -> (AP, ipiv, info)hptrf! computes the factorization of a complex Hermitian matrix $A$ stored in packed format AP using the Bunch-Kaufman diagonal pivoting method:
$A = U D U'$ or $A = L D L'$
where $U$ (or $L$) is a product of permutation and unit upper (lower) triangular matrices, and $D$ is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.sptrs! — Functionsptrs!(uplo, AP::AbstractVector, ipiv, B)
sptrs!(AP::SPMatrixUnscaled, ipiv, B)sptrs! solves a system of linear equations $A X = B$ with a real or complex symmetric matrix $A$ stored in packed format using the factorization $A = U D U^\top$ or $A = L D L^\top$ computed by sptrf!.
StandardPacked.hptrs! — Functionhptrs!(uplo, AP::AbstractVector, ipiv, B)
hptrs!(AP::SPMatrixUnscaled, ipiv, B)hptrs! solves a system of linear equations $A X = B$ with a complex Hermitian matrix $A$ stored in packed format using the factorization $A = U D U'$ or $A = L D L'$ computed by hptrf!.
StandardPacked.sptri! — Functionsptri!(uplo, AP::AbstractVector, ipiv, work=missing)
sptri!(AP::SPMatrixUnscaled, ipiv, work=missing)sptri! computes the inverse of a real or complex symmetric indefinite matrix $A$ in packed storage AP using the factorization $A = U D U^\top$ or $A = L D L^\top$ computed by sptrf!.
StandardPacked.hptri! — Functionhptri!(uplo, AP::AbstractVector, ipiv, work=missing)
hptri!(AP::SPMatrixUnscaled, ipiv, work=missing)hptri! computes the inverse of a complex Hermitian indefinite matrix $A$ in packed storage AP using the factorization $A = U D U'$ or $A = L D L'$ computed by hptrf!.
Eigenvalue
StandardPacked.spev! — Functionspev!('N', uplo, AP::AbstractVector, W=missing, work=missing[, rwork=missing]) -> W
spev!('N', AP::SPMatrix, W=missing, work=missing[, rwork=missing]) -> W
spev!('V', uplo, AP::AbstractVector, W=missing, Z=missing, work=missing
[, rwork=missing]) -> (W, Z)
spev!('V', AP::SPMatrix, W=missing, Z=missing, work=missing[, rwork=missing])
-> (W, Z)Finds the eigenvalues (first parameter 'N') or eigenvalues and eigenvectors (first parameter 'V') of a Hermitian matrix $A$ in packed storage AP. If uplo = 'U', AP is the upper triangle of $A$. If uplo = 'L', AP is the lower triangle of $A$.
The output vector W and matrix Z, as well as the temporaries work and rwork (in the complex case) may be specified to save allocations.
StandardPacked.spevx! — Functionspevx!('N', range, uplo, AP::AbstractVector, vl, vu, il, iu, abstol, W=missing,
work=missing[, rwork=missing], iwork=missing) -> view(W)
spevx!('N', range, AP::SPMatrix, vl, vu, il, iu, abstol, W=missing, work=missing
[, rwork=missing], iwork=missing) -> view(W)
spevx!('V', range, uplo, AP::AbstractVector, vl, vu, il, iu, abstol, W=missing,
Z=missing, work=missing[, rwork=missing], iwork=missing, ifail=missing)
-> (view(W), view(Z), info, view(ifail))
spevx!('V', range, AP::SPMatrix, vl, vu, il, iu, abstol, W=missing, Z=missing,
work=missing[, rwork=missing], iwork=missing, ifail=missing)
-> (view(W), view(Z), info, view(ifail))Finds the eigenvalues (first parameter 'N') or eigenvalues and eigenvectors (first parameter 'V') of a Hermitian matrix $A$ in packed storage AP. If uplo = 'U', AP is the upper triangle of $A$. If uplo = 'L', AP is the lower triangle of $A$. If range = 'A', all the eigenvalues are found. If range = 'V', the eigenvalues in the half-open interval (vl, vu] are found. If range = 'I', the eigenvalues with indices between il and iu are found. abstol can be set as a tolerance for convergence.
The eigenvalues are returned in W and the eigenvectors in Z. info is zero upon success or the number of eigenvalues that failed to converge; their indices are stored in ifail. The resulting views to the original data are sized according to the number of eigenvalues that fulfilled the bounds given by the range.
The output vector W, ifail and matrix Z, as well as the temporaries work, rwork (in the complex case), and iwork may be specified to save allocations.
StandardPacked.spevd! — Functionspevd!('N', uplo, AP::AbstractVector, W=missing, work=missing[, rwork=missing],
iwork=missing) -> W
spevd!('N', AP::SPMatrix, W=missing, work=missing[, rwork=missing],
iwork=missing) -> W
spevd!('V', uplo, AP::AbstractVector, W=missing, Z=missing, work=missing
[, rwork=missing], iwork=missing) -> (W, Z)
spevd!('V', AP::SPMatrix, W=missing, Z=missing, work=missing[, rwork=missing],
iwork=missing) -> (W, Z)Finds the eigenvalues (first parameter 'N') or eigenvalues and eigenvectors (first parameter 'V') of a Hermitian matrix $A$ in packed storage AP. If uplo = 'U', AP is the upper triangle of $A$. If uplo = 'L', AP is the lower triangle of $A$.
The output vector W and matrix Z, as well as the temporaries work, rwork (in the complex case), and iwork may be specified to save allocations.
If eigenvectors are desired, it uses a divide and conquer algorithm.
StandardPacked.spgv! — Functionspgv!(itype, 'N', uplo, AP::AbstractVector, BP::AbstractVector, W=missing,
work=missing[, rwork=missing]) -> (W, BP)
spgv!(itype, 'N', AP::SPMatrix, BP::SPMatrix, W=missing, work=missing
[, rwork=missing]) -> (W, BP)
spgv!(itype, 'V', uplo, AP::AbstractVector, BP::AbstractVector, W=missing, Z=missing,
work=missing[, rwork=missing]) -> (W, Z, BP)
spgv!(itype, 'V', AP::SPMatrix, BP::SPMatrix, W=missing, Z=missing,
work=missing[, rwork=missing]) -> (W, Z, BP)Finds the generalized eigenvalues (second parameter 'N') or eigenvalues and eigenvectors (second parameter 'V') of a Hermitian matrix $A$ in packed storage AP and Hermitian positive-definite matrix $B$ in packed storage BP. If uplo = 'U', AP and BP are the upper triangles of $A$ and $B$, respectively. If uplo = 'L', they are the lower triangles.
- If
itype = 1, the problem to solve is $A x = \lambda B x$. - If
itype = 2, the problem to solve is $A B x = \lambda x$. - If
itype = 3, the problem to solve is $B A x = \lambda x$.
On exit, $B$ is overwritten with the triangular factor $U$ or $L$ from the Cholesky factorization $B = U' U$ or $B = L L'$.
The output vector W and matrix Z, as well as the temporaries work and rwork (in the complex case) may be specified to save allocations.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.spgvx! — Functionspgvx!(itype, 'N', range, uplo, AP::AbstractVector, BP::AbstractVector, vl, vu,
il, iu, abstol, W=missing, work=missing[, rwork=missing], iwork=missing)
-> (view(W), BP)
spgvx!(itype, 'N', range, AP::SPMatrix, BP::SPMatrix, vl, vu, il, iu, abstol,
W=missing, work=missing[, rwork=missing], iwork=missing) -> (view(W), BP)
spgvx!(itype, 'V', range, uplo, AP::AbstractVector, BP::AbstractVector, vl, vu,
il, iu, abstol, W=missing, Z=missing, work=missing[, rwork=missing],
iwork=missing, ifail=missing) -> (view(W), view(Z), info, view(ifail), BP)
spgvx!(itype, 'V', range, AP::SPMatrix, BP::SPMatrix, vl, vu, il, iu, abstol,
W=missing, Z=missing, work=missing[, rwork=missing], iwork=missing, ifail=missing)
-> (view(W), view(Z), info, view(ifail), BP)Finds the generalized eigenvalues (second parameter 'N') or eigenvalues and eigenvectors (second parameter 'V') of a Hermitian matrix $A$ in packed storage AP and Hermitian positive-definite matrix $B$ in packed storage BP. If uplo = 'U', AP and BP are the upper triangles of $A$ and $B$, respectively. If uplo = 'L', they are the lower triangles. If range = 'A', all the eigenvalues are found. If range = 'V', the eigenvalues in the half-open interval (vl, vu] are found. If range = 'I', the eigenvalues with indices between il and iu are found. abstol can be set as a tolerance for convergence.
- If
itype = 1, the problem to solve is $A x = \lambda B x$. - If
itype = 2, the problem to solve is $A B x = \lambda x$. - If
itype = 3, the problem to solve is $B A x = \lambda x$.
On exit, $B$ is overwritten with the triangular factor $U$ or $L$ from the Cholesky factorization $B = U' U$ or $B = L L'$.
The eigenvalues are returned in W and the eigenvectors in Z. info is zero upon success or the number of eigenvalues that failed to converge; their indices are stored in ifail. The resulting views to the original data are sized according to the number of eigenvalues that fulfilled the bounds given by the range.
The output vector W, ifail and matrix Z, as well as the temporaries work, rwork (in the complex case), and iwork may be specified to save allocations.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.spgvd! — Functionspgvd!(itype, 'N', uplo, AP::AbstractVector, BP::AbstractVector, W=missing,
work=missing[, rwork=missing], iwork=missing) -> (W, BP)
spgvd!(itype, 'N', AP::SPMatrix, BP::SPMatrix, W=missing, work=missing
[, rwork=missing], iwork=missing) -> (W, BP)
spgvd!(itype, 'V', uplo, AP::AbstractVector, BP::AbstractVector, W=missing, Z=missing,
work=missing[, rwork=missing], iwork=missing) -> (W, Z, BP)
spgvd!(itype, 'V', AP::SPMatrix, BP::SPMatrix, W=missing, Z=missing,
work=missing[, rwork=missing], iwork=missing) -> (W, Z, BP)Finds the generalized eigenvalues (second parameter 'N') or eigenvalues and eigenvectors (second parameter 'V') of a Hermitian matrix $A$ in packed storage AP and Hermitian positive-definite matrix $B$ in packed storage BP. If uplo = 'U', AP and BP are the upper triangles of $A$, and $B$, respectively. If uplo = 'L', they are the lower triangles.
- If
itype = 1, the problem to solve is $A x = \lambda B x$. - If
itype = 2, the problem to solve is $A B x = \lambda x$. - If
itype = 3, the problem to solve is $B A x = \lambda x$.
On exit, $B$ is overwritten with the triangular factor $U$ or $L$ from the Cholesky factorization $B = U' U$ or $B = L L'$.
The output vector W and matrix Z, as well as the temporaries work, rwork (in the complex case), and iwork may be specified to save allocations.
If eigenvectors are desired, it uses a divide and conquer algorithm.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
Tridiagonal
StandardPacked.hptrd! — Functionhptrd!(uplo, AP::AbstractVector, D=missing, E=missing, τ=missing) -> (A, τ, D, E)
hptrd!(AP::SPMatrix, D=missing, E=missing, τ=missing) -> (A, τ, D, E)hptrd! reduces a Hermitian matrix $A$ stored in packed form AP to real-symmetric tridiagonal form $T$ by an orthogonal similarity transformation: $Q' A Q = T$. If uplo = 'U', AP is the upper triangle of $A$, else it is the lower triangle.
On exit, if uplo = 'U', the diagonal and first superdiagonal of $A$ are overwritten by the corresponding elements of the tridiagonal matrix $T$, and the elements above the first superdiagonal, with the array τ, represent the orthogonal matrix $Q$ as a product of elementary reflectors. If uplo = 'L' the diagonal and first subdiagonal of $A$ are over-written by the corresponding elements of the tridiagonal matrix $T$, and the elements below the first subdiagonal, with the array τ, represent the orthogonal matrix $Q$ as a product of elementary reflectors.
tau contains the scalar factors of the elementary reflectors, D contains the diagonal elements of $T$ and E contains the off-diagonal elements of $T$.
The output vectors D, E, and τ may be specified to save allocations.
The variant of this function that takes a SPMatrix also allows for scaled packed matrices. It will automatically call packed_unscale! on the matrix and return the unscaled result. Do not use the reference to the scaled matrix any more, only the result of this function!
StandardPacked.opgtr! — Functionopgtr!(uplo, AP::AbstractVector, τ, Q=missing, work=missing)
opgtr!(AP::SPMatrixUnscaled, τ, Q=missing, work=missing)Explicitly finds Q, the orthogonal/unitary matrix from hptrd!. uplo, AP, and τ must correspond to the input/output to hptrd!.
The output matrix Q and the temporary work may be specified to save allocations.
StandardPacked.opmtr! — Functionopmtr!(side, uplo, trans, AP::AbstractVector, τ, C, work=missing)
opmtr!(side, trans, AP::SPMatrixUnscaled, τ, C, work=missing)opmtr! overwrites the general $m \times n$ matrix C with
SIDE = 'L' | SIDE = 'R' | |
|---|---|---|
TRANS = 'N' | $Q C$ | $C Q$ |
TRANS = 'T' or 'C' (complex case) | $Q' C$ | $C Q'$ |
where $Q$ is a real orthogonal or complex unitary matrix of order $n_q$, with $n_q = m$ if SIDE = 'L' and $n_q = n$ if SIDE = 'R'. $Q$ is defined as the product of $n_q -1$ elementary reflectors, as returned by hptrd! using packed storage:
- if
UPLO = 'U', $Q = H_{n_q -1} \dotsm H_2 H_1$; - if
UPLO = 'L', $Q = H_1 H_2 \dotsm H_{n_q -1}$.
The temporary work may be specified to save allocations.