NoPivot

Pivoting is not performed. This is the default strategy forcholesky andqr factorizations. Note, however, that other matrix factorizations such as the LU factorization may fail without pivoting, and may also be numerically unstable for floating-point matrices in the face of roundoff error. In such cases, this pivot strategy is mainly useful for pedagogical purposes.

RowNonZero

First non-zero element in the remaining rows is chosen as the pivot element.

Beware that for floating-point matrices, the resulting LU algorithm is numerically unstable — this strategy is mainly useful for comparison to hand calculations (which typically use this strategy) or for other algebraic types (e.g. rational numbers) not susceptible to roundoff errors. Otherwise, the defaultRowMaximum pivoting strategy should be generally preferred in Gaussian elimination.

Note that theelement type of the matrix must admit aniszero method.

RowMaximum

A row (and potentially also column) pivot is chosen based on a maximum property. This is the default strategy for LU factorization and for pivoted Cholesky factorization (though [NoPivot] is the default forcholesky).

In the LU case, the maximum-magnitude element within the current column in the remaining rows is chosen as the pivot element. This is sometimes referred to as the "partial pivoting" algorithm. In this case, theelement type of the matrix must admit anabs method, whose result type must admit a< method.

In the Cholesky case, the maximal element among the remaining diagonal elements is chosen as the pivot element. This is sometimes referred to as the "diagonal pivoting" algorithm, and leads tocomplete pivoting (i.e., of both rows and columns by the same permutation). In this case, the (real part of the)element type of the matrix must admit a< method.

ColumnNorm

The column with the maximum norm is used for subsequent computation. This is used for pivoted QR factorization.

Note that theelement type of the matrix must admitnorm andabs methods, whose respective result types must admit a< method.

*(A::AbstractMatrix, B::AbstractMatrix)

Matrix multiplication.

Examples

julia> [1 1; 0 1] * [1 0; 1 1]2×2 Matrix{Int64}: 2  1 1  1

*(A, B::AbstractMatrix, C)A * B * C * D

Chained multiplication of 3 or 4 matrices is done in the most efficient sequence, based on the sizes of the arrays. That is, the number of scalar multiplications needed for(A * B) * C (with 3 dense matrices) is compared to that forA * (B * C) to choose which of these to execute.

If the last factor is a vector, or the first a transposed vector, then it is efficient to deal with these first. In particularx' * B * y means(x' * B) * y for an ordinary column-majorB::Matrix. Unlikedot(x, B, y), this allocates an intermediate array.

If the first or last factor is a number, this will be fused with the matrix multiplication, using 5-argmul!.

See alsomuladd,dot.

\(A, B)

Matrix division using a polyalgorithm. For input matricesA andB, the resultX is such thatA*X == B whenA is square. The solver that is used depends upon the structure ofA. IfA is upper or lower triangular (or diagonal), no factorization ofA is required and the system is solved with either forward or backward substitution. For non-triangular square matrices, 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.

WhenA 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.

See also:factorize,pinv.

Examples

julia> A = [1 0; 1 -2]; B = [32; -4];julia> X = A \ B2-element Vector{Float64}: 32.0 18.0julia> A * X == Btrue

A / B

Matrix right-division:A / B is equivalent to(B' \ A')' where\ is the left-division operator. For square matrices, the resultX is such thatA == X*B.

See also:rdiv!.

Examples

julia> A = Float64[1 4 5; 3 9 2]; B = Float64[1 4 2; 3 4 2; 8 7 1];julia> X = A / B2×3 Matrix{Float64}: -0.65   3.75  -1.2  3.25  -2.75   1.0julia> isapprox(A, X*B)truejulia> isapprox(X, A*pinv(B))true

SingularException

Exception thrown when the input matrix has one or more zero-valued eigenvalues, and is not invertible. A linear solve involving such a matrix cannot be computed. Theinfo field indicates the location of (one of) the singular value(s).

PosDefException

Exception thrown when the input matrix was notpositive definite. Some linear algebra functions and factorizations are only applicable to positive definite matrices. Theinfo field indicates the location of (one of) the eigenvalue(s) which is (are) less than/equal to 0.

ZeroPivotException <: Exception

Exception thrown when a matrix factorization/solve encounters a zero in a pivot (diagonal) position and cannot proceed. This maynot mean that the matrix is singular: it may be fruitful to switch to a different factorization such as pivoted LU that can re-order variables to eliminate spurious zero pivots. Theinfo field indicates the location of (one of) the zero pivot(s).

RankDeficientException

Exception thrown when the input matrix isrank deficient. Some linear algebra functions, such as the Cholesky decomposition, are only applicable to matrices that are not rank deficient. Theinfo field indicates the computed rank of the matrix.

LAPACKException

Generic LAPACK exception thrown either during direct calls to theLAPACK functions or during calls to other functions that use the LAPACK functions internally but lack specialized error handling. Theinfo field contains additional information on the underlying error and depends on the LAPACK function that was invoked.

dot(x, y)x ⋅ y

Compute the dot product between two vectors. For complex vectors, the first vector is conjugated.

dot also works on arbitrary iterable objects, including arrays of any dimension, as long asdot is defined on the elements.

dot is semantically equivalent tosum(dot(vx,vy) for (vx,vy) in zip(x, y)), with the added restriction that the arguments must have equal lengths.

x ⋅ y (where⋅ can be typed by tab-completing\cdot in the REPL) is a synonym fordot(x, y).

Examples

julia> dot([1; 1], [2; 3])5julia> dot([im; im], [1; 1])0 - 2imjulia> dot(1:5, 2:6)70julia> x = fill(2., (5,5));julia> y = fill(3., (5,5));julia> dot(x, y)150.0

dot(x, A, y)

Compute the generalized dot productdot(x, A*y) between two vectorsx andy, without storing the intermediate result ofA*y. As for the two-argumentdot(_,_), this acts recursively. Moreover, for complex vectors, the first vector is conjugated.

Examples

julia> dot([1; 1], [1 2; 3 4], [2; 3])26julia> dot(1:5, reshape(1:25, 5, 5), 2:6)4850julia> ⋅(1:5, reshape(1:25, 5, 5), 2:6) == dot(1:5, reshape(1:25, 5, 5), 2:6)true

cross(x, y)×(x,y)

Compute the cross product of two 3-vectors.

Examples

julia> a = [0;1;0]3-element Vector{Int64}: 0 1 0julia> b = [0;0;1]3-element Vector{Int64}: 0 0 1julia> cross(a,b)3-element Vector{Int64}: 1 0 0

axpy!(α, x::AbstractArray, y::AbstractArray)

Overwritey withx * α + y and returny. Ifx andy have the same axes, it's equivalent withy .+= x .* a.

Examples

julia> x = [1; 2; 3];julia> y = [4; 5; 6];julia> axpy!(2, x, y)3-element Vector{Int64}:  6  9 12

axpby!(α, x::AbstractArray, β, y::AbstractArray)

Overwritey withx * α + y * β and returny. Ifx andy have the same axes, it's equivalent withy .= x .* a .+ y .* β.

Examples

julia> x = [1; 2; 3];julia> y = [4; 5; 6];julia> axpby!(2, x, 2, y)3-element Vector{Int64}: 10 14 18

rotate!(x, y, c, s)

Overwritex withc*x + s*y andy with-conj(s)*x + c*y. Returnsx andy.

reflect!(x, y, c, s)

Overwritex withc*x + s*y andy withconj(s)*x - c*y. Returnsx andy.

factorize(A)

Compute a convenient factorization ofA, based upon the type of the input matrix. IfA is passed as a generic matrix,factorize checks to see if it is symmetric/triangular/etc. To this end,factorize may check every element ofA to verify/rule out each property. It will short-circuit as soon as it can rule out symmetry/triangular structure. The return value can be reused for efficient solving of multiple systems. For example:A=factorize(A); x=A\b; y=A\C.

Examples

julia> A = Array(Bidiagonal(fill(1.0, (5, 5)), :U))5×5 Matrix{Float64}: 1.0  1.0  0.0  0.0  0.0 0.0  1.0  1.0  0.0  0.0 0.0  0.0  1.0  1.0  0.0 0.0  0.0  0.0  1.0  1.0 0.0  0.0  0.0  0.0  1.0julia> factorize(A) # factorize will check to see that A is already factorized5×5 Bidiagonal{Float64, Vector{Float64}}: 1.0  1.0   ⋅    ⋅    ⋅  ⋅   1.0  1.0   ⋅    ⋅  ⋅    ⋅   1.0  1.0   ⋅  ⋅    ⋅    ⋅   1.0  1.0  ⋅    ⋅    ⋅    ⋅   1.0

This returns a5×5 Bidiagonal{Float64}, which can now be passed to other linear algebra functions (e.g. eigensolvers) which will use specialized methods forBidiagonal types.

Diagonal(V::AbstractVector)

Construct a lazy matrix withV as its diagonal.

See alsoUniformScaling for the lazy identity matrixI,diagm to make a dense matrix, anddiag to extract diagonal elements.

Examples

julia> d = Diagonal([1, 10, 100])3×3 Diagonal{Int64, Vector{Int64}}: 1   ⋅    ⋅ ⋅  10    ⋅ ⋅   ⋅  100julia> diagm([7, 13])2×2 Matrix{Int64}: 7   0 0  13julia> ans + I2×2 Matrix{Int64}: 8   0 0  14julia> I(2)2×2 Diagonal{Bool, Vector{Bool}}: 1  ⋅ ⋅  1

julia> A = transpose([7.0 13.0])2×1 transpose(::Matrix{Float64}) with eltype Float64:  7.0 13.0julia> Diagonal(A)1×1 Diagonal{Float64, Vector{Float64}}: 7.0

Diagonal(A::AbstractMatrix)

Construct a matrix from the principal diagonal ofA. The input matrixA may be rectangular, but the output will be square.

Examples

julia> A = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> D = Diagonal(A)2×2 Diagonal{Int64, Vector{Int64}}: 1  ⋅ ⋅  4julia> A = [1 2 3; 4 5 6]2×3 Matrix{Int64}: 1  2  3 4  5  6julia> Diagonal(A)2×2 Diagonal{Int64, Vector{Int64}}: 1  ⋅ ⋅  5

Diagonal{T}(undef, n)

Construct an uninitializedDiagonal{T} of lengthn. Seeundef.

Bidiagonal(dv::V, ev::V, uplo::Symbol) where V <: AbstractVector

Constructs an upper (uplo=:U) or lower (uplo=:L) 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 withconvert(Array, _) (orArray(_) for short). The length ofev must be one less than the length ofdv.

Examples

julia> dv = [1, 2, 3, 4]4-element Vector{Int64}: 1 2 3 4julia> ev = [7, 8, 9]3-element Vector{Int64}: 7 8 9julia> Bu = Bidiagonal(dv, ev, :U) # ev is on the first superdiagonal4×4 Bidiagonal{Int64, Vector{Int64}}: 1  7  ⋅  ⋅ ⋅  2  8  ⋅ ⋅  ⋅  3  9 ⋅  ⋅  ⋅  4julia> Bl = Bidiagonal(dv, ev, :L) # ev is on the first subdiagonal4×4 Bidiagonal{Int64, Vector{Int64}}: 1  ⋅  ⋅  ⋅ 7  2  ⋅  ⋅ ⋅  8  3  ⋅ ⋅  ⋅  9  4

Bidiagonal(A, uplo::Symbol)

Construct aBidiagonal matrix from the main diagonal ofA and its first super- (ifuplo=:U) or sub-diagonal (ifuplo=:L).

Examples

julia> A = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]4×4 Matrix{Int64}: 1  1  1  1 2  2  2  2 3  3  3  3 4  4  4  4julia> Bidiagonal(A, :U) # contains the main diagonal and first superdiagonal of A4×4 Bidiagonal{Int64, Vector{Int64}}: 1  1  ⋅  ⋅ ⋅  2  2  ⋅ ⋅  ⋅  3  3 ⋅  ⋅  ⋅  4julia> Bidiagonal(A, :L) # contains the main diagonal and first subdiagonal of A4×4 Bidiagonal{Int64, Vector{Int64}}: 1  ⋅  ⋅  ⋅ 2  2  ⋅  ⋅ ⋅  3  3  ⋅ ⋅  ⋅  4  4

SymTridiagonal(dv::V, ev::V) where V <: AbstractVector

Construct a symmetric tridiagonal matrix from the diagonal (dv) and first sub/super-diagonal (ev), respectively. The result is of typeSymTridiagonal and provides efficient specialized eigensolvers, but may be converted into a regular matrix withconvert(Array, _) (orArray(_) for short).

ForSymTridiagonal block matrices, the elements ofdv are symmetrized. The argumentev is interpreted as the superdiagonal. Blocks from the subdiagonal are (materialized) transpose of the corresponding superdiagonal blocks.

Examples

julia> dv = [1, 2, 3, 4]4-element Vector{Int64}: 1 2 3 4julia> ev = [7, 8, 9]3-element Vector{Int64}: 7 8 9julia> SymTridiagonal(dv, ev)4×4 SymTridiagonal{Int64, Vector{Int64}}: 1  7  ⋅  ⋅ 7  2  8  ⋅ ⋅  8  3  9 ⋅  ⋅  9  4julia> A = SymTridiagonal(fill([1 2; 3 4], 3), fill([1 2; 3 4], 2));julia> A[1,1]2×2 Symmetric{Int64, Matrix{Int64}}: 1  2 2  4julia> A[1,2]2×2 Matrix{Int64}: 1  2 3  4julia> A[2,1]2×2 Matrix{Int64}: 1  3 2  4

SymTridiagonal(A::AbstractMatrix)

Construct a symmetric tridiagonal matrix from the diagonal and first superdiagonal of the symmetric matrixA.

Examples

julia> A = [1 2 3; 2 4 5; 3 5 6]3×3 Matrix{Int64}: 1  2  3 2  4  5 3  5  6julia> SymTridiagonal(A)3×3 SymTridiagonal{Int64, Vector{Int64}}: 1  2  ⋅ 2  4  5 ⋅  5  6julia> B = reshape([[1 2; 2 3], [1 2; 3 4], [1 3; 2 4], [1 2; 2 3]], 2, 2);julia> SymTridiagonal(B)2×2 SymTridiagonal{Matrix{Int64}, Vector{Matrix{Int64}}}: [1 2; 2 3]  [1 3; 2 4] [1 2; 3 4]  [1 2; 2 3]

Tridiagonal(dl::V, d::V, du::V) where V <: AbstractVector

Construct a tridiagonal matrix from the first subdiagonal, diagonal, and first superdiagonal, respectively. The result is of typeTridiagonal and provides efficient specialized linear solvers, but may be converted into a regular matrix withconvert(Array, _) (orArray(_) for short). The lengths ofdl anddu must be one less than the length ofd.

Examples

julia> dl = [1, 2, 3];julia> du = [4, 5, 6];julia> d = [7, 8, 9, 0];julia> Tridiagonal(dl, d, du)4×4 Tridiagonal{Int64, Vector{Int64}}: 7  4  ⋅  ⋅ 1  8  5  ⋅ ⋅  2  9  6 ⋅  ⋅  3  0

Tridiagonal(A)

Construct a tridiagonal matrix from the first sub-diagonal, diagonal and first super-diagonal of the matrixA.

Examples

julia> A = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4]4×4 Matrix{Int64}: 1  2  3  4 1  2  3  4 1  2  3  4 1  2  3  4julia> Tridiagonal(A)4×4 Tridiagonal{Int64, Vector{Int64}}: 1  2  ⋅  ⋅ 1  2  3  ⋅ ⋅  2  3  4 ⋅  ⋅  3  4

Symmetric(A::AbstractMatrix, uplo::Symbol=:U)

Construct aSymmetric view of the upper (ifuplo = :U) or lower (ifuplo = :L) triangle of the matrixA.

Symmetric views are mainly useful for real-symmetric matrices, for which specialized algorithms (e.g. for eigenproblems) are enabled forSymmetric types. More generally, see alsoHermitian(A) for Hermitian matricesA == A', which is effectively equivalent toSymmetric for real matrices but is also useful for complex matrices. (Whereas complexSymmetric matrices are supported but have few if any specialized algorithms.)

To compute the symmetric part of a real matrix, or more generally the Hermitian part(A + A') / 2 of a real or complex matrixA, usehermitianpart.

Examples

julia> A = [1 2 3; 4 5 6; 7 8 9]3×3 Matrix{Int64}: 1  2  3 4  5  6 7  8  9julia> Supper = Symmetric(A)3×3 Symmetric{Int64, Matrix{Int64}}: 1  2  3 2  5  6 3  6  9julia> Slower = Symmetric(A, :L)3×3 Symmetric{Int64, Matrix{Int64}}: 1  4  7 4  5  8 7  8  9julia> hermitianpart(A)3×3 Hermitian{Float64, Matrix{Float64}}: 1.0  3.0  5.0 3.0  5.0  7.0 5.0  7.0  9.0

Note thatSupper will not be equal toSlower unlessA is itself symmetric (e.g. ifA == transpose(A)).

Hermitian(A::AbstractMatrix, uplo::Symbol=:U)

Construct aHermitian view of the upper (ifuplo = :U) or lower (ifuplo = :L) triangle of the matrixA.

To compute the Hermitian part ofA, usehermitianpart.

Examples

julia> A = [1 2+2im 3-3im; 4 5 6-6im; 7 8+8im 9]3×3 Matrix{Complex{Int64}}: 1+0im  2+2im  3-3im 4+0im  5+0im  6-6im 7+0im  8+8im  9+0imjulia> Hupper = Hermitian(A)3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: 1+0im  2+2im  3-3im 2-2im  5+0im  6-6im 3+3im  6+6im  9+0imjulia> Hlower = Hermitian(A, :L)3×3 Hermitian{Complex{Int64}, Matrix{Complex{Int64}}}: 1+0im  4+0im  7+0im 4+0im  5+0im  8-8im 7+0im  8+8im  9+0imjulia> hermitianpart(A)3×3 Hermitian{ComplexF64, Matrix{ComplexF64}}: 1.0+0.0im  3.0+1.0im  5.0-1.5im 3.0-1.0im  5.0+0.0im  7.0-7.0im 5.0+1.5im  7.0+7.0im  9.0+0.0im

Note thatHupper will not be equal toHlower unlessA is itself Hermitian (e.g. ifA == adjoint(A)).

All non-real parts of the diagonal will be ignored.

Hermitian(fill(complex(1,1), 1, 1)) == fill(1, 1, 1)

LowerTriangular(A::AbstractMatrix)

Construct aLowerTriangular view of the matrixA.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]3×3 Matrix{Float64}: 1.0  2.0  3.0 4.0  5.0  6.0 7.0  8.0  9.0julia> LowerTriangular(A)3×3 LowerTriangular{Float64, Matrix{Float64}}: 1.0   ⋅    ⋅ 4.0  5.0   ⋅ 7.0  8.0  9.0

UpperTriangular(A::AbstractMatrix)

Construct anUpperTriangular view of the matrixA.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]3×3 Matrix{Float64}: 1.0  2.0  3.0 4.0  5.0  6.0 7.0  8.0  9.0julia> UpperTriangular(A)3×3 UpperTriangular{Float64, Matrix{Float64}}: 1.0  2.0  3.0  ⋅   5.0  6.0  ⋅    ⋅   9.0

UnitLowerTriangular(A::AbstractMatrix)

Construct aUnitLowerTriangular view of the matrixA. Such a view has theoneunit of theeltype ofA on its diagonal.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]3×3 Matrix{Float64}: 1.0  2.0  3.0 4.0  5.0  6.0 7.0  8.0  9.0julia> UnitLowerTriangular(A)3×3 UnitLowerTriangular{Float64, Matrix{Float64}}: 1.0   ⋅    ⋅ 4.0  1.0   ⋅ 7.0  8.0  1.0

UnitUpperTriangular(A::AbstractMatrix)

Construct anUnitUpperTriangular view of the matrixA. Such a view has theoneunit of theeltype ofA on its diagonal.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]3×3 Matrix{Float64}: 1.0  2.0  3.0 4.0  5.0  6.0 7.0  8.0  9.0julia> UnitUpperTriangular(A)3×3 UnitUpperTriangular{Float64, Matrix{Float64}}: 1.0  2.0  3.0  ⋅   1.0  6.0  ⋅    ⋅   1.0

UpperHessenberg(A::AbstractMatrix)

Construct anUpperHessenberg view of the matrixA. Entries ofA below the first subdiagonal are ignored.

Efficient algorithms are implemented forH \ b,det(H), and similar.

See also thehessenberg function to factor any matrix into a similar upper-Hessenberg matrix.

IfF::Hessenberg 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 withconvert(Array, _) (orArray(_) for short).

Iterating the decomposition produces the factorsF.Q andF.H.

Examples

julia> A = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]4×4 Matrix{Int64}:  1   2   3   4  5   6   7   8  9  10  11  12 13  14  15  16julia> UpperHessenberg(A)4×4 UpperHessenberg{Int64, Matrix{Int64}}: 1   2   3   4 5   6   7   8 ⋅  10  11  12 ⋅   ⋅  15  16

UniformScaling{T<:Number}

Generically sized uniform scaling operator defined as a scalar times the identity operator,λ*I. Although without an explicitsize, it acts similarly to a matrix in many cases and includes support for some indexing. See alsoI.

Examples

julia> J = UniformScaling(2.)UniformScaling{Float64}2.0*Ijulia> A = [1. 2.; 3. 4.]2×2 Matrix{Float64}: 1.0  2.0 3.0  4.0julia> J*A2×2 Matrix{Float64}: 2.0  4.0 6.0  8.0julia> J[1:2, 1:2]2×2 Matrix{Float64}: 2.0  0.0 0.0  2.0

I

An object of typeUniformScaling, representing an identity matrix of any size.

Examples

julia> fill(1, (5,6)) * I == fill(1, (5,6))truejulia> [1 2im 3; 1im 2 3] * I2×3 Matrix{Complex{Int64}}: 1+0im  0+2im  3+0im 0+1im  2+0im  3+0im

(I::UniformScaling)(n::Integer)

Construct aDiagonal matrix from aUniformScaling.

Examples

julia> I(3)3×3 Diagonal{Bool, Vector{Bool}}: 1  ⋅  ⋅ ⋅  1  ⋅ ⋅  ⋅  1julia> (0.7*I)(3)3×3 Diagonal{Float64, Vector{Float64}}: 0.7   ⋅    ⋅  ⋅   0.7   ⋅  ⋅    ⋅   0.7

LinearAlgebra.Factorization

Abstract type formatrix factorizations a.k.a. matrix decompositions. Seeonline documentation for a list of available matrix factorizations.

LU <: Factorization

Matrix factorization type of theLU factorization of a square matrixA. This is the return type oflu, the corresponding matrix factorization function.

The individual components of the factorizationF::LU can be accessed viagetproperty:

Iterating the factorization produces the componentsF.L,F.U, andF.p.

Examples

julia> A = [4 3; 6 3]2×2 Matrix{Int64}: 4  3 6  3julia> F = lu(A)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0       0.0 0.666667  1.0U factor:2×2 Matrix{Float64}: 6.0  3.0 0.0  1.0julia> F.L * F.U == A[F.p, :]truejulia> l, u, p = lu(A); # destructuring via iterationjulia> l == F.L && u == F.U && p == F.ptrue

lu(A::AbstractSparseMatrixCSC; check = true, q = nothing, control = get_umfpack_control()) -> F::UmfpackLU

Compute the LU factorization of a sparse matrixA.

For sparseA with real or complex element type, the return type ofF isUmfpackLU{Tv, Ti}, withTv =Float64 orComplexF64 respectively andTi is an integer type (Int32 orInt64).

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

The permutationq can either be a permutation vector ornothing. If no permutation vector is provided orq isnothing, UMFPACK's default is used. If the permutation is not zero-based, a zero-based copy is made.

Thecontrol vector defaults to the Julia SparseArrays package's default configuration for UMFPACK (NB: this is modified from the UMFPACK defaults to disable iterative refinement), but can be changed by passing a vector of lengthUMFPACK_CONTROL, see the UMFPACK manual for possible configurations. For example to reenable iterative refinement:

umfpack_control = SparseArrays.UMFPACK.get_umfpack_control(Float64, Int64) # read Julia default configuration for a Float64 sparse matrixSparseArrays.UMFPACK.show_umf_ctrl(umfpack_control) # optional - display valuesumfpack_control[SparseArrays.UMFPACK.JL_UMFPACK_IRSTEP] = 2.0 # reenable iterative refinement (2 is UMFPACK default max iterative refinement steps)Alu = lu(A; control = umfpack_control)x = Alu \ b   # solve Ax = b, including UMFPACK iterative refinement

The individual components of the factorizationF can be accessed by indexing:

The relation betweenF andA is

F.L*F.U == (F.Rs .* A)[F.p, F.q]

F further supports the following functions:

• \
• det

See alsolu!

lu(A, pivot = RowMaximum(); check = true, allowsingular = false) -> F::LU

Compute the LU factorization ofA.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

By default, withcheck = true, an error is also thrown when the decomposition produces valid factors, but the upper-triangular factorU is rank-deficient. This may be changed by passingallowsingular = true.

In most cases, ifA is a subtypeS ofAbstractMatrix{T} with an element typeT supporting+,-,* and/, the return type isLU{T,S{T}}.

In general, LU factorization involves a permutation of the rows of the matrix (corresponding to theF.p output described below), known as "pivoting" (because it corresponds to choosing which row contains the "pivot", the diagonal entry ofF.U). One of the following pivoting strategies can be selected via the optionalpivot argument:

• RowMaximum() (default): the standard pivoting strategy; the pivot corresponds to the element of maximum absolute value among the remaining, to be factorized rows. This pivoting strategy requires the element type to also supportabs and<. (This is generally the only numerically stable option for floating-point matrices.)
• RowNonZero(): the pivot corresponds to the first non-zero element among the remaining, to be factorized rows. (This corresponds to the typical choice in hand calculations, and is also useful for more general algebraic number types that supportiszero but notabs or<.)
• NoPivot(): pivoting turned off (will fail if a zero entry is encountered in a pivot position, even whenallowsingular = true).

The individual components of the factorizationF can be accessed viagetproperty:

Iterating the factorization produces the componentsF.L,F.U, andF.p.

The relationship betweenF andA is

F.L*F.U == A[F.p, :]

F further supports the following functions:

Examples

julia> A = [4 3; 6 3]2×2 Matrix{Int64}: 4  3 6  3julia> F = lu(A)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0       0.0 0.666667  1.0U factor:2×2 Matrix{Float64}: 6.0  3.0 0.0  1.0julia> F.L * F.U == A[F.p, :]truejulia> l, u, p = lu(A); # destructuring via iterationjulia> l == F.L && u == F.U && p == F.ptruejulia> lu([1 2; 1 2], allowsingular = true)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0  0.0 1.0  1.0U factor (rank-deficient):2×2 Matrix{Float64}: 1.0  2.0 0.0  0.0

lu!(F::UmfpackLU, A::AbstractSparseMatrixCSC; check=true, reuse_symbolic=true, q=nothing) -> F::UmfpackLU

Compute the LU factorization of a sparse matrixA, reusing the symbolic factorization of an already existing LU factorization stored inF. Unlessreuse_symbolic is set to false, the sparse matrixA must have an identical nonzero pattern as the matrix used to create the LU factorizationF, otherwise an error is thrown. If the size ofA andF differ, all vectors will be resized accordingly.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

The permutationq can either be a permutation vector ornothing. If no permutation vector is provided orq isnothing, UMFPACK's default is used. If the permutation is not zero based, a zero based copy is made.

See alsolu

Examples

julia> A = sparse(Float64[1.0 2.0; 0.0 3.0]);julia> F = lu(A);julia> B = sparse(Float64[1.0 1.0; 0.0 1.0]);julia> lu!(F, B);julia> F \ ones(2)2-element Vector{Float64}: 0.0 1.0

lu!(A, pivot = RowMaximum(); check = true, allowsingular = false) -> LU

lu! is the same aslu, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.

Examples

julia> A = [4. 3.; 6. 3.]2×2 Matrix{Float64}: 4.0  3.0 6.0  3.0julia> F = lu!(A)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0       0.0 0.666667  1.0U factor:2×2 Matrix{Float64}: 6.0  3.0 0.0  1.0julia> iA = [4 3; 6 3]2×2 Matrix{Int64}: 4  3 6  3julia> lu!(iA)ERROR: InexactError: Int64(0.6666666666666666)Stacktrace:[...]

Cholesky <: Factorization

Matrix factorization type of the Cholesky factorization of a dense symmetric/Hermitian positive definite matrixA. This is the return type ofcholesky, the corresponding matrix factorization function.

The triangular Cholesky factor can be obtained from the factorizationF::Cholesky viaF.L andF.U, whereA ≈ F.U' * F.U ≈ F.L * F.L'.

The following functions are available forCholesky objects:size,\,inv,det,logdet andisposdef.

Iterating the decomposition produces the componentsL andU.

Examples

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]3×3 Matrix{Float64}:   4.0   12.0  -16.0  12.0   37.0  -43.0 -16.0  -43.0   98.0julia> C = cholesky(A)Cholesky{Float64, Matrix{Float64}}U factor:3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0  6.0  -8.0  ⋅   1.0   5.0  ⋅    ⋅    3.0julia> C.U3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0  6.0  -8.0  ⋅   1.0   5.0  ⋅    ⋅    3.0julia> C.L3×3 LowerTriangular{Float64, Matrix{Float64}}:  2.0   ⋅    ⋅  6.0  1.0   ⋅ -8.0  5.0  3.0julia> C.L * C.U == Atruejulia> l, u = C; # destructuring via iterationjulia> l == C.L && u == C.Utrue

CholeskyPivoted

Matrix factorization type of the pivoted Cholesky factorization of a dense symmetric/Hermitian positive semi-definite matrixA. This is the return type ofcholesky(_, ::RowMaximum), the corresponding matrix factorization function.

The triangular Cholesky factor can be obtained from the factorizationF::CholeskyPivoted viaF.L andF.U, and the permutation viaF.p, whereA[F.p, F.p] ≈ Ur' * Ur ≈ Lr * Lr' withUr = F.U[1:F.rank, :] andLr = F.L[:, 1:F.rank], or alternativelyA ≈ Up' * Up ≈ Lp * Lp' withUp = F.U[1:F.rank, invperm(F.p)] andLp = F.L[invperm(F.p), 1:F.rank].

The following functions are available forCholeskyPivoted objects:size,\,inv,det, andrank.

Iterating the decomposition produces the componentsL andU.

Examples

julia> X = [1.0, 2.0, 3.0, 4.0];julia> A = X * X';julia> C = cholesky(A, RowMaximum(), check = false)CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}U factor with rank 1:4×4 UpperTriangular{Float64, Matrix{Float64}}: 4.0  2.0  3.0  1.0  ⋅   0.0  6.0  2.0  ⋅    ⋅   9.0  3.0  ⋅    ⋅    ⋅   1.0permutation:4-element Vector{Int64}: 4 2 3 1julia> C.U[1:C.rank, :]' * C.U[1:C.rank, :] ≈ A[C.p, C.p]truejulia> l, u = C; # destructuring via iterationjulia> l == C.L && u == C.Utrue

cholesky(A, NoPivot(); check = true) -> Cholesky

Compute the Cholesky factorization of a dense symmetric positive definite matrixA and return aCholesky factorization. The matrixA can either be aSymmetric orHermitianAbstractMatrix or aperfectly symmetric or HermitianAbstractMatrix.

The triangular Cholesky factor can be obtained from the factorizationF viaF.L andF.U, whereA ≈ F.U' * F.U ≈ F.L * F.L'.

The following functions are available forCholesky objects:size,\,inv,det,logdet andisposdef.

If you have a matrixA that is slightly non-Hermitian due to roundoff errors in its construction, wrap it inHermitian(A) before passing it tocholesky in order to treat it as perfectly Hermitian.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

Examples

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]3×3 Matrix{Float64}:   4.0   12.0  -16.0  12.0   37.0  -43.0 -16.0  -43.0   98.0julia> C = cholesky(A)Cholesky{Float64, Matrix{Float64}}U factor:3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0  6.0  -8.0  ⋅   1.0   5.0  ⋅    ⋅    3.0julia> C.U3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0  6.0  -8.0  ⋅   1.0   5.0  ⋅    ⋅    3.0julia> C.L3×3 LowerTriangular{Float64, Matrix{Float64}}:  2.0   ⋅    ⋅  6.0  1.0   ⋅ -8.0  5.0  3.0julia> C.L * C.U == Atrue

cholesky(A, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivoted

Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrixA and return aCholeskyPivoted factorization. The matrixA can either be aSymmetric orHermitianAbstractMatrix or aperfectly symmetric or HermitianAbstractMatrix.

The triangular Cholesky factor can be obtained from the factorizationF viaF.L andF.U, and the permutation viaF.p, whereA[F.p, F.p] ≈ Ur' * Ur ≈ Lr * Lr' withUr = F.U[1:F.rank, :] andLr = F.L[:, 1:F.rank], or alternativelyA ≈ Up' * Up ≈ Lp * Lp' withUp = F.U[1:F.rank, invperm(F.p)] andLp = F.L[invperm(F.p), 1:F.rank].

The following functions are available forCholeskyPivoted objects:size,\,inv,det, andrank.

The argumenttol determines the tolerance for determining the rank. For negative values, the tolerance is equal toeps()*size(A,1)*maximum(diag(A)).

If you have a matrixA that is slightly non-Hermitian due to roundoff errors in its construction, wrap it inHermitian(A) before passing it tocholesky in order to treat it as perfectly Hermitian.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

Examples

julia> X = [1.0, 2.0, 3.0, 4.0];julia> A = X * X';julia> C = cholesky(A, RowMaximum(), check = false)CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}U factor with rank 1:4×4 UpperTriangular{Float64, Matrix{Float64}}: 4.0  2.0  3.0  1.0  ⋅   0.0  6.0  2.0  ⋅    ⋅   9.0  3.0  ⋅    ⋅    ⋅   1.0permutation:4-element Vector{Int64}: 4 2 3 1julia> C.U[1:C.rank, :]' * C.U[1:C.rank, :] ≈ A[C.p, C.p]truejulia> l, u = C; # destructuring via iterationjulia> l == C.L && u == C.Utrue

cholesky(A::SparseMatrixCSC; shift = 0.0, check = true, perm = nothing) -> CHOLMOD.Factor

Compute the Cholesky factorization of a sparse positive definite matrixA.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian. Ifperm is not given, a fill-reducing permutation is used.F = cholesky(A) is most frequently used to solve systems of equations withF\b, but also the methodsdiag,det, andlogdet 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 extract "combined" factors likePtL = F.PtL (the equivalent ofP'*L) andLtP = F.UP (the equivalent ofL'*P).

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

Setting the optionalshift keyword argument computes the factorization ofA+shift*I instead ofA. If theperm argument is provided, it should be a permutation of1:size(A,1) giving the ordering to use (instead of CHOLMOD's default AMD ordering).

See alsoldlt for a similar factorization that does not require positive definiteness, but can be significantly slower thancholesky.

Examples

In the following example, the fill-reducing permutation used is[3, 2, 1]. Ifperm is set to1:3 to enforce no permutation, the number of nonzero elements in the factor is 6.

julia> A = [2 1 1; 1 2 0; 1 0 2]3×3 Matrix{Int64}: 2  1  1 1  2  0 1  0  2julia> C = cholesky(sparse(A))SparseArrays.CHOLMOD.Factor{Float64, Int64}type:    LLtmethod:  simplicialmaxnnz:  5nnz:     5success: truejulia> C.p3-element Vector{Int64}: 3 2 1julia> L = sparse(C.L);julia> Matrix(L)3×3 Matrix{Float64}: 1.41421   0.0       0.0 0.0       1.41421   0.0 0.707107  0.707107  1.0julia> L * L' ≈ A[C.p, C.p]truejulia> P = sparse(1:3, C.p, ones(3))3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries:  ⋅    ⋅   1.0  ⋅   1.0   ⋅ 1.0   ⋅    ⋅julia> P' * L * L' * P ≈ Atruejulia> C = cholesky(sparse(A), perm=1:3)SparseArrays.CHOLMOD.Factor{Float64, Int64}type:    LLtmethod:  simplicialmaxnnz:  6nnz:     6success: truejulia> L = sparse(C.L);julia> Matrix(L)3×3 Matrix{Float64}: 1.41421    0.0       0.0 0.707107   1.22474   0.0 0.707107  -0.408248  1.1547julia> L * L' ≈ Atrue

cholesky!(A::AbstractMatrix, NoPivot(); check = true) -> Cholesky

The same ascholesky, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.

Examples

julia> A = [1 2; 2 50]2×2 Matrix{Int64}: 1   2 2  50julia> cholesky!(A)ERROR: InexactError: Int64(6.782329983125268)Stacktrace:[...]

cholesky!(A::AbstractMatrix, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivoted

The same ascholesky, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.

cholesky!(F::CHOLMOD.Factor, A::SparseMatrixCSC; shift = 0.0, check = true) -> CHOLMOD.Factor

Compute the Cholesky ($LL'$) factorization ofA, reusing the symbolic factorizationF.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian.

See alsocholesky.

lowrankupdate(C::Cholesky, v::AbstractVector) -> CC::Cholesky

Update a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U + v*v') but the computation ofCC only usesO(n^2) operations.

lowrankupdate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.Factor

Get anLDLt Factorization ofA + C*C' given anLDLt orLLt factorizationF ofA.

The returned factor is always anLDLt factorization.

See alsolowrankupdate!,lowrankdowndate,lowrankdowndate!.

lowrankdowndate(C::Cholesky, v::AbstractVector) -> CC::Cholesky

Downdate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U - v*v') but the computation ofCC only usesO(n^2) operations.

lowrankdowndate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.Factor

Get anLDLt Factorization ofA + C*C' given anLDLt orLLt factorizationF ofA.

The returned factor is always anLDLt factorization.

See alsolowrankdowndate!,lowrankupdate,lowrankupdate!.

lowrankupdate!(C::Cholesky, v::AbstractVector) -> CC::Cholesky

Update a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U + v*v') but the computation ofCC only usesO(n^2) operations. The input factorizationC is updated in place such that on exitC == CC. The vectorv is destroyed during the computation.

lowrankupdate!(F::CHOLMOD.Factor, C::AbstractArray)

Update anLDLt orLLt FactorizationF ofA to a factorization ofA + C*C'.

LLt factorizations are converted toLDLt.

See alsolowrankupdate,lowrankdowndate,lowrankdowndate!.

lowrankdowndate!(C::Cholesky, v::AbstractVector) -> CC::Cholesky

Downdate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U - v*v') but the computation ofCC only usesO(n^2) operations. The input factorizationC is updated in place such that on exitC == CC. The vectorv is destroyed during the computation.

lowrankdowndate!(F::CHOLMOD.Factor, C::AbstractArray)

Update anLDLt orLLt FactorizationF ofA to a factorization ofA - C*C'.

LLt factorizations are converted toLDLt.

See alsolowrankdowndate,lowrankupdate,lowrankupdate!.

LDLt <: Factorization

Matrix factorization type of theLDLt factorization of a realSymTridiagonal matrixS such thatS = L*Diagonal(d)*L', whereL is aUnitLowerTriangular matrix andd is a vector. The main use of anLDLt factorizationF = ldlt(S) is to solve the linear system of equationsSx = b withF\b. This is the return type ofldlt, the corresponding matrix factorization function.

The individual components of the factorizationF::LDLt can be accessed viagetproperty:

Examples

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 3.0  1.0   ⋅ 1.0  4.0  2.0  ⋅   2.0  5.0julia> F = ldlt(S)LDLt{Float64, SymTridiagonal{Float64, Vector{Float64}}}L factor:3×3 UnitLowerTriangular{Float64, SymTridiagonal{Float64, Vector{Float64}}}: 1.0        ⋅         ⋅ 0.333333  1.0        ⋅ 0.0       0.545455  1.0D factor:3×3 Diagonal{Float64, Vector{Float64}}: 3.0   ⋅        ⋅  ⋅   3.66667   ⋅  ⋅    ⋅       3.90909

ldlt(S::SymTridiagonal) -> LDLt

Compute anLDLt (i.e.,$LDL^T$) factorization of the real symmetric tridiagonal matrixS such thatS = L*Diagonal(d)*L' whereL is a unit lower triangular matrix andd is a vector. The main use of anLDLt factorizationF = ldlt(S) is to solve the linear system of equationsSx = b withF\b.

See alsobunchkaufman for a similar, but pivoted, factorization of arbitrary symmetric or Hermitian matrices.

Examples

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 3.0  1.0   ⋅ 1.0  4.0  2.0  ⋅   2.0  5.0julia> ldltS = ldlt(S);julia> b = [6., 7., 8.];julia> ldltS \ b3-element Vector{Float64}: 1.7906976744186047 0.627906976744186 1.3488372093023255julia> S \ b3-element Vector{Float64}: 1.7906976744186047 0.627906976744186 1.3488372093023255

ldlt(A::SparseMatrixCSC; shift = 0.0, check = true, perm=nothing) -> CHOLMOD.Factor

Compute the$LDL'$ factorization of a sparse matrixA.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian. A fill-reducing permutation is used.F = ldlt(A) is most frequently used to solve systems of equationsA*x = b withF\b. The returned factorization objectF also supports the methodsdiag,det,logdet, andinv. You can 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 is typically preferable to extract "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.

Unlike the related Cholesky factorization, the$LDL'$ factorization does not requireA to be positive definite. However, it still requires all leading principal minors to be well-conditioned and will fail if this is not satisfied.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

Setting the optionalshift keyword argument computes the factorization ofA+shift*I instead ofA. If theperm argument is provided, it should be a permutation of1:size(A,1) giving the ordering to use (instead of CHOLMOD's default AMD ordering).

See alsocholesky for a factorization that can be significantly faster thanldlt, but requiresA to be positive definite.

ldlt!(S::SymTridiagonal) -> LDLt

Same asldlt, but saves space by overwriting the inputS, instead of creating a copy.

Examples

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 3.0  1.0   ⋅ 1.0  4.0  2.0  ⋅   2.0  5.0julia> ldltS = ldlt!(S);julia> ldltS === Sfalsejulia> S3×3 SymTridiagonal{Float64, Vector{Float64}}: 3.0       0.333333   ⋅ 0.333333  3.66667   0.545455  ⋅        0.545455  3.90909

ldlt!(F::CHOLMOD.Factor, A::SparseMatrixCSC; shift = 0.0, check = true) -> CHOLMOD.Factor

Compute the$LDL'$ factorization ofA, reusing the symbolic factorizationF.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian.

See alsoldlt.

QR <: Factorization

A QR matrix factorization stored in a packed format, typically obtained fromqr. If$A$ is anm×n matrix, then

\[A = Q R\]

where$Q$ is an orthogonal/unitary matrix and$R$ is upper triangular. The matrix$Q$ is stored as a sequence of Householder reflectors$v_i$ and coefficients$\tau_i$ where:

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T).\]

Iterating the decomposition produces the componentsQ andR.

The object has two fields:

• factors is anm×n matrix.

  • The upper triangular part contains the elements of$R$, that isR = triu(F.factors) for aQR objectF.

  • The subdiagonal part contains the reflectors$v_i$ stored in a packed format where$v_i$ is the$i$th column of the matrixV = I + tril(F.factors, -1).

• τ is a vector of lengthmin(m,n) containing the coefficients$\tau_i$.

QRCompactWY <: Factorization

A QR matrix factorization stored in a compact blocked format, typically obtained fromqr. If$A$ is anm×n matrix, then

\[A = Q R\]

where$Q$ is an orthogonal/unitary matrix and$R$ is upper triangular. It is similar to theQR format except that the orthogonal/unitary matrix$Q$ is stored inCompact WY format^{[Schreiber1989]}. For the block size$n_b$, it is stored as am×n lower trapezoidal matrix$V$ and a matrix$T = (T_1 \; T_2 \; ... \; T_{b-1} \; T_b')$ composed of$b = \lceil \min(m,n) / n_b \rceil$ upper triangular matrices$T_j$ of size$n_b$×$n_b$ ($j = 1, ..., b-1$) and an upper trapezoidal$n_b$×$\min(m,n) - (b-1) n_b$ matrix$T_b'$ ($j=b$) whose upper square part denoted with$T_b$ satisfying

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T)= \prod_{j=1}^{b} (I - V_j T_j V_j^T)\]

such that$v_i$ is the$i$th column of$V$,$\tau_i$ is the$i$th element of[diag(T_1); diag(T_2); …; diag(T_b)], and$(V_1 \; V_2 \; ... \; V_b)$ is the leftm×min(m, n) block of$V$. When constructed usingqr, the block size is given by$n_b = \min(m, n, 36)$.

Iterating the decomposition produces the componentsQ andR.

The object has two fields:

• factors, as in theQR type, is anm×n matrix.

  • The upper triangular part contains the elements of$R$, that isR = triu(F.factors) for aQR objectF.

  • The subdiagonal part contains the reflectors$v_i$ stored in a packed format such thatV = I + tril(F.factors, -1).

• T is a$n_b$-by-$\min(m,n)$ matrix as described above. The subdiagonal elements for each triangular matrix$T_j$ are ignored.

QRPivoted <: Factorization

A QR matrix factorization with column pivoting in a packed format, typically obtained fromqr. If$A$ is anm×n matrix, then

\[A P = Q R\]

where$P$ is a permutation matrix,$Q$ is an orthogonal/unitary matrix and$R$ is upper triangular. The matrix$Q$ is stored as a sequence of Householder reflectors:

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T).\]

Iterating the decomposition produces the componentsQ,R, andp.

The object has three fields:

• factors is anm×n matrix.

  • The upper triangular part contains the elements of$R$, that isR = triu(F.factors) for aQR objectF.

  • The subdiagonal part contains the reflectors$v_i$ stored in a packed format where$v_i$ is the$i$th column of the matrixV = I + tril(F.factors, -1).

• τ is a vector of lengthmin(m,n) containing the coefficients$au_i$.

• jpvt is an integer vector of lengthn corresponding to the permutation$P$.

qr(A::SparseMatrixCSC; tol=_default_tol(A), ordering=ORDERING_DEFAULT) -> QRSparse

Compute theQR factorization of a sparse matrixA. Fill-reducing row and column permutations are used such thatF.R = F.Q'*A[F.prow,F.pcol]. The main application of this type is to solve least squares or underdetermined problems with\. The function calls the C library SPQR^[ACM933].

Examples

julia> A = sparse([1,2,3,4], [1,1,2,2], [1.0,1.0,1.0,1.0])4×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries: 1.0   ⋅ 1.0   ⋅  ⋅   1.0  ⋅   1.0julia> qr(A)SparseArrays.SPQR.QRSparse{Float64, Int64}Q factor:4×4 SparseArrays.SPQR.QRSparseQ{Float64, Int64}R factor:2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries: -1.41421    ⋅   ⋅       -1.41421Row permutation:4-element Vector{Int64}: 1 3 4 2Column permutation:2-element Vector{Int64}: 1 2

qr(A, pivot = NoPivot(); blocksize) -> F

Compute the QR factorization of the matrixA: an orthogonal (or unitary ifA is complex-valued) matrixQ, and an upper triangular matrixR such that

\[A = Q R\]

The returned objectF stores the factorization in a packed format:

• ifpivot == ColumnNorm() thenF is aQRPivoted object,

• otherwise if the element type ofA is a BLAS type (Float32,Float64,ComplexF32 orComplexF64), thenF is aQRCompactWY object,

• otherwiseF is aQR object.

The individual components of the decompositionF can be retrieved via property accessors:

• F.Q: the orthogonal/unitary matrixQ
• F.R: the upper triangular matrixR
• F.p: the permutation vector of the pivot (QRPivoted only)
• F.P: the permutation matrix of the pivot (QRPivoted only)

Iterating the decomposition produces the componentsQ,R, and if extantp.

The following functions are available for theQR objects:inv,size, and\. WhenA is rectangular,\ will return a least squares solution and if the solution is not unique, the one with smallest norm is returned. WhenA is not full rank, factorization with (column) pivoting is required to obtain a minimum norm solution.

Multiplication with respect to either full/square or non-full/squareQ is allowed, i.e. bothF.Q*F.R andF.Q*A are supported. AQ matrix can be converted into a regular matrix withMatrix. This operation returns the "thin" Q factor, i.e., ifA ism×n withm>=n, thenMatrix(F.Q) yields anm×n matrix with orthonormal columns. To retrieve the "full" Q factor, anm×m orthogonal matrix, useF.Q*I orcollect(F.Q). Ifm<=n, thenMatrix(F.Q) yields anm×m orthogonal matrix.

The block size for QR decomposition can be specified by keyword argumentblocksize :: Integer whenpivot == NoPivot() andA isa StridedMatrix{<:BlasFloat}. It is ignored whenblocksize > minimum(size(A)). SeeQRCompactWY.

Examples

julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0]3×2 Matrix{Float64}: 3.0  -6.0 4.0  -8.0 0.0   1.0julia> F = qr(A)LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}Q factor: 3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}R factor:2×2 Matrix{Float64}: -5.0  10.0  0.0  -1.0julia> F.Q * F.R == Atrue

qr!(A, pivot = NoPivot(); blocksize)

qr! is the same asqr whenA is a subtype ofAbstractMatrix, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.

Examples

julia> a = [1. 2.; 3. 4.]2×2 Matrix{Float64}: 1.0  2.0 3.0  4.0julia> qr!(a)LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}Q factor: 2×2 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}R factor:2×2 Matrix{Float64}: -3.16228  -4.42719  0.0      -0.632456julia> a = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> qr!(a)ERROR: InexactError: Int64(3.1622776601683795)Stacktrace:[...]

LQ <: Factorization

Matrix factorization type of theLQ factorization of a matrixA. TheLQ decomposition is theQR decomposition oftranspose(A). This is the return type oflq, the corresponding matrix factorization function.

IfS::LQ is the factorization object, the lower triangular component can be obtained viaS.L, and the orthogonal/unitary component viaS.Q, such thatA ≈ S.L*S.Q.

Iterating the decomposition produces the componentsS.L andS.Q.

Examples

julia> A = [5. 7.; -2. -4.]2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> S = lq(A)LQ{Float64, Matrix{Float64}, Vector{Float64}}L factor:2×2 Matrix{Float64}: -8.60233   0.0  4.41741  -0.697486Q factor: 2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}julia> S.L * S.Q2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> l, q = S; # destructuring via iterationjulia> l == S.L &&  q == S.Qtrue

lq(A) -> S::LQ

Compute the LQ decomposition ofA. The decomposition's lower triangular component can be obtained from theLQ objectS viaS.L, and the orthogonal/unitary component viaS.Q, such thatA ≈ S.L*S.Q.

Iterating the decomposition produces the componentsS.L andS.Q.

The LQ decomposition is the QR decomposition oftranspose(A), and it is useful in order to compute the minimum-norm solutionlq(A) \ b to an underdetermined system of equations (A has more columns than rows, but has full row rank).

Examples

julia> A = [5. 7.; -2. -4.]2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> S = lq(A)LQ{Float64, Matrix{Float64}, Vector{Float64}}L factor:2×2 Matrix{Float64}: -8.60233   0.0  4.41741  -0.697486Q factor: 2×2 LinearAlgebra.LQPackedQ{Float64, Matrix{Float64}, Vector{Float64}}julia> S.L * S.Q2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> l, q = S; # destructuring via iterationjulia> l == S.L &&  q == S.Qtrue

lq!(A) -> LQ

Compute theLQ factorization ofA, using the input matrix as a workspace. See alsolq.

BunchKaufman <: Factorization

Matrix factorization type of the Bunch-Kaufman factorization of a symmetric or Hermitian matrixA asP'UDU'P orP'LDL'P, depending on whether the upper (the default) or the lower triangle is stored inA. IfA is complex symmetric thenU' andL' denote the unconjugated transposes, i.e.transpose(U) andtranspose(L), respectively. This is the return type ofbunchkaufman, the corresponding matrix factorization function.

IfS::BunchKaufman is the factorization object, the components can be obtained viaS.D,S.U orS.L as appropriate givenS.uplo, andS.p.

Iterating the decomposition produces the componentsS.D,S.U orS.L as appropriate givenS.uplo, andS.p.

Examples

julia> A = Float64.([1 2; 2 3])2×2 Matrix{Float64}: 1.0  2.0 2.0  3.0julia> S = bunchkaufman(A) # A gets wrapped internally by Symmetric(A)BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}D factor:2×2 Tridiagonal{Float64, Vector{Float64}}: -0.333333  0.0  0.0       3.0U factor:2×2 UnitUpperTriangular{Float64, Matrix{Float64}}: 1.0  0.666667  ⋅   1.0permutation:2-element Vector{Int64}: 1 2julia> d, u, p = S; # destructuring via iterationjulia> d == S.D && u == S.U && p == S.ptruejulia> S = bunchkaufman(Symmetric(A, :L))BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}D factor:2×2 Tridiagonal{Float64, Vector{Float64}}: 3.0   0.0 0.0  -0.333333L factor:2×2 UnitLowerTriangular{Float64, Matrix{Float64}}: 1.0        ⋅ 0.666667  1.0permutation:2-element Vector{Int64}: 2 1

bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman

Compute the Bunch-Kaufman^[Bunch1977] factorization of a symmetric or Hermitian matrixA asP'*U*D*U'*P orP'*L*D*L'*P, depending on which triangle is stored inA, and return aBunchKaufman object. Note that ifA is complex symmetric thenU' andL' denote the unconjugated transposes, i.e.transpose(U) andtranspose(L).

Iterating the decomposition produces the componentsS.D,S.U orS.L as appropriate givenS.uplo, andS.p.

Ifrook istrue, rook pivoting is used. Ifrook is false, rook pivoting is not used.

Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.

The following functions are available forBunchKaufman objects:size,\,inv,issymmetric,ishermitian,getindex.

Examples

julia> A = Float64.([1 2; 2 3])2×2 Matrix{Float64}: 1.0  2.0 2.0  3.0julia> S = bunchkaufman(A) # A gets wrapped internally by Symmetric(A)BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}D factor:2×2 Tridiagonal{Float64, Vector{Float64}}: -0.333333  0.0  0.0       3.0U factor:2×2 UnitUpperTriangular{Float64, Matrix{Float64}}: 1.0  0.666667  ⋅   1.0permutation:2-element Vector{Int64}: 1 2julia> d, u, p = S; # destructuring via iterationjulia> d == S.D && u == S.U && p == S.ptruejulia> S.U*S.D*S.U' - S.P*A*S.P'2×2 Matrix{Float64}: 0.0  0.0 0.0  0.0julia> S = bunchkaufman(Symmetric(A, :L))BunchKaufman{Float64, Matrix{Float64}, Vector{Int64}}D factor:2×2 Tridiagonal{Float64, Vector{Float64}}: 3.0   0.0 0.0  -0.333333L factor:2×2 UnitLowerTriangular{Float64, Matrix{Float64}}: 1.0        ⋅ 0.666667  1.0permutation:2-element Vector{Int64}: 2 1julia> S.L*S.D*S.L' - A[S.p, S.p]2×2 Matrix{Float64}: 0.0  0.0 0.0  0.0

bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman

bunchkaufman! is the same asbunchkaufman, but saves space by overwriting the inputA, instead of creating a copy.

Eigen <: Factorization

Matrix factorization type of the eigenvalue/spectral decomposition of a square matrixA. This is the return type ofeigen, the corresponding matrix factorization function.

IfF::Eigen is the factorization object, the eigenvalues can be obtained viaF.values and the eigenvectors as the columns of the matrixF.vectors. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

Examples

julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}values:3-element Vector{Float64}:  1.0  3.0 18.0vectors:3×3 Matrix{Float64}: 1.0  0.0  0.0 0.0  1.0  0.0 0.0  0.0  1.0julia> F.values3-element Vector{Float64}:  1.0  3.0 18.0julia> F.vectors3×3 Matrix{Float64}: 1.0  0.0  0.0 0.0  1.0  0.0 0.0  0.0  1.0julia> vals, vecs = F; # destructuring via iterationjulia> vals == F.values && vecs == F.vectorstrue

GeneralizedEigen <: Factorization

Matrix factorization type of the generalized eigenvalue/spectral decomposition ofA andB. This is the return type ofeigen, the corresponding matrix factorization function, when called with two matrix arguments.

IfF::GeneralizedEigen is the factorization object, the eigenvalues can be obtained viaF.values and the eigenvectors as the columns of the matrixF.vectors. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

Examples

julia> A = [1 0; 0 -1]2×2 Matrix{Int64}: 1   0 0  -1julia> B = [0 1; 1 0]2×2 Matrix{Int64}: 0  1 1  0julia> F = eigen(A, B)GeneralizedEigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}values:2-element Vector{ComplexF64}: 0.0 - 1.0im 0.0 + 1.0imvectors:2×2 Matrix{ComplexF64}:  0.0+1.0im   0.0-1.0im -1.0+0.0im  -1.0-0.0imjulia> F.values2-element Vector{ComplexF64}: 0.0 - 1.0im 0.0 + 1.0imjulia> F.vectors2×2 Matrix{ComplexF64}:  0.0+1.0im   0.0-1.0im -1.0+0.0im  -1.0-0.0imjulia> vals, vecs = F; # destructuring via iterationjulia> vals == F.values && vecs == F.vectorstrue

eigvals(A; permute::Bool=true, scale::Bool=true, sortby) -> values

Return the eigenvalues ofA.

For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvalue calculation. Thepermute,scale, andsortby keywords are the same as foreigen.

Examples

julia> diag_matrix = [1 0; 0 4]2×2 Matrix{Int64}: 1  0 0  4julia> eigvals(diag_matrix)2-element Vector{Float64}: 1.0 4.0

For a scalar input,eigvals will return a scalar.

Examples

julia> eigvals(-2)-2

eigvals(A, B) -> values

Compute the generalized eigenvalues ofA andB.

Examples

julia> A = [1 0; 0 -1]2×2 Matrix{Int64}: 1   0 0  -1julia> B = [0 1; 1 0]2×2 Matrix{Int64}: 0  1 1  0julia> eigvals(A,B)2-element Vector{ComplexF64}: 0.0 - 1.0im 0.0 + 1.0im

eigvals(A::Union{Hermitian, Symmetric}; alg::Algorithm = default_eigen_alg(A))) -> values

Return the eigenvalues ofA.

alg specifies which algorithm and LAPACK method to use for eigenvalue decomposition:

• alg = DivideAndConquer(): CallsLAPACK.syevd!.
• alg = QRIteration(): CallsLAPACK.syev!.
• alg = RobustRepresentations() (default): Multiple relatively robust representations method, CallsLAPACK.syevr!.

See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for a comparison of the accuracy and performance of different methods.

The defaultalg used may change in the future.

eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values

Return the eigenvalues ofA. It is possible to calculate only a subset of the eigenvalues by specifying aUnitRangeirange covering indices of the sorted eigenvalues, e.g. the 2nd to 8th eigenvalues.

Examples

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 1.0  2.0   ⋅ 2.0  2.0  3.0  ⋅   3.0  1.0julia> eigvals(A, 2:2)1-element Vector{Float64}: 0.9999999999999996julia> eigvals(A)3-element Vector{Float64}: -2.1400549446402604  1.0000000000000002  5.140054944640259

eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> values

Return the eigenvalues ofA. It is possible to calculate only a subset of the eigenvalues by specifying a pairvl andvu for the lower and upper boundaries of the eigenvalues.

Examples

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 1.0  2.0   ⋅ 2.0  2.0  3.0  ⋅   3.0  1.0julia> eigvals(A, -1, 2)1-element Vector{Float64}: 1.0000000000000009julia> eigvals(A)3-element Vector{Float64}: -2.1400549446402604  1.0000000000000002  5.140054944640259

eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values

Same aseigvals, but saves space by overwriting the inputA, instead of creating a copy. Thepermute,scale, andsortby keywords are the same as foreigen.

Examples

julia> A = [1. 2.; 3. 4.]2×2 Matrix{Float64}: 1.0  2.0 3.0  4.0julia> eigvals!(A)2-element Vector{Float64}: -0.3722813232690143  5.372281323269014julia> A2×2 Matrix{Float64}: -0.372281  -1.0  0.0        5.37228

eigvals!(A, B; sortby) -> values

Same aseigvals, but saves space by overwriting the inputA (andB), instead of creating copies.

Examples

julia> A = [1. 0.; 0. -1.]2×2 Matrix{Float64}: 1.0   0.0 0.0  -1.0julia> B = [0. 1.; 1. 0.]2×2 Matrix{Float64}: 0.0  1.0 1.0  0.0julia> eigvals!(A, B)2-element Vector{ComplexF64}: 0.0 - 1.0im 0.0 + 1.0imjulia> A2×2 Matrix{Float64}: -0.0  -1.0  1.0  -0.0julia> B2×2 Matrix{Float64}: 1.0  0.0 0.0  1.0

eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values

Same aseigvals, but saves space by overwriting the inputA, instead of creating a copy.irange is a range of eigenvalueindices to search for - for instance, the 2nd to 8th eigenvalues.

eigvals!(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> values

Same aseigvals, but saves space by overwriting the inputA, instead of creating a copy.vl is the lower bound of the interval to search for eigenvalues, andvu is the upper bound.

eigmax(A; permute::Bool=true, scale::Bool=true)

Return the largest eigenvalue ofA. The optionpermute=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. Note that if the eigenvalues ofA are complex, this method will fail, since complex numbers cannot be sorted.

Examples

julia> A = [0 im; -im 0]2×2 Matrix{Complex{Int64}}: 0+0im  0+1im 0-1im  0+0imjulia> eigmax(A)1.0julia> A = [0 im; -1 0]2×2 Matrix{Complex{Int64}}:  0+0im  0+1im -1+0im  0+0imjulia> eigmax(A)ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:`A` cannot have complex eigenvalues.Stacktrace:[...]

eigmin(A; permute::Bool=true, scale::Bool=true)

Return the smallest eigenvalue ofA. The optionpermute=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. Note that if the eigenvalues ofA are complex, this method will fail, since complex numbers cannot be sorted.

Examples

julia> A = [0 im; -im 0]2×2 Matrix{Complex{Int64}}: 0+0im  0+1im 0-1im  0+0imjulia> eigmin(A)-1.0julia> A = [0 im; -1 0]2×2 Matrix{Complex{Int64}}:  0+0im  0+1im -1+0im  0+0imjulia> eigmin(A)ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:`A` cannot have complex eigenvalues.Stacktrace:[...]

eigvecs(A::SymTridiagonal[, eigvals]) -> Matrix

Return a matrixM whose columns are the eigenvectors ofA. (Thekth eigenvector can be obtained from the sliceM[:, k].)

If the optional vector of eigenvalueseigvals is specified,eigvecs returns the specific corresponding eigenvectors.

Examples

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 1.0  2.0   ⋅ 2.0  2.0  3.0  ⋅   3.0  1.0julia> eigvals(A)3-element Vector{Float64}: -2.1400549446402604  1.0000000000000002  5.140054944640259julia> eigvecs(A)3×3 Matrix{Float64}:  0.418304  -0.83205      0.364299 -0.656749  -7.39009e-16  0.754109  0.627457   0.5547       0.546448julia> eigvecs(A, [1.])3×1 Matrix{Float64}:  0.8320502943378438  4.263514128092366e-17 -0.5547001962252291

eigvecs(A; permute::Bool=true, scale::Bool=true, `sortby`) -> Matrix

Return a matrixM whose columns are the eigenvectors ofA. (Thekth eigenvector can be obtained from the sliceM[:, k].) Thepermute,scale, andsortby keywords are the same as foreigen.

Examples

julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])3×3 Matrix{Float64}: 1.0  0.0  0.0 0.0  1.0  0.0 0.0  0.0  1.0

eigvecs(A, B) -> Matrix

Return a matrixM whose columns are the generalized eigenvectors ofA andB. (Thekth eigenvector can be obtained from the sliceM[:, k].)

Examples

julia> A = [1 0; 0 -1]2×2 Matrix{Int64}: 1   0 0  -1julia> B = [0 1; 1 0]2×2 Matrix{Int64}: 0  1 1  0julia> eigvecs(A, B)2×2 Matrix{ComplexF64}:  0.0+1.0im   0.0-1.0im -1.0+0.0im  -1.0-0.0im

eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen

Compute the eigenvalue decomposition ofA, returning anEigen factorization objectF which contains the eigenvalues inF.values and the normalized eigenvectors in the columns of the matrixF.vectors. This corresponds to solving an eigenvalue problem of the formAx = λx, whereA is a matrix,x is an eigenvector, andλ is an eigenvalue. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

The following functions are available forEigen objects:inv,det, andisposdef.

For general nonsymmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The optionpermute=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.

By default, the eigenvalues and vectors are sorted lexicographically by(real(λ),imag(λ)). A different comparison functionby(λ) can be passed tosortby, or you can passsortby=nothing to leave the eigenvalues in an arbitrary order. Some special matrix types (e.g.Diagonal orSymTridiagonal) may implement their own sorting convention and not accept asortby keyword.

Examples

julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}values:3-element Vector{Float64}:  1.0  3.0 18.0vectors:3×3 Matrix{Float64}: 1.0  0.0  0.0 0.0  1.0  0.0 0.0  0.0  1.0julia> F.values3-element Vector{Float64}:  1.0  3.0 18.0julia> F.vectors3×3 Matrix{Float64}: 1.0  0.0  0.0 0.0  1.0  0.0 0.0  0.0  1.0julia> vals, vecs = F; # destructuring via iterationjulia> vals == F.values && vecs == F.vectorstrue

eigen(A, B; sortby) -> GeneralizedEigen

Compute the generalized eigenvalue decomposition ofA andB, returning aGeneralizedEigen factorization objectF which contains the generalized eigenvalues inF.values and the generalized eigenvectors in the columns of the matrixF.vectors. This corresponds to solving a generalized eigenvalue problem of the formAx = λBx, whereA, B are matrices,x is an eigenvector, andλ is an eigenvalue. (Thekth generalized eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

By default, the eigenvalues and vectors are sorted lexicographically by(real(λ),imag(λ)). A different comparison functionby(λ) can be passed tosortby, or you can passsortby=nothing to leave the eigenvalues in an arbitrary order.

Examples

julia> A = [1 0; 0 -1]2×2 Matrix{Int64}: 1   0 0  -1julia> B = [0 1; 1 0]2×2 Matrix{Int64}: 0  1 1  0julia> F = eigen(A, B);julia> F.values2-element Vector{ComplexF64}: 0.0 - 1.0im 0.0 + 1.0imjulia> F.vectors2×2 Matrix{ComplexF64}:  0.0+1.0im   0.0-1.0im -1.0+0.0im  -1.0-0.0imjulia> vals, vecs = F; # destructuring via iterationjulia> vals == F.values && vecs == F.vectorstrue

eigen(A::Union{Hermitian, Symmetric}; alg::LinearAlgebra.Algorithm = LinearAlgebra.default_eigen_alg(A)) -> Eigen

Compute the eigenvalue decomposition ofA, returning anEigen factorization objectF which contains the eigenvalues inF.values and the orthonormal eigenvectors in the columns of the matrixF.vectors. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

alg specifies which algorithm and LAPACK method to use for eigenvalue decomposition:

• alg = DivideAndConquer(): CallsLAPACK.syevd!.
• alg = QRIteration(): CallsLAPACK.syev!.
• alg = RobustRepresentations() (default): Multiple relatively robust representations method, CallsLAPACK.syevr!.

See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for a comparison of the accuracy and performance of different algorithms.

The defaultalg used may change in the future.

The following functions are available forEigen objects:inv,det, andisposdef.

eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen

Compute the eigenvalue decomposition ofA, returning anEigen factorization objectF which contains the eigenvalues inF.values and the orthonormal eigenvectors in the columns of the matrixF.vectors. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

The following functions are available forEigen objects:inv,det, andisposdef.

TheUnitRangeirange specifies indices of the sorted eigenvalues to search for.

eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen

Compute the eigenvalue decomposition ofA, returning anEigen factorization objectF which contains the eigenvalues inF.values and the orthonormal eigenvectors in the columns of the matrixF.vectors. (Thekth eigenvector can be obtained from the sliceF.vectors[:, k].)

Iterating the decomposition produces the componentsF.values andF.vectors.

The following functions are available forEigen objects:inv,det, andisposdef.

vl is the lower bound of the window of eigenvalues to search for, andvu is the upper bound.

eigen!(A; permute, scale, sortby)eigen!(A, B; sortby)

Same aseigen, but saves space by overwriting the inputA (andB), instead of creating a copy.

Hessenberg <: Factorization

AHessenberg object represents the Hessenberg factorizationQHQ' of a square matrix, or a shiftQ(H+μI)Q' thereof, which is produced by thehessenberg function.

hessenberg(A) -> Hessenberg

Compute the Hessenberg decomposition ofA and return aHessenberg object. IfF is the factorization object, the unitary matrix can be accessed withF.Q (of typeLinearAlgebra.HessenbergQ) and the Hessenberg matrix withF.H (of typeUpperHessenberg), either of which may be converted to a regular matrix withMatrix(F.H) orMatrix(F.Q).

IfA isHermitian or real-Symmetric, then the Hessenberg decomposition produces a real-symmetric tridiagonal matrix andF.H is of typeSymTridiagonal.

Note that the shifted factorizationA+μI = Q (H+μI) Q' can be constructed efficiently byF + μ*I using theUniformScaling objectI, which creates a newHessenberg object with shared storage and a modified shift. The shift of a givenF is obtained byF.μ. This is useful because multiple shifted solves(F + μ*I) \ b (for differentμ and/orb) can be performed efficiently onceF is created.

Iterating the decomposition produces the factorsF.Q, F.H, F.μ.

Examples

julia> A = [4. 9. 7.; 4. 4. 1.; 4. 3. 2.]3×3 Matrix{Float64}: 4.0  9.0  7.0 4.0  4.0  1.0 4.0  3.0  2.0julia> F = hessenberg(A)Hessenberg{Float64, UpperHessenberg{Float64, Matrix{Float64}}, Matrix{Float64}, Vector{Float64}, Bool}Q factor: 3×3 LinearAlgebra.HessenbergQ{Float64, Matrix{Float64}, Vector{Float64}, false}H factor:3×3 UpperHessenberg{Float64, Matrix{Float64}}:  4.0      -11.3137       -1.41421 -5.65685    5.0           2.0   ⋅        -8.88178e-16   1.0julia> F.Q * F.H * F.Q'3×3 Matrix{Float64}: 4.0  9.0  7.0 4.0  4.0  1.0 4.0  3.0  2.0julia> q, h = F; # destructuring via iterationjulia> q == F.Q && h == F.Htrue

hessenberg!(A) -> Hessenberg

hessenberg! is the same ashessenberg, but saves space by overwriting the inputA, instead of creating a copy.

Schur <: Factorization

Matrix factorization type of the Schur factorization of a matrixA. This is the return type ofschur(_), the corresponding matrix factorization function.

IfF::Schur is the factorization object, the (quasi) triangular Schur factor can be obtained via eitherF.Schur orF.T and the orthogonal/unitary Schur vectors viaF.vectors orF.Z such thatA = F.vectors * F.Schur * F.vectors'. The eigenvalues ofA can be obtained withF.values.

Iterating the decomposition produces the componentsF.T,F.Z, andF.values.

Examples

julia> A = [5. 7.; -2. -4.]2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> F = schur(A)Schur{Float64, Matrix{Float64}, Vector{Float64}}T factor:2×2 Matrix{Float64}: 3.0   9.0 0.0  -2.0Z factor:2×2 Matrix{Float64}:  0.961524  0.274721 -0.274721  0.961524eigenvalues:2-element Vector{Float64}:  3.0 -2.0julia> F.vectors * F.Schur * F.vectors'2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> t, z, vals = F; # destructuring via iterationjulia> t == F.T && z == F.Z && vals == F.valuestrue

GeneralizedSchur <: Factorization

Matrix factorization type of the generalized Schur factorization of two matricesA andB. This is the return type ofschur(_, _), the corresponding matrix factorization function.

IfF::GeneralizedSchur is the factorization object, the (quasi) triangular Schur factors can be obtained viaF.S andF.T, the left unitary/orthogonal Schur vectors viaF.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.α./F.β.

Iterating the decomposition produces the componentsF.S,F.T,F.Q,F.Z,F.α, andF.β.

schur(A) -> F::Schur

Computes the Schur factorization of the matrixA. The (quasi) triangular Schur factor can be obtained from theSchur objectF with eitherF.Schur orF.T and the orthogonal/unitary 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.

For realA, the Schur factorization is "quasitriangular", which means that it is upper-triangular except with 2×2 diagonal blocks for any conjugate pair of complex eigenvalues; this allows the factorization to be purely real even when there are complex eigenvalues. To obtain the (complex) purely upper-triangular Schur factorization from a real quasitriangular factorization, you can useSchur{Complex}(schur(A)).

Iterating the decomposition produces the componentsF.T,F.Z, andF.values.

Examples

julia> A = [5. 7.; -2. -4.]2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> F = schur(A)Schur{Float64, Matrix{Float64}, Vector{Float64}}T factor:2×2 Matrix{Float64}: 3.0   9.0 0.0  -2.0Z factor:2×2 Matrix{Float64}:  0.961524  0.274721 -0.274721  0.961524eigenvalues:2-element Vector{Float64}:  3.0 -2.0julia> F.vectors * F.Schur * F.vectors'2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> t, z, vals = F; # destructuring via iterationjulia> t == F.T && z == F.Z && vals == F.valuestrue

schur(A, B) -> F::GeneralizedSchur

Computes the Generalized Schur (or QZ) factorization of the matricesA 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.α./F.β.

Iterating the decomposition produces the componentsF.S,F.T,F.Q,F.Z,F.α, andF.β.

schur!(A) -> F::Schur

Same asschur but uses the input argumentA as workspace.

Examples

julia> A = [5. 7.; -2. -4.]2×2 Matrix{Float64}:  5.0   7.0 -2.0  -4.0julia> F = schur!(A)Schur{Float64, Matrix{Float64}, Vector{Float64}}T factor:2×2 Matrix{Float64}: 3.0   9.0 0.0  -2.0Z factor:2×2 Matrix{Float64}:  0.961524  0.274721 -0.274721  0.961524eigenvalues:2-element Vector{Float64}:  3.0 -2.0julia> A2×2 Matrix{Float64}: 3.0   9.0 0.0  -2.0

schur!(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

Same asschur but uses the input matricesA andB as workspace.

ordschur(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

Reorders the Schur factorizationF of a matrixA = Z*T*Z' according to the logical arrayselect returning the reordered factorizationF object. The selected eigenvalues appear in the leading diagonal ofF.Schur and the corresponding leading columns ofF.vectors form an orthogonal/unitary basis of the corresponding right invariant subspace. In the real case, a complex conjugate pair of eigenvalues must be either both included or both excluded viaselect.

ordschur(F::GeneralizedSchur, select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

Reorders the Generalized Schur factorizationF of a matrix pair(A, B) = (Q*S*Z', Q*T*Z') according to the logical arrayselect and returns a GeneralizedSchur objectF. The selected eigenvalues appear in the leading diagonal of bothF.S andF.T, and the left and right orthogonal/unitary Schur vectors are also reordered such that(A, B) = F.Q*(F.S, F.T)*F.Z' still holds and the generalized eigenvalues ofA andB can still be obtained withF.α./F.β.

ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

Same asordschur but overwrites the factorizationF.

ordschur!(F::GeneralizedSchur, select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

Same asordschur but overwrites the factorizationF.

SVD <: Factorization

Matrix factorization type of the singular value decomposition (SVD) of a matrixA. This is the return type ofsvd(_), the corresponding matrix factorization function.

IfF::SVD is the factorization object,U,S,V andVt can be obtained viaF.U,F.S,F.V andF.Vt, such thatA = U * Diagonal(S) * Vt. The singular values inS are sorted in descending order.

Iterating the decomposition produces the componentsU,S, andV.

Examples

julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]4×5 Matrix{Float64}: 1.0  0.0  0.0  0.0  2.0 0.0  0.0  3.0  0.0  0.0 0.0  0.0  0.0  0.0  0.0 0.0  2.0  0.0  0.0  0.0julia> F = svd(A)SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}U factor:4×4 Matrix{Float64}: 0.0  1.0   0.0  0.0 1.0  0.0   0.0  0.0 0.0  0.0   0.0  1.0 0.0  0.0  -1.0  0.0singular values:4-element Vector{Float64}: 3.0 2.23606797749979 2.0 0.0Vt factor:4×5 Matrix{Float64}: -0.0        0.0  1.0  -0.0  0.0  0.447214   0.0  0.0   0.0  0.894427  0.0       -1.0  0.0   0.0  0.0  0.0        0.0  0.0   1.0  0.0julia> F.U * Diagonal(F.S) * F.Vt4×5 Matrix{Float64}: 1.0  0.0  0.0  0.0  2.0 0.0  0.0  3.0  0.0  0.0 0.0  0.0  0.0  0.0  0.0 0.0  2.0  0.0  0.0  0.0julia> u, s, v = F; # destructuring via iterationjulia> u == F.U && s == F.S && v == F.Vtrue

GeneralizedSVD <: Factorization

Matrix factorization type of the generalized singular value decomposition (SVD) of two matricesA andB, such thatA = F.U*F.D1*F.R0*F.Q' andB = F.V*F.D2*F.R0*F.Q'. This is the return type ofsvd(_, _), the corresponding matrix factorization function.

For an M-by-N matrixA and P-by-N matrixB,

• U is a M-by-M orthogonal matrix,
• V is a P-by-P orthogonal matrix,
• Q is a N-by-N orthogonal matrix,
• D1 is a M-by-(K+L) diagonal matrix with 1s in the first K entries,
• D2 is a P-by-(K+L) matrix whose top right L-by-L block is diagonal,
• R0 is a (K+L)-by-N matrix whose rightmost (K+L)-by-(K+L) block is nonsingular upper block triangular,

K+L is the effective numerical rank of the matrix[A; B].

Iterating the decomposition produces the componentsU,V,Q,D1,D2, andR0.

The entries ofF.D1 andF.D2 are related, as explained in the LAPACK documentation for thegeneralized SVD and thexGGSVD3 routine which is called underneath (in LAPACK 3.6.0 and newer).

Examples

julia> A = [1. 0.; 0. -1.]2×2 Matrix{Float64}: 1.0   0.0 0.0  -1.0julia> B = [0. 1.; 1. 0.]2×2 Matrix{Float64}: 0.0  1.0 1.0  0.0julia> F = svd(A, B)GeneralizedSVD{Float64, Matrix{Float64}, Float64, Vector{Float64}}U factor:2×2 Matrix{Float64}: 1.0  0.0 0.0  1.0V factor:2×2 Matrix{Float64}: -0.0  -1.0  1.0   0.0Q factor:2×2 Matrix{Float64}: 1.0  0.0 0.0  1.0D1 factor:2×2 Matrix{Float64}: 0.707107  0.0 0.0       0.707107D2 factor:2×2 Matrix{Float64}: 0.707107  0.0 0.0       0.707107R0 factor:2×2 Matrix{Float64}: 1.41421   0.0 0.0      -1.41421julia> F.U*F.D1*F.R0*F.Q'2×2 Matrix{Float64}: 1.0   0.0 0.0  -1.0julia> F.V*F.D2*F.R0*F.Q'2×2 Matrix{Float64}: -0.0  1.0  1.0  0.0

svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD

Compute the singular value decomposition (SVD) ofA 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 * Diagonal(S) * Vt. The algorithm producesVt and henceVt is more efficient to extract thanV. The singular values inS are sorted in descending order.

Iterating the decomposition produces the componentsU,S, andV.

Iffull = false (default), a "thin" SVD is returned. For an$M \times N$ matrixA, in the full factorizationU is$M \times M$ andV is$N \times N$, while in the thin factorizationU is$M \times K$ andV is$N \times K$, where$K = \min(M,N)$ is the number of singular values.

alg specifies which algorithm and LAPACK method to use for SVD:

• alg = LinearAlgebra.DivideAndConquer() (default): CallsLAPACK.gesdd!.
• alg = LinearAlgebra.QRIteration(): CallsLAPACK.gesvd! (typically slower but more accurate) .

Examples

julia> A = rand(4,3);julia> F = svd(A); # Store the Factorization Objectjulia> A ≈ F.U * Diagonal(F.S) * F.Vttruejulia> U, S, V = F; # destructuring via iterationjulia> A ≈ U * Diagonal(S) * V'truejulia> Uonly, = svd(A); # Store U onlyjulia> Uonly == Utrue

svd(A, B) -> GeneralizedSVD

Compute the generalized SVD ofA andB, returning aGeneralizedSVD factorization objectF such that[A;B] = [F.U * F.D1; F.V * F.D2] * F.R0 * F.Q'

• U is a M-by-M orthogonal matrix,
• V is a P-by-P orthogonal matrix,
• Q is a N-by-N orthogonal matrix,
• D1 is a M-by-(K+L) diagonal matrix with 1s in the first K entries,
• D2 is a P-by-(K+L) matrix whose top right L-by-L block is diagonal,
• R0 is a (K+L)-by-N matrix whose rightmost (K+L)-by-(K+L) block is nonsingular upper block triangular,

K+L is the effective numerical rank of the matrix[A; B].

Iterating the decomposition produces the componentsU,V,Q,D1,D2, andR0.

The generalized SVD is used in applications such as when one wants to compare how much belongs toA vs. how much belongs toB, as in human vs yeast genome, or signal vs noise, or between clusters vs within clusters. (See Edelman and Wang for discussion: https://arxiv.org/abs/1901.00485)

It decomposes[A; B] into[UC; VS]H, where[UC; VS] is a natural orthogonal basis for the column space of[A; B], andH = RQ' is a natural non-orthogonal basis for the rowspace of[A;B], where the top rows are most closely attributed to theA matrix, and the bottom to theB matrix. The multi-cosine/sine matricesC andS provide a multi-measure of how muchA vs how muchB, andU andV provide directions in which these are measured.

Examples

julia> A = randn(3,2); B=randn(4,2);julia> F = svd(A, B);julia> U,V,Q,C,S,R = F;julia> H = R*Q';julia> [A; B] ≈ [U*C; V*S]*Htruejulia> [A; B] ≈ [F.U*F.D1; F.V*F.D2]*F.R0*F.Q'truejulia> Uonly, = svd(A,B);julia> U == Uonlytrue

svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD

svd! is the same assvd, but saves space by overwriting the inputA, instead of creating a copy. See documentation ofsvd for details.

svd!(A, B) -> GeneralizedSVD

svd! is the same assvd, but modifies the argumentsA andB in-place, instead of making copies. See documentation ofsvd for details.

svdvals(A)

Return the singular values ofA in descending order.

Examples

julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]4×5 Matrix{Float64}: 1.0  0.0  0.0  0.0  2.0 0.0  0.0  3.0  0.0  0.0 0.0  0.0  0.0  0.0  0.0 0.0  2.0  0.0  0.0  0.0julia> svdvals(A)4-element Vector{Float64}: 3.0 2.23606797749979 2.0 0.0

svdvals(A, B)

Return the generalized singular values from the generalized singular value decomposition ofA andB. See alsosvd.

Examples

julia> A = [1. 0.; 0. -1.]2×2 Matrix{Float64}: 1.0   0.0 0.0  -1.0julia> B = [0. 1.; 1. 0.]2×2 Matrix{Float64}: 0.0  1.0 1.0  0.0julia> svdvals(A, B)2-element Vector{Float64}: 1.0 1.0

svdvals!(A)

Return the singular values ofA, saving space by overwriting the input. See alsosvdvals andsvd.

svdvals!(A, B)

Return the generalized singular values from the generalized singular value decomposition ofA andB, saving space by overwritingA andB. See alsosvd andsvdvals.

LinearAlgebra.Givens(i1,i2,c,s) -> G

A Givens rotation linear operator. The fieldsc ands represent the cosine and sine of the rotation angle, respectively. 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'.

See alsogivens.

givens(f::T, g::T, i1::Integer, i2::Integer) where {T} -> (G::Givens, r::T)

Computes the Givens rotationG and scalarr such that for any vectorx where

x[i1] = fx[i2] = g

the result of the multiplication

y = G*x

has the property that

y[i1] = ry[i2] = 0

See alsoLinearAlgebra.Givens.

givens(A::AbstractArray, i1::Integer, i2::Integer, j::Integer) -> (G::Givens, r)

Computes the Givens rotationG and scalarr such that the result of the multiplication

B = G*A

has the property that

B[i1,j] = rB[i2,j] = 0

See alsoLinearAlgebra.Givens.

givens(x::AbstractVector, i1::Integer, i2::Integer) -> (G::Givens, r)

Computes the Givens rotationG and scalarr such that the result of the multiplication

B = G*x

has the property that

B[i1] = rB[i2] = 0

See alsoLinearAlgebra.Givens.

triu(M, k::Integer = 0)

Return the upper triangle ofM starting from thekth superdiagonal.

Examples

julia> a = fill(1.0, (4,4))4×4 Matrix{Float64}: 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0julia> triu(a,3)4×4 Matrix{Float64}: 0.0  0.0  0.0  1.0 0.0  0.0  0.0  0.0 0.0  0.0  0.0  0.0 0.0  0.0  0.0  0.0julia> triu(a,-3)4×4 Matrix{Float64}: 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0

triu!(M)

Upper triangle of a matrix, overwritingM in the process. See alsotriu.

triu!(M, k::Integer)

Return the upper triangle ofM starting from thekth superdiagonal, overwritingM in the process.

Examples

julia> M = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5]5×5 Matrix{Int64}: 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5julia> triu!(M, 1)5×5 Matrix{Int64}: 0  2  3  4  5 0  0  3  4  5 0  0  0  4  5 0  0  0  0  5 0  0  0  0  0

tril(M, k::Integer = 0)

Return the lower triangle ofM starting from thekth superdiagonal.

Examples

julia> a = fill(1.0, (4,4))4×4 Matrix{Float64}: 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0julia> tril(a,3)4×4 Matrix{Float64}: 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0 1.0  1.0  1.0  1.0julia> tril(a,-3)4×4 Matrix{Float64}: 0.0  0.0  0.0  0.0 0.0  0.0  0.0  0.0 0.0  0.0  0.0  0.0 1.0  0.0  0.0  0.0

tril!(M)

Lower triangle of a matrix, overwritingM in the process. See alsotril.

tril!(M, k::Integer)

Return the lower triangle ofM starting from thekth superdiagonal, overwritingM in the process.

Examples

julia> M = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5]5×5 Matrix{Int64}: 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5julia> tril!(M, 2)5×5 Matrix{Int64}: 1  2  3  0  0 1  2  3  4  0 1  2  3  4  5 1  2  3  4  5 1  2  3  4  5

diagind(M::AbstractMatrix, k::Integer = 0, indstyle::IndexStyle = IndexLinear())diagind(M::AbstractMatrix, indstyle::IndexStyle = IndexLinear())

AnAbstractRange giving the indices of thekth diagonal of the matrixM. Optionally, an index style may be specified which determines the type of the range returned. Ifindstyle isa IndexLinear (default), this returns anAbstractRange{Integer}. On the other hand, ifindstyle isa IndexCartesian, this returns anAbstractRange{CartesianIndex{2}}.

Ifk is not provided, it is assumed to be0 (corresponding to the main diagonal).

See also:diag,diagm,Diagonal.

Examples

julia> A = [1 2 3; 4 5 6; 7 8 9]3×3 Matrix{Int64}: 1  2  3 4  5  6 7  8  9julia> diagind(A, -1)2:4:6julia> diagind(A, IndexCartesian())StepRangeLen(CartesianIndex(1, 1), CartesianIndex(1, 1), 3)

diag(M, k::Integer=0)

Thekth diagonal of a matrix, as a vector.

See alsodiagm,diagind,Diagonal,isdiag.

Examples

julia> A = [1 2 3; 4 5 6; 7 8 9]3×3 Matrix{Int64}: 1  2  3 4  5  6 7  8  9julia> diag(A,1)2-element Vector{Int64}: 2 6

diagm(kv::Pair{<:Integer,<:AbstractVector}...)diagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)

Construct a matrix fromPairs of diagonals and vectors. Vectorkv.second will be placed on thekv.first diagonal. By default the matrix is square and its size is inferred fromkv, but a non-square sizem×n (padded with zeros as needed) can be specified by passingm,n as the first arguments. For repeated diagonal indiceskv.first the values in the corresponding vectorskv.second will be added.

diagm constructs a full matrix; if you want storage-efficient versions with fast arithmetic, seeDiagonal,BidiagonalTridiagonal andSymTridiagonal.

Examples

julia> diagm(1 => [1,2,3])4×4 Matrix{Int64}: 0  1  0  0 0  0  2  0 0  0  0  3 0  0  0  0julia> diagm(1 => [1,2,3], -1 => [4,5])4×4 Matrix{Int64}: 0  1  0  0 4  0  2  0 0  5  0  3 0  0  0  0julia> diagm(1 => [1,2,3], 1 => [1,2,3])4×4 Matrix{Int64}: 0  2  0  0 0  0  4  0 0  0  0  6 0  0  0  0

diagm(v::AbstractVector)diagm(m::Integer, n::Integer, v::AbstractVector)

Construct a matrix with elements of the vector as diagonal elements. By default, the matrix is square and its size is given bylength(v), but a non-square sizem×n can be specified by passingm,n as the first arguments. The diagonal will be zero-padded if necessary.

Examples

julia> diagm([1,2,3])3×3 Matrix{Int64}: 1  0  0 0  2  0 0  0  3julia> diagm(4, 5, [1,2,3])4×5 Matrix{Int64}: 1  0  0  0  0 0  2  0  0  0 0  0  3  0  0 0  0  0  0  0

rank(::QRSparse{Tv,Ti}) -> Ti

Return the rank of the QR factorization

rank(S::SparseMatrixCSC{Tv,Ti}; [tol::Real]) -> Ti

Calculate rank ofS by calculating its QR factorization. Values smaller thantol are considered as zero. See SPQR's manual.

rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)rank(A::AbstractMatrix, rtol::Real)

Compute the numerical rank of a matrix by counting how many outputs ofsvdvals(A) are greater thanmax(atol, rtol*σ₁) whereσ₁ isA's largest calculated singular value.atol andrtol are the absolute and relative tolerances, respectively. The default relative tolerance isn*ϵ, wheren is the size of the smallest dimension ofA, andϵ is theeps of the element type ofA.

Examples

julia> rank(Matrix(I, 3, 3))3julia> rank(diagm(0 => [1, 0, 2]))2julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.1)2julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.00001)3julia> rank(diagm(0 => [1, 0.001, 2]), atol=1.5)1

rank(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}

Compute the numerical rank of the SVD objectS by counting how many singular values are greater thanmax(atol, rtol*σ₁) whereσ₁ is the largest calculated singular value. This is equivalent to the defaultrank(::AbstractMatrix) method except that it re-uses an existing SVD factorization.atol andrtol are the absolute and relative tolerances, respectively. The default relative tolerance isn*ϵ, wheren is the size of the smallest dimension of UΣV' andϵ is theeps of the element type ofS.

rank(A::QRPivoted{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}

Compute the numerical rank of the QR factorizationA by counting how many diagonal entries ofA.factors are greater thanmax(atol, rtol*Δ₁) whereΔ₁ is the largest calculated such entry. This is similar to therank(::AbstractMatrix) method insofar as it counts the number of (numerically) nonzero coefficients from a matrix factorization, although the default method uses an SVD instead of a QR factorization. Likerank(::SVD), this method also re-uses an existing matrix factorization.

Using a QR factorization to compute rank should typically produce the same result as using SVD, although it may be more prone to overestimating the rank in pathological cases where the matrix is ill-conditioned. It is also worth noting that it is generally faster to compute a QR factorization than it is to compute an SVD, so this method may be preferred when performance is a concern.

atol andrtol are the absolute and relative tolerances, respectively. The default relative tolerance isn*ϵ, wheren is the size of the smallest dimension ofA andϵ is theeps of the element type ofA.

norm(A, p::Real=2)

For any iterable containerA (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.

Thep-norm is defined as

\[\|A\|_p = \left( \sum_{i=1}^n | a_i | ^p \right)^{1/p}\]

with$a_i$ the entries of$A$,$| a_i |$ thenorm of$a_i$, and$n$ the length of$A$. Since thep-norm is computed using thenorms of the entries ofA, thep-norm of a vector of vectors is not compatible with the interpretation of it as a block vector in general ifp != 2.

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. IfA is a matrix andp=2, then this is equivalent to the Frobenius norm.

The second argumentp is not necessarily a part of the interface fornorm, i.e. a custom type may only implementnorm(A) without second argument.

Useopnorm to compute the operator norm of a matrix.

Examples

julia> v = [3, -2, 6]3-element Vector{Int64}:  3 -2  6julia> norm(v)7.0julia> norm(v, 1)11.0julia> norm(v, Inf)6.0julia> norm([1 2 3; 4 5 6; 7 8 9])16.881943016134134julia> norm([1 2 3 4 5 6 7 8 9])16.881943016134134julia> norm(1:9)16.881943016134134julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) != norm([v,v], 1)truejulia> norm(hcat(v,v), 2) == norm(vcat(v,v), 2) == norm([v,v], 2)truejulia> norm(hcat(v,v), Inf) == norm(vcat(v,v), Inf) != norm([v,v], Inf)true

norm(x::Number, p::Real=2)

For numbers, return$\left( |x|^p \right)^{1/p}$.

Examples

julia> norm(2, 1)2.0julia> norm(-2, 1)2.0julia> norm(2, 2)2.0julia> norm(-2, 2)2.0julia> norm(2, Inf)2.0julia> norm(-2, Inf)2.0

opnorm(A::AbstractMatrix, p::Real=2)

Compute the operator norm (or matrix norm) induced by the vectorp-norm, where valid values ofp are1,2, orInf. (Note that for sparse matrices,p=2 is currently not implemented.) Usenorm to compute the Frobenius norm.

Whenp=1, the operator norm is the maximum absolute column sum ofA:

\[\|A\|_1 = \max_{1 ≤ j ≤ n} \sum_{i=1}^m | a_{ij} |\]

with$a_{ij}$ the entries of$A$, and$m$ and$n$ its dimensions.

Whenp=2, the operator norm is the spectral norm, equal to the largest singular value ofA.

Whenp=Inf, the operator norm is the maximum absolute row sum ofA:

\[\|A\|_\infty = \max_{1 ≤ i ≤ m} \sum _{j=1}^n | a_{ij} |\]

Examples

julia> A = [1 -2 -3; 2 3 -1]2×3 Matrix{Int64}: 1  -2  -3 2   3  -1julia> opnorm(A, Inf)6.0julia> opnorm(A, 1)5.0

opnorm(x::Number, p::Real=2)

For numbers, return$\left( |x|^p \right)^{1/p}$. This is equivalent tonorm.

opnorm(A::Adjoint{<:Any,<:AbstractVector}, q::Real=2)opnorm(A::Transpose{<:Any,<:AbstractVector}, q::Real=2)

For Adjoint/Transpose-wrapped vectors, return the operator$q$-norm ofA, which is equivalent to thep-norm with valuep = q/(q-1). They coincide atp = q = 2. Usenorm to compute thep norm ofA as a vector.

The difference in norm between a vector space and its dual arises to preserve the relationship between duality and the dot product, and the result is consistent with the operatorp-norm of a1 × n matrix.

Examples

julia> v = [1; im];julia> vc = v';julia> opnorm(vc, 1)1.0julia> norm(vc, 1)2.0julia> norm(v, 1)2.0julia> opnorm(vc, 2)1.4142135623730951julia> norm(vc, 2)1.4142135623730951julia> norm(v, 2)1.4142135623730951julia> opnorm(vc, Inf)2.0julia> norm(vc, Inf)1.0julia> norm(v, Inf)1.0

normalize!(a::AbstractArray, p::Real=2)

Normalize the arraya in-place so that itsp-norm equals unity, i.e.norm(a, p) == 1. See alsonormalize andnorm.

normalize(a, p::Real=2)

Normalizea so that itsp-norm equals unity, i.e.norm(a, p) == 1. For scalars, this is similar to sign(a), except normalize(0) = NaN. See alsonormalize!,norm, andsign.

Examples

julia> a = [1,2,4];julia> b = normalize(a)3-element Vector{Float64}: 0.2182178902359924 0.4364357804719848 0.8728715609439696julia> norm(b)1.0julia> c = normalize(a, 1)3-element Vector{Float64}: 0.14285714285714285 0.2857142857142857 0.5714285714285714julia> norm(c, 1)1.0julia> a = [1 2 4 ; 1 2 4]2×3 Matrix{Int64}: 1  2  4 1  2  4julia> norm(a)6.48074069840786julia> normalize(a)2×3 Matrix{Float64}: 0.154303  0.308607  0.617213 0.154303  0.308607  0.617213julia> normalize(3, 1)1.0julia> normalize(-8, 1)-1.0julia> normalize(0, 1)NaN

cond(M, p::Real=2)

Condition number of the matrixM, computed using the operatorp-norm. Valid values forp are1,2 (default), orInf.

condskeel(M, [x, p::Real=Inf])

\[\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) = \frac{\left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \left\vert x \right\vert \right\Vert_p}{\left \Vert x \right \Vert_p}\]

Skeel condition number$\kappa_S$ of the matrixM, optionally with respect to the vectorx, as computed using the operatorp-norm.$\left\vert M \right\vert$ denotes the matrix of (entry wise) absolute values of$M$;$\left\vert M \right\vert_{ij} = \left\vert M_{ij} \right\vert$. Valid values forp are1,2 andInf (default).

This quantity is also known in the literature as the Bauer condition number, relative condition number, or componentwise relative condition number.

tr(M)

Matrix trace. Sums the diagonal elements ofM.

Examples

julia> A = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> tr(A)5

det(M)

Matrix determinant.

See also:logdet andlogabsdet.

Examples

julia> M = [1 0; 2 2]2×2 Matrix{Int64}: 1  0 2  2julia> det(M)2.0

Note that, in general,det computes a floating-point approximation of the determinant, even for integer matrices, typically via Gaussian elimination. Julia includes an exact algorithm for integer determinants (the Bareiss algorithm), but only uses it by default forBigInt matrices (since determinants quickly overflow any fixed integer precision):

julia> det(BigInt[1 0; 2 2]) # exact integer determinant2

logdet(M)

Logarithm of matrix determinant. Equivalent tolog(det(M)), but may provide increased accuracy and avoids overflow/underflow.

Examples

julia> M = [1 0; 2 2]2×2 Matrix{Int64}: 1  0 2  2julia> logdet(M)0.6931471805599453julia> logdet(Matrix(I, 3, 3))0.0

logabsdet(M)

Log of absolute value of matrix determinant. Equivalent to(log(abs(det(M))), sign(det(M))), but may provide increased accuracy and/or speed.

Examples

julia> A = [-1. 0.; 0. 1.]2×2 Matrix{Float64}: -1.0  0.0  0.0  1.0julia> det(A)-1.0julia> logabsdet(A)(0.0, -1.0)julia> B = [2. 0.; 0. 1.]2×2 Matrix{Float64}: 2.0  0.0 0.0  1.0julia> det(B)2.0julia> logabsdet(B)(0.6931471805599453, 1.0)

inv(M)

Matrix inverse. Computes matrixN such thatM * N = I, whereI is the identity matrix. Computed by solving the left-divisionN = M \ I.

ASingularException is thrown ifM fails numerical inversion.

Examples

julia> M = [2 5; 1 3]2×2 Matrix{Int64}: 2  5 1  3julia> N = inv(M)2×2 Matrix{Float64}:  3.0  -5.0 -1.0   2.0julia> M*N == N*M == Matrix(I, 2, 2)true

pinv(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)pinv(M, rtol::Real) = pinv(M; rtol=rtol) # to be deprecated in Julia 2.0

Computes the Moore-Penrose pseudoinverse.

For matricesM with floating point elements, it is convenient to compute the pseudoinverse by inverting only singular values greater thanmax(atol, rtol*σ₁) whereσ₁ is the largest singular value ofM.

The optimal choice of absolute (atol) and relative tolerance (rtol) varies both with the value ofM and the intended application of the pseudoinverse. The default relative tolerance isn*ϵ, wheren is the size of the smallest dimension ofM, andϵ is theeps of the element type ofM.

For inverting dense ill-conditioned matrices in a least-squares sense,rtol = sqrt(eps(real(float(oneunit(eltype(M)))))) is recommended.

For more information, see^[issue8859],^[B96],^[S84],^[KY88].

Examples

julia> M = [1.5 1.3; 1.2 1.9]2×2 Matrix{Float64}: 1.5  1.3 1.2  1.9julia> N = pinv(M)2×2 Matrix{Float64}:  1.47287   -1.00775 -0.930233   1.16279julia> M * N2×2 Matrix{Float64}: 1.0          -2.22045e-16 4.44089e-16   1.0

nullspace(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)nullspace(M, rtol::Real) = nullspace(M; rtol=rtol) # to be deprecated in Julia 2.0

Computes a basis for the nullspace ofM by including the singular vectors ofM whose singular values have magnitudes smaller thanmax(atol, rtol*σ₁), whereσ₁ isM's largest singular value.

By default, the relative tolerancertol isn*ϵ, wheren is the size of the smallest dimension ofM, andϵ is theeps of the element type ofM.

Examples

julia> M = [1 0 0; 0 1 0; 0 0 0]3×3 Matrix{Int64}: 1  0  0 0  1  0 0  0  0julia> nullspace(M)3×1 Matrix{Float64}: 0.0 0.0 1.0julia> nullspace(M, rtol=3)3×3 Matrix{Float64}: 0.0  1.0  0.0 1.0  0.0  0.0 0.0  0.0  1.0julia> nullspace(M, atol=0.95)3×1 Matrix{Float64}: 0.0 0.0 1.0

kron(A, B)

Computes the Kronecker product of two vectors, matrices or numbers.

For real vectorsv andw, the Kronecker product is related to the outer product bykron(v,w) == vec(w * transpose(v)) orw * transpose(v) == reshape(kron(v,w), (length(w), length(v))). Note how the ordering ofv andw differs on the left and right of these expressions (due to column-major storage). For complex vectors, the outer productw * v' also differs by conjugation ofv.

Examples

julia> A = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> B = [im 1; 1 -im]2×2 Matrix{Complex{Int64}}: 0+1im  1+0im 1+0im  0-1imjulia> kron(A, B)4×4 Matrix{Complex{Int64}}: 0+1im  1+0im  0+2im  2+0im 1+0im  0-1im  2+0im  0-2im 0+3im  3+0im  0+4im  4+0im 3+0im  0-3im  4+0im  0-4imjulia> v = [1, 2]; w = [3, 4, 5];julia> w*transpose(v)3×2 Matrix{Int64}: 3   6 4   8 5  10julia> reshape(kron(v,w), (length(w), length(v)))3×2 Matrix{Int64}: 3   6 4   8 5  10

kron!(C, A, B)

Computes the Kronecker product ofA andB and stores the result inC, overwriting the existing content ofC. This is the in-place version ofkron.

exp(A::AbstractMatrix)

Compute the matrix exponential ofA, defined by

\[e^A = \sum_{n=0}^{\infty} \frac{A^n}{n!}.\]

For symmetric or HermitianA, an eigendecomposition (eigen) is used, otherwise the scaling and squaring algorithm (see^[H05]) is chosen.

Examples

julia> A = Matrix(1.0I, 2, 2)2×2 Matrix{Float64}: 1.0  0.0 0.0  1.0julia> exp(A)2×2 Matrix{Float64}: 2.71828  0.0 0.0      2.71828

cis(A::AbstractMatrix)

More efficient method forexp(im*A) of square matrixA (especially ifA isHermitian or real-Symmetric).

See alsocispi,sincos,exp.

Examples

julia> cis([π 0; 0 π]) ≈ -Itrue

^(A::AbstractMatrix, p::Number)

Matrix power, equivalent to$\exp(p\log(A))$

Examples

julia> [1 2; 0 3]^32×2 Matrix{Int64}: 1  26 0  27

^(b::Number, A::AbstractMatrix)

Matrix exponential, equivalent to$\exp(\log(b)A)$.

Examples

julia> 2^[1 2; 0 3]2×2 Matrix{Float64}: 2.0  6.0 0.0  8.0julia> ℯ^[1 2; 0 3]2×2 Matrix{Float64}: 2.71828  17.3673 0.0      20.0855

log(A::AbstractMatrix)

IfA 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 nonprincipal matrix function is returned whenever possible.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used, ifA is triangular an improved version of the inverse scaling and squaring method is employed (see^[AH12] and^[AHR13]). IfA is real with no negative eigenvalues, then the real Schur form is computed. Otherwise, the complex Schur form is computed. Then the upper (quasi-)triangular algorithm in^[AHR13] is used on the upper (quasi-)triangular factor.

Examples

julia> A = Matrix(2.7182818*I, 2, 2)2×2 Matrix{Float64}: 2.71828  0.0 0.0      2.71828julia> log(A)2×2 Matrix{Float64}: 1.0  0.0 0.0  1.0

sqrt(x)

Return$\sqrt{x}$.

Throw aDomainError for negativeReal arguments. UseComplex negative arguments instead to obtain aComplex result.

The prefix operator√ is equivalent tosqrt.

See also:hypot.

Examples

julia> sqrt(big(81))9.0julia> sqrt(big(-81))ERROR: DomainError with -81.0:NaN result for non-NaN input.Stacktrace: [1] sqrt(::BigFloat) at ./mpfr.jl:501[...]julia> sqrt(big(complex(-81)))0.0 + 9.0imjulia> sqrt(-81 - 0.0im)  # -0.0im is below the branch cut0.0 - 9.0imjulia> .√(1:4)4-element Vector{Float64}: 1.0 1.4142135623730951 1.7320508075688772 2.0

sqrt(A::AbstractMatrix)

IfA 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.

IfA is real-symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the square root. For such matrices, eigenvalues λ that appear to be slightly negative due to roundoff errors are treated as if they were zero. More precisely, matrices with all eigenvalues≥ -rtol*(max |λ|) are treated as semidefinite (yielding a Hermitian square root), with negative eigenvalues taken to be zero.rtol is a keyword argument tosqrt (in the Hermitian/real-symmetric case only) that defaults to machine precision scaled bysize(A,1).

Otherwise, the square root is determined by means of the Björck-Hammarling method^[BH83], which computes the complex Schur form (schur) and then the complex square root of the triangular factor. If a real square root exists, then an extension of this method^[H87] that computes the real Schur form and then the real square root of the quasi-triangular factor is instead used.

Examples

julia> A = [4 0; 0 4]2×2 Matrix{Int64}: 4  0 0  4julia> sqrt(A)2×2 Matrix{Float64}: 2.0  0.0 0.0  2.0

cbrt(A::AbstractMatrix{<:Real})

Computes the real-valued cube root of a real-valued matrixA. IfT = cbrt(A), then we haveT*T*T ≈ A, see example given below.

IfA is symmetric, i.e., of typeHermOrSym{<:Real}, then (eigen) is used to find the cube root. Otherwise, a specialized version of the p-th root algorithm^[S03] is utilized, which exploits the real-valued Schur decomposition (schur) to compute the cube root.

Examples

julia> A = [0.927524 -0.15857; -1.3677 -1.01172]2×2 Matrix{Float64}:  0.927524  -0.15857 -1.3677    -1.01172julia> T = cbrt(A)2×2 Matrix{Float64}:  0.910077  -0.151019 -1.30257   -0.936818julia> T*T*T ≈ Atrue

cos(A::AbstractMatrix)

Compute the matrix cosine of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the cosine. Otherwise, the cosine is determined by callingexp.

Examples

julia> cos(fill(1.0, (2,2)))2×2 Matrix{Float64}:  0.291927  -0.708073 -0.708073   0.291927

sin(A::AbstractMatrix)

Compute the matrix sine of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the sine. Otherwise, the sine is determined by callingexp.

Examples

julia> sin(fill(1.0, (2,2)))2×2 Matrix{Float64}: 0.454649  0.454649 0.454649  0.454649

sincos(A::AbstractMatrix)

Compute the matrix sine and cosine of a square matrixA.

Examples

julia> S, C = sincos(fill(1.0, (2,2)));julia> S2×2 Matrix{Float64}: 0.454649  0.454649 0.454649  0.454649julia> C2×2 Matrix{Float64}:  0.291927  -0.708073 -0.708073   0.291927

tan(A::AbstractMatrix)

Compute the matrix tangent of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the tangent. Otherwise, the tangent is determined by callingexp.

Examples

julia> tan(fill(1.0, (2,2)))2×2 Matrix{Float64}: -1.09252  -1.09252 -1.09252  -1.09252

sec(A::AbstractMatrix)

Compute the matrix secant of a square matrixA.

csc(A::AbstractMatrix)

Compute the matrix cosecant of a square matrixA.

cot(A::AbstractMatrix)

Compute the matrix cotangent of a square matrixA.

cosh(A::AbstractMatrix)

Compute the matrix hyperbolic cosine of a square matrixA.

sinh(A::AbstractMatrix)

Compute the matrix hyperbolic sine of a square matrixA.

tanh(A::AbstractMatrix)

Compute the matrix hyperbolic tangent of a square matrixA.

sech(A::AbstractMatrix)

Compute the matrix hyperbolic secant of square matrixA.

csch(A::AbstractMatrix)

Compute the matrix hyperbolic cosecant of square matrixA.

coth(A::AbstractMatrix)

Compute the matrix hyperbolic cotangent of square matrixA.

acos(A::AbstractMatrix)

Compute the inverse matrix cosine of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the inverse cosine. Otherwise, the inverse cosine is determined by usinglog andsqrt. For the theory and logarithmic formulas used to compute this function, see^[AH16_1].

Examples

julia> acos(cos([0.5 0.1; -0.2 0.3]))2×2 Matrix{ComplexF64}:  0.5-8.32667e-17im  0.1+0.0im -0.2+2.63678e-16im  0.3-3.46945e-16im

asin(A::AbstractMatrix)

Compute the inverse matrix sine of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the inverse sine. Otherwise, the inverse sine is determined by usinglog andsqrt. For the theory and logarithmic formulas used to compute this function, see^[AH16_2].

Examples

julia> asin(sin([0.5 0.1; -0.2 0.3]))2×2 Matrix{ComplexF64}:  0.5-4.16334e-17im  0.1-5.55112e-17im -0.2+9.71445e-17im  0.3-1.249e-16im

atan(A::AbstractMatrix)

Compute the inverse matrix tangent of a square matrixA.

IfA is symmetric or Hermitian, its eigendecomposition (eigen) is used to compute the inverse tangent. Otherwise, the inverse tangent is determined by usinglog. For the theory and logarithmic formulas used to compute this function, see^[AH16_3].

Examples

julia> atan(tan([0.5 0.1; -0.2 0.3]))2×2 Matrix{ComplexF64}:  0.5+1.38778e-17im  0.1-2.77556e-17im -0.2+6.93889e-17im  0.3-4.16334e-17im

asec(A::AbstractMatrix)

Compute the inverse matrix secant ofA.

acsc(A::AbstractMatrix)

Compute the inverse matrix cosecant ofA.

acot(A::AbstractMatrix)

Compute the inverse matrix cotangent ofA.

acosh(A::AbstractMatrix)

Compute the inverse hyperbolic matrix cosine of a square matrixA. For the theory and logarithmic formulas used to compute this function, see^[AH16_4].

asinh(A::AbstractMatrix)

Compute the inverse hyperbolic matrix sine of a square matrixA. For the theory and logarithmic formulas used to compute this function, see^[AH16_5].

atanh(A::AbstractMatrix)

Compute the inverse hyperbolic matrix tangent of a square matrixA. For the theory and logarithmic formulas used to compute this function, see^[AH16_6].

asech(A::AbstractMatrix)

Compute the inverse matrix hyperbolic secant ofA.

acsch(A::AbstractMatrix)

Compute the inverse matrix hyperbolic cosecant ofA.

acoth(A::AbstractMatrix)

Compute the inverse matrix hyperbolic cotangent ofA.

lyap(A, C)

Computes the solutionX 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.

Examples

julia> A = [3. 4.; 5. 6]2×2 Matrix{Float64}: 3.0  4.0 5.0  6.0julia> B = [1. 1.; 1. 2.]2×2 Matrix{Float64}: 1.0  1.0 1.0  2.0julia> X = lyap(A, B)2×2 Matrix{Float64}:  0.5  -0.5 -0.5   0.25julia> A*X + X*A' ≈ -Btrue

sylvester(A, B, C)

Computes the solutionX to the Sylvester equationAX + XB + C = 0, whereA,B andC have compatible dimensions andA and-B have no eigenvalues with equal real part.

Examples

julia> A = [3. 4.; 5. 6]2×2 Matrix{Float64}: 3.0  4.0 5.0  6.0julia> B = [1. 1.; 1. 2.]2×2 Matrix{Float64}: 1.0  1.0 1.0  2.0julia> C = [1. 2.; -2. 1]2×2 Matrix{Float64}:  1.0  2.0 -2.0  1.0julia> X = sylvester(A, B, C)2×2 Matrix{Float64}: -4.46667   1.93333  3.73333  -1.8julia> A*X + X*B ≈ -Ctrue

issuccess(F::Factorization)

Test that a factorization of a matrix succeeded.

Examples

julia> F = cholesky([1 0; 0 1]);julia> issuccess(F)true

issuccess(F::LU; allowsingular = false)

Test that the LU factorization of a matrix succeeded. By default a factorization that produces a valid but rank-deficient U factor is considered a failure. This can be changed by passingallowsingular = true.

Examples

julia> F = lu([1 2; 1 2], check = false);julia> issuccess(F)falsejulia> issuccess(F, allowsingular = true)true

issymmetric(A) -> Bool

Test whether a matrix is symmetric.

Examples

julia> a = [1 2; 2 -1]2×2 Matrix{Int64}: 1   2 2  -1julia> issymmetric(a)truejulia> b = [1 im; -im 1]2×2 Matrix{Complex{Int64}}: 1+0im  0+1im 0-1im  1+0imjulia> issymmetric(b)false

isposdef(A) -> Bool

Test whether a matrix is positive definite (and Hermitian) by trying to perform a Cholesky factorization ofA.

See alsoisposdef!,cholesky.

Examples

julia> A = [1 2; 2 50]2×2 Matrix{Int64}: 1   2 2  50julia> isposdef(A)true

isposdef!(A) -> Bool

Test whether a matrix is positive definite (and Hermitian) by trying to perform a Cholesky factorization ofA, overwritingA in the process. See alsoisposdef.

Examples

julia> A = [1. 2.; 2. 50.];julia> isposdef!(A)truejulia> A2×2 Matrix{Float64}: 1.0  2.0 2.0  6.78233

istril(A::AbstractMatrix, k::Integer = 0) -> Bool

Test whetherA is lower triangular starting from thekth superdiagonal.

Examples

julia> a = [1 2; 2 -1]2×2 Matrix{Int64}: 1   2 2  -1julia> istril(a)falsejulia> istril(a, 1)truejulia> c = [1 1 0; 1 1 1; 1 1 1]3×3 Matrix{Int64}: 1  1  0 1  1  1 1  1  1julia> istril(c)falsejulia> istril(c, 1)true

istriu(A::AbstractMatrix, k::Integer = 0) -> Bool

Test whetherA is upper triangular starting from thekth superdiagonal.

Examples

julia> a = [1 2; 2 -1]2×2 Matrix{Int64}: 1   2 2  -1julia> istriu(a)falsejulia> istriu(a, -1)truejulia> c = [1 1 1; 1 1 1; 0 1 1]3×3 Matrix{Int64}: 1  1  1 1  1  1 0  1  1julia> istriu(c)falsejulia> istriu(c, -1)true

isdiag(A) -> Bool

Test whether a matrix is diagonal in the sense thatiszero(A[i,j]) is true unlessi == j. Note that it is not necessary forA to be square; if you would also like to check that, you need to check thatsize(A, 1) == size(A, 2).

Examples

julia> a = [1 2; 2 -1]2×2 Matrix{Int64}: 1   2 2  -1julia> isdiag(a)falsejulia> b = [im 0; 0 -im]2×2 Matrix{Complex{Int64}}: 0+1im  0+0im 0+0im  0-1imjulia> isdiag(b)truejulia> c = [1 0 0; 0 2 0]2×3 Matrix{Int64}: 1  0  0 0  2  0julia> isdiag(c)truejulia> d = [1 0 0; 0 2 3]2×3 Matrix{Int64}: 1  0  0 0  2  3julia> isdiag(d)false

ishermitian(A) -> Bool

Test whether a matrix is Hermitian.

Examples

julia> a = [1 2; 2 -1]2×2 Matrix{Int64}: 1   2 2  -1julia> ishermitian(a)truejulia> b = [1 im; -im 1]2×2 Matrix{Complex{Int64}}: 1+0im  0+1im 0-1im  1+0imjulia> ishermitian(b)true

transpose(A)

Lazy transpose. Mutating the returned object should appropriately mutateA. Often, but not always, yieldsTranspose(A), whereTranspose is a lazy transpose wrapper. Note that this operation is recursive.

This operation is intended for linear algebra usage - for general data manipulation seepermutedims, which is non-recursive.

Examples

julia> A = [3 2; 0 0]2×2 Matrix{Int64}: 3  2 0  0julia> B = transpose(A)2×2 transpose(::Matrix{Int64}) with eltype Int64: 3  0 2  0julia> B isa Transposetruejulia> transpose(B) === A # the transpose of a transpose unwraps the parenttruejulia> Transpose(B) # however, the constructor always wraps its argument2×2 transpose(transpose(::Matrix{Int64})) with eltype Int64: 3  2 0  0julia> B[1,2] = 4; # modifying B will modify A automaticallyjulia> A2×2 Matrix{Int64}: 3  2 4  0

For complex matrices, theadjoint operation is equivalent to a conjugate-transpose.

julia> A = reshape([Complex(x, x) for x in 1:4], 2, 2)2×2 Matrix{Complex{Int64}}: 1+1im  3+3im 2+2im  4+4imjulia> adjoint(A) == conj(transpose(A))true

Thetranspose of anAbstractVector is a row-vector:

julia> v = [1,2,3]3-element Vector{Int64}: 1 2 3julia> transpose(v) # returns a row-vector1×3 transpose(::Vector{Int64}) with eltype Int64: 1  2  3julia> transpose(v) * v # compute the dot product14

For a matrix of matrices, the individual blocks are recursively operated on:

julia> C = [1 3; 2 4]2×2 Matrix{Int64}: 1  3 2  4julia> D = reshape([C, 2C, 3C, 4C], 2, 2) # construct a block matrix2×2 Matrix{Matrix{Int64}}: [1 3; 2 4]  [3 9; 6 12] [2 6; 4 8]  [4 12; 8 16]julia> transpose(D) # blocks are recursively transposed2×2 transpose(::Matrix{Matrix{Int64}}) with eltype Transpose{Int64, Matrix{Int64}}: [1 2; 3 4]   [2 4; 6 8] [3 6; 9 12]  [4 8; 12 16]

transpose(F::Factorization)

Lazy transpose of the factorizationF. By default, returns aTransposeFactorization, except forFactorizations with realeltype, in which case returns anAdjointFactorization.

transpose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

Transpose the matrixA and stores it in the matrixX.size(X) must be equal tosize(transpose(A)). No additional memory is allocated other than resizing the rowval and nzval ofX, if needed.

Seehalfperm!

transpose!(dest,src)

Transpose arraysrc 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.

Examples

julia> A = [3+2im 9+2im; 8+7im  4+6im]2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 8+7im  4+6imjulia> B = zeros(Complex{Int64}, 2, 2)2×2 Matrix{Complex{Int64}}: 0+0im  0+0im 0+0im  0+0imjulia> transpose!(B, A);julia> B2×2 Matrix{Complex{Int64}}: 3+2im  8+7im 9+2im  4+6imjulia> A2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 8+7im  4+6im

Transpose

Lazy wrapper type for a transpose view of the underlying linear algebra object, usually anAbstractVector/AbstractMatrix. Usually, theTranspose constructor should not be called directly, usetranspose instead. To materialize the view usecopy.

This type is intended for linear algebra usage - for general data manipulation seepermutedims.

Examples

julia> A = [2 3; 0 0]2×2 Matrix{Int64}: 2  3 0  0julia> Transpose(A)2×2 transpose(::Matrix{Int64}) with eltype Int64: 2  0 3  0

TransposeFactorization

Lazy wrapper type for the transpose of the underlyingFactorization object. Usually, theTransposeFactorization constructor should not be called directly, usetranspose(:: Factorization) instead.

A'adjoint(A)

Lazy adjoint (conjugate transposition). Note thatadjoint is applied recursively to elements.

For number types,adjoint returns the complex conjugate, and therefore it is equivalent to the identity function for real numbers.

This operation is intended for linear algebra usage - for general data manipulation seepermutedims.

Examples

julia> A = [3+2im 9+2im; 0  0]2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 0+0im  0+0imjulia> B = A' # equivalently adjoint(A)2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}: 3-2im  0+0im 9-2im  0+0imjulia> B isa Adjointtruejulia> adjoint(B) === A # the adjoint of an adjoint unwraps the parenttruejulia> Adjoint(B) # however, the constructor always wraps its argument2×2 adjoint(adjoint(::Matrix{Complex{Int64}})) with eltype Complex{Int64}: 3+2im  9+2im 0+0im  0+0imjulia> B[1,2] = 4 + 5im; # modifying B will modify A automaticallyjulia> A2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 4-5im  0+0im

For real matrices, theadjoint operation is equivalent to atranspose.

julia> A = reshape([x for x in 1:4], 2, 2)2×2 Matrix{Int64}: 1  3 2  4julia> A'2×2 adjoint(::Matrix{Int64}) with eltype Int64: 1  2 3  4julia> adjoint(A) == transpose(A)true

The adjoint of anAbstractVector is a row-vector:

julia> x = [3, 4im]2-element Vector{Complex{Int64}}: 3 + 0im 0 + 4imjulia> x'1×2 adjoint(::Vector{Complex{Int64}}) with eltype Complex{Int64}: 3+0im  0-4imjulia> x'x # compute the dot product, equivalently x' * x25 + 0im

For a matrix of matrices, the individual blocks are recursively operated on:

julia> A = reshape([x + im*x for x in 1:4], 2, 2)2×2 Matrix{Complex{Int64}}: 1+1im  3+3im 2+2im  4+4imjulia> C = reshape([A, 2A, 3A, 4A], 2, 2)2×2 Matrix{Matrix{Complex{Int64}}}: [1+1im 3+3im; 2+2im 4+4im]  [3+3im 9+9im; 6+6im 12+12im] [2+2im 6+6im; 4+4im 8+8im]  [4+4im 12+12im; 8+8im 16+16im]julia> C'2×2 adjoint(::Matrix{Matrix{Complex{Int64}}}) with eltype Adjoint{Complex{Int64}, Matrix{Complex{Int64}}}: [1-1im 2-2im; 3-3im 4-4im]    [2-2im 4-4im; 6-6im 8-8im] [3-3im 6-6im; 9-9im 12-12im]  [4-4im 8-8im; 12-12im 16-16im]

adjoint(F::Factorization)

Lazy adjoint of the factorizationF. By default, returns anAdjointFactorization wrapper.

adjoint!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}

Transpose the matrixA and stores the adjoint of the elements in the matrixX.size(X) must be equal tosize(transpose(A)). No additional memory is allocated other than resizing the rowval and nzval ofX, if needed.

Seehalfperm!

adjoint!(dest,src)

Conjugate transpose arraysrc 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.

Examples

julia> A = [3+2im 9+2im; 8+7im  4+6im]2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 8+7im  4+6imjulia> B = zeros(Complex{Int64}, 2, 2)2×2 Matrix{Complex{Int64}}: 0+0im  0+0im 0+0im  0+0imjulia> adjoint!(B, A);julia> B2×2 Matrix{Complex{Int64}}: 3-2im  8-7im 9-2im  4-6imjulia> A2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 8+7im  4+6im

Adjoint

Lazy wrapper type for an adjoint view of the underlying linear algebra object, usually anAbstractVector/AbstractMatrix. Usually, theAdjoint constructor should not be called directly, useadjoint instead. To materialize the view usecopy.

This type is intended for linear algebra usage - for general data manipulation seepermutedims.

Examples

julia> A = [3+2im 9+2im; 0 0]2×2 Matrix{Complex{Int64}}: 3+2im  9+2im 0+0im  0+0imjulia> Adjoint(A)2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}: 3-2im  0+0im 9-2im  0+0im

AdjointFactorization

Lazy wrapper type for the adjoint of the underlyingFactorization object. Usually, theAdjointFactorization constructor should not be called directly, useadjoint(:: Factorization) instead.

copy(A::Transpose)copy(A::Adjoint)

Eagerly evaluate the lazy matrix transpose/adjoint. Note that the transposition is applied recursively to elements.

This operation is intended for linear algebra usage - for general data manipulation seepermutedims, which is non-recursive.

Examples

julia> A = [1 2im; -3im 4]2×2 Matrix{Complex{Int64}}: 1+0im  0+2im 0-3im  4+0imjulia> T = transpose(A)2×2 transpose(::Matrix{Complex{Int64}}) with eltype Complex{Int64}: 1+0im  0-3im 0+2im  4+0imjulia> copy(T)2×2 Matrix{Complex{Int64}}: 1+0im  0-3im 0+2im  4+0im

stride1(A) -> Int

Return the distance between successive array elements in dimension 1 in units of element size.

Examples

julia> A = [1,2,3,4]4-element Vector{Int64}: 1 2 3 4julia> LinearAlgebra.stride1(A)1julia> B = view(A, 2:2:4)2-element view(::Vector{Int64}, 2:2:4) with eltype Int64: 2 4julia> LinearAlgebra.stride1(B)2

LinearAlgebra.checksquare(A)

Check that a matrix is square, then return its common dimension. For multiple arguments, return a vector.

Examples

julia> A = fill(1, (4,4)); B = fill(1, (5,5));julia> LinearAlgebra.checksquare(A, B)2-element Vector{Int64}: 4 5

LinearAlgebra.peakflops(n::Integer=4096; eltype::DataType=Float64, ntrials::Integer=3, parallel::Bool=false)

peakflops computes the peak flop rate of the computer by using double precisiongemm!. By default, if no arguments are specified, it multiplies twoFloat64 matrices of sizen x n, wheren = 4096. 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 argumenteltype is provided,peakflops will construct matrices with elements of typeeltype for calculating the peak flop rate.

By default,peakflops will use the best timing from 3 trials. If thentrials keyword argument is provided,peakflops will use those many trials for picking the best timing.

If the keyword argumentparallel 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.

hermitianpart(A::AbstractMatrix, uplo::Symbol=:U) -> Hermitian

Return the Hermitian part of the square matrixA, defined as(A + A') / 2, as aHermitian matrix. For real matricesA, this is also known as the symmetric part ofA; it is also sometimes called the "operator real part". The optional argumentuplo controls the corresponding argument of theHermitian view. For real matrices, the latter is equivalent to aSymmetric view.

See alsohermitianpart! for the corresponding in-place operation.

hermitianpart!(A::AbstractMatrix, uplo::Symbol=:U) -> Hermitian

Overwrite the square matrixA in-place with its Hermitian part(A + A') / 2, and returnHermitian(A, uplo). For real matricesA, this is also known as the symmetric part ofA.

See alsohermitianpart for the corresponding out-of-place operation.

copy_adjoint!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int},                A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B

Efficiently copy elements of matrixA toB with adjunction as follows:

B[ir_dest, jr_dest] = adjoint(A)[jr_src, ir_src]

The elementsB[ir_dest, jr_dest] are overwritten. Furthermore, the index range parameters must satisfylength(ir_dest) == length(jr_src) andlength(jr_dest) == length(ir_src).

copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int},                A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B

Efficiently copy elements of matrixA toB with transposition as follows:

B[ir_dest, jr_dest] = transpose(A)[jr_src, ir_src]

The elementsB[ir_dest, jr_dest] are overwritten. Furthermore, the index range parameters must satisfylength(ir_dest) == length(jr_src) andlength(jr_dest) == length(ir_src).

copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange, jr_dest::AbstractUnitRange,                tM::AbstractChar,                M::AbstractVecOrMat, ir_src::AbstractUnitRange, jr_src::AbstractUnitRange) -> B

Efficiently copy elements of matrixM toB conditioned on the character parametertM as follows:

The elementsB[ir_dest, jr_dest] are overwritten. Furthermore, the index range parameters must satisfylength(ir_dest) == length(jr_src) andlength(jr_dest) == length(ir_src).

See alsocopyto! andcopy_adjoint!.

mul!(Y, A, B) -> Y

Calculates the matrix-matrix or matrix-vector product$A B$ and stores the result inY, overwriting the existing value ofY. Note thatY must not be aliased with eitherA orB.

Examples

julia> A = [1.0 2.0; 3.0 4.0]; B = [1.0 1.0; 1.0 1.0]; Y = similar(B);julia> mul!(Y, A, B) === Ytruejulia> Y2×2 Matrix{Float64}: 3.0  3.0 7.0  7.0julia> Y == A * Btrue

Implementation

For custom matrix and vector types, it is recommended to implement 5-argumentmul! rather than implementing 3-argumentmul! directly if possible.

mul!(C, A, B, α, β) -> C

Combined inplace matrix-matrix or matrix-vector multiply-add$A B α + C β$. The result is stored inC by overwriting it. Note thatC must not be aliased with eitherA orB.

Examples

julia> A = [1.0 2.0; 3.0 4.0]; B = [1.0 1.0; 1.0 1.0]; C = [1.0 2.0; 3.0 4.0];julia> α, β = 100.0, 10.0;julia> mul!(C, A, B, α, β) === Ctruejulia> C2×2 Matrix{Float64}: 310.0  320.0 730.0  740.0julia> C_original = [1.0 2.0; 3.0 4.0]; # A copy of the original value of Cjulia> C == A * B * α + C_original * βtrue

lmul!(a::Number, B::AbstractArray)

Scale an arrayB by a scalara overwritingB in-place. Usermul! to multiply scalar from right. The scaling operation respects the semantics of the multiplication* betweena and an element ofB. In particular, this also applies to multiplication involving non-finite numbers such asNaN and±Inf.

Examples

julia> B = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> lmul!(2, B)2×2 Matrix{Int64}: 2  4 6  8julia> lmul!(0.0, [Inf])1-element Vector{Float64}: NaN

lmul!(A, B)

Calculate the matrix-matrix product$AB$, overwritingB, and return the result. Here,A must be of special matrix type, like, e.g.,Diagonal,UpperTriangular orLowerTriangular, or of some orthogonal type, seeQR.

Examples

julia> B = [0 1; 1 0];julia> A = UpperTriangular([1 2; 0 3]);julia> lmul!(A, B);julia> B2×2 Matrix{Int64}: 2  1 3  0julia> B = [1.0 2.0; 3.0 4.0];julia> F = qr([0 1; -1 0]);julia> lmul!(F.Q, B)2×2 Matrix{Float64}: 3.0  4.0 1.0  2.0

rmul!(A::AbstractArray, b::Number)

Scale an arrayA by a scalarb overwritingA in-place. Uselmul! to multiply scalar from left. The scaling operation respects the semantics of the multiplication* between an element ofA andb. In particular, this also applies to multiplication involving non-finite numbers such asNaN and±Inf.

Examples

julia> A = [1 2; 3 4]2×2 Matrix{Int64}: 1  2 3  4julia> rmul!(A, 2)2×2 Matrix{Int64}: 2  4 6  8julia> rmul!([NaN], 0.0)1-element Vector{Float64}: NaN

rmul!(A, B)

Calculate the matrix-matrix product$AB$, overwritingA, and return the result. Here,B must be of special matrix type, like, e.g.,Diagonal,UpperTriangular orLowerTriangular, or of some orthogonal type, seeQR.

Examples

julia> A = [0 1; 1 0];julia> B = UpperTriangular([1 2; 0 3]);julia> rmul!(A, B);julia> A2×2 Matrix{Int64}: 0  3 1  2julia> A = [1.0 2.0; 3.0 4.0];julia> F = qr([0 1; -1 0]);julia> rmul!(A, F.Q)2×2 Matrix{Float64}: 2.0  1.0 4.0  3.0

ldiv!(Y, A, B) -> Y

ComputeA \ B in-place and store the result inY, returning the result.

The argumentA shouldnot be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced byfactorize orcholesky). The reason for this is that factorization itself is both expensive and typically allocates memory (although it can also be done in-place via, e.g.,lu!), and performance-critical situations requiringldiv! usually also require fine-grained control over the factorization ofA.

Examples

julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];julia> B = [1, 2.5, 3];julia> Y = similar(B); # use similar since there is no need to read from itjulia> ldiv!(Y, qr(A), B); # you may also try qr!(A) to further reduce allocationjulia> Y ≈ A \ Btrue

ldiv!(A, B)

ComputeA \ B in-place and overwritingB to store the result.

The argumentA shouldnot be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced byfactorize orcholesky). The reason for this is that factorization itself is both expensive and typically allocates memory (although it can also be done in-place via, e.g.,lu!), and performance-critical situations requiringldiv! usually also require fine-grained control over the factorization ofA.

Examples

julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];julia> B = [1, 2.5, 3];julia> B0 = copy(B); # a backup copy to facilitate testingjulia> ldiv!(lu(A), B); # you may also try lu!(A) to further reduce allocationjulia> B ≈ A \ B0true

ldiv!(a::Number, B::AbstractArray)

Divide each entry in an arrayB by a scalara overwritingB in-place. Userdiv! to divide scalar from right.

Examples

julia> B = [1.0 2.0; 3.0 4.0]2×2 Matrix{Float64}: 1.0  2.0 3.0  4.0julia> ldiv!(2.0, B)2×2 Matrix{Float64}: 0.5  1.0 1.5  2.0

ldiv!(A::Tridiagonal, B::AbstractVecOrMat) -> B

ComputeA \ B in-place by Gaussian elimination with partial pivoting and store the result inB, returning the result. In the process, the diagonals ofA are overwritten as well.

rdiv!(A, B)

ComputeA / B in-place and overwritingA to store the result.

The argumentB shouldnot be a matrix. Rather, instead of matrices it should be a factorization object (e.g. produced byfactorize orcholesky). The reason for this is that factorization itself is both expensive and typically allocates memory (although it can also be done in-place via, e.g.,lu!), and performance-critical situations requiringrdiv! usually also require fine-grained control over the factorization ofB.

rdiv!(A::AbstractArray, b::Number)

Divide each entry in an arrayA by a scalarb overwritingA in-place. Useldiv! to divide scalar from left.

Examples

julia> A = [1.0 2.0; 3.0 4.0]2×2 Matrix{Float64}: 1.0  2.0 3.0  4.0julia> rdiv!(A, 2.0)2×2 Matrix{Float64}: 0.5  1.0 1.5  2.0

Interface to BLAS subroutines.

set_num_threads(n::Integer)set_num_threads(::Nothing)

Set the number of threads the BLAS library should use equal ton::Integer.

Also acceptsnothing, in which case julia tries to guess the default number of threads. Passingnothing is discouraged and mainly exists for historical reasons.

get_num_threads()

Get the number of threads the BLAS library is using.

rot!(n, X, incx, Y, incy, c, s)

OverwriteX withc*X + s*Y andY with-conj(s)*X + c*Y for the firstn elements of arrayX with strideincx and firstn elements of arrayY with strideincy. ReturnsX andY.

scal!(n, a, X, incx)scal!(a, X)

OverwriteX witha*X for the firstn elements of arrayX with strideincx. ReturnsX.

Ifn andincx are not provided,length(X) andstride(X,1) are used.

scal(n, a, X, incx)scal(a, X)

ReturnX scaled bya for the firstn elements of arrayX with strideincx.

Ifn andincx are not provided,length(X) andstride(X,1) are used.

blascopy!(n, X, incx, Y, incy)

Copyn elements of arrayX with strideincx to arrayY with strideincy. ReturnsY.

dot(n, X, incx, Y, incy)

Dot product of two vectors consisting ofn elements of arrayX with strideincx andn elements of arrayY with strideincy.

Examples

julia> BLAS.dot(10, fill(1.0, 10), 1, fill(1.0, 20), 2)10.0

dotu(n, X, incx, Y, incy)

Dot function for two complex vectors consisting ofn elements of arrayX with strideincx andn elements of arrayY with strideincy.

Examples

julia> BLAS.dotu(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)-10.0 + 10.0im

dotc(n, X, incx, U, incy)

Dot function for two complex vectors, consisting ofn elements of arrayX with strideincx andn elements of arrayU with strideincy, conjugating the first vector.

Examples

julia> BLAS.dotc(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)10.0 - 10.0im

nrm2(n, X, incx)

2-norm of a vector consisting ofn elements of arrayX with strideincx.

Examples

julia> BLAS.nrm2(4, fill(1.0, 8), 2)2.0julia> BLAS.nrm2(1, fill(1.0, 8), 2)1.0

asum(n, X, incx)

Sum of the magnitudes of the firstn elements of arrayX with strideincx.

For a real array, the magnitude is the absolute value. For a complex array, the magnitude is the sum of the absolute value of the real part and the absolute value of the imaginary part.

Examples

julia> BLAS.asum(5, fill(1.0im, 10), 2)5.0julia> BLAS.asum(2, fill(1.0im, 10), 5)2.0

iamax(n, dx, incx)iamax(dx)

Find the index of the element ofdx with the maximum absolute value.n is the length ofdx, andincx is the stride. Ifn andincx are not provided, they assume default values ofn=length(dx) andincx=stride1(dx).

gemv!(tA, alpha, A, x, beta, y)

Update the vectory asalpha*A*x + beta*y oralpha*A'x + beta*y according totA.alpha andbeta are scalars. Return the updatedy.

gemv(tA, alpha, A, x)

Returnalpha*A*x oralpha*A'x according totA.alpha is a scalar.

gemv(tA, A, x)

ReturnA*x orA'x according totA.

gbmv!(trans, m, kl, ku, alpha, A, x, beta, y)

Update vectory asalpha*A*x + beta*y oralpha*A'*x + beta*y according totrans. The matrixA is a general band matrix of dimensionm bysize(A,2) withkl sub-diagonals andku super-diagonals.alpha andbeta are scalars. Return the updatedy.

gbmv(trans, m, kl, ku, alpha, A, x)

Returnalpha*A*x oralpha*A'*x according totrans. The matrixA is a general band matrix of dimensionm bysize(A,2) withkl sub-diagonals andku super-diagonals, andalpha is a scalar.

hemv!(ul, alpha, A, x, beta, y)

Update the vectory asalpha*A*x + beta*y.A is assumed to be Hermitian. Only theul triangle ofA is used.alpha andbeta are scalars. Return the updatedy.

hemv(ul, alpha, A, x)

Returnalpha*A*x.A is assumed to be Hermitian. Only theul triangle ofA is used.alpha is a scalar.

hemv(ul, A, x)

ReturnA*x.A is assumed to be Hermitian. Only theul triangle ofA is used.

hpmv!(uplo, α, AP, x, β, y)

Update vectory asα*A*x + β*y, whereA is a Hermitian matrix provided in packed formatAP.

Withuplo = 'U', the array AP must contain the upper triangular part of the Hermitian matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[1, 2] andA[2, 2] respectively, and so on.

Withuplo = 'L', the array AP must contain the lower triangular part of the Hermitian matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[2, 1] andA[3, 1] respectively, and so on.

The scalar inputsα andβ must be complex or real numbers.

The array inputsx,y andAP must all be ofComplexF32 orComplexF64 type.

Return the updatedy.

symv!(ul, alpha, A, x, beta, y)

Update the vectory asalpha*A*x + beta*y.A is assumed to be symmetric. Only theul triangle ofA is used.alpha andbeta are scalars. Return the updatedy.

symv(ul, alpha, A, x)

Returnalpha*A*x.A is assumed to be symmetric. Only theul triangle ofA is used.alpha is a scalar.

symv(ul, A, x)

ReturnA*x.A is assumed to be symmetric. Only theul triangle ofA is used.

sbmv!(uplo, k, alpha, A, x, beta, y)

Update vectory asalpha*A*x + beta*y whereA is 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 athttps://www.netlib.org/lapack/explore-html/. Only theuplo triangle ofA is used.

Return the updatedy.

sbmv(uplo, k, alpha, A, x)

Returnalpha*A*x whereA is a symmetric band matrix of ordersize(A,2) withk super-diagonals stored in the argumentA. Only theuplo triangle ofA is used.

sbmv(uplo, k, A, x)

ReturnA*x whereA is a symmetric band matrix of ordersize(A,2) withk super-diagonals stored in the argumentA. Only theuplo triangle ofA is used.

spmv!(uplo, α, AP, x, β, y)

Update vectory asα*A*x + β*y, whereA is a symmetric matrix provided in packed formatAP.

Withuplo = 'U', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[1, 2] andA[2, 2] respectively, and so on.

Withuplo = 'L', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[2, 1] andA[3, 1] respectively, and so on.

The scalar inputsα andβ must be real.

The array inputsx,y andAP must all be ofFloat32 orFloat64 type.

Return the updatedy.

trmv!(ul, tA, dA, A, b)

Returnop(A)*b, whereop is determined bytA. Only theul triangle ofA is used.dA determines if the diagonal values are read or are assumed to be all ones. The multiplication occurs in-place onb.

trmv(ul, tA, dA, A, b)

Returnop(A)*b, whereop is determined bytA. Only theul triangle ofA is used.dA determines if the diagonal values are read or are assumed to be all ones.

trsv!(ul, tA, dA, A, b)

Overwriteb with the solution toA*x = b or one of the other two variants determined bytA andul.dA determines if the diagonal values are read or are assumed to be all ones. Return the updatedb.

trsv(ul, tA, dA, A, b)

Return the solution toA*x = b or one of the other two variants determined bytA andul.dA determines if the diagonal values are read or are assumed to be all ones.

ger!(alpha, x, y, A)

Rank-1 update of the matrixA with vectorsx andy asalpha*x*y' + A.

geru!(alpha, x, y, A)

Rank-1 update of the matrixA with vectorsx andy asalpha*x*transpose(y) + A.

her!(uplo, alpha, x, A)

Methods for complex arrays only. Rank-1 update of the Hermitian matrixA with vectorx asalpha*x*x' + A.uplo controls which triangle ofA is updated. ReturnsA.

syr!(uplo, alpha, x, A)

Rank-1 update of the symmetric matrixA with vectorx asalpha*x*transpose(x) + A.uplo controls which triangle ofA is updated. ReturnsA.

spr!(uplo, α, x, AP)

Update matrixA asA+α*x*x', whereA is a symmetric matrix provided in packed formatAP andx is a vector.

Withuplo = 'U', the array AP must contain the upper triangular part of the symmetric matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[1, 2] andA[2, 2] respectively, and so on.

Withuplo = 'L', the array AP must contain the lower triangular part of the symmetric matrix packed sequentially, column by column, so thatAP[1] containsA[1, 1],AP[2] andAP[3] containA[2, 1] andA[3, 1] respectively, and so on.

The scalar inputα must be real.

The array inputsx andAP must all be ofFloat32 orFloat64 type. Return the updatedAP.

gemmt!(uplo, tA, tB, alpha, A, B, beta, C)

Update the lower or upper triangular part specified byuplo ofC asalpha*A*B + beta*C or the other variants according totA andtB. Return the updatedC.

gemmt(uplo, tA, tB, alpha, A, B)

Return the lower or upper triangular part specified byuplo ofA*B or the other three variants according totA andtB.

gemmt(uplo, tA, tB, A, B)

Return the lower or upper triangular part specified byuplo ofA*B or the other three variants according totA andtB.

gemm!(tA, tB, alpha, A, B, beta, C)

UpdateC asalpha*A*B + beta*C or the other three variants according totA andtB. Return the updatedC.

gemm(tA, tB, alpha, A, B)

Returnalpha*A*B or the other three variants according totA andtB.

gemm(tA, tB, A, B)

ReturnA*B or the other three variants according totA andtB.

symm!(side, ul, alpha, A, B, beta, C)

UpdateC asalpha*A*B + beta*C oralpha*B*A + beta*C according toside.A is assumed to be symmetric. Only theul triangle ofA is used. Return the updatedC.

symm(side, ul, alpha, A, B)

Returnalpha*A*B oralpha*B*A according toside.A is assumed to be symmetric. Only theul triangle ofA is used.