function [ ap, det ] = zspdi ( ap, n, ipvt, job )

%*****************************************************************************80
%
%% ZSPDI sets the determinant and inverse of a complex symmetric packed matrix.
%
%  Discussion:
%
%    ZSPDI uses the factors from ZSPFA.
%
%    The matrix is stored in packed form.
%
%    A division by zero will occur if the inverse is requested and ZSPCO has
%    set RCOND to 0.0 or ZSPFA has set INFO nonzero.
%
%  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 matrix factors from ZSPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from ZSPFA.
%
%    Input, integer JOB, has the decimal expansion AB where
%    if B ~= 0, the inverse is computed,
%    if A ~= 0, the determinant is computed,
%    For example, JOB = 11 gives both.
%
%    Output, complex AP(N*(N+1)/2); if the inverse was requested, the upper
%    triangle of the inverse of the original matrix, stored in packed
%    form.  The columns of the upper triangle are stored sequentially
%    in a one-dimensional array.
%
%    Output, complex DET(2), the determinant of the original matrix.
%    Determinant = DET(1) * 10.0**DET(2) with 1.0 <= abs ( DET(1) ) < 10.0
%    or DET(1) = 0.0.  Also, DET(2) is strictly real.
%
  det = [];

  noinv = ( mod ( job, 10 ) == 0 );
  nodet = ( floor ( mod ( job, 100 ) / 10 ) == 0 );

  if ( ~nodet )

    det(1) = 1.0;
    det(2) = 0.0;
    t = 0.0;
    ik = 0;

    for k = 1 : n

      kk = ik + k;
      d = ap(kk);
%
%  2 by 2 block
%  Use det (D  T)  =  ( D / T * C - T ) * T
%          (T  C)
%  to avoid underflow/overflow troubles.
%  Take two passes through scaling.  Use T for flag.
%
      if ( ipvt(k) <= 0 )

        if ( zabs1 ( t ) == 0.0 )
          ikp1 = ik + k;
          kkp1 = ikp1 + k;
          t = ap(kkp1);
          d = ( d / t ) * ap(kkp1+1) - t;
        else
          d = t;
          t = 0.0;
        end

      end

      if ( ~nodet )

        det(1) = det(1) * d;

        if ( zabs1 ( det(1) ) ~= 0.0 )

          while ( zabs1 ( det(1) ) < 1.0 )
            det(1) = det(1) * 10.0;
            det(2) = det(2) - 1.0;
          end

          while ( 10.0 <= zabs1 ( det(1) ) )
            det(1) = det(1) / 10.0;
            det(2) = det(2) + 1.0;
          end

        end

      end

      ik = ik + k;

    end

  end
%
%  Compute inverse ( A ).
%
  if ( ~noinv )

    k = 1;
    ik = 0;

    while ( k <= n )

      km1 = k - 1;
      kk = ik + k;
      ikp1 = ik + k;

      if ( 0 <= ipvt(k) )
%
%  1 by 1
%
        ap(kk) = 1.0 / ap(kk);

        if ( 1 <= km1 )

          work(1:km1) = ap(ik+1:ik+km1);
          ij = 0;

          for j = 1 : km1
            jk = ik + j;
            ap(jk) = ap(ij+1:ij+j) * transpose ( work(1:j) );
            ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kk) = ap(kk) + work(1:km1) * transpose ( ap(ik+1:ik+km1) );

        end

        kstep = 1;
%
%  2 by 2
%
      else

        kkp1 = ikp1 + k;
        t = ap(kkp1);
        ak = ap(kk) / t;
        akp1 = ap(kkp1+1) / t;
        akkp1 = ap(kkp1) / t;
        d = t * ( ak * akp1 - 1.0 );
        ap(kk) = akp1 / d;
        ap(kkp1+1) = ak / d;
        ap(kkp1) = -akkp1 / d;

        if ( 1 <= km1 )

          work(1:km1) = ap(ikp1+1:ikp1+km1);
          ij = 0;

          for j = 1 : km1
            jkp1 = ikp1 + j;
            ap(jkp1) = ap(ij+1:ij+j) * work(1:j);
            ap(ikp1+1:ikp1+j-1) = ...
              ap(ikp1+1:ikp1+j-1) + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kkp1+1) = ap(kkp1+1) + work(1:km1) * transpose ( ap(ikp1+1:ikp1+km1) );

          ap(kkp1) = ap(kkp1) + ap(ik+1:ik+km1) * transpose ( ap(ikp1+1:ikp1+km1) );

          work(1:km1) = ap(ik+1:ik+km1);
          ij = 0;

          for j = 1 : km1
            jk = ik + j;
            ap(jk) = ap(ij+1:ij+j) * transpose ( work(1:j) );
            ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kk) = ap(kk) + work(1:km1) * ap(ik+1:ik+km1);

        end

        kstep = 2;

      end
%
%  Swap.
%
      ks = abs ( ipvt(k) );

      if ( ks ~= k )

        iks = ( ks * ( ks - 1 ) ) / 2;

        temp             = ap(iks+1:iks+ks);
        ap(iks+1:iks+ks) = ap(ik+1:ik+ks);
        ap(ik+1:ik+ks)   = temp;

        ksj = ik + ks;

        for jb = ks : k
          j = k + ks - jb;
          jk = ik + j;
          t       = ap(jk);
          ap(jk)  = ap(ksj);
          ap(ksj) = t;
          ksj = ksj - ( j - 1 );
        end

        if ( kstep ~= 1 )
          kskp1 = ikp1 + ks;
          t         = ap(kskp1);
          ap(kskp1) = ap(kkp1);
          ap(kkp1)  = t;
        end

      end

      ik = ik + k;

      if ( kstep == 2 )
        ik = ik + k + 1;
      end

      k = k + kstep;

    end

  end

  return
end
