Skip to content

Commit 048c1a4

Browse files
lkdvoslkdvos
committed
test: add tests for squareroot, logarithm and power
Mirrors the `exponential` test layout: BLAS floats against the LinearAlgebra reference, hermitian fast paths, domain-error and tolerance-clamping edge cases, `Diagonal` inputs over generic float types, and BigFloat coverage via the GenericSchur and GenericLinearAlgebra test directories.
1 parent bbd0c06 commit 048c1a4

9 files changed

Lines changed: 516 additions & 0 deletions

File tree

‎test/genericlinearalgebra/logarithm.jl‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericLinearAlgebra
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "logarithm! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
X = randn(rng, T, m, m)
16+
A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I
17+
D, V = @constinferred eigh_full(A)
18+
19+
algs = (MatrixFunctionViaEigh(GLA_QRIteration()),)
20+
@testset "algorithm $alg" for alg in algs
21+
logA = @constinferred logarithm!(copy(A); alg)
22+
logA2 = @constinferred logarithm(A; alg)
23+
@test logA2 ≈ logA
24+
25+
Dlog, Vlog = @constinferred eigh_full(logA)
26+
@test diagview(Dlog) ≈ log.(diagview(D))
27+
end
28+
end

‎test/genericlinearalgebra/power.jl‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericLinearAlgebra
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "power! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
X = randn(rng, T, m, m)
16+
A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I
17+
D, V = @constinferred eigh_full(A)
18+
19+
algs = (MatrixFunctionViaEigh(GLA_QRIteration()),)
20+
@testset "algorithm $alg" for alg in algs
21+
@test power(A, 3, alg) ≈ A * A * A
22+
@test power(A, -1, alg) * A ≈ LinearAlgebra.I
23+
@test power(A, big"0.5", alg) ≈ squareroot(A; alg)
24+
25+
powA = @constinferred power(A, big"0.5", alg)
26+
Dpow, Vpow = @constinferred eigh_full(powA)
27+
@test diagview(Dpow) ≈ sqrt.(diagview(D))
28+
end
29+
end

‎test/genericlinearalgebra/squareroot.jl‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericLinearAlgebra
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "squareroot! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
X = randn(rng, T, m, m)
16+
A = project_hermitian!(X * X') + one(real(T)) * LinearAlgebra.I
17+
D, V = @constinferred eigh_full(A)
18+
19+
algs = (MatrixFunctionViaEigh(GLA_QRIteration()),)
20+
@testset "algorithm $alg" for alg in algs
21+
sqrtA = @constinferred squareroot!(copy(A); alg)
22+
sqrtA2 = @constinferred squareroot(A; alg)
23+
@test sqrtA2 ≈ sqrtA
24+
@test sqrtA * sqrtA ≈ A
25+
@test LinearAlgebra.ishermitian(sqrtA)
26+
27+
Dsqrt, Vsqrt = @constinferred eigh_full(sqrtA)
28+
@test diagview(Dsqrt) ≈ sqrt.(diagview(D))
29+
end
30+
end

‎test/genericschur/logarithm.jl‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericSchur
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "logarithm! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
# spectrum inside a disk around 1, away from the negative real axis and zero
16+
A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m))
17+
D, V = @constinferred eig_full(A)
18+
19+
logA = @constinferred logarithm(A)
20+
@test eltype(logA) == T
21+
@test exponential(logA) ≈ A
22+
23+
algs = (MatrixFunctionViaEig(GS_QRIteration()),)
24+
@testset "algorithm $alg" for alg in algs
25+
logA2 = @constinferred logarithm(A; alg)
26+
@test logA2 ≈ logA
27+
28+
Dlog, Vlog = @constinferred eig_full(logA2)
29+
by = x -> (real(x), imag(x))
30+
@test sort(diagview(Dlog); by) ≈ sort(log.(diagview(D)); by)
31+
end
32+
end

‎test/genericschur/power.jl‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericSchur
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "power! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
# spectrum inside a disk around 1, away from the negative real axis and zero
16+
A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m))
17+
18+
powA = @constinferred power(A, 3)
19+
@test powA ≈ A * A * A
20+
@test eltype(powA) == T
21+
22+
powA = @constinferred power(A, big"0.5")
23+
@test powA ≈ squareroot(A)
24+
25+
algs = (MatrixFunctionViaEig(GS_QRIteration()),)
26+
@testset "algorithm $alg" for alg in algs
27+
@test power(A, 3, alg) ≈ A * A * A
28+
@test power(A, -1, alg) * A ≈ LinearAlgebra.I
29+
@test power(A, big"0.5", alg) ≈ squareroot(A)
30+
end
31+
end

‎test/genericschur/squareroot.jl‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using GenericSchur
8+
9+
GenericFloats = (BigFloat, Complex{BigFloat})
10+
11+
@testset "squareroot! for T = $T" for T in GenericFloats
12+
rng = StableRNG(123)
13+
m = 24
14+
15+
# spectrum inside a disk around 1, away from the negative real axis
16+
A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m))
17+
D, V = @constinferred eig_full(A)
18+
19+
sqrtA = @constinferred squareroot(A)
20+
@test sqrtA * sqrtA ≈ A
21+
@test eltype(sqrtA) == T
22+
23+
algs = (MatrixFunctionViaEig(GS_QRIteration()),)
24+
@testset "algorithm $alg" for alg in algs
25+
sqrtA2 = @constinferred squareroot(A; alg)
26+
@test sqrtA2 ≈ sqrtA
27+
28+
Dsqrt, Vsqrt = @constinferred eig_full(sqrtA2)
29+
by = x -> (real(x), imag(x))
30+
@test sort(diagview(Dsqrt); by) ≈ sort(sqrt.(diagview(D)); by)
31+
end
32+
end

‎test/logarithm.jl‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
using MatrixAlgebraKit
2+
using Test
3+
using TestExtras
4+
using StableRNGs
5+
using MatrixAlgebraKit: diagview
6+
using LinearAlgebra
7+
using LinearAlgebra: exp
8+
9+
BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
10+
GenericFloats = (Float16, ComplexF16, BigFloat, Complex{BigFloat})
11+
12+
@testset "logarithm! for T = $T" for T in BLASFloats
13+
rng = StableRNG(123)
14+
m = 54
15+
16+
# spectrum inside a disk around 1, away from the negative real axis and zero
17+
A = LinearAlgebra.I + LinearAlgebra.normalize!(randn(rng, T, m, m))
18+
Ac = copy(A)
19+
logA = LinearAlgebra.log(A)
20+
21+
logA2 = @constinferred logarithm(A)
22+
@test logA ≈ logA2
23+
@test exp(logA2) ≈ A
24+
@test A == Ac
25+
26+
algs = (MatrixFunctionViaLA(), MatrixFunctionViaEig(LAPACK_Simple()))
27+
@testset "algorithm $alg" for alg in algs
28+
logA2 = @constinferred logarithm(A, alg)
29+
@test logA ≈ logA2
30+
@test A == Ac
31+
end
32+
33+
@test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration()))
34+
35+
# roundtrip with exponential
36+
@test logarithm(exponential(logA)) ≈ logA
37+
end
38+
39+
@testset "logarithm! for hermitian T = $T" for T in BLASFloats
40+
rng = StableRNG(123)
41+
m = 54
42+
43+
X = randn(rng, T, m, m)
44+
A = Matrix(LinearAlgebra.Hermitian(X * X' + one(real(T)) * LinearAlgebra.I))
45+
Ac = copy(A)
46+
logA = LinearAlgebra.log(LinearAlgebra.Hermitian(A))
47+
48+
algs = (
49+
MatrixFunctionViaLA(), MatrixFunctionViaEigh(LAPACK_QRIteration()),
50+
MatrixFunctionViaEig(LAPACK_Simple()),
51+
)
52+
@testset "algorithm $alg" for alg in algs
53+
logA2 = @constinferred logarithm(A, alg)
54+
@test logA ≈ logA2
55+
@test A == Ac
56+
end
57+
end
58+
59+
@testset "logarithm! domain handling for T = $T" for T in (Float32, Float64)
60+
rng = StableRNG(123)
61+
m = 4
62+
63+
X = randn(rng, T, m, m)
64+
V = Matrix(LinearAlgebra.qr(X).Q)
65+
A = V * LinearAlgebra.Diagonal(T[-1, 1, 2, 3]) * V'
66+
A = (A + A') / 2
67+
68+
# negative eigenvalue: DomainError for real input, complex principal value otherwise
69+
@test_throws DomainError logarithm(A)
70+
@test_throws DomainError logarithm(A; alg = MatrixFunctionViaEigh(LAPACK_QRIteration()))
71+
@test_throws DomainError logarithm(A; alg = MatrixFunctionViaEig(LAPACK_Simple()))
72+
73+
logA = @constinferred logarithm(complex.(A))
74+
@test exp(logA) ≈ A
75+
76+
# (numerically) zero eigenvalue: no logarithm exists
77+
Asing = V * LinearAlgebra.Diagonal(T[0, 1, 2, 3]) * V'
78+
Asing = (Asing + Asing') / 2
79+
@test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEigh(LAPACK_QRIteration()))
80+
@test_throws DomainError logarithm(Asing; alg = MatrixFunctionViaEig(LAPACK_Simple()))
81+
@test_throws DomainError logarithm(complex.(Asing); alg = MatrixFunctionViaEig(LAPACK_Simple()))
82+
end
83+
84+
@testset "logarithm! for Diagonal{$T}" for T in (BLASFloats..., GenericFloats...)
85+
rng = StableRNG(123)
86+
m = 54
87+
88+
data = T <: Real ? (abs.(randn(rng, T, m)) .+ one(T)) : (randn(rng, T, m) .+ 4 * one(T))
89+
A = Diagonal(data)
90+
Ac = copy(A)
91+
logA = LinearAlgebra.log(A)
92+
93+
logA2 = @constinferred logarithm(A)
94+
@test logA2 isa Diagonal
95+
@test logA ≈ logA2
96+
@test A == Ac
97+
98+
if T <: Real
99+
@test_throws DomainError logarithm(Diagonal(-data))
100+
end
101+
@test_throws DomainError logarithm(Diagonal(zeros(T, m)))
102+
end

0 commit comments

Comments
 (0)