% Test scripts that illustrate various aspects of the Arnoldi iteration.
%
% Feb 1, 2022, P.G. Martinsson, UT-Austin

function main_arnoldi

%DRIVER1
%DRIVER2
%DRIVER3
DRIVER4

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script illustrates how things go south due to loss of orthogonality.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER1

rng(0)
n = 10;

%%% Generate A
V  = randn(n) + n*eye(n);
dd = 1 + 0.1*(1 - 2*rand(1,n)); % Vector of eigenvalues.
A  = V*diag(dd)*inv(V);

%%% Arnoldi iteration
b = randn(n,1);
Q = zeros(n);
H = zeros(n);
Q(:,1) = b/norm(b);
for j = 1:(n-1)
  v        = A*Q(:,j);
  H(1:j,j) = Q(:,1:j)'*v;
  v        = v - Q(:,1:j)*H(1:j,j);
  H(j+1,j) = norm(v);
  Q(:,j+1) = (1/H(j+1,j))*v;
end
H(:,n) = Q'*A*Q(:,n);

fprintf(1,'Max evec discrepancy = %10.3e\n',max(abs(sort(eig(A)) - sort(eig(H)))))
fprintf(1,'|A - Q*H*trans(Q)    = %10.3e\n',max(max(abs(A - Q*H*Q'))))
fprintf(1,'|I - Q*trans(Q)      = %10.3e\n',max(max(abs(eye(n) - Q*Q'))))

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here is fix the loss of orthogonality.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER2

rng(0)
n  = 10;

%%% Generate A
V  = randn(n) + n*eye(n);
dd = 1 + 0.1*(1 - 2*rand(1,n));
A  = V*diag(dd)*inv(V);

%%% Arnoldi iteration
b = randn(n,1);
Q = zeros(n);
H = zeros(n);
Q(:,1) = b/norm(b);
for j = 1:(n-1)
  v        = A*Q(:,j);
  H(1:j,j) = Q(:,1:j)'*v;
  v        = v - Q(:,1:j)*H(1:j,j);
  H(j+1,j) = norm(v);
  v        = v - Q(:,1:j)*(Q(:,1:j)'*v); % This is the extra step!!
  Q(:,j+1) = (1/norm(v))*v;
end
H(:,n) = Q'*A*Q(:,n);

fprintf(1,'Max evec discrepancy = %10.3e\n',max(abs(sort(eig(A)) - sort(eig(H)))))
fprintf(1,'|A - Q*H*trans(Q)    = %10.3e\n',max(max(abs(A - Q*H*Q'))))
fprintf(1,'|I - Q*trans(Q)      = %10.3e\n',max(max(abs(eye(n) - Q*Q'))))

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script plots the Ritz values at each step.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER3

rng(0)
n  = 15;

%%% Generate A
V  = randn(n) + n*eye(n);
dd = 1 + 0.1*(1 - 2*rand(1,n));
A  = V*diag(dd)*inv(V);
ee = eig(A);

%%% Arnoldi iteration
b = randn(n,1);
Q = zeros(n);
H = zeros(n);
Q(:,1) = b/norm(b);
for j = 1:(n-1)
  fprintf(1,'j=%d\n',j)
  v = A*Q(:,j);
  H(1:j,j) = Q(:,1:j)'*v;
  v        = v - Q(:,1:j)*H(1:j,j);
  H(j+1,j) = norm(v);
  v        = v - Q(:,1:j)*(Q(:,1:j)'*v);
  Q(:,j+1) = (1/norm(v))*v;
  figure(1)
  subplot(1,1,1)
  ff = eig(H(1:j,1:j));
  plot(real(ee),imag(ee),'b.',...
       real(ff),imag(ff),'ro',...
       'MarkerSize',25)
  axis equal
  keyboard
end
H(:,n) = Q'*A*Q(:,n);

fprintf(1,'Max evec discrepancy = %10.3e\n',max(abs(sort(eig(A)) - sort(eig(H)))))
fprintf(1,'|A - Q*H*trans(Q)    = %10.3e\n',max(max(abs(A - Q*H*Q'))))
fprintf(1,'|I - Q*trans(Q)      = %10.3e\n',max(max(abs(eye(n) - Q*Q'))))

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the Ritz values and produce a convergence graph.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER4

rng(0)
n  = 25;

%%% Generate A
V  = randn(n) + 1i*randn(n) + n*eye(n);
dd = 1 + 0.1*(1 - 2*rand(1,n)) + 0.1*1i*(1 - 2*rand(1,n));
dd(n) = 0.4 + 1i*0.3;
%dd = 1 + 0.2*round(2*rand(1,n)) + 0.1*1i*round(2*rand(1,n));
%[V,~] = qr(randn(n));
%dd = 1 + 0.1*(1 - 2*rand(1,n));
A  = V*diag(dd)*inv(V);
ee = eig(A);

%%% Arnoldi iteration
b = randn(n,1);
Q = zeros(n);
H = zeros(n);
Q(:,1) = b/norm(b);
ERR = zeros(n);
for j = 1:(n-1)
  fprintf(1,'j=%d\n',j)
  v        = A*Q(:,j);
  H(1:j,j) = Q(:,1:j)'*v;
  v        = v - Q(:,1:j)*H(1:j,j);
  H(j+1,j) = norm(v);
  v        = v - Q(:,1:j)*(Q(:,1:j)'*v);
  Q(:,j+1) = (1/norm(v))*v;
  figure(1)
  subplot(1,1,1)
  ff = eig(H(1:j,1:j));
  DIST = abs(ee*ones(1,j) - ones(n,1)*ff.');
  ERR(j,1:j) = sort(min(DIST));
  han = plot(real(ee),imag(ee),'b.',...
             real(ff),imag(ff),'ro',...
             'MarkerSize',15);
  set(han(2),'MarkerSize',7)
  axis equal
  keyboard
end
H(:,n) = Q'*A*Q(:,n);

fprintf(1,'Max evec discrepancy = %10.3e\n',max(abs(sort(eig(A)) - sort(eig(H)))))
fprintf(1,'|A - Q*H*trans(Q)    = %10.3e\n',max(max(abs(A - Q*H*Q'))))
fprintf(1,'|I - Q*trans(Q)      = %10.3e\n',max(max(abs(eye(n) - Q*Q'))))

figure(2)
semilogy(ERR)
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
