julia-common 0.4.7-6 / usr / share / julia / test / linalg / dense.jl

# This file is a part of Julia. License is MIT: http://julialang.org/license

debug = false
using Base.Test

# Check that non-floats are correctly promoted
@test_approx_eq [1 0 0; 0 1 0]\[1,1] [1;1;0]

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

a = rand(n,n)
for elty in (Float32, Float64, Complex64, Complex128)
    a = convert(Matrix{elty}, a)
    # cond
    @test_approx_eq_eps cond(a, 1) 4.837320054554436e+02 0.01
    @test_approx_eq_eps cond(a, 2) 1.960057871514615e+02 0.01
    @test_approx_eq_eps cond(a, Inf) 3.757017682707787e+02 0.01
    @test_approx_eq_eps cond(a[:,1:5]) 10.233059337453463 0.01
    @test_throws ArgumentError cond(a,3)
end

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

for eltya in (Float32, Float64, Complex64, Complex128, Int)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(a2real, a2img) : a2real)
    asym = a'+a                  # symmetric indefinite
    apd  = a'*a                 # symmetric positive-definite
    ε = εa = eps(abs(float(one(eltya))))

    @test isposdef(apd,:U)
    for eltyb in (Float32, Float64, Complex64, Complex128, Int)
        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")

debug && println("Solve square general system of equations")
    κ = cond(a,1)
    x = a \ b
    @test_throws DimensionMismatch b'\b
    @test_throws DimensionMismatch b\b'
    @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit!
    @test zeros(eltya,n)\ones(eltya,n) ≈ zeros(eltya,n,1)\ones(eltya,n,1)


debug && println("Test nullspace")
        a15null = nullspace(a[:,1:n1]')
        @test rank([a[:,1:n1] a15null]) == 10
        @test_approx_eq_eps norm(a[:,1:n1]'a15null, Inf) zero(eltya) 300ε
        @test_approx_eq_eps norm(a15null'a[:,1:n1], Inf) zero(eltya) 400ε
        @test size(nullspace(b), 2) == 0
        @test nullspace(zeros(eltya,n)) == eye(eltya,1)
    end # for eltyb

debug && println("Test pinv")
    pinva15 = pinv(a[:,1:n1])
    @test_approx_eq a[:,1:n1]*pinva15*a[:,1:n1] a[:,1:n1]
    @test_approx_eq pinva15*a[:,1:n1]*pinva15 pinva15

    @test size(pinv(ones(eltya,0,0))) == (0,0)

debug && println("Lyapunov/Sylvester")
    let
        x = lyap(a, a2)
        @test_approx_eq -a2 a*x + x*a'
        x2 = sylvester(a[1:3, 1:3], a[4:n, 4:n], a2[1:3,4:n])
        @test_approx_eq -a2[1:3, 4:n] a[1:3, 1:3]*x2 + x2*a[4:n, 4:n]
    end

debug && println("Matrix square root")
    asq = sqrtm(a)
    @test_approx_eq asq*asq a
    asymsq = sqrtm(asym)
    @test_approx_eq asymsq*asymsq asym

debug && println("Numbers")
    α = rand(eltya)
    A = zeros(eltya,1,1)
    A[1,1] = α
    @test diagm(α) == A
    @test expm(α) == exp(α)

debug && println("Powers")
    if eltya <: AbstractFloat
        z = zero(eltya)
        t = convert(eltya,2)
        r = convert(eltya,2.5)
        @test a^z ≈ eye(a)
        @test a^t ≈ a^2
        @test eye(eltya,n,n)^r ≈ eye(a)
    end

debug && println("Factorize")
    d = rand(eltya,n)
    e = rand(eltya,n-1)
    e2 = rand(eltya,n-1)
    f = rand(eltya,n-2)
    A = diagm(d)
    @test factorize(A) == Diagonal(d)
    A += diagm(e,-1)
    @test factorize(A) == Bidiagonal(d,e,false)
    A += diagm(f,-2)
    @test factorize(A) == LowerTriangular(A)
    A = diagm(d) + diagm(e,1)
    @test factorize(A) == Bidiagonal(d,e,true)
    if eltya <: Real
        A = diagm(d) + diagm(e,1) + diagm(e,-1)
        @test full(factorize(A)) ≈ full(factorize(SymTridiagonal(d,e)))
        A = diagm(d) + diagm(e,1) + diagm(e,-1) + diagm(f,2) + diagm(f,-2)
        @test inv(factorize(A)) ≈ inv(factorize(Symmetric(A)))
    end
    A = diagm(d) + diagm(e,1) + diagm(e2,-1)
    @test full(factorize(A)) ≈ full(factorize(Tridiagonal(e2,d,e)))
    A = diagm(d) + diagm(e,1) + diagm(f,2)
    @test factorize(A) == UpperTriangular(A)
end # for eltya

# test triu/tril bounds checking
A = rand(5,7)
@test_throws(ArgumentError,triu(A,8))
@test_throws(ArgumentError,triu(A,-6))
@test_throws(ArgumentError,tril(A,8))
@test_throws(ArgumentError,tril(A,-6))

# Test gradient
for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128)
    if elty <: Real
        x = convert(Vector{elty}, [1:3;])
        g = ones(elty, 3)
    else
        x = convert(Vector{elty}, complex([1:3;], [1:3;]))
        g = convert(Vector{elty}, complex(ones(3), ones(3)))
    end
    @test_approx_eq gradient(x) g
    @test gradient(ones(elty,1)) == zeros(elty,1)
end

# Tests norms
nnorm = 10
mmat = 10
nmat = 8
for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
    debug && println(elty)

    ## Vector
    x = ones(elty,10)
    xs = sub(x,1:2:10)
    @test_approx_eq norm(x, -Inf) 1
    @test_approx_eq norm(x, -1) 1/10
    @test_approx_eq norm(x, 0) 10
    @test_approx_eq norm(x, 1) 10
    @test_approx_eq norm(x, 2) sqrt(10)
    @test_approx_eq norm(x, 3) cbrt(10)
    @test_approx_eq norm(x, Inf) 1
    if elty <: Base.LinAlg.BlasFloat
        @test_approx_eq norm(x, 1:4) 2
        @test_throws BoundsError norm(x,-1:4)
        @test_throws BoundsError norm(x,1:11)
    end
    @test_approx_eq norm(xs, -Inf) 1
    @test_approx_eq norm(xs, -1) 1/5
    @test_approx_eq norm(xs, 0) 5
    @test_approx_eq norm(xs, 1) 5
    @test_approx_eq norm(xs, 2) sqrt(5)
    @test_approx_eq norm(xs, 3) cbrt(5)
    @test_approx_eq norm(xs, Inf) 1

    # Issue #12552:
    if real(elty) <: AbstractFloat
        for p in [-Inf,-1,1,2,3,Inf]
            @test isnan(norm(elty[0,NaN],p))
            @test isnan(norm(elty[NaN,0],p))
        end
    end

    ## Number
    norm(x[1:1]) === norm(x[1], -Inf)
    norm(x[1:1]) === norm(x[1], 0)
    norm(x[1:1]) === norm(x[1], 1)
    norm(x[1:1]) === norm(x[1], 2)
    norm(x[1:1]) === norm(x[1], Inf)

    for i = 1:10
        x = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) :
            elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) :
            convert(Vector{elty}, randn(nnorm))
        xs = sub(x,1:2:nnorm)
        y = elty <: Integer ? convert(Vector{elty}, rand(1:10, nnorm)) :
            elty <: Complex ? convert(Vector{elty}, complex(randn(nnorm), randn(nnorm))) :
            convert(Vector{elty}, randn(nnorm))
        ys = sub(y,1:2:nnorm)
        α = elty <: Integer ? randn() :
            elty <: Complex ? convert(elty, complex(randn(),randn())) :
            convert(elty, randn())
        # Absolute homogeneity
        @test_approx_eq norm(α*x,-Inf) abs(α)*norm(x,-Inf)
        @test_approx_eq norm(α*x,-1) abs(α)*norm(x,-1)
        @test_approx_eq norm(α*x,1) abs(α)*norm(x,1)
        @test_approx_eq norm(α*x) abs(α)*norm(x) # two is default
        @test_approx_eq norm(α*x,3) abs(α)*norm(x,3)
        @test_approx_eq norm(α*x,Inf) abs(α)*norm(x,Inf)

        @test_approx_eq norm(α*xs,-Inf) abs(α)*norm(xs,-Inf)
        @test_approx_eq norm(α*xs,-1) abs(α)*norm(xs,-1)
        @test_approx_eq norm(α*xs,1) abs(α)*norm(xs,1)
        @test_approx_eq norm(α*xs) abs(α)*norm(xs) # two is default
        @test_approx_eq norm(α*xs,3) abs(α)*norm(xs,3)
        @test_approx_eq norm(α*xs,Inf) abs(α)*norm(xs,Inf)

        # Triangle inequality
        @test norm(x + y,1) <= norm(x,1) + norm(y,1)
        @test norm(x + y) <= norm(x) + norm(y) # two is default
        @test norm(x + y,3) <= norm(x,3) + norm(y,3)
        @test norm(x + y,Inf) <= norm(x,Inf) + norm(y,Inf)

        @test norm(xs + ys,1) <= norm(xs,1) + norm(ys,1)
        @test norm(xs + ys) <= norm(xs) + norm(ys) # two is default
        @test norm(xs + ys,3) <= norm(xs,3) + norm(ys,3)
        @test norm(xs + ys,Inf) <= norm(xs,Inf) + norm(ys,Inf)

        # Against vectorized versions
        @test_approx_eq norm(x,-Inf) minimum(abs(x))
        @test_approx_eq norm(x,-1) inv(sum(1./abs(x)))
        @test_approx_eq norm(x,0) sum(x .!= 0)
        @test_approx_eq norm(x,1) sum(abs(x))
        @test_approx_eq norm(x) sqrt(sum(abs2(x)))
        @test_approx_eq norm(x,3) cbrt(sum(abs(x).^3.))
        @test_approx_eq norm(x,Inf) maximum(abs(x))
    end
    ## Matrix (Operator)
        A = ones(elty,10,10)
        As = sub(A,1:5,1:5)
        @test_approx_eq norm(A, 1) 10
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(A, 2) 10
        @test_approx_eq norm(A, Inf) 10
        @test_approx_eq norm(As, 1) 5
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(As, 2) 5
        @test_approx_eq norm(As, Inf) 5

    for i = 1:10
        A = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) :
            elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) :
            convert(Matrix{elty}, randn(mmat, nmat))
        As = sub(A,1:nmat,1:nmat)
        B = elty <: Integer ? convert(Matrix{elty}, rand(1:10, mmat, nmat)) :
            elty <: Complex ? convert(Matrix{elty}, complex(randn(mmat, nmat), randn(mmat, nmat))) :
            convert(Matrix{elty}, randn(mmat, nmat))
        Bs = sub(B,1:nmat,1:nmat)
        α = elty <: Integer ? randn() :
            elty <: Complex ? convert(elty, complex(randn(),randn())) :
            convert(elty, randn())

        # Absolute homogeneity
        @test_approx_eq norm(α*A,1) abs(α)*norm(A,1)
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*A) abs(α)*norm(A) # two is default
        @test_approx_eq norm(α*A,Inf) abs(α)*norm(A,Inf)

        @test_approx_eq norm(α*As,1) abs(α)*norm(As,1)
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test_approx_eq norm(α*As) abs(α)*norm(As) # two is default
        @test_approx_eq norm(α*As,Inf) abs(α)*norm(As,Inf)

        # Triangle inequality
        @test norm(A + B,1) <= norm(A,1) + norm(B,1)
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(A + B) <= norm(A) + norm(B) # two is default
        @test norm(A + B,Inf) <= norm(A,Inf) + norm(B,Inf)

        @test norm(As + Bs,1) <= norm(As,1) + norm(Bs,1)
        elty <: Union{BigFloat,Complex{BigFloat},BigInt} || @test norm(As + Bs) <= norm(As) + norm(Bs) # two is default
        @test norm(As + Bs,Inf) <= norm(As,Inf) + norm(Bs,Inf)

        # vecnorm:
        for p = -2:3
            @test norm(reshape(A, length(A)), p) == vecnorm(A, p)
        end

        # issue #10234
        if elty <: AbstractFloat || elty <: Complex
            let z = zeros(elty, 100)
                z[1] = -Inf
                for p in [-2,-1.5,-1,-0.5,0.5,1,1.5,2,Inf]
                    @test norm(z, p) == (p < 0 ? 0 : Inf)
                    @test norm(elty[Inf],p) == Inf
                end
            end
        end
    end
end

# issue #10234
@test norm(Any[Inf],-2) == norm(Any[Inf],-1) == norm(Any[Inf],1) == norm(Any[Inf],1.5) == norm(Any[Inf],2) == norm(Any[Inf],Inf) == Inf

# overflow/underflow in norms:
@test_approx_eq norm(Float64[1e-300, 1], -3)*1e300 1
@test_approx_eq norm(Float64[1e300, 1], 3)*1e-300 1

## Issue related tests
# issue #1447
let
    A = [1.+0.im 0; 0 1]
    B = pinv(A)
    for i = 1:4
        @test_approx_eq A[i] B[i]
    end
end

# issue #2246
let
    A = [1 2 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]
    Asq = sqrtm(A)
    @test_approx_eq Asq*Asq A
    A2 = sub(A, 1:2, 1:2)
    A2sq = sqrtm(A2)
    @test_approx_eq A2sq*A2sq A2
end

let
    N = 3
    @test_approx_eq log(det(eye(N))) logdet(eye(N))
end

# issue #2637
let
  a = [1, 2, 3]
  b = [4, 5, 6]
  @test kron(eye(2),eye(2)) == eye(4)
  @test kron(a,b) == [4,5,6,8,10,12,12,15,18]
  @test kron(a',b') == [4 5 6 8 10 12 12 15 18]
  @test kron(a,b')  == [4 5 6; 8 10 12; 12 15 18]
  @test kron(a',b)  == [4 8 12; 5 10 15; 6 12 18]
  @test kron(a,eye(2)) == [1 0; 0 1; 2 0; 0 2; 3 0; 0 3]
  @test kron(eye(2),a) == [ 1 0; 2 0; 3 0; 0 1; 0 2; 0 3]
  @test kron(eye(2),2) == 2*eye(2)
  @test kron(3,eye(3)) == 3*eye(3)
  @test kron(a,2) == [2, 4, 6]
  @test kron(b',2) == [8 10 12]
end

# issue #4796
let
    dim=2
    S=zeros(Complex,dim,dim)
    T=zeros(Complex,dim,dim)
    T[:] = 1
    z = 2.5 + 1.5im
    S[1] = z
    @test S*T == [z z; 0 0]
end

# Matrix exponential
for elty in (Float32, Float64, Complex64, Complex128)
    A1  = convert(Matrix{elty}, [4 2 0; 1 4 1; 1 1 4])
    eA1 = convert(Matrix{elty}, [147.866622446369 127.781085523181  127.781085523182;
                                 183.765138646367 183.765138646366  163.679601723179;
                                 71.797032399996  91.8825693231832 111.968106246371]')
    @test_approx_eq expm(A1) eA1

    A2  = convert(Matrix{elty},
                  [29.87942128909879    0.7815750847907159 -2.289519314033932;
                   0.7815750847907159 25.72656945571064    8.680737820540137;
                   -2.289519314033932   8.680737820540137  34.39400925519054])
    eA2 = convert(Matrix{elty},
                  [  5496313853692458.0 -18231880972009236.0 -30475770808580460.0;
                   -18231880972009252.0  60605228702221920.0 101291842930249760.0;
                   -30475770808580480.0 101291842930249728.0 169294411240851968.0])
    @test_approx_eq expm(A2) eA2

    A3  = convert(Matrix{elty}, [-131 19 18;-390 56 54;-387 57 52])
    eA3 = convert(Matrix{elty}, [-1.50964415879218 -5.6325707998812  -4.934938326092;
                                 0.367879439109187 1.47151775849686  1.10363831732856;
                                 0.135335281175235 0.406005843524598 0.541341126763207]')
    @test_approx_eq expm(A3) eA3

    A4 = convert(Matrix{elty}, [0.25 0.25; 0 0])
    eA4 = convert(Matrix{elty}, [1.2840254166877416 0.2840254166877415; 0 1])
    @test expm(A4) ≈ eA4

    A5 = convert(Matrix{elty}, [0 0.02; 0 0])
    eA5 = convert(Matrix{elty}, [1 0.02; 0 1])
    @test expm(A5) ≈ eA5

    # Hessenberg
    @test_approx_eq hessfact(A1)[:H] convert(Matrix{elty},
                                             [4.000000000000000  -1.414213562373094  -1.414213562373095
                                              -1.414213562373095   4.999999999999996  -0.000000000000000
                                              0  -0.000000000000002   3.000000000000000])
end

for elty in (Float64, Complex{Float64})

    A4  = convert(Matrix{elty}, [1/2 1/3 1/4 1/5+eps();
                                 1/3 1/4 1/5 1/6;
                                 1/4 1/5 1/6 1/7;
                                 1/5 1/6 1/7 1/8])
    @test_approx_eq expm(logm(A4)) A4

    OLD_STDERR = STDERR
    rd,wr = redirect_stderr()
    A5  = convert(Matrix{elty}, [1 1 0 1; 0 1 1 0; 0 0 1 1; 1 0 0 1])
    @test_approx_eq expm(logm(A5)) A5
    s = readline(rd)
    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.")

    A6  = convert(Matrix{elty}, [-5 2 0 0 ; 1/2 -7 3 0; 0 1/3 -9 4; 0 0 1/4 -11])
    @test_approx_eq expm(logm(A6)) A6
    s = readline(rd)
    @test contains(s, "WARNING: Matrix with nonpositive real eigenvalues, a nonprincipal matrix logarithm will be returned.")
    redirect_stderr(OLD_STDERR)

    A7  = convert(Matrix{elty}, [1 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1])
    @test_approx_eq expm(logm(A7)) A7

end

A8 = 100 * [-1+1im 0 0 1e-8; 0 1 0 0; 0 0 1 0; 0 0 0 1];
@test_approx_eq expm(logm(A8)) A8

# issue 5116
A9  = [0 10 0 0; -1 0 0 0; 0 0 0 0; -2 0 0 0]
eA9 = [-0.999786072879326  -0.065407069689389   0.0   0.0
       0.006540706968939  -0.999786072879326   0.0   0.0
       0.0                 0.0                 1.0   0.0
       0.013081413937878  -3.999572145758650   0.0   1.0]
@test_approx_eq expm(A9) eA9

# issue 5116
A10  = [ 0. 0. 0. 0. ; 0. 0. -im 0.; 0. im 0. 0.; 0. 0. 0. 0.]
eA10 = [ 1.0+0.0im   0.0+0.0im                 0.0+0.0im                0.0+0.0im
        0.0+0.0im   1.543080634815244+0.0im   0.0-1.175201193643801im  0.0+0.0im
        0.0+0.0im   0.0+1.175201193643801im   1.543080634815243+0.0im  0.0+0.0im
        0.0+0.0im   0.0+0.0im                 0.0+0.0im                1.0+0.0im]
@test_approx_eq expm(A10) eA10

# issue #7181
A = [ 1  5  9
      2  6 10
      3  7 11
      4  8 12 ]
@test_throws ArgumentError diag(A, -5)
@test diag(A,-4) == []
@test diag(A,-3) == [4]
@test diag(A,-2) == [3,8]
@test diag(A,-1) == [2,7,12]
@test diag(A, 0) == [1,6,11]
@test diag(A, 1) == [5,10]
@test diag(A, 2) == [9]
@test diag(A, 3) == []
@test_throws ArgumentError diag(A, 4)

@test diag(zeros(0,0)) == []
@test_throws ArgumentError diag(zeros(0,0),1)
@test_throws ArgumentError diag(zeros(0,0),-1)

@test diag(zeros(1,0)) == []
@test diag(zeros(1,0),-1) == []
@test_throws ArgumentError diag(zeros(1,0),1)
@test_throws ArgumentError diag(zeros(1,0),-2)

@test diag(zeros(0,1)) == []
@test diag(zeros(0,1),1) == []
@test_throws ArgumentError diag(zeros(0,1),-1)
@test_throws ArgumentError diag(zeros(0,1),2)

## Least squares solutions
a = [ones(20) 1:20 1:20]
b = reshape(eye(8, 5), 20, 2)
for elty in (Float32, Float64, Complex64, Complex128)
    a = convert(Matrix{elty}, a)
    b = convert(Matrix{elty}, b)

    # Vector rhs
    x = a[:,1:2]\b[:,1]
    @test_approx_eq ((a[:,1:2]*x-b[:,1])'*(a[:,1:2]*x-b[:,1]))[1] convert(elty, 2.546616541353384)

    # Matrix rhs
    x = a[:,1:2]\b
    @test_approx_eq det((a[:,1:2]*x-b)'*(a[:,1:2]*x-b)) convert(elty, 4.437969924812031)

    # Rank deficient
    x = a\b
    @test_approx_eq det((a*x-b)'*(a*x-b)) convert(elty, 4.437969924812031)

    # Underdetermined minimum norm
    x = convert(Matrix{elty}, [1 0 0; 0 1 -1]) \ convert(Vector{elty}, [1,1])
    @test_approx_eq x convert(Vector{elty}, [1, 0.5, -0.5])

    # symmetric, positive definite
    @test_approx_eq inv(convert(Matrix{elty}, [6. 2; 2 1])) convert(Matrix{elty}, [0.5 -1; -1 3])

    # symmetric, indefinite
    @test_approx_eq inv(convert(Matrix{elty}, [1. 2; 2 1])) convert(Matrix{elty}, [-1. 2; 2 -1]/3)
end

# test ops on Numbers
for elty in [Float32,Float64,Complex64,Complex128]
    a = rand(elty)
    @test expm(a) == exp(a)
    @test isposdef(one(elty))
    @test sqrtm(a) == sqrt(a)
    @test logm(a) ≈ log(a)
    @test lyap(one(elty),a) == -a/2
end

# stride1
a = rand(10)
b = slice(a,2:2:10)
@test Base.LinAlg.stride1(b) == 2