Added logabspfaffian(::AbstractMatrix{<:Complex}) (#139)

* Reduced allocations in pfaffian(::AbstractMatrix{<:Complex})

* Added logabspfaffian(::AbstractMatrix{<:Complex})

I realised that the algorithm used to compute
`pfaffian(::AbstractMatrix{<:Complex})` can be used to compute
the `logabspfaffian` with little modifications and decided to
add that too for completeness sake

* Added logabspfaffian(::AbstractMatrix{<:Complex})

I realised that the algorithm used to compute
`pfaffian(::AbstractMatrix{<:Complex})` can be used to compute
the `logabspfaffian` with little modifications and decided to
add that too for completeness sake

* Update src/pfaffian.jl

Co-authored-by: Steven G. Johnson <[email protected]>

* Revert "Added logabspfaffian(::AbstractMatrix{<:Complex})"

This reverts commit 0b49126.

* Removed superfluous @inbounds

* Added argmaxabs() for compatibility with julia 1.6

* Merged orign/main in and added argmaxabs() to _logabspfaffian()

---------

Co-authored-by: Steven G. Johnson <[email protected]>

Author: jlbosse

docs/src/pfaffian.md

@@ -49,7 +49,7 @@ julia> pfaffian!(BigInt[0 2 -7; -2 0 -8; 7 8 0])
 ```
 (Note that the Pfaffian is *always zero* for any *odd* size skew-symmetric matrix.)
 
-Since the computation of the pfaffian can easily overflow/underflow the maximum/minimum representable floating-point value, we also provide a function `logabspfaffian` for real skew-symmetric matrices (along with an in-place variant `logabspfaffian!`) that returns a tuple `(logpf, sign)` such
+Since the computation of the pfaffian can easily overflow/underflow the maximum/minimum representable floating-point value, we also provide a function `logabspfaffian` (along with an in-place variant `logabspfaffian!`) that returns a tuple `(logpf, sign)` such
 that the Pfaffian is `sign * exp(logpf)`.   (This is similar to the [`logabsdet`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.logabsdet) function in Julia's `LinearAlgebra` library to compute the log of the determinant.)
 
 ```jl

src/pfaffian.jl

@@ -106,21 +106,21 @@ function _pfaffian!(A::AbstractMatrix{<:Complex})
         # Pivot if neccessary
         @views kp = k + argmaxabs(A[k+1:end, k])
         if kp != k + 1
-            @inbounds @simd for l in k:n
+            @simd for l in k:n
                 A[k+1,l], A[kp,l] = A[kp,l], A[k+1,l]
             end
-            @inbounds @simd for l in k:n
+            @simd for l in k:n
                 A[l,k+1], A[l,kp] = A[l,kp], A[l,k+1]
             end
             pf *= -1
         end
 
         # Apply Gauss transformation and update `pf`
-        @inbounds if A[k+1,k] != zero(eltype(A))
-            @inbounds @views tauk .= A[k,k+2:end] ./ A[k,k+1]
-            pf *= @inbounds A[k,k+1]
+        if A[k+1,k] != zero(eltype(A))
+            @views tauk .= A[k,k+2:end] ./ A[k,k+1]
+            pf *= A[k,k+1]
             if k + 2 <= n
-                @inbounds for l1 in eachindex(tauk)
+                for l1 in eachindex(tauk)
                     @simd for l2 in eachindex(tauk)
                         @fastmath A[k+1+l2,k+1+l1] +=
                             tauk[l2] * A[k+1+l1,k+1] -
@@ -177,6 +177,7 @@ function _logabspfaffian!(A::SkewHermitian{<:Real})
     end
     return logpf, sgn*sign(det(H.Q))
 end
+
 function _logabspfaffian!(A::SkewHermTridiagonal{<:Real})
     n = size(A, 1)
     isodd(n) && return convert(eltype(A.ev), -Inf), zero(eltype(A.ev))
@@ -188,18 +189,61 @@ function _logabspfaffian!(A::SkewHermTridiagonal{<:Real})
     end
     return logpf, sgn
 end
+
+function _logabspfaffian!(A::AbstractMatrix{<:Complex})
+    n = size(A, 1)
+    isodd(n) && return convert(real(eltype(A)), -Inf), zero(eltype(A))
+    tau = Array{eltype(A)}(undef, n - 2)
+    logpf = phase = zero(real(eltype(A)))
+    @inbounds for k in 1:2:n-1
+        tauk = @view tau[k:end]
+
+        # Pivot if neccessary
+        @views kp = k + argmaxabs(A[k+1:end, k])
+        if kp != k + 1
+            @simd for l in k:n
+                A[k+1,l], A[kp,l] = A[kp,l], A[k+1,l]
+            end
+            @simd for l in k:n
+                A[l,k+1], A[l,kp] = A[l,kp], A[l,k+1]
+            end
+            phase += π
+        end
+
+        # Apply Gauss transformation and update `pf`
+        if A[k+1,k] != zero(eltype(A))
+            @views tauk .= A[k,k+2:end] ./ A[k,k+1]
+            logpf += log(abs(A[k,k+1]))
+            phase += angle(A[k,k+1])
+            if k + 2 <= n
+                for l1 in eachindex(tauk)
+                    @simd for l2 in eachindex(tauk)
+                        @fastmath A[k+1+l2,k+1+l1] +=
+                            tauk[l2] * A[k+1+l1,k+1] -
+                            tauk[l1] * A[k+1+l2,k+1]
+                    end
+                end
+            end
+        else
+            return convert(real(eltype(A)), -Inf), zero(eltype(A))
+        end
+    end
+
+    return logpf, cis(phase)
+end
+
 logabspfaffian!(A::Union{SkewHermitian{<:Real},SkewHermTridiagonal{<:Real}})= _logabspfaffian!(A)
 logabspfaffian(A::Union{SkewHermitian{<:Real},SkewHermTridiagonal{<:Real}})= logabspfaffian!(copyeigtype(A))
 
 """
     logabspfaffian(A)
 
-Returns a tuple `(log|Pf A|, sign)`, with the log of the absolute value of the pfaffian of `A` as first output
-and the sign (±1) of the pfaffian as second output. A must be a real skew-Hermitian matrix.
-If `A` is not of type `SkewHermitian{<:Real}`, then `isskewhermitian(A)`
-is checked to ensure that `A == -A'`
+Returns a tuple `(log|pf(A)|, sign)`, with the log of the absolute value of the pfaffian of
+`A` as first output and the sign of the pfaffian as second output such that
+`pfaffian(A) ≈ sign * exp(log|pf(A)|)`. `A` must be a skew-symmetric matrix. If `A` is not of type
+`SkewHermitian{<:Real}`, then `isskewsymmetric(A)` is checked to ensure that `A == -transpose(A)`
 """
-logabspfaffian(A::AbstractMatrix{<:Real}) = logabspfaffian!(copyeigtype(A))
+logabspfaffian(A::AbstractMatrix) = logabspfaffian!(copyeigtype(A))
 
 
 """
@@ -211,3 +255,9 @@ function logabspfaffian!(A::AbstractMatrix{<:Real})
     isskewhermitian(A) || throw(ArgumentError("Pfaffian requires a skew-Hermitian matrix"))
     return _logabspfaffian!(SkewHermitian(A))
 end
+
+function logabspfaffian!(A::AbstractMatrix{<:Complex})
+    LinearAlgebra.require_one_based_indexing(A)
+    isskewsymmetric(A) || throw(ArgumentError("Pfaffian requires a skew-symmetric matrix"))
+    return _logabspfaffian!(A)
+end

test/runtests.jl

@@ -414,6 +414,9 @@ end
         if VERSION ≥ v"1.7" # for exact det of BigInt matrices
             @test SkewLinearAlgebra.exactpfaffian(Abig)^2 == det(Abig)
         end
+
+        logpf, sign = logabspfaffian(A)
+        @test pfaffian(A) ≈ sign * exp(logpf)
     end
 
     # issue #49