function [ a, ipvt, info ] = zsifa ( a, lda, n )

%*****************************************************************************80
%
%% ZSIFA factors a complex symmetric matrix.
%
%  Discussion:
%
%    The factorization is accomplished by elimination with symmetric pivoting.
%
%    To solve A*X = B, follow ZSIFA by ZSISL.
%
%    To compute inverse(A)*C, follow ZSIFA by ZSISL.
%
%    To compute determinant(A), follow ZSIFA by ZSIDI.
%
%    To compute inverse(A), follow ZSIFA by ZSIDI.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 May 2007
%
%  Author:
%
%    This version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex A(LDA,N); the symmetric matrix to be factored.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Output, complex A(LDA,N); a block diagonal matrix and the multipliers which
%    were used to obtain it.  The factorization can be written A = U*D*U'
%    where U is a product of permutation and unit upper triangular matrices,
%    U' is the transpose of U, and D is block diagonal with 1 by 1 and 2 by 2
%    blocks.  Only the diagonal and upper triangle are used.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, integer INFO.
%    0, normal value.
%    K, if the K-th pivot block is singular.  This is not an error condition
%    for this subroutine, but it does indicate that ZSISL or ZSIDI may
%    divide by zero if called.
%

%
%  Initialize.
%
%  ALPHA is used in choosing pivot block size.
%
  alpha = ( 1.0 + sqrt ( 17.0 ) ) / 8.0;

  info = 0;
%
%  Main loop on K, which goes from N to 1.
%
  k = n;

  while ( true )
%
%  Leave the loop if K = 0 or K = 1.
%
    if ( k == 0 )
      break
    end

    if ( k == 1 )
      ipvt(1) = 1;
      if ( zabs1 ( a(1,1) ) == 0.0 )
        info = 1;
      end
      break
    end
%
%  This section of code determines the kind of
%  elimination to be performed.  When it is completed,
%  KSTEP will be set to the size of the pivot block, and
%  SWAP will be set to TRUE if an interchange is
%  required.
%
    km1 = k - 1;
    absakk = zabs1 ( a(k,k) );
%
%  Determine the largest off-diagonal element in column K.
%
    imax = izamax ( k-1, a(1:k-1,k), 1 );
    colmax = zabs1 ( a(imax,k) );

    if ( alpha * colmax < absakk )

      kstep = 1;
      swap = 0;
%
%  Determine the largest off-diagonal element in row IMAX.
%
    else

      rowmax = 0.0;
      for j = imax + 1 : k
        rowmax = max ( rowmax, zabs1 ( a(imax,j) ) );
      end

      if ( imax ~= 1 )
        jmax = izamax ( imax-1, a(1:imax-1,imax), 1 );
        rowmax = max ( rowmax, zabs1 ( a(jmax,imax) ) );
      end

      if ( alpha * rowmax <= zabs1 ( a(imax,imax) ) )
        kstep = 1;
        swap = 1;
      elseif ( alpha * colmax * ( colmax / rowmax ) <= absakk )
        kstep = 1;
        swap = 0;
      else
        kstep = 2;
        swap = ( imax ~= km1 );
      end

    end
%
%  Column K is zero.  Set INFO and iterate the loop.
%
    if ( max ( absakk, colmax ) == 0.0 )
      ipvt(k) = k;
      info = k;
      k = k - kstep;
      continue
    end

    if ( kstep ~= 2 )
%
%  1 x 1 pivot block.
%
%  Perform an interchange.
%
      if ( swap )

        temp           = a(1:imax,imax);
        a(1:imax,imax) = a(1:imax,k);
        a(1:imax,k)    = temp;

        for jj = imax : k
          j = k + imax - jj;
          t         = a(j,k);
          a(j,k)    = a(imax,j);
          a(imax,j) = t;
        end

      end
%
%  Perform the elimination.
%
      for jj = 1 : km1
        j = k - jj;
        mulk = -a(j,k) / a(k,k);
        t = mulk;
        a(1:j,j) = a(1:j,j) + t * a(1:j,k);
        a(j,k) = mulk;
      end
%
%  Set the pivot array.
%
      if ( swap )
        ipvt(k) = imax;
      else
        ipvt(k) = k;
      end
%
%  2 x 2 pivot block.
%
    else

      if ( swap )

        temp           = a(1:imax,imax);
        a(1:imax,imax) = a(1:imax,k-1);
        a(1:imax,k-1)  = temp;

        for jj = imax : km1
          j = km1 + imax - jj;
          t         = a(j,k-1);
          a(j,k-1)  = a(imax,j);
          a(imax,j) = t;
        end

        t         = a(k-1,k);
        a(k-1,k)  = a(imax,k);
        a(imax,k) = t;

      end
%
%  Perform the elimination.
%
      km2 = k - 2;

      if ( km2 ~= 0 )

        ak = a(k,k) / a(k-1,k);
        akm1 = a(k-1,k-1) / a(k-1,k);
        denom = 1.0 - ak * akm1;

        for jj = 1 : km2
          j = km1 - jj;
          bk = a(j,k) / a(k-1,k);
          bkm1 = a(j,k-1) / a(k-1,k);
          mulk = ( akm1 * bk - bkm1 ) / denom;
          mulkm1 = ( ak * bkm1 - bk ) / denom;
          t = mulk;
          a(1:j,j) = a(1:j,j) + t * a(1:j,k);
          t = mulkm1;
          a(1:j,j) = a(1:j,j) + t * a(1:j,k-1);
          a(j,k) = mulk;
          a(j,k-1) = mulkm1;
        end

      end
%
%  Set the pivot array.
%
      if ( swap )
        ipvt(k) = -imax;
      else
        ipvt(k) = 1 - k;
      end

      ipvt(k-1) = ipvt(k);

    end

    k = k - kstep;

  end

  return
end
