% Usage: [E,psi,x,Vs] = schrodinger_galerkin_sparse(N, rho, V, m)
%
% Return the m lowest eigenvalues E and eigenfunctions psi of the Schrodinger
% operator for a potential function V(x), on [-1,1] with periodic boundaries.
%
% Uses a linear-element Galerkin scheme with N points whose spacing
% is proportional to the function rho(x).  (See makegrid.m).  This
% function differs from schrodinger_galerkin.m in that it uses
% sparse-matrix iterative solvers in order to require O(Nm) storage
% and O(Nm^2) time, instead of O(N^2) and O(N^3) for dense matrix methods.
%
% E(i) is an array of eigenvalues, and the columns psi(:,i) are the
% eigenfunctions corresponding to E(i) at the points x(:).  Vs(:)
% is an array giving the average potential V(x) at the x(:) elements.

function [E,psi,x,Vs] = schrodinger_galerkin_sparse(N, rho, V, m)
  x = makegrid(N, rho);       % x_1 .. x_N
  xm = [ x(N)-2, x(1:N-1) ];  % x_0 .. x_{N-1}
  xp = [x(2:N), x(1)+2];      % x_2 .. x_{N+1}
  dx = xp - x(1:N);
  dxm = x(1:N) - xm;

  B = spdiags([dx'/6 (dxm+dx)'/3 dxm'/6], -1:1, N, N);
  B(1,N) = dxm(1)/6;
  B(N,1) = dx(N)/6;

  % check: B should be Hermitian
  if max(max(abs(B - B'))) > 1e-8
    error('B is not Hermitian!');
  end

  % Force B to be exactly Hermitian (in case of numerical errors in makegrid):
  B = (B + B')/2;

  K = spdiags([-1./dx' (1./dxm+1./dx)' -1./dxm'], -1:1, N, N);
  K(1,N) = -1/dxm(1);
  K(N,1) = -1/dx(N);

  V0 = zeros(N,1);
  for n = 1:N
    P = quadl(@(y) V(y + xm(n)) .* y.^2, 0, dxm(n), 1e-8) / dxm(n)^2;
    M = quadl(@(y) V(y + x(n)) .* (dx(n) - y).^2, 0, dx(n), 1e-8) / dx(n)^2;
    V0(n) = P + M;
  end
  V1 = zeros(N,1);
  for n = 1:N-1
    V1(n) = quadl(@(y) V(y + x(n)) .* y .* (dx(n)-y), 0, dx(n), 1e-8) /dx(n)^2;
  end
  Vm = spdiags([V1 V0 [0;V1(1:end-1)]], -1:1, N, N);
  Vm(N,1) = quadl(@(y) V(y + x(N)) .* y .* (dx(N)-y), 0, dx(N), 1e-8) /dx(N)^2;
  Vm(1,N) = Vm(N,1);

  % hack: we'll shift eigenvalues by +40 to make them all positive
  % (since we know smallest eig is around -35), and then we can
  % use 'sm' (smallest-magnitude) Arnoldi, which seems better behaved
  % than 'sa' for some reason I don't understand.
  A = K + Vm + 40*B;
  
  % check: A should be Hermitian
  if max(max(abs(A - A'))) > 1e-8
    error('A is not Hermitian!');
  end

  opts.disp = 0;
  opts.maxit = 1000;
  [psis,Es] = eigs(A,B, m, 'sm', opts);
  E = diag(Es) - 40;
  [E,Ei] = sort(E);
  E = E(1:m);
  psi = psis(:,Ei(1:m));

  Vs = V0 ./ diag(B); % potential values on grid