function [x,lambda] = fastfiedler(G,k);
% FASTFIEDLER : Fiedler vector of a graph.
%
% x = fastfiedler(G)  returns the eigenvector corresponding to 
%                 the second-smallest eigenvalue of the laplacian matrix 
%                 of the graph whose adjacency matrix is G|G'.
%                 (If G is connected this is the smallest positive eval.)
% [x,lambda] = fastfiedler(G)  also returns the eigenvalue.
% [x,lambda] = fastfiedler(G,k) returns the k smallest nonzero eigenvalues
%                 and their eigenvectors.
%
% If G is full or smaller than 100 vertices we use Matlab's full "eig";
% if G is sparse we use an iterative method.
%
% See also LAPLACIAN, SPECPARTITION, FIEDLER.
%
% John Gilbert  2 February 1994
% Copyright (c) 1993-1994 by Xerox Corporation.  All rights reserved.

if nargin < 2
    k = 1;
end;

[n,n] = size(G);
if n < 100, 
    G = full(G); 
end;
L = laplacian(G);

%L(1,1) = L(1,1) - 2;
%L(3,3) = L(3,3) +2;

%for i = 1:length(L),
% for j = 1: i-1,
%   L(i,j) = - L(i,j);
% end;
% end;

%L(8,8) = L(8,8) - 10
%L(1,1) = L(1,1) -2
%L(2,2) = L(2,2) -2

%for i = 1: length(L),
%   L(i,i) = 0;
%end


if issparse(G)
    % Sparse graph -- use iteration.
    x = ones(n,1)/sqrt(n);
    lambda = 1;
    shift = -eps^.25;
    I = speye(n);
    p = symmmd(L);
    L = L(p,p);
    R = chol(L-shift*I);
    Rt = R';
    for j = 1:k
        vector = j;
        % Append the j+1'st eigenvector to x.
    
        % Use a few power iterations with a small negative shift first.
        y = rand(n,1);
        y = y - x*(x'*y);
        for i = 1:20*j
             preiteration = i;
             y = R\(Rt\y);
             estimate = 1/norm(y) + shift;
             y = y - x*(x'*y);
             y = y/norm(y);
        end;

        % Now use Rayleigh quotient iteration with a variable shift.
        tol = eps^.25;
        maxiter = 10;
        iteration = 0;
        done = 0;
        while ~done,
            iteration = iteration+1;
            % Get shift as Rayleigh quotient
            estimate = y'*L*y;
            % Get new y by inverse iteration with Rayleigh shift
            newy = (L - estimate*I) \ y;
            newy = newy/norm(newy);
            % Test for convergence
            measure = norm(newy-y)*norm(newy+y);
            done = measure < tol;
            y = newy;
            % Project y to be orthogonal to earlier eigenvectors.
            y = y - x*(x'*y);
            y = y/norm(y);
            if iteration == maxiter
                done = 1;
%                disp('Rayleigh iteration reached maxiter.');
            end
        end;
        x = [x y];
        lambda = [lambda estimate];
    end;
    x = x(:,2:k+1);
    x(p,:) = x;
    lambda = lambda(2:k+1);

else

    % Full graph -- just use eig.
    [V,D] = eig(L);
    [d,i] = sort(diag(D));
    lambda = d(2:k+1);
    x = V(:,i(2:k+1));

end;