function [Ta,T] = prim(A,s)
% function [Ta,T] = prim(A,s)
%
% generate Prim's minimum spanning tree in the Graph starting at s
% Ta is the tree in adj matrix format.
%
% T is the parent pointer vector

n = length(A);

if (nargin < 2),
   s = 1;
end

W = max(max(A));


%initialize the min-heap

H = W*ones(n+1,2) + 1;
H(1,:) = -1;
H(2,1) = 0;
H(:,2) = [-1,1:n]';
H(2,2) = s;
H(s+1,2) = 1;


hindex = [2:n+1];
hindex(1) = s+1;
hindex(s) = 2;

ppi = -1*ones(n,1);

while (length(H(:,1)) > 1)

   [u,H,hindex] = extractMin(H,hindex);
   nbrs = find(A(:,u));
   for k = 1:length(nbrs),
      v = nbrs(k);
      if (hindex(v) > 0) & (A(u,v) < H(hindex(v),1)),
           ppi(v) = u;
          [H,hindex] = heapUpdate(H,hindex,v,A(u,v));

      end
   end
end

T = ppi;

Ta = sparse(n,n);

for i = 1:n,
  if (T(i) ~= -1)
    Ta(i,T(i)) = 1;
  end
end

Ta = Ta + Ta';



function [H,hindex] = heapUpdate(H,hindex,v,w)
   p = hindex(v);
   H(p,1) = w;
   H(p,2) = v;
   while (p > 1) & (H(p,1) < H(ceil(p/2),1))
     hindex(H(p,2)) = ceil(p/2);
     hindex(H(ceil(p/2),2)) = p;	
     temp = H(p,:);
     H(p,:) = H(ceil(p/2),:);
     H(ceil(p/2),:) = temp;
     p = ceil(p/2);
   end



function [u,H,hindex] = extractMin(H,hindex)
   m = length(H(:,1));
   if (m > 1),
     u = H(2,2);
     hindex(u) = -1;
     H(2,:) = H(m,:);
     hindex(H(m,2))=2;
     H = H(1:m-1,:);
     m = m -1;
     p = 2;
     done = 0;
     while (2*p - 1 < m) & (done == 0)
        c = 2*p -1;
        if c < m
          if H(c+1,1) < H(c,1), c = c+1; end
        end
        if H(p,1) > H(c,1),
             hindex(H(p,2)) = c;
             hindex(H(c,2)) = p;	
             temp = H(p,:);
             H(p,:) = H(c,:);
             H(c,:) = temp;
             p = c;
        else done = 1;
        end
      end;
   end



