function [H, DH, X0] = dks(n, k, alpha)
%DKS: Discrete Kohn-Sham equation with a single particle Hamiltonian.
%
% 	H = dks(n, k, alpha) produces a matrix-valued function in the form of
%  		H(X) := L + alpha * Diag( L^-1 * diag(X*X') ),
% 	where L = tridiag(-1,2,-1) is a discrete 1D Laplacian of size n,
% 	Diag is a diagonal matrix, diag extracts the diagonal vector of a
% 	matrix, and alpha>0 is a controlling parameter. Here the variable X
% 	is a real matrix of size n-by-k.
%
% 	The default n = 30 and alpha = 0.5.  
%	Sparse matrices are used if n >= 3000.
%
%	[H, DH, X0] = dks(...) also generates the directional derivative
%		DH(X,Y) = lim_{t to 0} ( H(X+tY) - H(X) ) / t,
%	where Y is of size n-by-k. 
%
%	X0 is an approximate eigenvector to be used for iterative methods such as SCF.
%
% -------
%	Usage:
%	[Hfun, DHfun, X0] = dks(30, 5, .5); 
%	Hx0 = Hfun(X0); % size 30*30
%
% ------
% 	Note: The discrete Kohn-Sham equation is in the form of 
%		H(X)*X=X*Lam 
% 	where Lam is a k-by-k matrix containing the k smallest eigenvalues
% 	of H(X).
%
% 	NEPv features: Real symmetric H(X); Block variable X of size n-by-k;
% 	Unitarily invariant H(X); Smallest eigenvalues;
%
% 	References:
%	[1]	C. Yang, W. Gao, and J. C. Meza, On the convergence of the
%		self-consistent field iteration for a class of nonlinear
%		eigenvalue problems, SIAM J. Matrix Anal. Appl., 30 (2009), pp.
%		1773–1788.
%	[2]	Z. Bai, R.-C. Li, and D. Lu. Sharp estimation of convergence 
%		rate for self-consistent field iteration to solve
%		eigenvector-dependent nonlinear eigenvalue problems, SIAM J.
%		Matrix Anal. Appl., 43 (2022), pp.301-327.
% 
%	Ding.Lu@uky.edu. Version date July-20-2022.
%

if nargin < 1, n = 30; end
if nargin < 2, k = min([n, 5]); end
if nargin < 3, alpha = .5; end

rho = @(X) vecnorm(X,2,2).^2; % rho(X) = diag(X*X');

if n < 3000 % use dense matrices 
	L = 2*diag(ones(n,1)) - diag(ones(n-1,1),1) - diag(ones(n-1,1),-1);
	LL = chol(L); 	
	H = @(X) L + alpha * diag( LL\(LL'\rho(X)) );

	% Generate directional derivative 
	if nargout > 1
		DH = @(X,Y) 2 * alpha * diag( LL\(LL'\sum(Y.*X,2)) );
		vv = [1:1:n]'/(n+1)*pi;
		X0 = sin(vv * [1:1:k])*sqrt(2/(n+1));
	end

else % use sparse matrices
	e = ones(n,1);
	L = spdiags([-e 2*e -e], -1:1, n, n);
	LL = chol(L); 	
	H = @(X) L + alpha * spdiags( LL\(LL'\rho(X)), 0, n, n);

	% Generate directional derivative 
	if nargout > 1
		DH = @(X,Y) 2 * alpha * spdiags( LL\(LL'\sum(Y.*X,2)), 0, n, n);
		vv = [1:1:n]'/(n+1)*pi;
		X0 = sin(vv * [1:1:k])*sqrt(2/(n+1));
	end
end

return;
