function P3()
clc
randn('seed',0)
% Example 1
% A small example to check correctness of your method

% 
m = 100; n = 10; A = randn(m,n); B = randn(n,n); B = B'*B;
G = chol(B,'lower');
[U,S,V] = svd((G\A')');
xDense = G'\V(:,1);
rhoMaxDense = S(1,1);
[x,rhoMax] = BigX(sparse(A),sparse(B));
disp('Example 1.')
rhoMaxCheck = [rhoMaxDense rhoMax]
xCheck = [xDense x]

% Example 2

idx = 67; % Identifies a sym pos def matrix from the UF collection
Problem = UFget(67);
B = Problem.A;
% Generate a sparse m-by-n matrix with m = 5*n...
[n,n] = size(B);
A = sprandn(5*n,n,.01);
fprintf('\n\n\nExample 2.\n\n m = %5d    n = %5d\n\n',5*n,n)
tic
[xCVL,rhoMaxCVL] = BigXCVL(A,B);
tCVL = toc;
fprintf('tCVL = %5.3f    rhoMaxCVL = %10.5f\n',tCVL,rhoMaxCVL)
tic
[xYOU,rhoMaxYOU] = BigX(A,B);
tYOU = toc;
fprintf('tYOU = %5.3f    rhoMaxYOU = %10.5f\n',tYOU,full(rhoMaxYOU))
cosine = abs((xCVL/norm(xCVL))'*(xYOU/norm(xYOU)));
fprintf('\nCosine of angle between your x and my x = %6.3f\n',cosine)


 function [x,rhoMax] = BigXCVL(A,B)
% minimizes norm(Ax) subject to the constraint that  x'*B*x = 1
% where B is n-by-n symmetric positive definite and in sparse format and
% A is m-by-n in sparse format.

% First compute PAP' = GG' where P = I(:,p)
[G,v,p] = chol(B,'lower','vector');

% Through a change of variable we have a max singular value problem.
% Since B = (G'P)'(G'P), we set y = G'Px and A1 = AP'*inv(G') observe that
%       1 = x'*B*x = (G'Px)'(G'Px) = y'*y
% and
%     norm(Ax) = norm(A*inv(G'P)y) = norm(A1*y)
%
% It follows that rhoMax is the largest singular value of A1 and 
% x = inv(G'P)*y1 where y1 is a right singular vector of A1 associated
% with rhoMax.
% 
% Notice that A1 is dense but data-sparse. We cannot use svds because
% that function requires the matrix to be specified in sparse format.
% So we go after the largest eigenvalue and corresponding eigenvactor 
% of A1'*A1. That largest eigenvalue  is rhoMax^2 and the eigenvector is
% y1. 

[n,n] = size(B);
options = struct('issym',1,'disp',0);
[z,LambdaMax] = eigs( @(w) MatProd(w,A,p,G), n, [], 1, 'lm', options );
z = z/norm(z); 
% Solve (G'P)x = y1
x = G'\z;
x(p)  = x;
rhoMax = sqrt(LambdaMax);



     function z = MatProd(w,A,p,G)
     % z = A1'*A1*w 
     %   = (A*P'*inv(G'))'*(A*P'*inv(G')*w
     %   = inv(G)*P*A'*A*P'*inv(G')*w
     
     % Six steps:
     z    = G'\w;
     z(p) = z;
     z    = A*z;
     z    = A'*z;
     z    = z(p);
     z    = G\z;
