Skip to content

1-update syntax julia1.6.2 #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

1-update syntax julia1.6.2 #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fa8d5cf
updated syntax, updated tests, still not working with matrices made b…
gdelvecchiodelvecchio Aug 28, 2021
c021c71
updated syntax julia 1.6.2, updated tests
gdelvecchiodelvecchio Aug 28, 2021
3acdb23
Merge branch '1-update-syntax-julia1.6.2' of https://github.com/gdelv…
gdelvecchiodelvecchio Aug 28, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 19 additions & 17 deletions Pfaffian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)
Expand All @@ -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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


n = size(A)[1]

Expand All @@ -58,16 +60,16 @@ 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
v::Vector{T}, tau::Float64, alpha::T = Householder(A[i+1:end, i])
@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)
Expand Down Expand Up @@ -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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


T = eltype(A)
n = size(A)[1]
Expand All @@ -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])
Expand Down Expand Up @@ -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
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


T = eltype(A)
n = size(A)[1]
Expand All @@ -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]))
Expand All @@ -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")
49 changes: 25 additions & 24 deletions test-Pfaffian.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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

Expand All @@ -21,31 +20,31 @@ 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)
@show abs(alpha - alpha_w)
println()

A = rand(Float64, (N, N))
A = A - A.'
A = A - A'
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


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()
Expand All @@ -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'
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


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()
Expand All @@ -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'
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change to transpose.


T, Q = skew_tridiagonalize(A)
T_w, Q_w = pfaffian[:skew_tridiagonalize](A)
T_w, Q_w = pf.skew_tridiagonalize(A)
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This currently throws

  File "/Users/mishmash/opt/anaconda3/envs/pfaffiantests37/lib/python3.7/site-packages/pfapack/pfaffian.py", line 93, in skew_tridiagonalize
    assert abs((A + A.T).max()) < 1e-14

since A' is conjugate transpose. Should probably use transpose throughout to make the antisymmetric matrices.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And also in the corresponding antisymmetric-checking asserts in the Julia code.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mishmash do you mean in the corresponding Julia code right? If so yes, you should use transpose(A) and not the ' as it does the conjugation too.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, in the Julia code we need to be careful about transpose vs. conjugate transpose throughout.


@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()