Linear Algebra¶
Standard Functions¶
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 a pivoted QR factorization ofA
and a rank estimate ofA
based on the R factor.When
A
is sparse, a similar polyalgorithm is used. For indefinite matrices, theLDLt
factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices.
vecdot
(x, y)¶For any iterable containers
x
andy
(including arrays of any dimension) of numbers (or any element type for whichdot
is defined), compute the Euclidean dot product (the sum ofdot(x[i],y[i])
) as if they were vectors.
factorize
(A)¶Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, LowerTriangular, UpperTriangular) 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
.
full
(F)¶Reconstruct the matrix
A
from the factorizationF=factorize(A)
.
lu
(A) → L, U, p¶Compute the LU factorization of
A
, such thatA[p,:]=L*U
.
lufact
(A[, pivot=Val{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}}
F[:L]*F[:U]==A[F[:p],:]
SparseMatrixCSC()
UmfpackLU
F[:L]*F[:U]==(F[: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
✓ ✓ ✓ logdet
✓ ✓ logabsdet
✓ ✓ size
✓ ✓
lufact!
(A) → LU¶lufact!
is the same aslufact()
, but saves space by overwriting the inputA
, 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
isVal{:U}
(Upper),F
is of typeUpperTriangular
andA=F'*F
. IfLU
isVal{:L}
(Lower),F
is of typeLowerTriangular
andA=F*F'
.LU
defaults toVal{:U}
.
cholfact
(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) → Cholesky¶Compute the Cholesky factorization of a dense symmetric positive (semi)definite matrix
A
and return either aCholesky
ifpivot==Val{false}
orCholeskyPivoted
ifpivot==Val{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==Val{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; shift=0, perm=Int[]) → CHOLMOD.FactorCompute the Cholesky factorization of a sparse positive definite matrix
A
. A fill-reducing permutation is used.F=cholfact(A)
is most frequently used to solve systems of equations withF\b
, but also the methodsdiag
,det
,logdet
are defined forF
. You can also extract individual factors fromF
, usingF[:L]
. However, since pivoting is on by default, the factorization is internally represented asA==P'*L*L'*P
with a permutation matrixP
; using justL
without accounting forP
will give incorrect answers. To include the effects of permutation, it’s typically preferable to extact “combined” factors likePtL=F[:PtL]
(the equivalent ofP'*L
) andLtP=F[:UP]
(the equivalent ofL'*P
).Setting optional
shift
keyword argument computes the factorization ofA+shift*I
instead ofA
. If theperm
argument is nonempty, it should be a permutation of1:size(A,1)
giving the ordering to use (instead of CHOLMOD’s default AMD ordering).The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
cholfact!
(A [,LU=:U [,pivot=Val{false}]][;tol=-1.0]) → Cholesky¶cholfact!
is the same ascholfact()
, but saves space by overwriting the inputA
, instead of creating a copy.cholfact!
can also reuse the symbolic factorization from a different matrixF
with the same structure when used as:cholfact!(F::CholmodFactor,A)
.
ldltfact
(::SymTridiagonal) → LDLt¶Compute an
LDLt
factorization of a real symmetric tridiagonal matrix such thatA=L*Diagonal(d)*L'
whereL
is a unit lower triangular matrix andd
is a vector. The main use of anLDLt
factorizationF=ldltfact(A)
is to solve the linear system of equationsAx=b
withF\b
.
ldltfact
(::Union{SparseMatrixCSC, Symmetric{Float64, SparseMatrixCSC{Flaot64, SuiteSparse_long}}, Hermitian{Complex{Float64}, SparseMatrixCSC{Complex{Float64}, SuiteSparse_long}}}; shift=0, perm=Int[]) → CHOLMOD.FactorCompute the
LDLt
factorization of a sparse symmetric or Hermitian matrix. A fill-reducing permutation is used.F=ldltfact(A)
is most frequently used to solve systems of equationsA*x=b
withF\b
, but also the methodsdiag
,det
,logdet
are defined forF
. You can also extract individual factors fromF
, usingF[:L]
. However, since pivoting is on by default, the factorization is internally represented asA==P'*L*D*L'*P
with a permutation matrixP
; using justL
without accounting forP
will give incorrect answers. To include the effects of permutation, it’s typically preferable to extact “combined” factors likePtL=F[:PtL]
(the equivalent ofP'*L
) andLtP=F[:UP]
(the equivalent ofL'*P
). The complete list of supported factors is:L,:PtL,:D,:UP,:U,:LD,:DU,:PtLD,:DUP
.Setting optional
shift
keyword argument computes the factorization ofA+shift*I
instead ofA
. If theperm
argument is nonempty, it should be a permutation of1:size(A,1)
giving the ordering to use (instead of CHOLMOD’s default AMD ordering).The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
ldltfact!
(::SymTridiagonal) → LDLt¶Same as
ldltfact
, but saves space by overwriting the inputA
, instead of creating a copy.
qr
(A[, pivot=Val{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=Val{false}]) → F¶Computes the QR factorization of
A
. The return type ofF
depends on the element type ofA
and whether pivoting is specified (withpivot==Val{true}
).Return type eltype(A)
pivot
Relationship between F
andA
QR
not BlasFloat
either A==F[:Q]*F[:R]
QRCompactWY
BlasFloat
Val{false}
A==F[:Q]*F[:R]
QRPivoted
BlasFloat
Val{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
.Note
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] (1, 2) 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) → SPQR.FactorizationCompute the QR factorization of a sparse matrix
A
. A fill-reducing permutation is used. The main application of this type is to solve least squares problems with\
. The function calls the C library SPQR and a few additional functions from the library are wrapped but not exported.
qrfact!
(A[, pivot=Val{false}])¶qrfact!
is the same asqrfact()
whenA
is a subtype ofStridedMatrix
, but saves space by overwriting the inputA
, instead of creating a copy.
full
(QRCompactWYQ[, thin=true]) → MatrixConverts an orthogonal or unitary matrix stored as a
QRCompactWYQ
object, i.e. in the compact WY format [Bischof1987], to a dense matrix.Optionally takes a
thin
Boolean argument, which iftrue
omits the columns that span the rows ofR
in the QR factorization that are zero. The resulting matrix is theQ
in a thin QR factorization (sometimes called the reduced QR factorization). Iffalse
, returns aQ
that spans all rows ofR
in its corresponding QR factorization.
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.
eig
(A,[irange,][vl,][vu,][permute=true,][scale=true]) → D, V¶Computes eigenvalues and eigenvectors of
A
. Seeeigfact()
for details on thebalance
keyword argument.julia>eig([1.00.00.0;0.03.00.0;0.00.018.0])([1.0,3.0,18.0],3x3Array{Float64,2}:1.00.00.00.01.00.00.00.01.0)
eig
is a wrapper aroundeigfact()
, extracting all parts of the factorization to a tuple; where possible, usingeigfact()
is recommended.
eig
(A, B) → D, VComputes 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]) → values¶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 moreequal in norm. The default istrue
for both options.
eigvals!
(A,[irange,][vl,][vu]) → values¶Same as
eigvals
, but saves space by overwriting the inputA
(andB
), instead of creating a copy.
eigmax
(A)¶Returns the largest eigenvalue of
A
.
eigmin
(A)¶Returns the smallest eigenvalue of
A
.
eigvecs
(A, [eigvals,][permute=true,][scale=true]) → Matrix¶Returns a matrix
M
whose columns are the eigenvectors ofA
. (Thek
th eigenvector can be obtained from the sliceM[:,k]
.) 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,[irange,][vl,][vu,][permute=true,][scale=true]) → Eigen¶Computes the eigenvalue decomposition of
A
, returning anEigen
factorization objectF
which contains the eigenvalues inF[:values]
and the eigenvectors in the columns of the matrixF[:vectors]
. (Thek
th eigenvector can be obtained from the sliceF[:vectors][:,k]
.)The following functions are available for
Eigen
objects:inv
,det
.If
A
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 nonsymmetric 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) → GeneralizedEigenComputes the generalized eigenvalue decomposition of
A
andB
, returning aGeneralizedEigen
factorization objectF
which contains the generalized eigenvalues inF[:values]
and the generalized eigenvectors in the columns of the matrixF[:vectors]
. (Thek
th generalized eigenvector can be obtained from the sliceF[:vectors][:,k]
.)
eigfact!
(A[, B])¶Same as
eigfact()
, but saves space by overwriting the inputA
(andB
), 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 inputA
, 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)¶Computes the Schur factorization of
A
, overwritingA
in the process. Seeschurfact()
schur
(A) → Schur[:T], Schur[:Z], Schur[:values]¶See
schurfact()
ordschur
(Q, T, select) → Schur¶Reorders the Schur factorization of a real matrix
A=Q*T*Q'
according to the logical arrayselect
returning a Schur objectF
. The selected eigenvalues appear in the leading diagonal ofF[:Schur]
and the corresponding leading columns ofF[:vectors]
form an orthonormal basis of the corresponding right invariant subspace. A complex conjugate pair of eigenvalues must be either both included or excluded viaselect
.
ordschur!
(Q, T, select) → Schur¶Reorders the Schur factorization of a real matrix
A=Q*T*Q'
, overwritingQ
andT
in the process. Seeordschur()
ordschur
(S, select) → SchurReorders the Schur factorization
S
of typeSchur
.
ordschur!
(S, select) → SchurReorders the Schur factorization
S
of typeSchur
, overwritingS
in the process. Seeordschur()
schurfact
(A, B) → GeneralizedSchurComputes 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()
ordschur
(S, T, Q, Z, select) → GeneralizedSchurReorders the Generalized Schur factorization of a matrix
(A,B)=(Q*S*Z^{H},Q*T*Z^{H})
according to the logical arrayselect
and returns a GeneralizedSchur objectGS
. The selected eigenvalues appear in the leading diagonal of both(GS[:S],GS[:T])
and the left and right unitary/orthogonal Schur vectors are also reordered such that(A,B)=GS[:Q]*(GS[:S],GS[:T])*GS[:Z]^{H}
still holds and the generalized eigenvalues ofA
andB
can still be obtained withGS[:alpha]./GS[:beta]
.
ordschur!
(S, T, Q, Z, select) → GeneralizedSchurReorders the Generalized Schur factorization of a matrix by overwriting the matrices
(S,T,Q,Z)
in the process. Seeordschur()
.
ordschur
(GS, select) → GeneralizedSchurReorders the Generalized Schur factorization of a Generalized Schur object. See
ordschur()
.
ordschur!
(GS, select) → GeneralizedSchurReorders the Generalized Schur factorization of a Generalized Schur object by overwriting the object with the new factorization. See
ordschur()
.
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 inputA
, 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 ofA
, 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) → GeneralizedSVDCompute 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, R0Wrapper 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
.
givens{T}
(::T, ::T, ::Integer, ::Integer) → {Givens, T}¶Computes the tuple
(G,r)=givens(f,g,i1,i2)
whereG
is a Givens rotation andr
is a scalar such thatG*x=y
withx[i1]=f
,x[i2]=g
,y[i1]=r
, andy[i2]=0
. The cosine and sine of the rotation angle can be extracted from theGivens
type withG.c
andG.s
respectively. The argumentsf
andg
can be eitherFloat32
,Float64
,Complex{Float32}
, orComplex{Float64}
. TheGivens
type supports left multiplicationG*A
and conjugated transpose right multiplicationA*G'
. The type doesn’t have asize
and can therefore be multiplied with matrices of arbitrary size as long asi2<=size(A,2)
forG*A
ori2<=size(A,1)
forA*G'
.
givens{T}
(::AbstractArray{T}, ::Integer, ::Integer, ::Integer) → {Givens, T}Computes the tuple
(G,r)=givens(A,i1,i2,col)
whereG
is Givens rotation andr
is a scalar such thatG*A[:,col]=y
withy[i1]=r
, andy[i2]=0
. The cosine and sine of the rotation angle can be extracted from theGivens
type withG.c
andG.s
respectively. The element type ofA
can be eitherFloat32
,Float64
,Complex{Float32}
, orComplex{Float64}
. TheGivens
type supports left multiplicationG*A
and conjugated transpose right multiplicationA*G'
. The type doesn’t have asize
and can therefore be multiplied with matrices of arbitrary size as long asi2<=size(A,2)
forG*A
ori2<=size(A,1)
forA*G'
.
triu
(M)¶Upper triangle of a matrix.
triu
(M, k)Returns the upper triangle of
M
starting from thek
th superdiagonal.
triu!
(M)¶Upper triangle of a matrix, overwriting
M
in the process.
triu!
(M, k)Returns the upper triangle of
M
starting from thek
th superdiagonal, overwritingM
in the process.
tril
(M)¶Lower triangle of a matrix.
tril
(M, k)Returns the lower triangle of
M
starting from thek
th superdiagonal.
tril!
(M)¶Lower triangle of a matrix, overwriting
M
in the process.
tril!
(M, k)Returns the lower triangle of
M
starting from thek
th superdiagonal, overwritingM
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()
.
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, the matrix norm induced by the vector
p
-norm is used, where valid values ofp
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 (or any element type for whichnorm
is defined), 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.
logabsdet
(M)¶Log of absolute value of determinant of real matrix. Equivalent to
(log(abs(det(M))),sign(det(M)))
, but may provide increased accuracy and/or speed.
inv
(M)¶Matrix inverse
pinv
(M[, tol])¶Computes the Moore-Penrose pseudoinverse.
For matrices
M
with floating point elements, it is convenient to compute the pseudoinverse by inverting only singular values above a given threshold,tol
.The optimal choice of
tol
varies both with the value ofM
and the intended application of the pseudoinverse. The default value oftol
iseps(real(float(one(eltype(M)))))*maximum(size(A))
, which is essentially machine epsilon for the real part of a matrix element multiplied by the larger matrix dimension. For inverting dense ill-conditioned matrices in a least-squares sense,tol=sqrt(eps(real(float(one(eltype(M))))))
is recommended.For more information, see [issue8859], [B96], [S84], [KY88].
[issue8859] Issue 8859, “Fix least squares”, https://github.com/JuliaLang/julia/pull/8859 [B96] Åke Björck, “Numerical Methods for Least Squares Problems”, SIAM Press, Philadelphia, 1996, “Other Titles in Applied Mathematics”, Vol. 51. doi:10.1137/1.9781611971484 [S84] G. W. Stewart, “Rank Degeneracy”, SIAM Journal on Scientific and Statistical Computing, 5(2), 1984, 403-413. doi:10.1137/0905030 [KY88] Konstantinos Konstantinides and Kung Yao, “Statistical analysis of effective singular values in matrix rank determination”, IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. doi:10.1109/29.1585
nullspace
(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¶Perform linear regression. Returns
a
andb
such thata+b*x
is the closest straight line to the given points(x,y)
, i.e., such that the squared error betweeny
anda+b*x
is minimized.Example:
usingPyPlotx=[1.0:12.0;]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 regressionplot(x,y,"o")# Plot (x, y) pointsplot(x,[a+b*iforiinx])# Plot line determined by linear regression
linreg
(x, y, w)Weighted least-squares linear regression.
expm
(A)¶Compute the matrix exponential of
A
, defined by\[e^A = \sum_{n=0}^{\infty} \frac{A^n}{n!}.\]For symmetric or Hermitian
A
, an eigendecomposition (eigfact()
) is used, otherwise the scaling and squaring algorithm (see [H05]) is chosen.[H05] Nicholas J. Higham, “The squaring and scaling method for the matrix exponential revisited”, SIAM Journal on Matrix Analysis and Applications, 26(4), 2005, 1179-1193. doi:10.1137/090768539
logm
(A)¶If
A
has no negative real eigenvalue, compute the principal matrix logarithm ofA
, i.e. the unique matrix \(X\) such that \(e^X = A\) and \(-\pi < Im(\lambda) < \pi\) for all the eigenvalues \(\lambda\) of \(X\). IfA
has nonpositive eigenvalues, a warning is printed and whenever possible a nonprincipal matrix function is returned.If
A
is symmetric or Hermitian, its eigendecomposition (eigfact()
) is used, ifA
is triangular an improved version of the inverse scaling and squaring method is employed (see [AH12] and [AHR13]). For general matrices, the complex Schur form (schur()
) is computed and the triangular algorithm is used on the triangular factor.[AH12] Awad H. Al-Mohy and Nicholas J. Higham, “Improved inverse scaling and squaring algorithms for the matrix logarithm”, SIAM Journal on Scientific Computing, 34(4), 2012, C153-C169. doi:10.1137/110852553 [AHR13] Awad H. Al-Mohy, Nicholas J. Higham and Samuel D. Relton, “Computing the Fréchet derivative of the matrix logarithm and estimating the condition number”, SIAM Journal on Scientific Computing, 35(4), 2013, C394-C410. doi:10.1137/120885991
sqrtm
(A)¶If
A
has no negative real eigenvalues, compute the principal matrix square root ofA
, that is the unique matrix \(X\) with eigenvalues having positive real part such that \(X^2 = A\). Otherwise, a nonprincipal square root is returned.If
A
is symmetric or Hermitian, its eigendecomposition (eigfact()
) is used to compute the square root. Otherwise, the square root is determined by means of the Björck-Hammarling method, which computes the complex Schur form (schur()
) and then the complex square root of the triangular factor.[BH83] Åke Björck and Sven Hammarling, “A Schur method for the square root of a matrix”, Linear Algebra and its Applications, 52-53, 1983, 127-140. doi:10.1016/0024-3795(83)80010-X
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.
isdiag
(A) → Bool¶Test whether a matrix is diagonal.
ishermitian
(A) → Bool¶Test whether a matrix is Hermitian.
transpose
(A)¶The transposition operator (
.'
).
transpose!
(dest, src)¶Transpose array
src
and store the result in the preallocated arraydest
, which should have a size corresponding to(size(src,2),size(src,1))
. No in-place transposition is supported and unexpected results will happen ifsrc
anddest
have overlapping memory regions.
ctranspose
(A)¶The conjugate transposition operator (
'
).
ctranspose!
(dest, src)¶Conjugate transpose array
src
and store the result in the preallocated arraydest
, which should have a size corresponding to(size(src,2),size(src,1))
. No in-place transposition is supported and unexpected results will happen ifsrc
anddest
have overlapping memory regions.
eigs
(A; nev=6, ncv=max(20, 2*nev+1), which="LM", tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0, ))) -> (d, [v, ]nconv, niter, nmult, resid)¶Computes eigenvalues
d
ofA
using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively.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, wheren
is the size of the input matrixA
. The default isncv=max(20,2*nev+1)
. Note that these restrictions limit the input matrixA
to be of dimension at least 2.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
.Note
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}\)
eigs
(A, B; nev=6, ncv=max(20, 2*nev+1), which="LM", tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0, ))) -> (d, [v, ]nconv, niter, nmult, resid)Computes generalized eigenvalues
d
ofA
andB
using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively.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, wheren
is the size of the input matricesA
andB
. The default isncv=max(20,2*nev+1)
. Note that these restrictions limit the input matrixA
to be of dimension at least 2.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
.Note
The
sigma
andwhich
keywords interact: the description of eigenvalues searched for bywhich
do _not_ necessarily refer to the eigenvalue problem \(Av = Bv\lambda\), but rather the linear operator constructed by the specification of the iteration mode implied bysigma
.sigma
iteration mode which
refers to the problemnothing
ordinary (forward) \(Av = Bv\lambda\) real or complex inverse with level shift sigma
\((A - \sigma B )^{-1}B = v\nu\)
svds
(A; nsv=6, ritzvec=true, tol=0.0, maxiter=1000) -> (left_sv, s, right_sv, nconv, niter, nmult, resid)¶svds
computes largest singular valuess
ofA
using Lanczos or Arnoldi iterations. Useseigs()
underneath.Inputs are:
A
: Linear operator. It can either subtype ofAbstractArray
(e.g., sparse matrix) or duck typed. For duck typingA
has to supportsize(A)
,eltype(A)
,A*vector
andA'*vector
.nsv
: Number of singular values.ritzvec
: Whether to return the left and right singular vectorsleft_sv
andright_sv
, default istrue
. Iffalse
the singular vectors are omitted from the output.tol
: tolerance, seeeigs()
.maxiter
: Maximum number of iterations, seeeigs()
.
Example:
X=sprand(10,5,0.2)svds(X,nsv=2)
peakflops
(n; parallel=false)¶peakflops
computes the peak flop rate of the computer by using double precisionBase.LinAlg.BLAS.gemm!()
. By default, if no arguments are specified, it multiplies a matrix of sizenxn
, 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¶
Base.LinAlg.BLAS
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!
(a, X, Y)¶Overwrite
Y
witha*X+Y
. ReturnsY
.
scal!
(n, a, X, incx)¶Overwrite
X
witha*X
. ReturnsX
.
scal
(n, a, X, incx)¶Returns
a*X
.
ger!
(alpha, x, y, A)¶Rank-1 update of the matrix
A
with vectorsx
andy
asalpha*x*y'+A
.
syr!
(uplo, alpha, x, A)¶Rank-1 update of the symmetric matrix
A
with vectorx
asalpha*x*x.'+A
. Whenuplo
is ‘U’ the upper triangle ofA
is updated (‘L’ for lower triangle). ReturnsA
.
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’).
her!
(uplo, alpha, x, A)¶Methods for complex arrays only. Rank-1 update of the Hermitian matrix
A
with vectorx
asalpha*x*x'+A
. Whenuplo
is ‘U’ the upper triangle ofA
is updated (‘L’ for lower triangle). ReturnsA
.
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*y
oralpha*A'x+beta*y
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*x+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
(transposeA
). 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
(transposeA
). 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
(transposeA
). 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
(transposeA
). 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
(transposeA
). 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
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
trsv!
(ul, tA, dA, A, b)¶Overwrite
b
with the solution toA*x=b
or one of the other two variants determined bytA
(transposeA
) andul
(triangle ofA
used).dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedb
.
trsv
(ul, tA, dA, A, b)¶Returns the solution to
A*x=b
or one of the other two variants determined bytA
(transposeA
) andul
(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.
I
¶An object of type
UniformScaling
, representing an identity matrix of any size.
LAPACK Functions¶
Base.LinAlg.LAPACK
provides wrappers for some of the LAPACK functions for
linear algebra. Those 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.
Note that the LAPACK API provided by Julia can and will change in the future. Since this API is not user-facing, there is no commitment to support/deprecate this specific set of functions in future releases.
gbtrf!
(kl, ku, m, AB) -> (AB, ipiv)¶Compute the LU factorization of a banded matrix
AB
.kl
is the first subdiagonal containing a nonzero band,ku
is the last superdiagonal containing one, andm
is the first dimension of the matrixAB
. Returns the LU factorization in-place andipiv
, the vector of pivots used.
gbtrs!
(trans, kl, ku, m, AB, ipiv, B)¶Solve the equation
AB*X=B
.trans
determines the orientation ofAB
. It may beN
(no transpose),T
(transpose), orC
(conjugate transpose).kl
is the first subdiagonal containing a nonzero band,ku
is the last superdiagonal containing one, andm
is the first dimension of the matrixAB
.ipiv
is the vector of pivots returned fromgbtrf!
. Returns the vector or matrixX
, overwritingB
in-place.
gebal!
(job, A) -> (ilo, ihi, scale)¶Balance the matrix
A
before computing its eigensystem or Schur factorization.job
can be one ofN
(A
will not be permuted or scaled),P
(A
will only be permuted),S
(A
will only be scaled), orB
(A
will be both permuted and scaled). ModifiesA
in-place and returnsilo
,ihi
, andscale
. If permuting was turned on,A[i,j]=0
ifj>i
and1<j<ilo
orj>ihi
.scale
contains information about the scaling/permutations performed.
gebak!
(job, side, ilo, ihi, scale, V)¶Transform the eigenvectors
V
of a matrix balanced usinggebal!
to the unscaled/unpermuted eigenvectors of the original matrix. ModifiesV
in-place.side
can beL
(left eigenvectors are transformed) orR
(right eigenvectors are transformed).
gebrd!
(A) -> (A, d, e, tauq, taup)¶Reduce
A
in-place to bidiagonal formA=QBP'
. ReturnsA
, containing the bidiagonal matrixB
;d
, containing the diagonal elements ofB
;e
, containing the off-diagonal elements ofB
;tauq
, containing the elementary reflectors representingQ
; andtaup
, containing the elementary reflectors representingP
.
gelqf!
(A, tau)¶Compute the
LQ
factorization ofA
,A=LQ
.tau
contains scalars which parameterize the elementary reflectors of the factorization.tau
must have length greater than or equal to the smallest dimension ofA
.Returns
A
andtau
modified in-place.
gelqf!
(A) -> (A, tau)Compute the
LQ
factorization ofA
,A=LQ
.Returns
A
, modified in-place, andtau
, which contains scalars which parameterize the elementary reflectors of the factorization.
geqlf!
(A, tau)¶Compute the
QL
factorization ofA
,A=QL
.tau
contains scalars which parameterize the elementary reflectors of the factorization.tau
must have length greater than or equal to the smallest dimension ofA
.Returns
A
andtau
modified in-place.
geqlf!
(A) -> (A, tau)Compute the
QL
factorization ofA
,A=QL
.Returns
A
, modified in-place, andtau
, which contains scalars which parameterize the elementary reflectors of the factorization.
geqrf!
(A, tau)¶Compute the
QR
factorization ofA
,A=QR
.tau
contains scalars which parameterize the elementary reflectors of the factorization.tau
must have length greater than or equal to the smallest dimension ofA
.Returns
A
andtau
modified in-place.
geqrf!
(A) -> (A, tau)Compute the
QR
factorization ofA
,A=QR
.Returns
A
, modified in-place, andtau
, which contains scalars which parameterize the elementary reflectors of the factorization.
geqp3!
(A, jpvt, tau)¶Compute the pivoted
QR
factorization ofA
,AP=QR
using BLAS level 3.P
is a pivoting matrix, represented byjpvt
.tau
stores the elementary reflectors.jpvt
must have length length greater than or equal ton
ifA
is an(mxn)
matrix.tau
must have length greater than or equal to the smallest dimension ofA
.A
,jpvt
, andtau
are modified in-place.
geqp3!
(A, jpvt) -> (A, jpvt, tau)Compute the pivoted
QR
factorization ofA
,AP=QR
using BLAS level 3.P
is a pivoting matrix, represented byjpvt
.jpvt
must have length greater than or equal ton
ifA
is an(mxn)
matrix.Returns
A
andjpvt
, modified in-place, andtau
, which stores the elementary reflectors.
geqp3!
(A) -> (A, jpvt, tau)Compute the pivoted
QR
factorization ofA
,AP=QR
using BLAS level 3.Returns
A
, modified in-place,jpvt
, which represents the pivoting matrixP
, andtau
, which stores the elementary reflectors.
gerqf!
(A, tau)¶Compute the
RQ
factorization ofA
,A=RQ
.tau
contains scalars which parameterize the elementary reflectors of the factorization.tau
must have length greater than or equal to the smallest dimension ofA
.Returns
A
andtau
modified in-place.
gerqf!
(A) -> (A, tau)Compute the
RQ
factorization ofA
,A=RQ
.Returns
A
, modified in-place, andtau
, which contains scalars which parameterize the elementary reflectors of the factorization.
geqrt!
(A, T)¶Compute the blocked
QR
factorization ofA
,A=QR
.T
contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization. The first dimension ofT
sets the block size and it must be between 1 andn
. The second dimension ofT
must equal the smallest dimension ofA
.Returns
A
andT
modified in-place.
geqrt!
(A, nb) -> (A, T)Compute the blocked
QR
factorization ofA
,A=QR
.nb
sets the block size and it must be between 1 andn
, the second dimension ofA
.Returns
A
, modified in-place, andT
, which contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization.
geqrt3!
(A, T)¶Recursively computes the blocked
QR
factorization ofA
,A=QR
.T
contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization. The first dimension ofT
sets the block size and it must be between 1 andn
. The second dimension ofT
must equal the smallest dimension ofA
.Returns
A
andT
modified in-place.
geqrt3!
(A) -> (A, T)Recursively computes the blocked
QR
factorization ofA
,A=QR
.Returns
A
, modified in-place, andT
, which contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization.
getrf!
(A) -> (A, ipiv, info)¶Compute the pivoted
LU
factorization ofA
,A=LU
.Returns
A
, modified in-place,ipiv
, the pivoting information, and aninfo
code which indicates success (info=0
), a singular value inU
(info=i
, in which caseU[i,i]
is singular), or an error code (info<0
).
tzrzf!
(A) -> (A, tau)¶Transforms the upper trapezoidal matrix
A
to upper triangular form in-place. ReturnsA
andtau
, the scalar parameters for the elementary reflectors of the transformation.
ormrz!
(side, trans, A, tau, C)¶Multiplies the matrix
C
byQ
from the transformation supplied bytzrzf!
. Depending onside
ortrans
the multiplication can be left-sided (side=L,Q*C
) or right-sided (side=R,C*Q
) andQ
can be unmodified (trans=N
), transposed (trans=T
), or conjugate transposed (trans=C
). Returns matrixC
which is modified in-place with the result of the multiplication.
gels!
(trans, A, B) -> (F, B, ssr)¶Solves the linear equation
A*X=B
,A.'*X=B
, orA'*X=B
using a QR or LQ factorization. Modifies the matrix/vectorB
in place with the solution.A
is overwritten with itsQR
orLQ
factorization.trans
may be one ofN
(no modification),T
(transpose), orC
(conjugate transpose).gels!
searches for the minimum norm/least squares solution.A
may be under or over determined. The solution is returned inB
.
gesv!
(A, B) -> (B, A, ipiv)¶Solves the linear equation
A*X=B
whereA
is a square matrix using theLU
factorization ofA
.A
is overwritten with itsLU
factorization andB
is overwritten with the solutionX
.ipiv
contains the pivoting information for theLU
factorization ofA
.
getrs!
(trans, A, ipiv, B)¶Solves the linear equation
A*X=B
,A.'*X=B
, orA'*X=B
for squareA
. Modifies the matrix/vectorB
in place with the solution.A
is theLU
factorization fromgetrf!
, withipiv
the pivoting information.trans
may be one ofN
(no modification),T
(transpose), orC
(conjugate transpose).
getri!
(A, ipiv)¶Computes the inverse of
A
, using itsLU
factorization found bygetrf!
.ipiv
is the pivot information output andA
contains theLU
factorization ofgetrf!
.A
is overwritten with its inverse.
gesvx!
(fact, trans, A, AF, ipiv, equed, R, C, B) -> (X, equed, R, C, B, rcond, ferr, berr, work)¶Solves the linear equation
A*X=B
(trans=N
),A.'*X=B
(trans=T
), orA'*X=B
(trans=C
) using theLU
factorization ofA
.fact
may beE
, in which caseA
will be equilibrated and copied toAF
;F
, in which caseAF
andipiv
from a previousLU
factorization are inputs; orN
, in which caseA
will be copied toAF
and then factored. Iffact=F
,equed
may beN
, meaningA
has not been equilibrated;R
, meaningA
was multiplied bydiagm(R)
from the left;C
, meaningA
was multiplied bydiagm(C)
from the right; orB
, meaningA
was multiplied bydiagm(R)
from the left anddiagm(C)
from the right. Iffact=F
andequed=R
orB
the elements ofR
must all be positive. Iffact=F
andequed=C
orB
the elements ofC
must all be positive.Returns the solution
X
;equed
, which is an output iffact
is notN
, and describes the equilibration that was performed;R
, the row equilibration diagonal;C
, the column equilibration diagonal;B
, which may be overwritten with its equilibrated formdiagm(R)*B
(iftrans=N
andequed=R,B
) ordiagm(C)*B
(iftrans=T,C
andequed=C,B
);rcond
, the reciprocal condition number ofA
after equilbrating;ferr
, the forward error bound for each solution vector inX
;berr
, the forward error bound for each solution vector inX
; andwork
, the reciprocal pivot growth factor.
gesvx!
(A, B)The no-equilibration, no-transpose simplification of
gesvx!
.
gelsd!
(A, B, rcond) -> (B, rnk)¶Computes the least norm solution of
A*X=B
by finding theSVD
factorization ofA
, then dividing-and-conquering the problem.B
is overwritten with the solutionX
. Singular values belowrcond
will be treated as zero. Returns the solution inB
and the effective rank ofA
inrnk
.
gelsy!
(A, B, rcond) -> (B, rnk)¶Computes the least norm solution of
A*X=B
by finding the fullQR
factorization ofA
, then dividing-and-conquering the problem.B
is overwritten with the solutionX
. Singular values belowrcond
will be treated as zero. Returns the solution inB
and the effective rank ofA
inrnk
.
gglse!
(A, c, B, d) -> (X, res)¶Solves the equation
A*x=c
wherex
is subject to the equality constraintB*x=d
. Uses the formula||c-A*x||^2=0
to solve. ReturnsX
and the residual sum-of-squares.
geev!
(jobvl, jobvr, A) -> (W, VL, VR)¶Finds the eigensystem of
A
. Ifjobvl=N
, the left eigenvectors ofA
aren’t computed. Ifjobvr=N
, the right eigenvectors ofA
aren’t computed. Ifjobvl=V
orjobvr=V
, the corresponding eigenvectors are computed. Returns the eigenvalues inW
, the right eigenvectors inVR
, and the left eigenvectors inVL
.
gesdd!
(job, A) -> (U, S, VT)¶Finds the singular value decomposition of
A
,A=U*S*V'
, using a divide and conquer approach. Ifjob=A
, all the columns ofU
and the rows ofV'
are computed. Ifjob=N
, no columns ofU
or rows ofV'
are computed. Ifjob=O
,A
is overwritten with the columns of (thin)U
and the rows of (thin)V'
. Ifjob=S
, the columns of (thin)U
and the rows of (thin)V'
are computed and returned separately.
gesvd!
(jobu, jobvt, A) -> (U, S, VT)¶Finds the singular value decomposition of
A
,A=U*S*V'
. Ifjobu=A
, all the columns ofU
are computed. Ifjobvt=A
all the rows ofV'
are computed. Ifjobu=N
, no columns ofU
are computed. Ifjobvt=N
no rows ofV'
are computed. Ifjobu=O
,A
is overwritten with the columns of (thin)U
. Ifjobvt=O
,A
is overwritten with the rows of (thin)V'
. Ifjobu=S
, the columns of (thin)U
are computed and returned separately. Ifjobvt=S
the rows of (thin)V'
are computed and returned separately.jobu
andjobvt
can’t both beO
.Returns
U
,S
, andVt
, whereS
are the singular values ofA
.
ggsvd!
(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)¶Finds the generalized singular value decomposition of
A
andB
,U'*A*Q=D1*R
andV'*B*Q=D2*R
.D1
hasalpha
on its diagonal andD2
hasbeta
on its diagonal. Ifjobu=U
, the orthogonal/unitary matrixU
is computed. Ifjobv=V
the orthogonal/unitary matrixV
is computed. Ifjobq=Q
, the orthogonal/unitary matrixQ
is computed. Ifjobu
,jobv
orjobq
isN
, that matrix is not computed. This function is only available in LAPACK versions prior to 3.6.0.
ggsvd3!
(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)¶Finds the generalized singular value decomposition of
A
andB
,U'*A*Q=D1*R
andV'*B*Q=D2*R
.D1
hasalpha
on its diagonal andD2
hasbeta
on its diagonal. Ifjobu=U
, the orthogonal/unitary matrixU
is computed. Ifjobv=V
the orthogonal/unitary matrixV
is computed. Ifjobq=Q
, the orthogonal/unitary matrixQ
is computed. Ifjobu
,jobv
, orjobq
isN
, that matrix is not computed. This function requires LAPACK 3.6.0.
geevx!
(balanc, jobvl, jobvr, sense, A) -> (A, w, VL, VR, ilo, ihi, scale, abnrm, rconde, rcondv)¶Finds the eigensystem of
A
with matrix balancing. Ifjobvl=N
, the left eigenvectors ofA
aren’t computed. Ifjobvr=N
, the right eigenvectors ofA
aren’t computed. Ifjobvl=V
orjobvr=V
, the corresponding eigenvectors are computed. Ifbalanc=N
, no balancing is performed. Ifbalanc=P
,A
is permuted but not scaled. Ifbalanc=S
,A
is scaled but not permuted. Ifbalanc=B
,A
is permuted and scaled. Ifsense=N
, no reciprocal condition numbers are computed. Ifsense=E
, reciprocal condition numbers are computed for the eigenvalues only. Ifsense=V
, reciprocal condition numbers are computed for the right eigenvectors only. Ifsense=B
, reciprocal condition numbers are computed for the right eigenvectors and the eigenvectors. Ifsense=E,B
, the right and left eigenvectors must be computed.
ggev!
(jobvl, jobvr, A, B) -> (alpha, beta, vl, vr)¶Finds the generalized eigendecomposition of
A
andB
. Ifjobvl=N
, the left eigenvectors aren’t computed. Ifjobvr=N
, the right eigenvectors aren’t computed. Ifjobvl=V
orjobvr=V
, the corresponding eigenvectors are computed.
gtsv!
(dl, d, du, B)¶Solves the equation
A*X=B
whereA
is a tridiagonal matrix withdl
on the subdiagonal,d
on the diagonal, anddu
on the superdiagonal.Overwrites
B
with the solutionX
and returns it.
gttrf!
(dl, d, du) -> (dl, d, du, du2, ipiv)¶Finds the
LU
factorization of a tridiagonal matrix withdl
on the subdiagonal,d
on the diagonal, anddu
on the superdiagonal.Modifies
dl
,d
, anddu
in-place and returns them and the second superdiagonaldu2
and the pivoting vectoripiv
.
gttrs!
(trans, dl, d, du, du2, ipiv, B)¶Solves the equation
A*X=B
(trans=N
),A.'*X=B
(trans=T
), orA'*X=B
(trans=C
) using theLU
factorization computed bygttrf!
.B
is overwritten with the solutionX
.
orglq!
(A, tau, k = length(tau))¶Explicitly finds the matrix
Q
of aLQ
factorization after callinggelqf!
onA
. Uses the output ofgelqf!
.A
is overwritten byQ
.
orgqr!
(A, tau, k = length(tau))¶Explicitly finds the matrix
Q
of aQR
factorization after callinggeqrf!
onA
. Uses the output ofgeqrf!
.A
is overwritten byQ
.
ormlq!
(side, trans, A, tau, C)¶Computes
Q*C
(trans=N
),Q.'*C
(trans=T
),Q'*C
(trans=C
) forside=L
or the equivalent right-sided multiplication forside=R
usingQ
from aLQ
factorization ofA
computed usinggelqf!
.C
is overwritten.
ormqr!
(side, trans, A, tau, C)¶Computes
Q*C
(trans=N
),Q.'*C
(trans=T
),Q'*C
(trans=C
) forside=L
or the equivalent right-sided multiplication forside=R
usingQ
from aQR
factorization ofA
computed usinggeqrf!
.C
is overwritten.
gemqrt!
(side, trans, V, T, C)¶Computes
Q*C
(trans=N
),Q.'*C
(trans=T
),Q'*C
(trans=C
) forside=L
or the equivalent right-sided multiplication forside=R
usingQ
from aQR
factorization ofA
computed usinggeqrt!
.C
is overwritten.
posv!
(uplo, A, B) -> (A, B)¶Finds the solution to
A*X=B
whereA
is a symmetric or Hermitian positive definite matrix. Ifuplo=U
the upper Cholesky decomposition ofA
is computed. Ifuplo=L
the lower Cholesky decomposition ofA
is computed.A
is overwritten by its Cholesky decomposition.B
is overwritten with the solutionX
.
potrf!
(uplo, A)¶Computes the Cholesky (upper if
uplo=U
, lower ifuplo=L
) decomposition of positive-definite matrixA
.A
is overwritten and returned with an info code.
potri!
(uplo, A)¶Computes the inverse of positive-definite matrix
A
after callingpotrf!
to find its (upper ifuplo=U
, lower ifuplo=L
) Cholesky decomposition.A
is overwritten by its inverse and returned.
potrs!
(uplo, A, B)¶Finds the solution to
A*X=B
whereA
is a symmetric or Hermitian positive definite matrix whose Cholesky decomposition was computed bypotrf!
. Ifuplo=U
the upper Cholesky decomposition ofA
was computed. Ifuplo=L
the lower Cholesky decomposition ofA
was computed.B
is overwritten with the solutionX
.
pstrf!
(uplo, A, tol) -> (A, piv, rank, info)¶Computes the (upper if
uplo=U
, lower ifuplo=L
) pivoted Cholesky decomposition of positive-definite matrixA
with a user-set tolerancetol
.A
is overwritten by its Cholesky decomposition.Returns
A
, the pivotspiv
, the rank ofA
, and aninfo
code. Ifinfo=0
, the factorization succeeded. Ifinfo=i>0`,then`A
is indefinite or rank-deficient.
ptsv!
(D, E, B)¶Solves
A*X=B
for positive-definite tridiagonalA
.D
is the diagonal ofA
andE
is the off-diagonal.B
is overwritten with the solutionX
and returned.
pttrf!
(D, E)¶Computes the LDLt factorization of a positive-definite tridiagonal matrix with
D
as diagonal andE
as off-diagonal.D
andE
are overwritten and returned.
pttrs!
(D, E, B)¶Solves
A*X=B
for positive-definite tridiagonalA
with diagonalD
and off-diagonalE
after computingA
‘s LDLt factorization usingpttrf!
.B
is overwritten with the solutionX
.
trtri!
(uplo, diag, A)¶Finds the inverse of (upper if
uplo=U
, lower ifuplo=L
) triangular matrixA
. Ifdiag=N
,A
has non-unit diagonal elements. Ifdiag=U
, all diagonal elements ofA
are one.A
is overwritten with its inverse.
trtrs!
(uplo, trans, diag, A, B)¶Solves
A*X=B
(trans=N
),A.'*X=B
(trans=T
), orA'*X=B
(trans=C
) for (upper ifuplo=U
, lower ifuplo=L
) triangular matrixA
. Ifdiag=N
,A
has non-unit diagonal elements. Ifdiag=U
, all diagonal elements ofA
are one.B
is overwritten with the solutionX
.
trcon!
(norm, uplo, diag, A)¶Finds the reciprocal condition number of (upper if
uplo=U
, lower ifuplo=L
) triangular matrixA
. Ifdiag=N
,A
has non-unit diagonal elements. Ifdiag=U
, all diagonal elements ofA
are one. Ifnorm=I
, the condition number is found in the infinity norm. Ifnorm=O
or1
, the condition number is found in the one norm.
trevc!
(side, howmny, select, T, VL = similar(T), VR = similar(T))¶Finds the eigensystem of an upper triangular matrix
T
. Ifside=R
, the right eigenvectors are computed. Ifside=L
, the left eigenvectors are computed. Ifside=B
, both sets are computed. Ifhowmny=A
, all eigenvectors are found. Ifhowmny=B
, all eigenvectors are found and backtransformed usingVL
andVR
. Ifhowmny=S
, only the eigenvectors corresponding to the values inselect
are computed.
trrfs!
(uplo, trans, diag, A, B, X, Ferr, Berr) -> (Ferr, Berr)¶Estimates the error in the solution to
A*X=B
(trans=N
),A.'*X=B
(trans=T
),A'*X=B
(trans=C
) forside=L
, or the equivalent equations a right-handedside=R
X*A
after computingX
usingtrtrs!
. Ifuplo=U
,A
is upper triangular. Ifuplo=L
,A
is lower triangular. Ifdiag=N
,A
has non-unit diagonal elements. Ifdiag=U
, all diagonal elements ofA
are one.Ferr
andBerr
are optional inputs.Ferr
is the forward error andBerr
is the backward error, each component-wise.
stev!
(job, dv, ev) -> (dv, Zmat)¶Computes the eigensystem for a symmetric tridiagonal matrix with
dv
as diagonal andev
as off-diagonal. Ifjob=N
only the eigenvalues are found and returned indv
. Ifjob=V
then the eigenvectors are also found and returned inZmat
.
stebz!
(range, order, vl, vu, il, iu, abstol, dv, ev) -> (dv, iblock, isplit)¶Computes the eigenvalues for a symmetric tridiagonal matrix with
dv
as diagonal andev
as off-diagonal. Ifrange=A
, all the eigenvalues are found. Ifrange=V
, the eigenvalues in the half-open interval(vl,vu]
are found. Ifrange=I
, the eigenvalues with indices betweenil
andiu
are found. Iforder=B
, eigvalues are ordered within a block. Iforder=E
, they are ordered across all the blocks.abstol
can be set as a tolerance for convergence.
stegr!
(jobz, range, dv, ev, vl, vu, il, iu) -> (w, Z)¶Computes the eigenvalues (
jobz=N
) or eigenvalues and eigenvectors (jobz=V
) for a symmetric tridiagonal matrix withdv
as diagonal andev
as off-diagonal. Ifrange=A
, all the eigenvalues are found. Ifrange=V
, the eigenvalues in the half-open interval(vl,vu]
are found. Ifrange=I
, the eigenvalues with indices betweenil
andiu
are found. The eigenvalues are returned inw
and the eigenvectors inZ
.
stein!
(dv, ev_in, w_in, iblock_in, isplit_in)¶Computes the eigenvectors for a symmetric tridiagonal matrix with
dv
as diagonal andev_in
as off-diagonal.w_in
specifies the input eigenvalues for which to find corresponding eigenvectors.iblock_in
specifies the submatrices corresponding to the eigenvalues inw_in
.isplit_in
specifies the splitting points between the submatrix blocks.
syconv!
(uplo, A, ipiv) -> (A, work)¶Converts a symmetric matrix
A
(which has been factorized into a triangular matrix) into two matricesL
andD
. Ifuplo=U
,A
is upper triangular. Ifuplo=L
, it is lower triangular.ipiv
is the pivot vector from the triangular factorization.A
is overwritten byL
andD
.
sysv!
(uplo, A, B) -> (B, A, ipiv)¶Finds the solution to
A*X=B
for symmetric matrixA
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.B
is overwritten by the solutionX
.A
is overwritten by its Bunch-Kaufman factorization.ipiv
contains pivoting information about the factorization.
sytrf!
(uplo, A) -> (A, ipiv)¶Computes the Bunch-Kaufman factorization of a symmetric matrix
A
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.Returns
A
, overwritten by the factorization, and a pivot vectoripiv
.
sytri!
(uplo, A, ipiv)¶Computes the inverse of a symmetric matrix
A
using the results ofsytrf!
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.A
is overwritten by its inverse.
sytrs!
(uplo, A, ipiv, B)¶Solves the equation
A*X=B
for a symmetric matrixA
using the results ofsytrf!
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.B
is overwritten by the solutionX
.
hesv!
(uplo, A, B) -> (B, A, ipiv)¶Finds the solution to
A*X=B
for Hermitian matrixA
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.B
is overwritten by the solutionX
.A
is overwritten by its Bunch-Kaufman factorization.ipiv
contains pivoting information about the factorization.
hetrf!
(uplo, A) -> (A, ipiv)¶Computes the Bunch-Kaufman factorization of a Hermitian matrix
A
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.Returns
A
, overwritten by the factorization, and a pivot vector.
hetri!
(uplo, A, ipiv)¶Computes the inverse of a Hermitian matrix
A
using the results ofsytrf!
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.A
is overwritten by its inverse.
hetrs!
(uplo, A, ipiv, B)¶Solves the equation
A*X=B
for a Hermitian matrixA
using the results ofsytrf!
. Ifuplo=U
, the upper half ofA
is stored. Ifuplo=L
, the lower half is stored.B
is overwritten by the solutionX
.
syev!
(jobz, uplo, A)¶Finds the eigenvalues (
jobz=N
) or eigenvalues and eigenvectors (jobz=V
) of a symmetric matrixA
. Ifuplo=U
, the upper triangle ofA
is used. Ifuplo=L
, the lower triangle ofA
is used.
syevr!
(jobz, range, uplo, A, vl, vu, il, iu, abstol) -> (W, Z)¶Finds the eigenvalues (
jobz=N
) or eigenvalues and eigenvectors (jobz=V
) of a symmetric matrixA
. Ifuplo=U
, the upper triangle ofA
is used. Ifuplo=L
, the lower triangle ofA
is used. Ifrange=A
, all the eigenvalues are found. Ifrange=V
, the eigenvalues in the half-open interval(vl,vu]
are found. Ifrange=I
, the eigenvalues with indices betweenil
andiu
are found.abstol
can be set as a tolerance for convergence.The eigenvalues are returned in
W
and the eigenvectors inZ
.
sygvd!
(jobz, range, uplo, A, vl, vu, il, iu, abstol) -> (w, A, B)¶Finds the generalized eigenvalues (
jobz=N
) or eigenvalues and eigenvectors (jobz=V
) of a symmetric matrixA
and symmetric positive-definite matrixB
. Ifuplo=U
, the upper triangles ofA
andB
are used. Ifuplo=L
, the lower triangles ofA
andB
are used. Ifitype=1
, the problem to solve isA*x=lambda*B*x
. Ifitype=2
, the problem to solve isA*B*x=lambda*x
. Ifitype=3
, the problem to solve isB*A*x=lambda*x
.
bdsqr!
(uplo, d, e_, Vt, U, C) -> (d, Vt, U, C)¶Computes the singular value decomposition of a bidiagonal matrix with
d
on the diagonal ande_
on the off-diagonal. Ifuplo=U
,e_
is the superdiagonal. Ifuplo=L
,e_
is the subdiagonal. Can optionally also compute the productQ'*C
.Returns the singular values in
d
, and the matrixC
overwritten withQ'*C
.
bdsdc!
(uplo, compq, d, e_) -> (d, e, u, vt, q, iq)¶Computes the singular value decomposition of a bidiagonal matrix with
d
on the diagonal ande_
on the off-diagonal using a divide and conqueq method. Ifuplo=U
,e_
is the superdiagonal. Ifuplo=L
,e_
is the subdiagonal. Ifcompq=N
, only the singular values are found. Ifcompq=I
, the singular values and vectors are found. Ifcompq=P
, the singular values and vectors are found in compact form. Only works for real types.Returns the singular values in
d
, and ifcompq=P
, the compact singular vectors iniq
.
gecon!
(normtype, A, anorm)¶Finds the reciprocal condition number of matrix
A
. Ifnormtype=I
, the condition number is found in the infinity norm. Ifnormtype=O
or1
, the condition number is found in the one norm.A
must be the result ofgetrf!
andanorm
is the norm ofA
in the relevant norm.
gehrd!
(ilo, ihi, A) -> (A, tau)¶Converts a matrix
A
to Hessenberg form. IfA
is balanced withgebal!
thenilo
andihi
are the outputs ofgebal!
. Otherwise they should beilo=1
andihi=size(A,2)
.tau
contains the elementary reflectors of the factorization.
orghr!
(ilo, ihi, A, tau)¶Explicitly finds
Q
, the orthogonal/unitary matrix fromgehrd!
.ilo
,ihi
,A
, andtau
must correspond to the input/output togehrd!
.
gees!
(jobvs, A) -> (A, vs, w)¶Computes the eigenvalues (
jobvs=N
) or the eigenvalues and Schur vectors (jobvs=V
) of matrixA
.A
is overwritten by its Schur form.Returns
A
,vs
containing the Schur vectors, andw
, containing the eigenvalues.
gges!
(jobvsl, jobvsr, A, B) -> (A, B, alpha, beta, vsl, vsr)¶Computes the generalized eigenvalues, generalized Schur form, left Schur vectors (
jobsvl=V
), or right Schur vectors (jobvsr=V
) ofA
andB
.The generalized eigenvalues are returned in
alpha
andbeta
. The left Schur vectors are returned invsl
and the right Schur vectors are returned invsr
.
trsen!
(compq, job, select, T, Q) -> (T, Q, w)¶Reorder the Schur factorization of a matrix and optionally finds reciprocal condition numbers. If
job=N
, no condition numbers are found. Ifjob=E
, only the condition number for this cluster of eigenvalues is found. Ifjob=V
, only the condition number for the invariant subspace is found. Ifjob=B
then the condition numbers for the cluster and subspace are found. Ifcompq=V
the Schur vectorsQ
are updated. Ifcompq=N
the Schur vectors are not modified.select
determines which eigenvalues are in the cluster.Returns
T
,Q
, and reordered eigenvalues inw
.
tgsen!
(select, S, T, Q, Z) -> (S, T, alpha, beta, Q, Z)¶Reorders the vectors of a generalized Schur decomposition.
select
specifices the eigenvalues in each cluster.
trsyl!
(transa, transb, A, B, C, isgn=1) -> (C, scale)¶Solves the Sylvester matrix equation
A*X+/-X*B=scale*C
whereA
andB
are both quasi-upper triangular. Iftransa=N
,A
is not modified. Iftransa=T
,A
is transposed. Iftransa=C
,A
is conjugate transposed. Similarly fortransb
andB
. Ifisgn=1
, the equationA*X+X*B=scale*C
is solved. Ifisgn=-1
, the equationA*X-X*B=scale*C
is solved.Returns
X
(overwritingC
) andscale
.