function P3()
% Checks out prodEIGS
randn('seed',0)
rand('seed',0)
clc
n = 1000; m = 4; p = .03; A = cell(m,1); neigs = 5;
for k=1:m
A{k} = randns(n,n,p);
end
opts = struct('disp',0,'tol',.000001);
[V,d] = prodEIGS(A,neigs,opts);
disp(' k Real(d(k)) Imag(d(k)) norm(A*V(:,k) - d(k)*V(:,k))')
disp('-------------------------------------------------------------')
for k=1:5
z = V(:,k);
for j=m:-1:1
z = A{j}*z;
end
Re = real(d(k));
Im = imag(d(k));
resid = norm(z - d(k)*V(:,k));
fprintf('%2d %10.5f %10.5f %10.3e\n',k,Re,Im,resid)
end

function A = randns(m,n,p)
% Generates an m-by-n random sparse matrix A
% The probability that A(i,j) is nonzero is p, 0<=p<=1.
% If A(i,j) is nonzero, its value is chosen from the N(0,1) distribution.
A = sparse(m,n);
for i=1:m
for j=1:n
if rand(1,1)<p
A(i,j) = randn(1,1);
end
end
end