function [K,U,V] = assembleFermiMatrix(p,N,Ufun,Vfun,bc)
%ASSEMBLEFERMIMATRIX Assembly of discretized many-particle Hamiltonian.
%   [K,U,V] = ASSEMBLEFERMIMATRIX(P,N,UFUN,VFUN,BC) assembles the stiffness
%   and potential matrices for a spectral discretization of the Schroedinger
%   equation on [-1,1] for N interaction fermions with P basis functions per
%   particle. UFUN and VFUN are function handles describing the confinement
%   and interaction potential, respectively. BC is a string specifying the
%   boundary conditions; supported are 'Dirichlet' (homogeneous) and 
%   'periodic'. Returned are the stiffness matrix K, confinement potential
%   matrix U and interaction potential matrix V.
%
%   For example, this computes the first eigenvalues (by increasing 
%   magnitude) of the stiffness matrix for non-interacting electrons:
%      K = assembleFermiMatrix(15,4,@(x) 0*x,@(x) 0*x(:,1),'dirichlet');
%      lambda = eigs(K,3,'SM')
%   (the exact eigenvalues are 30*pi^2/4, 39*pi^2/4, 46*pi^2/4).
%
%   This computes the two lowest energy levels of the helium atom:
%     [K,U,V] = assembleFermiMatrix(8,3,@(x) 0*x,@(x) cos(x(:,1)-x(:,2)),'periodic');
%     lambda = eigs(K+U+V,2,'SM')
%   (no exact values are known, but the result should be close to 0.9642
%   and 3.9642).

%   February 23, 2012
%   Christian Clason (christian.clason@uni-graz.at)
%   Greg von Winckel (gregory.von-winckel@uni-graz.at)

%% Precompute
[xi,K1d] = computeLGL(p,bc); % 1D Legendre-Gauss-Lobatto quadrature points
diagK1d  = diag(K1d);        % diagonal elements of 1D stiffness matrix

d  = nchoosek(1:N,2);      % all pairs of interacting particles
T  = nchoosek(1:p,N);      % matrix of nD basis elements tuples
Np = size(T,1);            % number of nD basis elements
Nk = Np*N/2*(p-N);         % number of nonzero off-diagonal entries in K
Ta = [T ones(Np,1)*(p+1)]; % precompute: augment T to account for G_N
G  = diff(Ta,1,2)-1;       % precompute: length of gaps in alpha
vN = 1:Np;                 % precompute: vector (1,...,Np)
rInd = triu(ones(N,N-1))+ones(N,1)*(1:(N-1)); % indices after removal of m
hash = (p-N).^(N-1:-1:0)'; % assign unique id to each basis element
dict = T * hash;           % dictionary for fast lookup of basis elements
lookup = sparse(dict,ones(Np,1),vN,dict(Np),1); % lookup via sparse index

%% Compute potential terms
Vc = arrayfun(@(j)Vfun(xi(T(:,d(j,:)))),1:size(d,1),'UniformOutput',false);
V  = sparse(vN,vN,sum(cat(2,Vc{:}),2),Np,Np);  % interaction potential
U  = sparse(vN,vN,sum(Ufun(xi(T)),2),Np,Np);   % confinement potential

%% Compute stiffness matrix
row = zeros(Nk,1);               % preallocate vector of row indices of K
col = zeros(Nk,1);               % preallocate vector of column indices
val = zeros(Nk,1);               % preallocate vector of entries
indstart = 0;                    % running index for above vectors

for j = 1:(Np-1)                 % loop over all tuples: j = tau(alpha)
    for n = vN(G(j,:)>0)         % insert at all indices n in G_alpha
        bLength = 1:G(j,n);      % number of possible values for b (vector)
        b = T(j,n) + bLength;    % values to be inserted at beta_n
        hash1 = hash(rInd(n,:)); % common part of hash for beta_i, i != n
        hash2 = hash(n)*b;       % beta_n = b
        sgn = (-1).^(n+1:-1:2);  % sign of K_jk: (-1)^(m+n)
        for m = 1:n              % remove at all indices m <= n
            ind = indstart + bLength;           % index of entries to set
            dex = T(j,rInd(m,:))*hash1 + hash2; % id of beta
            row(ind) = j;                       % j = tau(alpha)
            col(ind) = lookup(dex);             % k = tau(beta)
            val(ind) = sgn(m)*K1d(T(j,m),b);    % K_jk
            indstart = ind(end);                % increment running index
        end
    end
end

Ko = sparse(row,col,val,Np,Np);                 % off-diagonal entries of K
Kd = sparse(vN,vN,sum(diagK1d(T),2),Np,Np);     % diagonal entries of K
K  = Kd + Ko + Ko';                             % stiffness matrix
% end function fermion_setup

function [xi,K] = computeLGL(p,bc)
switch lower(bc)
    case 'periodic'
        xi = 2*pi*(0:p-1)'/p;
        K  = fft(ifft(diag([0:floor(p/2) floor((p-1)/2:-1:1)].^2)).');
    case 'dirichlet'
        p2 = p+2;
        b  = sqrt([(1:p).^2./(4*(1:p).^2-1) (p+1)/(2*p+1)])';
        [V, X] = eig(diag(b,1) + diag(b,-1));
        [x, ord] = sort(diag(X));
        V  = V(:,ord)';
        w  = 2*V(:,1).^2;
        xi = x(2:end-1);
        X  = x*ones(1,p2);
        Xdiff = X-X'+eye(p2);
        W = (1./prod(Xdiff,2))*ones(1,p2);
        D = W./(W'.*Xdiff);
        D(1:(p2+1):(p2^2)) = 1-sum(D);
        D  = -D'*diag(1./sqrt(w));
        Di = D(:,2:end-1);
        K  = Di'*diag(w)*Di;
    otherwise
        error(['Boundary condition "' bc '" is not implemented.'])
end
% end function computeLGL
