function [H, DH, x0] = dgp2d(N, beta, mode)
%DGP2D: Discrete 2D Gross-Pitaevskii equation.
%
% 	H = dgp2d(N, beta) produces a matrix-valued function in the form of
%  		H(X) := Af + beta * Diag(|X|)^2,
% 	where Af is a positive definite Hermitian matrix of size n = N^2,
% 	Diag is a diagonal matrix, and beta>0 is a controlling parameter.
% 	Here the variable X is a length-n complex vector and |X| represents
% 	absolute values of the vector.
%
%	The matrix Af is dependent of a potential function f(x,y). 
%	If mode == 1, then a symmetric potential f = (x.^2+y.^2)/2 is used;
%	If mode == 2, then a non-symmetric potential f = (x.^2+100*y.^2)/2
%	is used;
%
% 	The default N = 10, beta = 2, and mode == 1. 
%	Sparse matrices are used if N >= 100, i.e., matrix size n = N^2 >= 10000.
%
%	[H, DH, x0] = dgp2d(...) 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] = dgp2d(10, 2.0, 1); 
%	Hx0 = Hfun(X0); % size 100*100
%
% -------
% 	Note: The discrete Gross-Pitaevskii equation is in the form of 
%		H(X)*X=X*lam, 
% 	where lam is the smallest eigenvalues of H(X).
%
% 	NEPv features: Complex Hermitian H(X); Single vector X; Unitarily
% 	invariant; Smallest eigenvalues;
%
% 	References:
%	[1]	E. Jarlebring, S. Kvaal, and W. Michiels, An inverse iteration
%		method for eigenvalue problems with eigenvector nonlinearities, 
%		SIAM J. Sci. Comput., 36 (2014), pp. A1978– A2001.
%	[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 < 3, mode = 1; end;
if nargin < 2, beta = 2; end;
if nargin < 1, N = 10; end;

% set potential function
if mode == 1
   Vf = @(x,y) (x.^2+y.^2)/2;      	% symmetric potential 
elseif mode == 2
   Vf = @(x,y) (x.^2+100*y.^2)/2;	% non-symmetric potential 
end
L = 1;          % domain [-L, L] * [-L, L]
h = 2*L/(N+1);  % stepsize
omega = 0.85;

DN = (diag(ones(N-1,1),1) - diag(ones(N-1,1),-1))/2/h; % central difference for 1st deriv
D2N = (diag(ones(N-1,1),1) -diag(ones(N,1))*2 +diag(ones(N-1,1),-1))/h/h; % for 2nd deriv
LN = kron(D2N,speye(N)) + kron(speye(N), D2N);
xx = linspace(-L, L, N+2); xx = xx(2:end-1)';
LFN = kron(spdiags(xx, 0, N, N), DN) - kron(DN, spdiags(xx, 0, N, N));
[XX, YY] = meshgrid(xx,xx);
Vtilde = arrayfun(Vf, XX,YY);
Vtilde = Vtilde(:); % vectorize
A0tilde_real = -0.5*LN + spdiags(Vtilde, 0, N^2, N^2);
A0tilde_imag = -1*omega*LFN;
A0 = h^2*(A0tilde_real + 1i*A0tilde_imag);

if N < 100 	% use dense matrices if N < 100 (i.e., n < 10000)
	A0 = full(A0);
	Bv = @(X) diag(abs(X).^2);
	H = @(X)  A0 + Bv(X)*beta;

	% Generate directional derivative 
	if nargout > 1
		DH = @(X,Y) 2 * beta * diag( real(conj(X).*Y) );
		[V0, D0] = eig(A0);
		[~, idx] = sort(diag(D0), 'ascend');
		x0 = V0(:,idx(1));
	end


else % use dense matrices
	Bv = @(X) spdiags(abs(X).^2, 0, N^2, N^2);
	H = @(X)  A0 + Bv(X)*beta;

	% Generate directional derivative 
	if nargout > 1
		DH = @(X,Y) 2 * beta * spdiags( real(conj(X).*Y), 0, N^2, N^2);
		[V0, D0] = eigs(A0, 20, 'sr');
		[~, idx] = sort(diag(D0), 'ascend');
		x0 = V0(:,idx(1));
	end
end

return;
