function P1()
% Checks Out FieldOfValues
% Superimposes the boundary of my W(A) on your plot of W(A).

randn('seed',0)
close all

% Example 1 (Small random example.)
n = 3; i = sqrt(-1); A = randn(n,n) + i*randn(n,n);
figure
FieldOfValues(A)
axis equal
hold on
FieldOfValuesCVL(A)
title('Example 1')

% Example 2 (Larger random example.)
n = 20; i = sqrt(-1); A = randn(n,n) + i*randn(n,n);
figure
% We benchmark the "big" example just to guide the grading. The idea is not
% to "beat" solution code but to help identify serious inefficiencies.
tic
FieldOfValues(A)
t1 = toc;
axis equal
hold on
tic;
FieldOfValuesCVL(A)
s1 = toc;
%title('Example 2')
title(sprintf('Example 2. YourTime/MyTime = %5.3f',t1/s1))

% Example 3 (A normal matrix example.)
n = 20; i = sqrt(-1);
[Q,R] = qr(randn(n,n) + i*randn(n,n));
A = 3*Q - 7*Q^2 + 2*Q^3;
figure
FieldOfValues(A)
axis equal
hold on
FieldOfValuesCVL(A)
title('Example 3')

% Example 4 (A stochastic example.)
n = 10;
A = rand(n,n);
for k=1:n
A(:,k) = A(:,k)/sum(A(:,k));
end
figure
FieldOfValues(A)
axis equal
hold on
FieldOfValuesCVL(A)
title('Example 4')

% Example 5 (An Hermitian example.)
n = 5; A = randn(n,n) + i*randn(n,n); A = A + A';
figure
FieldOfValues(A)
axis equal
hold on
FieldOfValuesCVL(A)
title('Example 5')

% Example 6 (A positive definite example.)
n = 5; B = randn(n,n) + i*randn(n,n); B = B'*B;
C = randn(n,n)+i*randn(n,n); C = C-C'; A = B+C;
figure
FieldOfValues(A)
axis equal
hold on
FieldOfValuesCVL(A)
title('Example 6')

function FieldOfValuesCVL(A)
% Displays the field of values of an n-by-n complex matrix A
[n,n] = size(A);
% The number of boundary points (even)
m = 50;
% The termination tolerance for Jacobi...
tol = .001*norm(A,1);

% Note that the largest and smallest eigenvalue of B(theta) is the smallest
% and largest eigenvalue of B(theta+pi). This observation halves the amount
% of work
delTheta = 2*pi/m;
theta = 0;
BoundaryPt = zeros(m,1);
w = zeros(m,1);
Q = eye(n,n);
for k=1:m/2
% Generate the Hermitian matrix...
Btheta = A*exp(-i*theta); Btheta = (Btheta+Btheta')/2;
% Compute its extremal eigenvectors..
%[z1,zn] = JacHerm1(Btheta);
[z1,zn,Q] = HermJac(Btheta,Q,tol);
% and the boundary points...
w(k) = z1'*A*z1;
w(k+m/2) = zn'*A*zn;
theta = theta + delTheta;
end
% There is no need to worry about the counterclockwise order of the w(k).
% But it is a harmless overhead to sort w so that the boundary points have
% that property.
% fill(real(z),imag(z),'y')
plot(real([w;w(1)]),imag([w;w(1)]),'m','linewidth',2)
hold on
e = eig(A);
plot(real(e),imag(e),'ok')
axis equal
shg

function [z1,zn,Q1] = HermJac(B,Q,tol)
% B is an n-by-n Hermitian matrix.
% Q1 is unitary such that off(Q1'*B*Q1)<=tol
% Q is unitary and shuold be regarded as an approximation to Q
% z1 is a column unit vector and approximates an eigenvector of B
% associated with B's largest eigenvalue.
% zn is a column unit vector and approximates an eigenvector of B
% associated with B's smallest eigenvalue.
[n,n] = size(B);
tol1 = norm(B,'fro')*tol;
B = Q'*B*Q;
sweep = 0;
while off(B)>tol1 && sweep<=10
for p=1:n-1
for q=p+1:n
% Solve the (p,q) subproblem...
U = Two(B([p q],[p q]));
Q(:,[p q]) = Q(:,[p q])*U;
B(:,[p q]) = B(:,[p q])*U;
B([p q],:) = U'*B([p q],:);
end
end
sweep = sweep + 1;
end
[lambda1,j] = max(real(diag(B)));
z1 = Q(:,j);
[lambdan,j] = min(real(diag(B)));
zn = Q(:,j);
Q1 = Q;

function U = Two(B)
% B is 2-by-2 Hermitian
% U is 2-by-2 unitary so that U'8B*U is diagonal.

r = abs(B(2,1)); z = B(2,1)/r;
% Note that if V = [1 0 ; 0 z] then V is unitary and
% V'*B*V = [ B(1,1) r; r B(2,2)] = Btilde is real symmetric.
% Using GVL page 428 we compute c and s so that
% [c s;-s c]'*Btilde*[c s; -s c] is diagonal...
if r ~=0
tau = (B(2,2)-B(1,1))/(2*r);
if tau>=0
t = 1/(tau + sqrt(1+tau^2));
else
t = -1/(-tau + sqrt(1 + tau^2));
end
c = 1/sqrt(1 + t^2); s = t*c;
else
c = 1; s = 0;
end
U = [1 0 ; 0 z]*[c s; -s c];
%[U,D] = schur(B);

function tau = off(A)
tau = sqrt(norm(A,'fro')^2 - sum(diag(A).^2));