From fa8d5cfe52ccf608202d298f98198db960a6d0ad Mon Sep 17 00:00:00 2001 From: Giuseppe Del Vecchio Del Vecchio Date: Sat, 28 Aug 2021 15:18:58 +0200 Subject: [PATCH 1/2] updated syntax, updated tests, still not working with matrices made by integers --- Pfaffian.jl | 36 ++++++++++++++++++----------------- test-Pfaffian.jl | 49 ++++++++++++++++++++++++------------------------ 2 files changed, 44 insertions(+), 41 deletions(-) diff --git a/Pfaffian.jl b/Pfaffian.jl index 251df5f..128be32 100644 --- a/Pfaffian.jl +++ b/Pfaffian.jl @@ -5,6 +5,8 @@ # (file, filename, data) = imp.find_module("pfaffian", [path]); # pfaffian = imp.load_module(name, file, filename, data) +using LinearAlgebra + function Householder(x::Vector{Float64})::Tuple{Vector{Float64}, Float64, Float64} @assert length(x) > 0 @@ -29,17 +31,17 @@ end Householder(x::Vector{Int}) = Householder(convert(Vector{Float64}, x)) -function Householder(x::Vector{Complex128})::Tuple{Vector{Complex128}, Float64, Complex128} +function Householder(x::Vector{ComplexF64})::Tuple{Vector{ComplexF64}, Float64, ComplexF64} @assert length(x) > 1 sigma = dot(x[2:end], x[2:end]) if sigma == 0 - return (zeros(Complex128, length(x)), 0., x[1]) + return (zeros(ComplexF64, length(x)), 0., x[1]) else norm_x = sqrt(abs2(x[1]) + sigma) v = copy(x) - phase = exp(im * atan2(imag(x[1]), real(x[1]))) + phase = exp(im * atan(imag(x[1]), real(x[1]))) v[1] += phase * norm_x normalize!(v) return (v, 2., -phase * norm_x) @@ -48,7 +50,7 @@ end function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where {T <: Number} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < 1e-12 + @assert maximum(abs, A + A') < 1e-12 n = size(A)[1] @@ -58,7 +60,7 @@ function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where end if calc_Q - Q = eye(T, n) + Q = diagm(ones(T, n)) end for i in 1:n-2 @@ -66,8 +68,8 @@ function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where @inbounds begin A[i+1, i] = alpha A[i, i+1] = -alpha - A[i+2:end, i] = 0. - A[i, i+2:end] = 0. + A[i+2:end, i] .= 0. + A[i, i+2:end] .= 0. end @inbounds w = @view(A[i+1:end, i+1:end]) * conj(v) @@ -96,9 +98,9 @@ function skew_tridiagonalize(A::Matrix{Int}; overwrite_A=false, calc_Q=true) return skew_tridiagonalize(convert(Matrix{Float64}, A), overwrite_A=overwrite_A, calc_Q=calc_Q) end -function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} +function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{ComplexF64}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < skewsymtol + @assert maximum(abs, A + A') < skewsymtol T = eltype(A) n = size(A)[1] @@ -115,7 +117,7 @@ function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A pf_val = T(1.) for k in 1:2:n-1 - kp = k + indmax(abs.(A[k+1:end, k])) + kp = k + argmax(abs.(A[k+1:end, k])) if kp != k+1 temp = copy(A[k+1, k:end]) @@ -151,9 +153,9 @@ function Pfaffian_LTL(A::Matrix{Int}; overwrite_A=false, skewsymtol=1e-6)::Float return Pfaffian_LTL(convert(Matrix{Float64}, A), overwrite_A=overwrite_A, skewsymtol=skewsymtol) end -function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} +function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{ComplexF64}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < skewsymtol + @assert maximum(abs, A + A') < skewsymtol T = eltype(A) n = size(A)[1] @@ -173,8 +175,8 @@ function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{Complex128}}; ove v, tau, alpha = Householder(A[i+1:end, i]) A[i+1, i] = alpha A[i, i+1] = -alpha - A[i+2:end, i] = 0. - A[i, i+2:end] = 0. + A[i+2:end, i] .= 0. + A[i, i+2:end] .= 0. w = @view(A[i+1:end, i+1:end]) * conj(v) BLAS.ger!(T(tau), v, conj(w), @view(A[i+1:end, i+1:end])) @@ -198,8 +200,8 @@ function Pfaffian_Householder(A::Matrix{Int}; overwrite_A=false, skewsymtol=1e-6 end # default to LTL for speed -Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = Pfaffian_LTL(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) -# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = Pfaffian_Householder(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) +Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = Pfaffian_LTL(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) +# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = Pfaffian_Householder(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) # Wimmer's version -# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = pfaffian[:pfaffian](A, overwrite_a=overwrite_A, method="P") +# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = pfaffian[:pfaffian](A, overwrite_a=overwrite_A, method="P") diff --git a/test-Pfaffian.jl b/test-Pfaffian.jl index a7ae21d..7cb225d 100644 --- a/test-Pfaffian.jl +++ b/test-Pfaffian.jl @@ -1,14 +1,13 @@ # test against Wimmer's PFAPACK Python implementation: https://michaelwimmer.org/downloads.html # (PFAPACK must be installed in the same directory as Pfaffian.jl and this script) - +# run the present script in a REPL opened in the same folder as pfapack folder downloaded above +# note that the Julia Python installation requires scipy to be installed. Do this using Conda.jl before +# running the test. include("Pfaffian.jl") using PyCall -@pyimport imp -path = dirname("pfapack/python/pfaffian.py") -name = basename("pfapack/python/pfaffian.py") -(file, filename, data) = imp.find_module("pfaffian", [path]); -pfaffian = imp.load_module(name, file, filename, data) +pf = pyimport("pfapack.pfaffian") + N = 100 @@ -21,7 +20,7 @@ println() x = rand(Float64, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_real](x) +v_w, tau_w, alpha_w = pf.householder_real(x) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @@ -29,23 +28,23 @@ v_w, tau_w, alpha_w = pfaffian[:householder_real](x) println() A = rand(Float64, (N, N)) -A = A - A.' +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](A) +T_w, Q_w = pf.skew_tridiagonalize(A) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println() @@ -57,31 +56,33 @@ println() x = rand(-3:3, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_real](convert(Vector{Float64}, x)) +v_w, tau_w, alpha_w = pf.householder_real(convert(Vector{Float64}, x)) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @show abs(alpha - alpha_w) println() -A = rand(-3:3, (N, N)) -A = A - A.' + + +A = rand(-10:10, (N, N)) +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](convert(Matrix{Float64}, A)) +T_w, Q_w = pf.skew_tridiagonalize(convert(Matrix{Float64}, A)) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println() @@ -90,34 +91,34 @@ println() println("Testing complex methods:") println() -x = rand(Complex128, N) +x = rand(ComplexF64, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_complex](x) +v_w, tau_w, alpha_w = pf.householder_complex(x) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @show abs(alpha - alpha_w) println() -A = rand(Complex128, (N, N)) -A = A - A.' +A = rand(ComplexF64, (N, N)) +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](A) +T_w, Q_w = pf.skew_tridiagonalize(A) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println() From c021c71c5f0e3394dafb26fa7b734c1939e12a41 Mon Sep 17 00:00:00 2001 From: Giuseppe Del Vecchio Del Vecchio Date: Sat, 28 Aug 2021 17:50:57 +0200 Subject: [PATCH 2/2] updated syntax julia 1.6.2, updated tests --- Pfaffian.jl | 36 ++++++++++++++++++----------------- test-Pfaffian.jl | 49 ++++++++++++++++++++++++------------------------ 2 files changed, 44 insertions(+), 41 deletions(-) diff --git a/Pfaffian.jl b/Pfaffian.jl index 251df5f..128be32 100644 --- a/Pfaffian.jl +++ b/Pfaffian.jl @@ -5,6 +5,8 @@ # (file, filename, data) = imp.find_module("pfaffian", [path]); # pfaffian = imp.load_module(name, file, filename, data) +using LinearAlgebra + function Householder(x::Vector{Float64})::Tuple{Vector{Float64}, Float64, Float64} @assert length(x) > 0 @@ -29,17 +31,17 @@ end Householder(x::Vector{Int}) = Householder(convert(Vector{Float64}, x)) -function Householder(x::Vector{Complex128})::Tuple{Vector{Complex128}, Float64, Complex128} +function Householder(x::Vector{ComplexF64})::Tuple{Vector{ComplexF64}, Float64, ComplexF64} @assert length(x) > 1 sigma = dot(x[2:end], x[2:end]) if sigma == 0 - return (zeros(Complex128, length(x)), 0., x[1]) + return (zeros(ComplexF64, length(x)), 0., x[1]) else norm_x = sqrt(abs2(x[1]) + sigma) v = copy(x) - phase = exp(im * atan2(imag(x[1]), real(x[1]))) + phase = exp(im * atan(imag(x[1]), real(x[1]))) v[1] += phase * norm_x normalize!(v) return (v, 2., -phase * norm_x) @@ -48,7 +50,7 @@ end function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where {T <: Number} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < 1e-12 + @assert maximum(abs, A + A') < 1e-12 n = size(A)[1] @@ -58,7 +60,7 @@ function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where end if calc_Q - Q = eye(T, n) + Q = diagm(ones(T, n)) end for i in 1:n-2 @@ -66,8 +68,8 @@ function skew_tridiagonalize(A::Matrix{T}; overwrite_A=false, calc_Q=true) where @inbounds begin A[i+1, i] = alpha A[i, i+1] = -alpha - A[i+2:end, i] = 0. - A[i, i+2:end] = 0. + A[i+2:end, i] .= 0. + A[i, i+2:end] .= 0. end @inbounds w = @view(A[i+1:end, i+1:end]) * conj(v) @@ -96,9 +98,9 @@ function skew_tridiagonalize(A::Matrix{Int}; overwrite_A=false, calc_Q=true) return skew_tridiagonalize(convert(Matrix{Float64}, A), overwrite_A=overwrite_A, calc_Q=calc_Q) end -function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} +function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{ComplexF64}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < skewsymtol + @assert maximum(abs, A + A') < skewsymtol T = eltype(A) n = size(A)[1] @@ -115,7 +117,7 @@ function Pfaffian_LTL(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A pf_val = T(1.) for k in 1:2:n-1 - kp = k + indmax(abs.(A[k+1:end, k])) + kp = k + argmax(abs.(A[k+1:end, k])) if kp != k+1 temp = copy(A[k+1, k:end]) @@ -151,9 +153,9 @@ function Pfaffian_LTL(A::Matrix{Int}; overwrite_A=false, skewsymtol=1e-6)::Float return Pfaffian_LTL(convert(Matrix{Float64}, A), overwrite_A=overwrite_A, skewsymtol=skewsymtol) end -function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{Complex128}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} +function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{ComplexF64}}; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} @assert size(A)[1] == size(A)[2] > 0 - @assert maximum(abs, A + A.') < skewsymtol + @assert maximum(abs, A + A') < skewsymtol T = eltype(A) n = size(A)[1] @@ -173,8 +175,8 @@ function Pfaffian_Householder(A::Union{Matrix{Float64}, Matrix{Complex128}}; ove v, tau, alpha = Householder(A[i+1:end, i]) A[i+1, i] = alpha A[i, i+1] = -alpha - A[i+2:end, i] = 0. - A[i, i+2:end] = 0. + A[i+2:end, i] .= 0. + A[i, i+2:end] .= 0. w = @view(A[i+1:end, i+1:end]) * conj(v) BLAS.ger!(T(tau), v, conj(w), @view(A[i+1:end, i+1:end])) @@ -198,8 +200,8 @@ function Pfaffian_Householder(A::Matrix{Int}; overwrite_A=false, skewsymtol=1e-6 end # default to LTL for speed -Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = Pfaffian_LTL(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) -# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = Pfaffian_Householder(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) +Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = Pfaffian_LTL(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) +# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = Pfaffian_Householder(A, overwrite_A=overwrite_A, skewsymtol=skewsymtol) # Wimmer's version -# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, Complex128} = pfaffian[:pfaffian](A, overwrite_a=overwrite_A, method="P") +# Pfaffian(A::Matrix; overwrite_A=false, skewsymtol=1e-6)::Union{Float64, ComplexF64} = pfaffian[:pfaffian](A, overwrite_a=overwrite_A, method="P") diff --git a/test-Pfaffian.jl b/test-Pfaffian.jl index a7ae21d..7cb225d 100644 --- a/test-Pfaffian.jl +++ b/test-Pfaffian.jl @@ -1,14 +1,13 @@ # test against Wimmer's PFAPACK Python implementation: https://michaelwimmer.org/downloads.html # (PFAPACK must be installed in the same directory as Pfaffian.jl and this script) - +# run the present script in a REPL opened in the same folder as pfapack folder downloaded above +# note that the Julia Python installation requires scipy to be installed. Do this using Conda.jl before +# running the test. include("Pfaffian.jl") using PyCall -@pyimport imp -path = dirname("pfapack/python/pfaffian.py") -name = basename("pfapack/python/pfaffian.py") -(file, filename, data) = imp.find_module("pfaffian", [path]); -pfaffian = imp.load_module(name, file, filename, data) +pf = pyimport("pfapack.pfaffian") + N = 100 @@ -21,7 +20,7 @@ println() x = rand(Float64, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_real](x) +v_w, tau_w, alpha_w = pf.householder_real(x) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @@ -29,23 +28,23 @@ v_w, tau_w, alpha_w = pfaffian[:householder_real](x) println() A = rand(Float64, (N, N)) -A = A - A.' +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](A) +T_w, Q_w = pf.skew_tridiagonalize(A) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println() @@ -57,31 +56,33 @@ println() x = rand(-3:3, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_real](convert(Vector{Float64}, x)) +v_w, tau_w, alpha_w = pf.householder_real(convert(Vector{Float64}, x)) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @show abs(alpha - alpha_w) println() -A = rand(-3:3, (N, N)) -A = A - A.' + + +A = rand(-10:10, (N, N)) +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](convert(Matrix{Float64}, A)) +T_w, Q_w = pf.skew_tridiagonalize(convert(Matrix{Float64}, A)) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println() @@ -90,34 +91,34 @@ println() println("Testing complex methods:") println() -x = rand(Complex128, N) +x = rand(ComplexF64, N) v, tau, alpha = Householder(x) -v_w, tau_w, alpha_w = pfaffian[:householder_complex](x) +v_w, tau_w, alpha_w = pf.householder_complex(x) @show maximum(abs, v - v_w) @show abs(tau - tau_w) @show abs(alpha - alpha_w) println() -A = rand(Complex128, (N, N)) -A = A - A.' +A = rand(ComplexF64, (N, N)) +A = A - A' T, Q = skew_tridiagonalize(A) -T_w, Q_w = pfaffian[:skew_tridiagonalize](A) +T_w, Q_w = pf.skew_tridiagonalize(A) @show maximum(abs, T - T_w) @show maximum(abs, Q - Q_w) println() pf_ltl = Pfaffian_LTL(A) -pf_ltl_w = pfaffian[:pfaffian](A, method="P") +pf_ltl_w = pf.pfaffian(A, method="P") @show abs((pf_ltl - pf_ltl_w)/pf_ltl_w) println() pf_h = Pfaffian_Householder(A) -pf_h_w = pfaffian[:pfaffian](A, method="H") +pf_h_w = pf.pfaffian(A, method="H") @show abs((pf_h - pf_h_w)/pf_h_w) println()