From 8f47d59ebb293c95908398fffce03238c9928141 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Thu, 8 Aug 2019 16:32:29 +0200 Subject: [PATCH 01/21] Add `dot` method for symmetric matrices --- stdlib/LinearAlgebra/src/symmetric.jl | 37 ++++++++++++++++++++++++++ stdlib/LinearAlgebra/test/symmetric.jl | 8 ++++++ 2 files changed, 45 insertions(+) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 1f3829d3727bf..840bbe8f26585 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -414,6 +414,43 @@ function triu(A::Symmetric, k::Integer=0) end end +function dot(A::Symmetric{T, Matrix{T}}, B::Symmetric{T, Matrix{T}}) where T + dotprod = zero(T) + m, n = size(A) + mB, nB = size(B) + + if m != mB || n != nB + throw(DimensionMismatch("A has dimensions ($m,$n) but B has dimensions ($mB,$nB)")) + end + + @inbounds if A.uplo == 'U' && B.uplo == 'U' + for j in 1:n + for i in 1:(j - 1) + dotprod += 2 * A.data[i, j] * conj(B.data[i, j]) + end + dotprod += A.data[j, j] * conj(B.data[j, j]) + end + elseif A.uplo == 'L' && B.uplo == 'L' + for j in 1:n + dotprod += A.data[j, j] * conj(B.data[j, j]) + for i in (j + 1):m + dotprod += 2 * A.data[i, j] * conj(B.data[i, j]) + end + end + else + if A.uplo == 'L' && B.uplo == 'U' + A, B = B, A + end + for j in 1:n + for i in 1:(j - 1) + dotprod += 2 * A.data[i, j] * conj(B.data[j, i]) + end + dotprod += A.data[j, j] * conj(B.data[j, j]) + end + end + return dotprod +end + (-)(A::Symmetric) = Symmetric(-A.data, sym_uplo(A.uplo)) (-)(A::Hermitian) = Hermitian(-A.data, sym_uplo(A.uplo)) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index ebdfc9c93c72f..c5f61a1186964 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -376,6 +376,14 @@ end @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), x) end end + + @testset "dot product of symmetric matrices" begin + rhs = dot(asym, asym) + @test dot(Symmetric(asym), Symmetric(asym)) ≈ rhs + @test dot(Symmetric(asym, :L), Symmetric(asym, :L)) ≈ rhs + @test dot(Symmetric(asym), Symmetric(asym, :L)) ≈ rhs + @test dot(Symmetric(asym, :L), Symmetric(asym)) ≈ rhs + @test_throws DimensionMismatch dot(Symmetric(asym), Symmetric(Matrix{eltya}(undef, n-1, n-1))) end end end From 7d6ba61b17fe41d3ca55c83f0c5d414d85e0c516 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Fri, 9 Aug 2019 17:08:24 +0200 Subject: [PATCH 02/21] Relax 'same type' requirement --- stdlib/LinearAlgebra/src/symmetric.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 840bbe8f26585..e008ad604408b 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -414,8 +414,8 @@ function triu(A::Symmetric, k::Integer=0) end end -function dot(A::Symmetric{T, Matrix{T}}, B::Symmetric{T, Matrix{T}}) where T - dotprod = zero(T) +function dot(A::Symmetric, B::Symmetric) + dotprod = zero(dot(first(A), first(B))) m, n = size(A) mB, nB = size(B) From 8aa2bb22908dbf83afc918ff4e64e98c17711dd4 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Fri, 9 Aug 2019 17:14:51 +0200 Subject: [PATCH 03/21] Avoid to store unnecessary matrix dimensions --- stdlib/LinearAlgebra/src/symmetric.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index e008ad604408b..e20122e465a70 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -416,11 +416,10 @@ end function dot(A::Symmetric, B::Symmetric) dotprod = zero(dot(first(A), first(B))) - m, n = size(A) - mB, nB = size(B) - if m != mB || n != nB - throw(DimensionMismatch("A has dimensions ($m,$n) but B has dimensions ($mB,$nB)")) + n = size(A, 2) + if n != size(B, 2) + throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) end @inbounds if A.uplo == 'U' && B.uplo == 'U' @@ -433,7 +432,7 @@ function dot(A::Symmetric, B::Symmetric) elseif A.uplo == 'L' && B.uplo == 'L' for j in 1:n dotprod += A.data[j, j] * conj(B.data[j, j]) - for i in (j + 1):m + for i in (j + 1):n dotprod += 2 * A.data[i, j] * conj(B.data[i, j]) end end From fa2e275fe7035ab9c0b08a9b122b7f645a408716 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Fri, 9 Aug 2019 17:20:30 +0200 Subject: [PATCH 04/21] Make use of inner dot product --- stdlib/LinearAlgebra/src/symmetric.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index e20122e465a70..b62b49ac74b05 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -425,15 +425,15 @@ function dot(A::Symmetric, B::Symmetric) @inbounds if A.uplo == 'U' && B.uplo == 'U' for j in 1:n for i in 1:(j - 1) - dotprod += 2 * A.data[i, j] * conj(B.data[i, j]) + dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end - dotprod += A.data[j, j] * conj(B.data[j, j]) + dotprod += dot(A.data[j, j], B.data[j, j]) end elseif A.uplo == 'L' && B.uplo == 'L' for j in 1:n - dotprod += A.data[j, j] * conj(B.data[j, j]) + dotprod += dot(A.data[j, j], B.data[j, j]) for i in (j + 1):n - dotprod += 2 * A.data[i, j] * conj(B.data[i, j]) + dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end end else @@ -442,9 +442,9 @@ function dot(A::Symmetric, B::Symmetric) end for j in 1:n for i in 1:(j - 1) - dotprod += 2 * A.data[i, j] * conj(B.data[j, i]) + dotprod += 2 * dot(A.data[i, j], B.data[j, i]) end - dotprod += A.data[j, j] * conj(B.data[j, j]) + dotprod += dot(A.data[j, j], B.data[j, j]) end end return dotprod From 181faeb1e95c018ba6a0d3b6ed207b1029a27ed5 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Fri, 9 Aug 2019 22:40:25 +0200 Subject: [PATCH 05/21] Fix test (don't use `undef` entries) --- stdlib/LinearAlgebra/test/symmetric.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index c5f61a1186964..43968dbdbd115 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -383,7 +383,7 @@ end @test dot(Symmetric(asym, :L), Symmetric(asym, :L)) ≈ rhs @test dot(Symmetric(asym), Symmetric(asym, :L)) ≈ rhs @test dot(Symmetric(asym, :L), Symmetric(asym)) ≈ rhs - @test_throws DimensionMismatch dot(Symmetric(asym), Symmetric(Matrix{eltya}(undef, n-1, n-1))) + @test_throws DimensionMismatch dot(Symmetric(asym), Symmetric(zeros(eltya, n-1, n-1))) end end end From 2751b404a41d5d36aea7aabcd27f30b507279d2d Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Fri, 9 Aug 2019 23:48:06 +0200 Subject: [PATCH 06/21] Initialize return value after dimension check --- stdlib/LinearAlgebra/src/symmetric.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index b62b49ac74b05..02aeb5d508dc8 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -415,13 +415,12 @@ function triu(A::Symmetric, k::Integer=0) end function dot(A::Symmetric, B::Symmetric) - dotprod = zero(dot(first(A), first(B))) - n = size(A, 2) if n != size(B, 2) throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) end + dotprod = zero(dot(first(A), first(B))) @inbounds if A.uplo == 'U' && B.uplo == 'U' for j in 1:n for i in 1:(j - 1) From bd2afe306ba9cb23796830c4bf918c3467c53a5e Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Sun, 11 Aug 2019 04:27:56 +0200 Subject: [PATCH 07/21] Avoid accesing A.data and B.data directly My previous code produced wrong results when A and B were symmetric block matrices, because they accessed the raw data in a way that only works if A and B are symmetric matrices of numbers. --- stdlib/LinearAlgebra/src/symmetric.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 02aeb5d508dc8..367e3a872695a 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -424,15 +424,15 @@ function dot(A::Symmetric, B::Symmetric) @inbounds if A.uplo == 'U' && B.uplo == 'U' for j in 1:n for i in 1:(j - 1) - dotprod += 2 * dot(A.data[i, j], B.data[i, j]) + dotprod += 2 * dot(A[i, j], B[i, j]) end - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) end elseif A.uplo == 'L' && B.uplo == 'L' for j in 1:n - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) for i in (j + 1):n - dotprod += 2 * dot(A.data[i, j], B.data[i, j]) + dotprod += 2 * dot(A[i, j], B[i, j]) end end else @@ -441,9 +441,9 @@ function dot(A::Symmetric, B::Symmetric) end for j in 1:n for i in 1:(j - 1) - dotprod += 2 * dot(A.data[i, j], B.data[j, i]) + dotprod += 2 * dot(A[i, j], B[j, i]) end - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) end end return dotprod From d3d1523b64b90752883a4e7faad5818a86f40b06 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Mon, 12 Aug 2019 06:23:22 +0200 Subject: [PATCH 08/21] Correct code for case `A.uplo=='L' && B.uplo=='U'` I wronlgy assumed that the dot product is commutative. It is not for complex matrices. Since we do not access `A.data` or `B.data`, we can use the same code for the cases where `A.uplo` and `B.uplo` do not coincide. --- stdlib/LinearAlgebra/src/symmetric.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 367e3a872695a..dab7158478b5f 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -436,9 +436,6 @@ function dot(A::Symmetric, B::Symmetric) end end else - if A.uplo == 'L' && B.uplo == 'U' - A, B = B, A - end for j in 1:n for i in 1:(j - 1) dotprod += 2 * dot(A[i, j], B[j, i]) From 69ab753e998ab1248b71eedddb0b16e257ae9e56 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Mon, 12 Aug 2019 05:15:26 +0200 Subject: [PATCH 09/21] Extend basic tests --- stdlib/LinearAlgebra/test/symmetric.jl | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 43968dbdbd115..e2179e82f38a4 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -378,12 +378,24 @@ end end @testset "dot product of symmetric matrices" begin - rhs = dot(asym, asym) - @test dot(Symmetric(asym), Symmetric(asym)) ≈ rhs - @test dot(Symmetric(asym, :L), Symmetric(asym, :L)) ≈ rhs - @test dot(Symmetric(asym), Symmetric(asym, :L)) ≈ rhs - @test dot(Symmetric(asym, :L), Symmetric(asym)) ≈ rhs - @test_throws DimensionMismatch dot(Symmetric(asym), Symmetric(zeros(eltya, n-1, n-1))) + symau = Symmetric(a, :U) + symal = Symmetric(a, :L) + msymau = Matrix(symau) + msymal = Matrix(symal) + @test_throws DimensionMismatch dot(symau, Symmetric(zeros(eltya, n-1, n-1))) + for eltyc in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) + creal = randn(n, n)/2 + cimag = randn(n, n)/2 + c = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(creal, cimag) : creal) + symcu = Symmetric(c, :U) + symcl = Symmetric(c, :L) + msymcu = Matrix(symcu) + msymcl = Matrix(symcl) + @test dot(symau, symcu) ≈ dot(msymau, msymcu) + @test dot(symau, symcl) ≈ dot(msymau, msymcl) + @test dot(symal, symcu) ≈ dot(msymal, msymcu) + @test dot(symal, symcl) ≈ dot(msymal, msymcl) + end end end end From 2943bf19f9bacc3a4b1e0f059ad9e25ed459d2ab Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Mon, 12 Aug 2019 07:05:17 +0200 Subject: [PATCH 10/21] Add tests for block matrices --- stdlib/LinearAlgebra/test/symmetric.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index e2179e82f38a4..2ebac20348cbe 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -396,6 +396,17 @@ end @test dot(symal, symcu) ≈ dot(msymal, msymcu) @test dot(symal, symcl) ≈ dot(msymal, msymcl) end + + # block matrices + blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3] + symblockmu = Symmetric(blockm, :U) + symblockml = Symmetric(blockm, :L) + msymblockmu = Matrix(symblockmu) + msymblockml = Matrix(symblockml) + @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu) + @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml) + @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu) + @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml) end end end From 1421c28c6f5e0a942b384f66387fc17af3e22408 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 13 Aug 2019 12:12:06 +0200 Subject: [PATCH 11/21] Restrict dot method to Symmetric matrices of numbers --- stdlib/LinearAlgebra/src/symmetric.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index dab7158478b5f..38a0f8911d8f1 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -414,7 +414,7 @@ function triu(A::Symmetric, k::Integer=0) end end -function dot(A::Symmetric, B::Symmetric) +function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) where {Ta,Tb<:Number} n = size(A, 2) if n != size(B, 2) throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) From e0b8574bb1080dc293abe82af73c35eb1c8f3486 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 13 Aug 2019 12:12:59 +0200 Subject: [PATCH 12/21] Remove tests for block matrices --- stdlib/LinearAlgebra/test/symmetric.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 2ebac20348cbe..e2179e82f38a4 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -396,17 +396,6 @@ end @test dot(symal, symcu) ≈ dot(msymal, msymcu) @test dot(symal, symcl) ≈ dot(msymal, msymcl) end - - # block matrices - blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3] - symblockmu = Symmetric(blockm, :U) - symblockml = Symmetric(blockm, :L) - msymblockmu = Matrix(symblockmu) - msymblockml = Matrix(symblockml) - @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu) - @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml) - @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu) - @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml) end end end From 367f3136b6c6dfe09ed7180a14fe704d638ad783 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 13 Aug 2019 12:42:01 +0200 Subject: [PATCH 13/21] Access data directly --- stdlib/LinearAlgebra/src/symmetric.jl | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 38a0f8911d8f1..519b067848146 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -424,23 +424,30 @@ function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) @inbounds if A.uplo == 'U' && B.uplo == 'U' for j in 1:n for i in 1:(j - 1) - dotprod += 2 * dot(A[i, j], B[i, j]) + dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end - dotprod += dot(A[j, j], B[j, j]) + dotprod += dot(A.data[j, j], B.data[j, j]) end elseif A.uplo == 'L' && B.uplo == 'L' for j in 1:n - dotprod += dot(A[j, j], B[j, j]) + dotprod += dot(A.data[j, j], B.data[j, j]) for i in (j + 1):n - dotprod += 2 * dot(A[i, j], B[i, j]) + dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end end - else + elseif A.uplo == 'L' && B.uplo == 'U' for j in 1:n for i in 1:(j - 1) - dotprod += 2 * dot(A[i, j], B[j, i]) + dotprod += 2 * dot(A.data[i, j], B.data[j, i]) + end + dotprod += dot(A.data[j, j], B.data[j, j]) + end + elseif A.uplo == 'U' && B.uplo == 'L' + for j in 1:n + dotprod += dot(A.data[j, j], B.data[j, j]) + for i in (j + 1):n + dotprod += 2 * dot(A.data[i,j], B.data[j, i]) end - dotprod += dot(A[j, j], B[j, j]) end end return dotprod From a93ce2a620ed7d2e0f19bc3e1c22dc8738da73a5 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 13 Aug 2019 13:06:50 +0200 Subject: [PATCH 14/21] Correct mixed cases --- stdlib/LinearAlgebra/src/symmetric.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 519b067848146..6543901375f83 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -435,14 +435,14 @@ function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end end - elseif A.uplo == 'L' && B.uplo == 'U' + elseif A.uplo == 'U' && B.uplo == 'L' for j in 1:n for i in 1:(j - 1) dotprod += 2 * dot(A.data[i, j], B.data[j, i]) end dotprod += dot(A.data[j, j], B.data[j, j]) end - elseif A.uplo == 'U' && B.uplo == 'L' + elseif A.uplo == 'L' && B.uplo == 'U' for j in 1:n dotprod += dot(A.data[j, j], B.data[j, j]) for i in (j + 1):n From 73171b1d09072f4a36b9cc2d5d1aeccfc39428ce Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 13 Aug 2019 14:18:21 +0200 Subject: [PATCH 15/21] Add missing space --- stdlib/LinearAlgebra/src/symmetric.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 6543901375f83..42d5d4de71a4f 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -446,7 +446,7 @@ function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) for j in 1:n dotprod += dot(A.data[j, j], B.data[j, j]) for i in (j + 1):n - dotprod += 2 * dot(A.data[i,j], B.data[j, i]) + dotprod += 2 * dot(A.data[i, j], B.data[j, i]) end end end From 4dde9fc4fd411ab3137939abbdd2fad10cf0a26f Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Sat, 17 Aug 2019 18:27:23 +0200 Subject: [PATCH 16/21] Make method generic again --- stdlib/LinearAlgebra/src/symmetric.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 42d5d4de71a4f..133607a788055 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -414,7 +414,7 @@ function triu(A::Symmetric, k::Integer=0) end end -function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) where {Ta,Tb<:Number} +function dot(A::Symmetric, B::Symmetric) n = size(A, 2) if n != size(B, 2) throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) @@ -426,11 +426,11 @@ function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) for i in 1:(j - 1) dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) end elseif A.uplo == 'L' && B.uplo == 'L' for j in 1:n - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) for i in (j + 1):n dotprod += 2 * dot(A.data[i, j], B.data[i, j]) end @@ -438,15 +438,15 @@ function dot(A::Symmetric{Ta,<:AbstractArray}, B::Symmetric{Tb,<:AbstractArray}) elseif A.uplo == 'U' && B.uplo == 'L' for j in 1:n for i in 1:(j - 1) - dotprod += 2 * dot(A.data[i, j], B.data[j, i]) + dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i])) end - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) end elseif A.uplo == 'L' && B.uplo == 'U' for j in 1:n - dotprod += dot(A.data[j, j], B.data[j, j]) + dotprod += dot(A[j, j], B[j, j]) for i in (j + 1):n - dotprod += 2 * dot(A.data[i, j], B.data[j, i]) + dotprod += 2 * dot(A.data[i, j], transpose(B.data[j, i])) end end end From bc771d4d47afd19978323f5a62f3d9f79c056743 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Sat, 17 Aug 2019 18:30:43 +0200 Subject: [PATCH 17/21] Add tests for block matrices --- stdlib/LinearAlgebra/test/symmetric.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index e2179e82f38a4..2ebac20348cbe 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -396,6 +396,17 @@ end @test dot(symal, symcu) ≈ dot(msymal, msymcu) @test dot(symal, symcl) ≈ dot(msymal, msymcl) end + + # block matrices + blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3] + symblockmu = Symmetric(blockm, :U) + symblockml = Symmetric(blockm, :L) + msymblockmu = Matrix(symblockmu) + msymblockml = Matrix(symblockml) + @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu) + @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml) + @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu) + @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml) end end end From 1c3f1f96949f390d55e360dec608ecaea8e03929 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Tue, 20 Aug 2019 19:14:26 +0200 Subject: [PATCH 18/21] Add `dot` method for Hermitian matrices --- stdlib/LinearAlgebra/src/symmetric.jl | 39 +++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 133607a788055..794cb810ca75d 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -453,6 +453,45 @@ function dot(A::Symmetric, B::Symmetric) return dotprod end +function dot(A::Hermitian, B::Hermitian) + n = size(A, 2) + if n != size(B, 2) + throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))")) + end + + dotprod = zero(dot(first(A), first(B))) + @inbounds if A.uplo == 'U' && B.uplo == 'U' + for j in 1:n + for i in 1:(j - 1) + dotprod += 2 * real(dot(A.data[i, j], B.data[i, j])) + end + dotprod += dot(A[j, j], B[j, j]) + end + elseif A.uplo == 'L' && B.uplo == 'L' + for j in 1:n + dotprod += dot(A[j, j], B[j, j]) + for i in (j + 1):n + dotprod += 2 * real(dot(A.data[i, j], B.data[i, j])) + end + end + elseif A.uplo == 'U' && B.uplo == 'L' + for j in 1:n + for i in 1:(j - 1) + dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i]))) + end + dotprod += dot(A[j, j], B[j, j]) + end + elseif A.uplo == 'L' && B.uplo == 'U' + for j in 1:n + dotprod += dot(A[j, j], B[j, j]) + for i in (j + 1):n + dotprod += 2 * real(dot(A.data[i, j], adjoint(B.data[j, i]))) + end + end + end + return dotprod +end + (-)(A::Symmetric) = Symmetric(-A.data, sym_uplo(A.uplo)) (-)(A::Hermitian) = Hermitian(-A.data, sym_uplo(A.uplo)) From 34fb1d81298d0cbd75c83dee1de8068815934ad8 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Wed, 21 Aug 2019 02:57:58 +0200 Subject: [PATCH 19/21] Add tests for Hermitian matrices --- stdlib/LinearAlgebra/test/symmetric.jl | 60 +++++++++++++------------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 2ebac20348cbe..cd09e72f5fa60 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -377,36 +377,38 @@ end end end - @testset "dot product of symmetric matrices" begin - symau = Symmetric(a, :U) - symal = Symmetric(a, :L) - msymau = Matrix(symau) - msymal = Matrix(symal) - @test_throws DimensionMismatch dot(symau, Symmetric(zeros(eltya, n-1, n-1))) - for eltyc in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) - creal = randn(n, n)/2 - cimag = randn(n, n)/2 - c = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(creal, cimag) : creal) - symcu = Symmetric(c, :U) - symcl = Symmetric(c, :L) - msymcu = Matrix(symcu) - msymcl = Matrix(symcl) - @test dot(symau, symcu) ≈ dot(msymau, msymcu) - @test dot(symau, symcl) ≈ dot(msymau, msymcl) - @test dot(symal, symcu) ≈ dot(msymal, msymcu) - @test dot(symal, symcl) ≈ dot(msymal, msymcl) - end + @testset "dot product of symmetric and Hermitian matrices" begin + for mtype in (Symmetric, Hermitian) + symau = mtype(a, :U) + symal = mtype(a, :L) + msymau = Matrix(symau) + msymal = Matrix(symal) + @test_throws DimensionMismatch dot(symau, mtype(zeros(eltya, n-1, n-1))) + for eltyc in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) + creal = randn(n, n)/2 + cimag = randn(n, n)/2 + c = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(creal, cimag) : creal) + symcu = mtype(c, :U) + symcl = mtype(c, :L) + msymcu = Matrix(symcu) + msymcl = Matrix(symcl) + @test dot(symau, symcu) ≈ dot(msymau, msymcu) + @test dot(symau, symcl) ≈ dot(msymau, msymcl) + @test dot(symal, symcu) ≈ dot(msymal, msymcu) + @test dot(symal, symcl) ≈ dot(msymal, msymcl) + end - # block matrices - blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3] - symblockmu = Symmetric(blockm, :U) - symblockml = Symmetric(blockm, :L) - msymblockmu = Matrix(symblockmu) - msymblockml = Matrix(symblockml) - @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu) - @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml) - @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu) - @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml) + # block matrices + blockm = [eltya == Int ? rand(1:7, 3, 3) : convert(Matrix{eltya}, eltya <: Complex ? complex.(randn(3, 3)/2, randn(3, 3)/2) : randn(3, 3)/2) for _ in 1:3, _ in 1:3] + symblockmu = mtype(blockm, :U) + symblockml = mtype(blockm, :L) + msymblockmu = Matrix(symblockmu) + msymblockml = Matrix(symblockml) + @test dot(symblockmu, symblockmu) ≈ dot(msymblockmu, msymblockmu) + @test dot(symblockmu, symblockml) ≈ dot(msymblockmu, msymblockml) + @test dot(symblockml, symblockmu) ≈ dot(msymblockml, msymblockmu) + @test dot(symblockml, symblockml) ≈ dot(msymblockml, msymblockml) + end end end end From 35328ca1419ef46e86b71f18601c5cc9d78b0096 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Wed, 21 Aug 2019 21:14:30 +0200 Subject: [PATCH 20/21] Replace `elseif` by `else` --- stdlib/LinearAlgebra/src/symmetric.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 794cb810ca75d..ca401f3c1c441 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -442,7 +442,7 @@ function dot(A::Symmetric, B::Symmetric) end dotprod += dot(A[j, j], B[j, j]) end - elseif A.uplo == 'L' && B.uplo == 'U' + else for j in 1:n dotprod += dot(A[j, j], B[j, j]) for i in (j + 1):n @@ -481,7 +481,7 @@ function dot(A::Hermitian, B::Hermitian) end dotprod += dot(A[j, j], B[j, j]) end - elseif A.uplo == 'L' && B.uplo == 'U' + else for j in 1:n dotprod += dot(A[j, j], B[j, j]) for i in (j + 1):n From d9a44372dc83c553da5ea0706a4d13ee377ee190 Mon Sep 17 00:00:00 2001 From: Alexander Seiler Date: Mon, 9 Sep 2019 13:47:22 +0200 Subject: [PATCH 21/21] Fix syntax error (introduced by resolving merge conflict) --- stdlib/LinearAlgebra/test/symmetric.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index cd09e72f5fa60..c8984be219fd3 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -376,6 +376,7 @@ end @test dot(x, Symmetric(aherm, uplo), x) ≈ dot(x, Symmetric(aherm, uplo)*x) ≈ dot(x, Matrix(Symmetric(aherm, uplo)), x) end end + end @testset "dot product of symmetric and Hermitian matrices" begin for mtype in (Symmetric, Hermitian)