function b = zspsl ( ap, n, ipvt, b )

%*****************************************************************************80
%
%% zspsl() solves a complex symmetric system factored by ZSPFA.
%
%  Discussion:
%
%    A division by zero may occur if ZSPCO has set RCOND == 0.0
%    or ZSPFA has set INFO ~= 0.
%
%    To compute inverse(A) * C where C is a matrix with P columns:
%
%      call zspfa(ap,n,ipvt,info)
%
%      if (info ~= 0)
%        error
%      end
%
%      do j = 1, p
%        call zspsl(ap,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 AP(N*(N+1)/2), the output from ZSPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from ZSPFA.
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%
  ap = ap(:);
%
%  Loop backward applying the transformations and d inverse to b.
%
  k = n;
  ik = ( n * ( n - 1 ) ) / 2;

  while ( 0 < k )

    kk = ik + k;
    if ( 0 <= ipvt(k) )
%
%  1 x 1 pivot block.
%
      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) * ap(ik+1:ik+k-1,1);

      end
%
%  Apply D inverse.
%
      b(k,1) = b(k,1) / ap(kk,1);
      k = k - 1;
      ik = ik - k;
%
%  2 x 2 pivot block.
%
    else

      ikm1 = ik - ( k - 1 );

      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)   * ap(ik+1:ik+k-2,1);
        b(1:k-2,1) = b(1:k-2,1) + b(k-1,1) * ap(ikm1+1:ikm1+k-2,1);

      end
%
%  Apply D inverse.
%
      km1k = ik + k - 1;
      kk = ik + k;
      ak = ap(kk,1) / ap(km1k,1);
      km1km1 = ikm1 + k - 1;
      akm1 = ap(km1km1,1) / ap(km1k,1);
      bk = b(k,1) / ap(km1k,1);
      bkm1 = b(k-1,1) / ap(km1k,1);
      denom = ak * akm1 - 1.0;
      b(k,1) = ( akm1 * bk - bkm1 ) / denom;
      b(k-1,1) = ( ak * bkm1 - bk ) / denom;
      k = k - 2;
      ik = ik - ( k + 1 ) - k;

    end

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

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

      if ( k ~= 1 )
        b(k,1) = b(k,1) + ap(ik+1:ik+k-1,1)' * b(1:k-1,1);
        kp = ipvt(k);
        if ( kp ~= k )
          t       = b(k,1);
          b(k,1)  = b(kp,1);
          b(kp,1) = t;
        end
      end

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

      if ( k ~= 1 )

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

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

      end

      ik = ik + k + k + 1;
      k = k + 2;

    end

  end

  return
end
