function b = zsisl ( a, lda, n, ipvt, b )

%*****************************************************************************80
%
%% zsisl() solves a complex symmetric system that was factored by ZSIFA.
%
%  Discussion:
%
%    A division by zero may occur if ZSICO has set RCOND == 0.0
%    or ZSIFA has set INFO ~= 0.
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call zsifa(a,lda,n,ipvt,info)
%
%      if ( info ~= 0 )
%        error
%      end
%
%      do j = 1, p
%        call zsisl(a,lda,n,ipvt,c(1,j))
%      end
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 May 2006
%
%  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 output from ZSICO or ZSIFA.
%
%    Input, integer LDA, the leading dimension of A.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from ZSICO or ZSIFA.
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%

%
%  Loop backward applying the transformations and D inverse to B.
%
  k = n;

  while ( 0 < k )
%
%  1 x 1 pivot block.
%
    if ( 0 <= ipvt(k) )

      if ( k ~= 1 )

        kp = ipvt(k);

        if ( kp ~= k )
          t       = b(k,1);
          b(k,1)  = b(kp,1);
          b(kp,1) = t;
        end

        b(1:k-1,1) = b(1:k-1,1) + b(k,1) * a(1:k-1,k);

      end

      b(k,1) = b(k,1) / a(k,k);
      k = k - 1;
%
%  2 x 2 pivot block.
%
    else

      if ( k ~= 2 )

        kp = abs ( ipvt(k) );

        if ( kp ~= k - 1 )
          t        = b(k-1,1);
          b(k-1,1) = b(kp,1);
          b(kp,1)  = t;
        end

        b(1:k-2,1) = b(1:k-2,1) + b(k,1)   * a(1:k-2,k);
        b(1:k-2,1) = b(1:k-2,1) + b(k-1,1) * a(1:k-2,k-1);

      end

      ak = a(k,k) / a(k-1,k);
      akm1 = a(k-1,k-1) / a(k-1,k);
      bk = b(k,1) / a(k-1,k);
      bkm1 = b(k-1,1) / a(k-1,k);
      denom = ak * akm1 - 1.0;
      b(k,1) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1,1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;

    end

  end
%
%  Loop forward applying the transformations.
%
  k = 1;

  while ( k <= n )

    if ( 0 <= ipvt(k) )
%
%  1 x 1 pivot block.
%
      if ( k ~= 1 )

        b(k,1) = b(k,1) + b(1:k-1,1)' * a(1:k-1,k);
        kp = ipvt(k);

        if ( kp ~= k )
          t       = b(k,1);
          b(k,1)  = b(kp,1);
          b(kp,1) = t;
        end

      end

      k = k + 1;
%
%  2 x 2 pivot block.
%
    else

      if ( k ~= 1 )

        b(k,1)   = b(k,1)   + b(1:k-1,1)' * a(1:k-1,k,1);
        b(k+1,1) = b(k+1,1) + b(1:k-1,1)' * a(1:k-1,k+1);
        kp = abs ( ipvt(k) );

        if ( kp ~= k )
          t       = b(k,1);
          b(k,1)  = b(kp,1);
          b(kp,1) = t;
        end

      end

      k = k + 2;

    end

  end

  return
end
