Linear Algebra¶
Linear algebra functions in Julia are largely implemented by calling functions from LAPACK. Sparse factorizations call functions from SuiteSparse.
-
*
(A, B) Matrix multiplication
-
\
(A, B) Matrix division using a polyalgorithm. For input matrices
A
andB
, the resultX
is such thatA*X == B
whenA
is square. The solver that is used depends upon the structure ofA
. A direct solver is used for upper- or lower triangularA
. For HermitianA
(equivalent to symmetricA
for non-complexA
) theBunchKaufman
factorization is used. Otherwise an LU factorization is used. For rectangularA
the result is the minimum-norm least squares solution computed by reducingA
to bidiagonal form and solving the bidiagonal least squares problem. For sparse, squareA
the LU factorization (from UMFPACK) is used.
-
dot
(x, y)¶ Compute the dot product. For complex vectors, the first vector is conjugated.
-
cross
(x, y)¶ Compute the cross product of two 3-vectors.
-
rref
(A)¶ Compute the reduced row echelon form of the matrix A.
-
factorize
(A)¶ Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, Triangular) of A, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example:
A=factorize(A); x=A\\b; y=A\\C
.
-
factorize!
(A)¶ factorize!
is the same asfactorize()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
lu
(A) → L, U, p¶ Compute the LU factorization of
A
, such thatA[p,:] = L*U
.
-
lufact
(A[, pivot=true]) → F¶ Compute the LU factorization of
A
. The return type ofF
depends on the type ofA
. In most cases, ifA
is a subtypeS
of AbstractMatrix with an element typeT`
supporting+
,-
,*
and/
the return type isLU{T,S{T}}
. If pivoting is chosen (default) the element type should also supportabs
and<
. WhenA
is sparse and have element of typeFloat32
,Float64
,Complex{Float32}
, orComplex{Float64}
the return type isUmfpackLU
. Some examples are shown in the table below.Type of input A
Type of output F
Relationship between F
andA
Matrix()
LU
F[:L]*F[:U] == A[F[:p], :]
Tridiagonal()
LU{T,Tridiagonal{T}}
N/A SparseMatrixCSC()
UmfpackLU
F[:L]*F[:U] == Rs .* A[F[:p], F[:q]]
The individual components of the factorization
F
can be accessed by indexing:Component Description LU
LU{T,Tridiagonal{T}}
UmfpackLU
F[:L]
L
(lower triangular) part ofLU
✓ ✓ F[:U]
U
(upper triangular) part ofLU
✓ ✓ F[:p]
(right) permutation Vector
✓ ✓ F[:P]
(right) permutation Matrix
✓ F[:q]
left permutation Vector
✓ F[:Rs]
Vector
of scaling factors✓ F[:(:)]
(L,U,p,q,Rs)
components✓ Supported function LU
LU{T,Tridiagonal{T}}
UmfpackLU
/
✓ \
✓ ✓ ✓ cond
✓ ✓ det
✓ ✓ ✓ size
✓ ✓
-
lufact!
(A) → LU¶ lufact!
is the same aslufact()
, but saves space by overwriting the input A, instead of creating a copy. For sparseA
thenzval
field is not overwritten but the index fields,colptr
androwval
are decremented in place, converting from 1-based indices to 0-based indices.
-
chol
(A[, LU]) → F¶ Compute the Cholesky factorization of a symmetric positive definite matrix
A
and return the matrixF
. IfLU
is:L
(Lower),A = L*L'
. IfLU
is:U
(Upper),A = R'*R
.
-
cholfact
(A, [LU,][pivot=false,][tol=-1.0]) → Cholesky¶ Compute the Cholesky factorization of a dense symmetric positive (semi)definite matrix
A
and return either aCholesky
ifpivot=false
orCholeskyPivoted
ifpivot=true
.LU
may be:L
for using the lower part or:U
for the upper part. The default is to use:U
. The triangular matrix can be obtained from the factorizationF
with:F[:L]
andF[:U]
. The following functions are available forCholesky
objects:size
,\
,inv
,det
. ForCholeskyPivoted
there is also defined arank
. Ifpivot=false
aPosDefException
exception is thrown in case the matrix is not positive definite. The argumenttol
determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision.
-
cholfact
(A[, ll]) → CholmodFactor Compute the sparse Cholesky factorization of a sparse matrix
A
. IfA
is Hermitian its Cholesky factor is determined. IfA
is not Hermitian the Cholesky factor ofA*A'
is determined. A fill-reducing permutation is used. Methods forsize
,solve
,\
,findn_nzs
,diag
,det
andlogdet
. One of the solve methods includes an integer argument that can be used to solve systems involving parts of the factorization only. The optional boolean argument,ll
determines whether the factorization returned is of theA[p,p] = L*L'
form, whereL
is lower triangular orA[p,p] = L*Diagonal(D)*L'
form whereL
is unit lower triangular andD
is a non-negative vector. The default is LDL.
-
cholfact!
(A, [LU,][pivot=false,][tol=-1.0]) → Cholesky¶ cholfact!
is the same ascholfact()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
ldltfact
(A) → LDLtFactorization¶ Compute a factorization of a positive definite matrix
A
such thatA=L*Diagonal(d)*L'
whereL
is a unit lower triangular matrix andd
is a vector with non-negative elements.
-
qr
(A, [pivot=false,][thin=true]) → Q, R, [p]¶ Compute the (pivoted) QR factorization of
A
such that eitherA = Q*R
orA[:,p] = Q*R
. Also seeqrfact
. The default is to compute a thin factorization. Note thatR
is not extended with zeros when the fullQ
is requested.
-
qrfact
(A[, pivot=false]) → F¶ Computes the QR factorization of
A
. The return type ofF
depends on the element type ofA
and whether pivoting is specified (withpivot=true
).Return type eltype(A)
pivot
Relationship between F
andA
QR
not BlasFloat
either A==F[:Q]*F[:R]
QRCompactWY
BlasFloat
false
A==F[:Q]*F[:R]
QRPivoted
BlasFloat
true
A[:,F[:p]]==F[:Q]*F[:R]
BlasFloat
refers to any of:Float32
,Float64
,Complex64
orComplex128
.The individual components of the factorization
F
can be accessed by indexing:Component Description QR
QRCompactWY
QRPivoted
F[:Q]
Q
(orthogonal/unitary) part ofQR
✓ ( QRPackedQ
)✓ ( QRCompactWYQ
)✓ ( QRPackedQ
)F[:R]
R
(upper right triangular) part ofQR
✓ ✓ ✓ F[:p]
pivot Vector
✓ F[:P]
(pivot) permutation Matrix
✓ The following functions are available for the
QR
objects:size
,\
. WhenA
is rectangular,\
will return a least squares solution and if the solution is not unique, the one with smallest norm is returned.Multiplication with respect to either thin or full
Q
is allowed, i.e. bothF[:Q]*F[:R]
andF[:Q]*A
are supported. AQ
matrix can be converted into a regular matrix withfull()
which has a named argumentthin
.注解
qrfact
returns multiple types because LAPACK uses several representations that minimize the memory storage requirements of products of Householder elementary reflectors, so that theQ
andR
matrices can be stored compactly rather as two separate dense matrices.The data contained in
QR
orQRPivoted
can be used to construct theQRPackedQ
type, which is a compact representation of the rotation matrix:\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T)\]where \(\tau_i\) is the scale factor and \(v_i\) is the projection vector associated with the \(i^{th}\) Householder elementary reflector.
The data contained in
QRCompactWY
can be used to construct theQRCompactWYQ
type, which is a compact representation of the rotation matrix\[Q = I + Y T Y^T\]where
Y
is \(m \times r\) lower trapezoidal andT
is \(r \times r\) upper triangular. The compact WY representation [Schreiber1989] is not to be confused with the older, WY representation [Bischof1987]. (The LAPACK documentation usesV
in lieu ofY
.)[Bischof1987] C Bischof and C Van Loan, The WY representation for products of Householder matrices, SIAM J Sci Stat Comput 8 (1987), s2-s13. doi:10.1137/0908009 [Schreiber1989] R Schreiber and C Van Loan, A storage-efficient WY representation for products of Householder transformations, SIAM J Sci Stat Comput 10 (1989), 53-57. doi:10.1137/0910005
-
qrfact!
(A[, pivot=false])¶ qrfact!
is the same asqrfact()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
bkfact
(A) → BunchKaufman¶ Compute the Bunch-Kaufman [Bunch1977] factorization of a real symmetric or complex Hermitian matrix
A
and return aBunchKaufman
object. The following functions are available forBunchKaufman
objects:size
,\
,inv
,issym
,ishermitian
.
[Bunch1977] | J R Bunch and L Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation 31:137 (1977), 163-179. url. |
-
bkfact!
(A) → BunchKaufman¶ bkfact!
is the same asbkfact()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
sqrtm
(A)¶ Compute the matrix square root of
A
. IfB = sqrtm(A)
, thenB*B == A
within roundoff error.sqrtm
uses a polyalgorithm, computing the matrix square root using Schur factorizations (schurfact()
) unless it detects the matrix to be Hermitian or real symmetric, in which case it computes the matrix square root from an eigendecomposition (eigfact()
). In the latter situation for positive definite matrices, the matrix square root hasReal
elements, otherwise it hasComplex
elements.
-
eig
(A,[irange,][vl,][vu,][permute=true,][scale=true]) → D, V¶ Compute eigenvalues and eigenvectors of
A
. Seeeigfact()
for details on thebalance
keyword argument.julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]) ([1.0,3.0,18.0], 3x3 Array{Float64,2}: 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0)
eig
is a wrapper aroundeigfact()
, extracting all parts of the factorization to a tuple; where possible, usingeigfact()
is recommended.
-
eig
(A, B) → D, V Computes generalized eigenvalues and vectors of
A
with respect toB
.eig
is a wrapper aroundeigfact()
, extracting all parts of the factorization to a tuple; where possible, usingeigfact()
is recommended.
-
eigvals
(A,[irange,][vl,][vu])¶ Returns the eigenvalues of
A
. IfA
isSymmetric()
,Hermitian()
orSymTridiagonal()
, it is possible to calculate only a subset of the eigenvalues by specifying either aUnitRange()
irange
covering indices of the sorted eigenvalues, or a pairvl
andvu
for the lower and upper boundaries of the eigenvalues.For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option
permute=true
permutes the matrix to become closer to upper triangular, andscale=true
scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default istrue
for both options.
-
eigmax
(A)¶ Returns the largest eigenvalue of
A
.
-
eigmin
(A)¶ Returns the smallest eigenvalue of
A
.
-
eigvecs
(A, [eigvals,][permute=true,][scale=true])¶ Returns the eigenvectors of
A
. Thepermute
andscale
keywords are the same as foreigfact()
.For
SymTridiagonal()
matrices, if the optional vector of eigenvalueseigvals
is specified, returns the specific corresponding eigenvectors.
-
eigfact
(A,[il,][iu,][vl,][vu,][permute=true,][scale=true])¶ Compute the eigenvalue decomposition of
A
and return anEigen
object. IfF
is the factorization object, the eigenvalues can be accessed withF[:values]
and the eigenvectors withF[:vectors]
. The following functions are available forEigen
objects:inv
,det
.If
A
isSymmetric
,Hermitian
orSymTridiagonal
, it is possible to calculate only a subset of the eigenvalues by specifying either a UnitRange`irange
covering indices of the sorted eigenvalues or a pairvl
andvu
for the lower and upper boundaries of the eigenvalues.For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option
permute=true
permutes the matrix to become closer to upper triangular, andscale=true
scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default istrue
for both options.
-
eigfact
(A, B) Compute the generalized eigenvalue decomposition of
A
andB
and return anGeneralizedEigen
object. IfF
is the factorization object, the eigenvalues can be accessed withF[:values]
and the eigenvectors withF[:vectors]
.
-
eigfact!
(A[, B])¶ eigfact!
is the same aseigfact()
, but saves space by overwriting the input A (and B), instead of creating a copy.
-
hessfact
(A)¶ Compute the Hessenberg decomposition of
A
and return aHessenberg
object. IfF
is the factorization object, the unitary matrix can be accessed withF[:Q]
and the Hessenberg matrix withF[:H]
. WhenQ
is extracted, the resulting type is theHessenbergQ
object, and may be converted to a regular matrix withfull()
.
-
hessfact!
(A)¶ hessfact!
is the same ashessfact()
, but saves space by overwriting the input A, instead of creating a copy.
-
schurfact
(A) → Schur¶ Computes the Schur factorization of the matrix
A
. The (quasi) triangular Schur factor can be obtained from theSchur
objectF
with eitherF[:Schur]
orF[:T]
and the unitary/orthogonal Schur vectors can be obtained withF[:vectors]
orF[:Z]
such thatA=F[:vectors]*F[:Schur]*F[:vectors]'
. The eigenvalues ofA
can be obtained withF[:values]
.
-
schurfact!
(A)¶ Computer the Schur factorization of
A
, overwritingA
in the process. Seeschurfact()
-
schur
(A) → Schur[:T], Schur[:Z], Schur[:values]¶ See
schurfact()
-
schurfact
(A, B) → GeneralizedSchur Computes the Generalized Schur (or QZ) factorization of the matrices
A
andB
. The (quasi) triangular Schur factors can be obtained from theSchur
objectF
withF[:S]
andF[:T]
, the left unitary/orthogonal Schur vectors can be obtained withF[:left]
orF[:Q]
and the right unitary/orthogonal Schur vectors can be obtained withF[:right]
orF[:Z]
such thatA=F[:left]*F[:S]*F[:right]'
andB=F[:left]*F[:T]*F[:right]'
. The generalized eigenvalues ofA
andB
can be obtained withF[:alpha]./F[:beta]
.
-
schur
(A, B) → GeneralizedSchur[:S], GeneralizedSchur[:T], GeneralizedSchur[:Q], GeneralizedSchur[:Z] See
schurfact()
-
svdfact
(A[, thin=true]) → SVD¶ Compute the Singular Value Decomposition (SVD) of
A
and return anSVD
object.U
,S
,V
andVt
can be obtained from the factorizationF
withF[:U]
,F[:S]
,F[:V]
andF[:Vt]
, such thatA = U*diagm(S)*Vt
. Ifthin
istrue
, an economy mode decomposition is returned. The algorithm producesVt
and henceVt
is more efficient to extract thanV
. The default is to produce a thin decomposition.
-
svdfact!
(A[, thin=true]) → SVD¶ svdfact!
is the same assvdfact()
, but saves space by overwriting the input A, instead of creating a copy. Ifthin
istrue
, an economy mode decomposition is returned. The default is to produce a thin decomposition.
-
svd
(A[, thin=true]) → U, S, V¶ Wrapper around
svdfact
extracting all parts the factorization to a tuple. Direct use ofsvdfact
is therefore generally more efficient. Computes the SVD of A, returningU
, vectorS
, andV
such thatA == U*diagm(S)*V'
. Ifthin
istrue
, an economy mode decomposition is returned. The default is to produce a thin decomposition.
-
svdvals
(A)¶ Returns the singular values of
A
.
-
svdvals!
(A)¶ Returns the singular values of
A
, while saving space by overwriting the input.
-
svdfact
(A, B) → GeneralizedSVD Compute the generalized SVD of
A
andB
, returning aGeneralizedSVD
Factorization objectF
, such thatA = F[:U]*F[:D1]*F[:R0]*F[:Q]'
andB = F[:V]*F[:D2]*F[:R0]*F[:Q]'
.
-
svd
(A, B) → U, V, Q, D1, D2, R0 Wrapper around
svdfact
extracting all parts the factorization to a tuple. Direct use ofsvdfact
is therefore generally more efficient. The function returns the generalized SVD ofA
andB
, returningU
,V
,Q
,D1
,D2
, andR0
such thatA = U*D1*R0*Q'
andB = V*D2*R0*Q'
.
-
svdvals
(A, B) Return only the singular values from the generalized singular value decomposition of
A
andB
.
-
triu
(M)¶ Upper triangle of a matrix.
-
triu!
(M)¶ Upper triangle of a matrix, overwriting
M
in the process.
-
tril
(M)¶ Lower triangle of a matrix.
-
tril!
(M)¶ Lower triangle of a matrix, overwriting
M
in the process.
-
diagind
(M[, k])¶ A
Range
giving the indices of thek
-th diagonal of the matrixM
.
-
diag
(M[, k])¶ The
k
-th diagonal of a matrix, as a vector. Usediagm
to construct a diagonal matrix.
-
diagm
(v[, k])¶ Construct a diagonal matrix and place
v
on thek
-th diagonal.
-
scale
(A, b)¶
-
scale
(b, A) Scale an array
A
by a scalarb
, returning a new array.If
A
is a matrix andb
is a vector, thenscale(A,b)
scales each columni
ofA
byb[i]
(similar toA*diagm(b)
), whilescale(b,A)
scales each rowi
ofA
byb[i]
(similar todiagm(b)*A
), returning a new array.Note: for large
A
,scale
can be much faster thanA .* b
orb .* A
, due to the use of BLAS.
-
scale!
(A, b)¶
-
scale!
(b, A) Scale an array
A
by a scalarb
, similar toscale()
but overwritingA
in-place.If
A
is a matrix andb
is a vector, thenscale!(A,b)
scales each columni
ofA
byb[i]
(similar toA*diagm(b)
), whilescale!(b,A)
scales each rowi
ofA
byb[i]
(similar todiagm(b)*A
), again operating in-place onA
.
-
Tridiagonal
(dl, d, du)¶ Construct a tridiagonal matrix from the lower diagonal, diagonal, and upper diagonal, respectively. The result is of type
Tridiagonal
and provides efficient specialized linear solvers, but may be converted into a regular matrix withfull()
.
-
Bidiagonal
(dv, ev, isupper)¶ Constructs an upper (
isupper=true
) or lower (isupper=false
) bidiagonal matrix using the given diagonal (dv
) and off-diagonal (ev
) vectors. The result is of typeBidiagonal
and provides efficient specialized linear solvers, but may be converted into a regular matrix withfull()
.
-
SymTridiagonal
(d, du)¶ Construct a real symmetric tridiagonal matrix from the diagonal and upper diagonal, respectively. The result is of type
SymTridiagonal
and provides efficient specialized eigensolvers, but may be converted into a regular matrix withfull()
.
-
Woodbury
(A, U, C, V)¶ Construct a matrix in a form suitable for applying the Woodbury matrix identity.
-
rank
(M)¶ Compute the rank of a matrix.
-
norm
(A[, p])¶ Compute the
p
-norm of a vector or the operator norm of a matrixA
, defaulting to thep=2
-norm.For vectors,
p
can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular,norm(A, Inf)
returns the largest value inabs(A)
, whereasnorm(A, -Inf)
returns the smallest.For matrices, valid values of
p
are1
,2
, orInf
. (Note that for sparse matrices,p=2
is currently not implemented.) Usevecnorm()
to compute the Frobenius norm.
-
vecnorm
(A[, p])¶ For any iterable container
A
(including arrays of any dimension) of numbers, compute thep
-norm (defaulting top=2
) as ifA
were a vector of the corresponding length.For example, if
A
is a matrix andp=2
, then this is equivalent to the Frobenius norm.
-
cond
(M[, p])¶ Condition number of the matrix
M
, computed using the operatorp
-norm. Valid values forp
are1
,2
(default), orInf
.
-
condskeel
(M[, x, p])¶ - \[\begin{split}\kappa_S(M, p) & = \left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \right\Vert_p \\ \kappa_S(M, x, p) & = \left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \left\vert x \right\vert \right\Vert_p\end{split}\]
Skeel condition number \(\kappa_S\) of the matrix
M
, optionally with respect to the vectorx
, as computed using the operatorp
-norm.p
isInf
by default, if not provided. Valid values forp
are1
,2
, orInf
.This quantity is also known in the literature as the Bauer condition number, relative condition number, or componentwise relative condition number.
-
trace
(M)¶ Matrix trace
-
det
(M)¶ Matrix determinant
-
logdet
(M)¶ Log of matrix determinant. Equivalent to
log(det(M))
, but may provide increased accuracy and/or speed.
-
inv
(M)¶ Matrix inverse
-
pinv
(M)¶ Moore-Penrose pseudoinverse
-
null
(M)¶ Basis for nullspace of
M
.
-
repmat
(A, n, m)¶ Construct a matrix by repeating the given matrix
n
times in dimension 1 andm
times in dimension 2.
-
repeat
(A, inner = Int[], outer = Int[])¶ Construct an array by repeating the entries of
A
. The i-th element ofinner
specifies the number of times that the individual entries of the i-th dimension ofA
should be repeated. The i-th element ofouter
specifies the number of times that a slice along the i-th dimension ofA
should be repeated.
-
kron
(A, B)¶ Kronecker tensor product of two vectors or two matrices.
-
blkdiag
(A...)¶ Concatenate matrices block-diagonally. Currently only implemented for sparse matrices.
-
linreg
(x, y) → [a; b]¶ Linear Regression. Returns
a
andb
such thata+b*x
is the closest line to the given points(x,y)
. In other words, this function determines parameters[a, b]
that minimize the squared error betweeny
anda+b*x
.Example:
using PyPlot; x = float([1:12]) y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99] a, b = linreg(x,y) # Linear regression plot(x, y, "o") # Plot (x,y) points plot(x, [a+b*i for i in x]) # Plot the line determined by the linear regression
-
linreg
(x, y, w) Weighted least-squares linear regression.
-
expm
(A)¶ Matrix exponential.
-
lyap
(A, C)¶ Computes the solution
X
to the continuous Lyapunov equationAX + XA' + C = 0
, where no eigenvalue ofA
has a zero real part and no two eigenvalues are negative complex conjugates of each other.
-
sylvester
(A, B, C)¶ Computes the solution
X
to the Sylvester equationAX + XB + C = 0
, whereA
,B
andC
have compatible dimensions andA
and-B
have no eigenvalues with equal real part.
-
issym
(A) → Bool¶ Test whether a matrix is symmetric.
-
isposdef
(A) → Bool¶ Test whether a matrix is positive definite.
-
isposdef!
(A) → Bool¶ Test whether a matrix is positive definite, overwriting
A
in the processes.
-
istril
(A) → Bool¶ Test whether a matrix is lower triangular.
-
istriu
(A) → Bool¶ Test whether a matrix is upper triangular.
-
ishermitian
(A) → Bool¶ Test whether a matrix is Hermitian.
-
transpose
(A)¶ The transposition operator (
.'
).
-
ctranspose
(A)¶ The conjugate transposition operator (
'
).
-
eigs
(A, [B, ]; nev=6, which="LM", tol=0.0, maxiter=1000, sigma=nothing, ritzvec=true, v0=zeros((0, ))) -> (d, [v, ]nconv, niter, nmult, resid)¶ eigs
computes eigenvaluesd
ofA
using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. IfB
is provided, the generalized eigen-problem is solved. The following keyword arguments are supported:nev
: Number of eigenvaluesncv
: Number of Krylov vectors used in the computation; should satisfynev+1 <= ncv <= n
for real symmetric problems andnev+2 <= ncv <= n
for other problems; default isncv = max(20,2*nev+1)
.which
: type of eigenvalues to compute. See the note below.which
type of eigenvalues :LM
eigenvalues of largest magnitude (default) :SM
eigenvalues of smallest magnitude :LR
eigenvalues of largest real part :SR
eigenvalues of smallest real part :LI
eigenvalues of largest imaginary part (nonsymmetric or complex A
only):SI
eigenvalues of smallest imaginary part (nonsymmetric or complex A
only):BE
compute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric A
only)tol
: tolerance (\(tol \le 0.0\) defaults toDLAMCH('EPS')
)maxiter
: Maximum number of iterations (default = 300)sigma
: Specifies the level shift used in inverse iteration. Ifnothing
(default), defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close tosigma
using shift and invert iterations.ritzvec
: Returns the Ritz vectorsv
(eigenvectors) iftrue
v0
: starting vector from which to start the iterations
eigs
returns thenev
requested eigenvalues ind
, the corresponding Ritz vectorsv
(only ifritzvec=true
), the number of converged eigenvaluesnconv
, the number of iterationsniter
and the number of matrix vector multiplicationsnmult
, as well as the final residual vectorresid
.注解
The
sigma
andwhich
keywords interact: the description of eigenvalues searched for bywhich
do _not_ necessarily refer to the eigenvalues ofA
, but rather the linear operator constructed by the specification of the iteration mode implied bysigma
.sigma
iteration mode which
refers to eigenvalues ofnothing
ordinary (forward) \(A\) real or complex inverse with level shift sigma
\((A - \sigma I )^{-1}\)
-
peakflops
(n; parallel=false)¶ peakflops
computes the peak flop rate of the computer by using BLAS double precisiongemm!()
. By default, if no arguments are specified, it multiplies a matrix of sizen x n
, wheren = 2000
. If the underlying BLAS is using multiple threads, higher flop rates are realized. The number of BLAS threads can be set withblas_set_num_threads(n)
.If the keyword argument
parallel
is set totrue
,peakflops
is run in parallel on all the worker processors. The flop rate of the entire parallel computer is returned. When running in parallel, only 1 BLAS thread is used. The argumentn
still refers to the size of the problem that is solved on each processor.
BLAS Functions¶
This module provides wrappers for some of the BLAS functions for
linear algebra. Those BLAS functions that overwrite one of the input
arrays have names ending in '!'
.
Usually a function has 4 methods defined, one each for Float64
,
Float32
, Complex128
and Complex64
arrays.
-
dot
(n, X, incx, Y, incy)¶ Dot product of two vectors consisting of
n
elements of arrayX
with strideincx
andn
elements of arrayY
with strideincy
.
-
dotu
(n, X, incx, Y, incy)¶ Dot function for two complex vectors.
-
dotc
(n, X, incx, U, incy)¶ Dot function for two complex vectors conjugating the first vector.
-
blascopy!
(n, X, incx, Y, incy)¶ Copy
n
elements of arrayX
with strideincx
to arrayY
with strideincy
. ReturnsY
.
-
nrm2
(n, X, incx)¶ 2-norm of a vector consisting of
n
elements of arrayX
with strideincx
.
-
asum
(n, X, incx)¶ sum of the absolute values of the first
n
elements of arrayX
with strideincx
.
-
axpy!
(n, a, X, incx, Y, incy)¶ Overwrite
Y
witha*X + Y
. ReturnsY
.
-
scal!
(n, a, X, incx)¶ Overwrite
X
witha*X
. ReturnsX
.
-
scal
(n, a, X, incx)¶ Returns
a*X
.
-
syrk!
(uplo, trans, alpha, A, beta, C)¶ Rank-k update of the symmetric matrix
C
asalpha*A*A.' + beta*C
oralpha*A.'*A + beta*C
according to whethertrans
is ‘N’ or ‘T’. Whenuplo
is ‘U’ the upper triangle ofC
is updated (‘L’ for lower triangle). ReturnsC
.
-
syrk
(uplo, trans, alpha, A)¶ Returns either the upper triangle or the lower triangle, according to
uplo
(‘U’ or ‘L’), ofalpha*A*A.'
oralpha*A.'*A
, according totrans
(‘N’ or ‘T’).
-
herk!
(uplo, trans, alpha, A, beta, C)¶ Methods for complex arrays only. Rank-k update of the Hermitian matrix
C
asalpha*A*A' + beta*C
oralpha*A'*A + beta*C
according to whethertrans
is ‘N’ or ‘T’. Whenuplo
is ‘U’ the upper triangle ofC
is updated (‘L’ for lower triangle). ReturnsC
.
-
herk
(uplo, trans, alpha, A)¶ Methods for complex arrays only. Returns either the upper triangle or the lower triangle, according to
uplo
(‘U’ or ‘L’), ofalpha*A*A'
oralpha*A'*A
, according totrans
(‘N’ or ‘T’).
-
gbmv!
(trans, m, kl, ku, alpha, A, x, beta, y)¶ Update vector
y
asalpha*A*x + beta*y
oralpha*A'*x + beta*y
according totrans
(‘N’ or ‘T’). The matrixA
is a general band matrix of dimensionm
bysize(A,2)
withkl
sub-diagonals andku
super-diagonals. Returns the updatedy
.
-
gbmv
(trans, m, kl, ku, alpha, A, x, beta, y)¶ Returns
alpha*A*x
oralpha*A'*x
according totrans
(‘N’ or ‘T’). The matrixA
is a general band matrix of dimensionm
bysize(A,2)
withkl
sub-diagonals andku
super-diagonals.
-
sbmv!
(uplo, k, alpha, A, x, beta, y)¶ Update vector
y
asalpha*A*x + beta*y
whereA
is a a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
. The storage layout forA
is described the reference BLAS module, level-2 BLAS at http://www.netlib.org/lapack/explore-html/.Returns the updated
y
.
-
sbmv
(uplo, k, alpha, A, x)¶ Returns
alpha*A*x
whereA
is a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
.
-
sbmv
(uplo, k, A, x) Returns
A*x
whereA
is a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
.
-
gemm!
(tA, tB, alpha, A, B, beta, C)¶ Update
C
asalpha*A*B + beta*C
or the other three variants according totA
(transposeA
) andtB
. Returns the updatedC
.
-
gemm
(tA, tB, alpha, A, B)¶ Returns
alpha*A*B
or the other three variants according totA
(transposeA
) andtB
.
-
gemm
(tA, tB, A, B) Returns
A*B
or the other three variants according totA
(transposeA
) andtB
.
-
gemv!
(tA, alpha, A, x, beta, y)¶ Update the vector
y
asalpha*A*x + beta*x
oralpha*A'x + beta*x
according totA
(transposeA
). Returns the updatedy
.
-
gemv
(tA, alpha, A, x)¶ Returns
alpha*A*x
oralpha*A'x
according totA
(transposeA
).
-
gemv
(tA, A, x) Returns
A*x
orA'x
according totA
(transposeA
).
-
symm!
(side, ul, alpha, A, B, beta, C)¶ Update
C
asalpha*A*B + beta*C
oralpha*B*A + beta*C
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used. Returns the updatedC
.
-
symm
(side, ul, alpha, A, B)¶ Returns
alpha*A*B
oralpha*B*A
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symm
(side, ul, A, B) Returns
A*B
orB*A
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symm
(tA, tB, alpha, A, B) Returns
alpha*A*B
or the other three variants according totA
(transposeA
) andtB
.
-
symv!
(ul, alpha, A, x, beta, y)¶ Update the vector
y
asalpha*A*y + beta*y
.A
is assumed to be symmetric. Only theul
triangle ofA
is used. Returns the updatedy
.
-
symv
(ul, alpha, A, x)¶ Returns
alpha*A*x
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symv
(ul, A, x) Returns
A*x
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
trmm!
(side, ul, tA, dA, alpha, A, B)¶ Update
B
asalpha*A*B
or one of the other three variants determined byside
(A on left or right) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedB
.
-
trmm
(side, ul, tA, dA, alpha, A, B)¶ Returns
alpha*A*B
or one of the other three variants determined byside
(A on left or right) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trsm!
(side, ul, tA, dA, alpha, A, B)¶ Overwrite
B
with the solution toA*X = alpha*B
or one of the other three variants determined byside
(A on left or right ofX
) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedB
.
-
trsm
(side, ul, tA, dA, alpha, A, B)¶ Returns the solution to
A*X = alpha*B
or one of the other three variants determined byside
(A on left or right ofX
) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trmv!
(side, ul, tA, dA, alpha, A, b)¶ Update
b
asalpha*A*b
or one of the other three variants determined byside
(A on left or right) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedb
.
-
trmv
(side, ul, tA, dA, alpha, A, b)¶ Returns
alpha*A*b
or one of the other three variants determined byside
(A on left or right) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trsv!
(side, ul, tA, dA, alpha, A, b)¶ Overwrite
b
with the solution toA*X = alpha*b
or one of the other three variants determined byside
(A on left or right ofX
) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedb
.
-
trsv
(side, ul, tA, dA, alpha, A, b)¶ Returns the solution to
A*X = alpha*b
or one of the other three variants determined byside
(A on left or right ofX
) andtA
(transpose A). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
blas_set_num_threads
(n)¶ Set the number of threads the BLAS library should use.