function [Q,d] = DataSparseEig(A,B,U,k)
% A is nxn, symmetric, and in sparse format.
% B is rxr, symmetric positive definite and in full format.
% U is nxr in sparse format.
% k is a positive integer
% d is a column k-vector that consists of the k most positive eigenvalues
% of M = A + U*inv(B)*U'.
% The columns of Q (nxr) are orthonormal eigenvectors associated with
% d(1),...,d(k). In particular, M*Q(:,j) = d(j)*Q(:,j), j=1:k

[n,r] = size(U);
G = chol(B,'lower');
OPTS.issym = 1;
[Q,D] = eigs(@(x) MProd(x,A,G,U),n,k,'LA',OPTS);
d = diag(D);

function y = MProd(x,A,G,U)

y = A*x + U*(G'\(G\(U'*x)));

% Grading:
%
% 2 for solving it somehow, using eigs
% 1 for passing a function to eigs
% 2 for that function being as efficient as possible