function [ pf ] = pfaffian_ltl ( A )
%*****************************************************************************80
%
%% pfaffian_ltl computes the Pfaffian of a skew-symmetric matrix
%
%  Discussion:
%
%    pf = pfaffian_LTL(A) computes the Pfaffian of the skew-symmetric matrix A
%    using the Parlett-Reid algorithm
%
%    The determinant of A is the square of the Pfaffian.
%
%  Modified:
%
%    12 May 2020
%
%  Author:
%
%    Michael Wimmer
%
%  Reference:
%
%    Michael Wimmer,
%    Algorithm 923: Efficient Numerical Computation of the Pfaffian for
%    Dense and Banded Skew-Symmetric Matrices,
%    ACM Transactions on Mathematical Software,
%    Volume 38, Number 4, Article 30, August 2012.
%
%  Input:
%
%    real A(N,N): a skew symmetric matrix.
%
%  Output:
%
%    real PF: the value of the Pfaffian.
%
  assert ( ndims ( A ) == 2, 'Argument is not two-dimensional.' )
  assert ( size ( A, 1 ) == size ( A, 2 ), 'Argument is not square.' )
%
%  Make sure input is skew-symmetric.
%
  assert ( norm ( A + A.') < 1e-14 * size ( A, 1 ), 'Argument does not seem skew-symmetric.' )

  N = size ( A, 1 );
    
  if ( mod ( N, 2 ) == 1 )
    pf = 0.0;
    return
  end
    
  pf = 1.0;

  for k = 1 : 2 : N - 2
%
%  Find the largest entry in A(k+1:N,k) and permute it to A(k+1,k).
%
    [ c, kp ] = max ( A(k+1:N,k) );
    kp = kp + k;
%
%  Check if we need to pivot.
%
    if( kp ~= k + 1 )
%
%  interchange rows k+1 and kp.
%
      temp = A(k+1,k:N);
      A(k+1,k:N) = A(kp,k:N);
      A(kp,k:N) = temp;
%
%  Then interchange columns k+1 and kp.
%
      temp = A(k:N,k+1);
      A(k:N,k+1) = A(k:N,kp);
      A(k:N,kp) = temp;
%
%  Every interchange gives a "-" for the determinant of the permutation.
%
      pf = - pf;
    end
        
    pf = pf * A(k,k+1);
%        
%  Form the Gauss vector.
%
    if ( A(k+1,k) ~= 0.0 )
      tau = A(k+2:N,k) / A(k+1,k);
%
%  Update the matrix block A(k+2:,k+2).
%
      A(k+2:N,k+2:N) = A(k+2:N,k+2:N) + tau * A(k+2:N,k+1).' - A(k+2:N,k+1) * tau.';
    end

  end

  pf = pf * A(N-1,N);

  return
end

