function testspa
% Test script to test subspacea.m 
% Test results are stored in testspa.txt
% Rico Argentati, 9/7/00
%
% Requires MATLAB 5.3 or newer 
% Previous Revision 1.7 Beta by MEA, 2001/6/19
% Revision 1.9 Beta by AK, 2001/6/23
% Revision 1.10 Beta by MEA, 2001/7/5
% Revision 1.11 Beta by MEA, 2001/7/20
% Revision 1.12 Beta by MEA, 2001/8/28
%
% Need to create the following function separately 
% in order to test passing A as a function
%
% function x=Afunc(v)
% global A
% x=A*v;
% return

global A
testvectors=1;
status_pass='TEST PASS'
status_warning='WARNING'
status_error='ERROR'
numb_pass=0;
numb_warning=0;
numb_error=0;


%Create output file
if exist('testspa.txt')
  delete 'testspa.txt';
end
diary('testspa.txt');

'Tests for subspacea.m'
date
ver


'Test4 - Look for very specific known result for small/medium values'
D=[1e-25 1e-10 1 10 1000];
k=size(D,2);
F = [eye(k) zeros(k)]';
G = [eye(k) diag(D)]';
[theta2,U,V]=subspacea(F,G);
sintheta1=(D.*(1./(sqrt(1+D.*D))))';
sintheta2=sin(theta2);
sintheta1
sintheta2
costheta1=(1./(sqrt(1+D.*D)))';
costheta2=cos(theta2);
costheta1
costheta2
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.062
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test4a try with A as identity'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,A);
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.072
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test4b try with A as an identity and function'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,'Afunc');
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.072
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test5 - Look for very specific known result for medium values'
D=[1e-5 1 10 1000 1e8];
k=size(D,2);
F = [eye(k) zeros(k)]';
G = [eye(k) diag(D)]';
[theta2,U,V]=subspacea(F,G);
sintheta1=(D.*(1./(sqrt(1+D.*D))))';
sintheta2=sin(theta2);
sintheta1
sintheta2
costheta1=(1./(sqrt(1+D.*D)))';
costheta2=cos(theta2);
costheta1
costheta2
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.383
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test5a try with A as identity'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,A);
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.391
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test5b try with A as an identity and function'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,'Afunc');
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.391
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test6 - Look for very specific known result for medium/large values'
'Note this fails if largest value in D is >1e14'
D=[1 10 1000 1e10 1e14];
k=size(D,2);
F = [eye(k) zeros(k)]';
G = [eye(k) diag(D)]';
[theta2,U,V]=subspacea(F,G);
sintheta1=(D.*(1./(sqrt(1+D.*D))))';
sintheta2=sin(theta2);
sintheta1
sintheta2
costheta1=(1./(sqrt(1+D.*D)))';
costheta2=cos(theta2);
costheta1
costheta2
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.762
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test6a try with A as identity'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,A);
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.771
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test6b try with A as an identity and function'
A=eye(size(F,1));
[theta2,U,V]=subspacea(F,G,'Afunc');
sintheta2=sin(theta2); costheta2=cos(theta2);
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.771
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 7: Check singular vectors'
F=rand(50,20);
G=rand(50,10);
[theta1,U1,V1]=subspacea(F,G);
temp=diag(U1'*V1);
diff=norm(cos(theta1)-temp)
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.9
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 8: Check singular vectors with A as identity'
F=rand(50,20);
G=rand(50,10);
A=eye(50);
[theta1,U1,V1]=subspacea(F,G,A);
temp=diag(U1'*V1);
diff=norm(cos(theta1)-temp)
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.9
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 9: Check singular vectors with A scalar product'
F=rand(50,20);
G=rand(50,10);
A=rand(50);
A=A'*A;
[theta1,U1,V1]=subspacea(F,G,A);
temp=diag(U1'*A*V1);
diff=norm(cos(theta1)-temp)
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.9
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 10: Check singular vectors with A as a function'
F=rand(50,20);
G=rand(50,10);
A=rand(50);
A=A'*A;
[theta1,U1,V1]=subspacea(F,G,'Afunc');
temp=diag(U1'*A*V1);
diff=norm(cos(theta1)-temp)
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.9
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 11: Test perturbation'
x=[1 1.0e-2 1e-5 1e-10 1e-15 1e-20 1e-25];
for i=1:size(x,2)
  c1=x(i);
  F=rand(50,20);
  G=rand(50,20);
  F1=F+c1*rand(50,20);
  G1=G+c1*rand(50,20);
  temp=[rank(F) == 20 rank(F1)==20 rank(G)==20 rank(G1)==20];
  if min(temp)==1
    [theta1,U1,V1]=subspacea(F,G);
    [theta2,U1,V1]=subspacea(F1,G1);
    maxdiffcos=max(abs(cos(theta2)-cos(theta1)))
    maxdiffsin=max(abs(sin(theta2)-sin(theta1)))
    gap=sin(max(subspacea(F,F1)))+sin(max(subspacea(G,G1)))
    if gap < (maxdiffcos+maxdiffsin)
      status_warning
      numb_warning=numb_warning+1;
    end
  end
end


'Test12 - test scaling'
'This test fails for smallest value of D<1e-15'
D=[1 1e-15];
k=size(D,2);
F = [eye(k) zeros(k)]';
G = [eye(k) diag(D)]';
[theta2,U,V]=subspacea(F,G);
ans_subspace=subspace(F,G)
sintheta1=flipud((D.*(1./(sqrt(1+D.*D))))');
sintheta2=sin(theta2);
sintheta1
sintheta2
costheta1=flipud((1./(sqrt(1+D.*D)))');
costheta2=cos(theta2);
costheta1
costheta2
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^.916
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

% rescale G
G(2,2)=1e15;
G(4,2)=1;
[theta2,U,V]=subspacea(F,G);
ans_subspace=subspace(F,G)
sintheta1=flipud((D.*(1./(sqrt(1+D.*D))))');
sintheta2=sin(theta2);
sintheta1
sintheta2
costheta1=flipud((1./(sqrt(1+D.*D)))');
costheta2=cos(theta2);
costheta1
costheta2
diff=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
%Adjustment makes tol consistent with experimental results
tol=(norm(F)+norm(G))*eps^1.85
if diff>tol
  status_warning
  numb_warning=numb_warning+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 13 - F,G and A in different combinations'
' and check orthonormality of principal vectors'
F=rand(50,30);
G=rand(50,10);
A=rand(50);
A=A'*A;

'Test 13a - permute F and G, no A'
[theta1,U1,V1]=subspacea(F,G);
[theta2,U2,V2]=subspacea(G,F);
diff_theta=norm(theta2-theta1)
diff_orth=norm(U1'*U1-eye(size(U1,2)));
diff_orth=diff_orth+norm(U2'*U2-eye(size(U2,2)));
diff_orth=diff_orth+norm(V1'*V1-eye(size(V1,2)));
diff_orth=diff_orth+norm(V2'*V2-eye(size(V2,2)))
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.852
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

tol=max(size(F,2),size(G,2))*max(size(F,2),size(G,2))*eps
if diff_orth>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_orth
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 13b - permute F and G, with A as a matrix'
[theta1,U1,V1]=subspacea(F,G,A);
[theta2,U2,V2]=subspacea(G,F,A);
diff_theta=norm(theta2-theta1)
diff_orth=norm(U1'*A*U1-eye(size(U1,2)));
diff_orth=diff_orth+norm(U2'*A*U2-eye(size(U2,2)));
diff_orth=diff_orth+norm(V1'*A*V1-eye(size(V1,2)));
diff_orth=diff_orth+norm(V2'*A*V2-eye(size(V2,2)))
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.828
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end
%Adjustment makes tol consistent with experimental results
tol=size(A,1)*size(A,1)*eps^.936
if diff_orth>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_orth
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 13c - permute F and G, with A as a function'
[theta1,U1,V1]=subspacea(F,G,'Afunc');
[theta2,U2,V2]=subspacea(G,F,'Afunc');
diff_theta=norm(theta2-theta1)
diff_orth=norm(U1'*A*U1-eye(size(U1,2)));
diff_orth=diff_orth+norm(U2'*A*U2-eye(size(U2,2)));
diff_orth=diff_orth+norm(V1'*A*V1-eye(size(V1,2)));
diff_orth=diff_orth+norm(V2'*A*V2-eye(size(V2,2)))
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.889
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

%Adjustment makes tol consistent with experimental results
tol=size(A,1)*size(A,1)*eps^.942
if diff_orth>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_orth
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 14 with Vandermond matrix'
F=gallery('vander',[1:20]);
F=F(:,12:20);
G=eye(20);
G=G(:,1:10);
A=eye(20);

'Test 14a - permute F and G, no A'
[theta1,U1,V1]=subspacea(F,G);
[theta2,U2,V2]=subspacea(G,F);
sintheta1=sin(theta1);
sintheta2=sin(theta2);
costheta1=cos(theta1);
costheta2=cos(theta2);
diff_theta=norm(sintheta2-sintheta1) + norm(costheta2-costheta1)
diff_orth=norm(U1'*U1-eye(size(U1,2)));
diff_orth=diff_orth+norm(U2'*U2-eye(size(U2,2)));
diff_orth=diff_orth+norm(V1'*V1-eye(size(V1,2)));
diff_orth=diff_orth+norm(V2'*V2-eye(size(V2,2)))
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.889
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*max(size(F,2),size(G,2))*eps^.901
if diff_orth>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_orth
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 14b - with A as a matrix'
[theta1A,U1,V1]=subspacea(F,G,A);
[theta2A,U2,V2]=subspacea(G,F,A);
sintheta1A=sin(theta1A);
sintheta2A=sin(theta2A);
costheta1A=cos(theta1A);
costheta2A=cos(theta2A);
diff_theta=norm(sintheta1A-sintheta1) ...
   + norm(costheta1A-costheta1) ...
   + norm(sintheta2A-sintheta2) ...
   + norm(costheta2A-costheta2) 
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.872
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 14c - with A as a function'
[theta1AA,U1,V1]=subspacea(F,G,'Afunc');
[theta2AA,U2,V2]=subspacea(G,F,'Afunc');
sintheta1AA=sin(theta1AA);
sintheta2AA=sin(theta2AA);
costheta1AA=cos(theta1AA);
costheta2AA=cos(theta2AA);
diff_theta=norm(sintheta1AA-sintheta1) ...
   + norm(costheta1AA-costheta1) ...
   + norm(sintheta2AA-sintheta2) ...
   + norm(costheta2AA-costheta2) 
%Adjustment makes tol consistent with experimental results
tol=max(size(F,2),size(G,2))*eps^.872
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 15, from the paper'
D = [ 1.0e10 1.0e8 2 1 0.5 1.0e-10 1.0e-20 ]; %used in the paper
%values in D must be in the decreasing order
exsin = fliplr(D./sqrt(1 + D.*D))';
excos = fliplr(1./sqrt(1 + D.*D))';
D = diag(D);
QF = [eye(size(D)) zeros(size(D))]';
QG = [eye(size(D)) D]';
QF=fliplr(QF);
QG=fliplr(QG);

'Data for Tab. 7.2.'
[thetaQFQG,UQFQG,VQFQG] = subspacea(QF,QG);
sinQFQG = sin(thetaQFQG);
cosQFQG = cos(thetaQFQG);

absersinQFQG = sinQFQG - exsin;
%relersinQFQG = absersinQFQG./exsin
absercosQFQG = cosQFQG - excos;
%relercosQFQG = absercosQFQG./excos

diff_theta=norm(absersinQFQG) ...
   + norm(absercosQFQG)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.918
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

if testvectors
  for i=1:size(D,1),
    diff_vectorsUQFQG(i) = subspacea(UQFQG(:,i),QF(:,i));
    diff_vectorsVQFQG(i) = subspacea(VQFQG(:,i),QG(:,i));
  end
  val=norm(diff_vectorsUQFQG)+norm(diff_vectorsVQFQG)
  %Adjustment makes tol consistent with experimental results
  tol=max(size(QF,2),size(QG,2))*eps^.918
  if val>tol;
    status_warning
    numb_warning=numb_warning+1;
    diff_vectorsUQFQG'
    diff_vectorsVQFQG'
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end

'Data for Tab. 7.3.'
n = size(QF,1);
Q=orth(randn(n,n));

F = Q*QF; G = Q*QG;

[thetaFG,UFG,VFG] = subspacea(F,G);

sinFG = sin(thetaFG);
cosFG = cos(thetaFG);

absersinFG = sinFG - exsin;
%relersinFG = absersinFG./exsin

absercosFG = cosFG - excos;
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFG) ...
   + norm(absercosFG)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.874
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

if testvectors
  for i=1:size(D,1),
    diff_vectorsUFG(i) = subspacea(UFG(:,i),F(:,i));
    diff_vectorsVFG(i) = subspacea(VFG(:,i),G(:,i));
  end
  
  val=norm(diff_vectorsUFG)+norm(diff_vectorsVFG)
  %Adjustment makes tol consistent with experimental results
  tol=max(size(QF,2),size(QG,2))*eps^.304
  if val>tol
    status_warning
    numb_warning=numb_warning+1;
    diff_vectorsUFG'
    diff_vectorsVFG'
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end

'Data for Tab. 7.4.'
p = size(QF,2); q = size(QG,2);
PF = randn(p,p); PG = randn(q,q);
cond(PF)
cond(PG)

FP = F*PF; GP = G*PG;

[thetaFPGP,UFPGP,VFPGP] = subspacea(FP,GP);
sinFPGP = sin(thetaFPGP);
cosFPGP = cos(thetaFPGP);

absersinFPGP = sinFPGP - exsin;
%relersinFG = absersinFG./exsin
absercosFPGP = cosFPGP - excos;
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFPGP) ...
   + norm(absercosFPGP)
 
%Adjustment makes tol consistent with experimental results
tol=max(cond(FP),cond(GP))*eps^.934
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

%Skip this test
if testvectors*0
  for i=1:size(D,1),
    diff_vectorsUFPGP(i) = subspacea(UFPGP(:,i),F(:,i));
    diff_vectorsVFPGP(i) = subspacea(VFPGP(:,i),G(:,i));
  end
  
  val=norm(diff_vectorsUFPGP)+norm(diff_vectorsVFPGP)
  tol=max(size(QF,2),size(QG,2))*eps
  if val>tol 
    status_warning
    numb_warning=numb_warning+1;
    diff_vectorsUFPGP'
    diff_vectorsVFPGP'
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end

'Test 15b with A=I, similar to 14b, but based on 15'
A = eye(size(QF,1));
'Test 15b.1 with A=I, simple QF,QG' 
[thetaQFQGA,UQFQGA,VQFQGA] = subspacea(QF,QG,A);
sinQFQGA = sin(thetaQFQGA);
cosQFQGA = cos(thetaQFQGA);

absersinQFQGA = sinQFQGA - exsin;
%relersinQFQG = absersinQFQG./exsin
absercosQFQGA = cosQFQGA - excos;
%relercosQFQG = absercosQFQG./excos

diff_theta=norm(absersinQFQGA) ...
   + norm(absercosQFQGA)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.915
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

if testvectors
  for i=1:size(D,1),
    diff_vectorsUQFQGA(i) = subspacea(UQFQGA(:,i),QF(:,i));
    diff_vectorsVQFQGA(i) = subspacea(VQFQGA(:,i),QG(:,i));
  end
  
  val=norm(diff_vectorsUQFQGA)+norm(diff_vectorsVQFQGA)
  %Adjustment makes tol consistent with experimental results
  tol=max(size(QF,2),size(QG,2))*eps^.918
  if val>tol
    status_warning
    numb_warning=numb_warning+1;
    diff_vectorsUQFQGA'
    diff_vectorsVQFQGA'
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end

'Test 15b.2 with A=I, general orthogonal F,G' 
[thetaFGA,UFGA,VFGA] = subspacea(F,G,A);
sinFGA = sin(thetaFGA);
cosFGA = cos(thetaFGA);

%absersinFGA = sinFGA - sinFG;
absersinFGA = sinFGA - exsin;
%relersinFG = absersinFG./exsin
%absercosFGA = cosFGA - cosFG;
absercosFGA = cosFGA - excos;
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFGA) ...
   + norm(absercosFGA)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.881
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

if testvectors
  for i=1:size(D,1),
    diff_vectorsUFGA(i) = subspacea(UFGA(:,i),F(:,i));
    diff_vectorsVFGA(i) = subspacea(VFGA(:,i),G(:,i));
  end
  
  val=norm(diff_vectorsUFGA)+norm(diff_vectorsVFGA)
  %Adjustment makes tol consistent with experimental results
  tol=max(size(QF,2),size(QG,2))*eps^.283
  if val>tol
    status_warning
    numb_warning=numb_warning+1;
    diff_vectorsUFGA'
    diff_vectorsVFGA'
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end

'Test 15b.3 with A=I, general FP,GP' 
thetaFPGPA = subspacea(FP,GP,A);
sinFPGPA = sin(thetaFPGPA);
cosFPGPA = cos(thetaFPGPA);

absersinFPGPA = sinFPGPA - sinFPGP;
%relersinFG = absersinFG./exsin
absercosFPGPA = cosFPGPA - cosFPGP;     
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFPGPA) ...
   + norm(absercosFPGPA)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*max(cond(FP),cond(GP))*eps^1.56
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end


'Test 16, general A'
%D = [ 1.0e10 1.0e8 2 1 0.5 1.0e-10 1.0e-14 ]; %used in the paper
D = [ 1.0e1 1.1 1.0001 1 1.1e-10 1.0e-12 1.0e-16 ];
%D = [ 1 2 3 4 5 6 7] %for testing
%testsize = 1000
%D = rand(1,testsize);

%Adiag1 = [ 1 1 1 1 1 1 1]; %should be the same size as D
%Adiag2 = [ 1 1 1 1 1 1 1]; 
Adiag1 = [ 1 2 4 5 2 1 2]; %should be the same size as D
Adiag2 = [ 9 7 2 2 5 1 4]; 
%Adiag1 = rand(1,testsize); Adiag2 = rand(1,testsize);

A = diag([Adiag1 Adiag2]);
    
exsin = sort(D.*sqrt(Adiag2)./sqrt(Adiag1 + Adiag2.*D.*D))'
excos = fliplr(sort( sqrt(Adiag1)./sqrt(Adiag1 + Adiag2.*D.*D)))'

D = diag(D);
QF = [eye(size(D)) zeros(size(D))]';
QG = [eye(size(D)) D]';

'Test 16.1, Diagonal A, simple QF, QG'

thetaQFQG = subspacea(QF,QG,A);
sinQFQG = sin(thetaQFQG);
cosQFQG = cos(thetaQFQG);

absersinQFQG = sinQFQG - exsin
%relersinQFQG = absersinQFQG./exsin
absercosQFQG = cosQFQG - excos
%relercosQFQG = absercosQFQG./excos

diff_theta=norm(absersinQFQG) ...
   + norm(absercosQFQG)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.902
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 16.2, Diagonal A, A-orthogonal QF, QG'
n = size(QF,1);
QA=ortha(A,randn(n,n));

F = QA*sqrt(A)*QF; G = QA*sqrt(A)*QG;

thetaFG = subspacea(F,G,A);
sinFG = sin(thetaFG)
cosFG = cos(thetaFG)

absersinFG = sinFG - exsin
%relersinFG = absersinFG./exsin

absercosFG = cosFG - excos
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFG) ...
   + norm(absercosFG)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.88
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 16.3, Diagonal A, general FP, GP'
p = size(QF,2); q = size(QG,2);
PF = randn(p,p); PG = randn(q,q);
cond(PF)
cond(PG)

FP = QF*PF; GP = QG*PG;

thetaFPGP = subspacea(FP,GP,A);
sinFPGP = sin(thetaFPGP)
cosFPGP = cos(thetaFPGP)

absersinFPGP = sinFPGP - exsin
%relersinFG = absersinFG./exsin
absercosFPGP = cosFPGP - excos
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFPGP) ...
   + norm(absercosFPGP)
 
%Adjustment makes tol consistent with experimental results
tol=max(cond(FP),cond(GP))*eps^.9
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 16.4, general A'
n = size(QF,1);
Q=orth(randn(n,n));
p = size(QF,2); q = size(QG,2);
PF = randn(p,p); PG = randn(q,q);
cond(PF)
cond(PG)

A = Q'*A*Q; 
F = Q'*QF; G = Q'*QG;

thetaFG = subspacea(F,G,A);
sinFG = sin(thetaFG)
cosFG = cos(thetaFG)

absersinFG = sinFG - exsin
%relersinFG = absersinFG./exsin

absercosFG = cosFG - excos
%relercosFG = absercosFG./excos

diff_theta=norm(absersinFG) ...
   + norm(absercosFG)
 
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.872
if diff_theta>tol
  status_warning
  numb_warning=numb_warning+1;
  diff_theta
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 16.5, general A, using sqrt(A)'
for i=1:40;
  condA=10^i;
  D=rand(n,1);
  minD=min(D); maxD=max(D);
  D=1+(D-minD)/(maxD-minD)*(condA-1);
  D=spdiags(D, 0, n, n); 
  QQ=orth(randn(n,n));
  A=QQ'*D*QQ;
  sqA=QQ'*sqrt(D)*QQ;
  err(:,i)=norm(subspacea(F,G,A)-subspacea(sqA*F,sqA*G));
  %Adjustment makes tol consistent with experimental results
  tol=eps^.808;
  if norm(err(i)) > tol 
    status_error
    numb_error=numb_error+1;
    diff
  else
    status_pass
    numb_pass=numb_pass+1;
  end
end
semilogy(err')
title('Error Versus Condition Number','fontsize',16)
print -depsc2 sqrta.eps

'Test 17, zero column in input'
F=rand(20,10);
G=rand(20,7);
theta1=subspacea(F,G);
F=[F zeros(20,1)];
theta2=subspacea(F,G);
G=[G zeros(20,2)];
theta3=subspacea(F,G);
diff=norm(theta1-theta2)+norm(theta1-theta3)

%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.936
if diff>tol
  status_error
  numb_error=numb_error+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

F=rand(20,10);
A=eye(size(F,1));
G=rand(20,7);
theta1=subspacea(F,G,A);
F=[F zeros(20,1)];
theta2=subspacea(F,G,A);
G=[G zeros(20,2)];
theta3=subspacea(F,G,A);
diff=norm(theta1-theta2)+norm(theta1-theta3)

%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.936
if diff>tol
  status_error
  numb_error=numb_error+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

F=rand(20,10);
A=eye(size(F,1));
G=rand(20,7);
theta1=subspacea(F,G,'Afunc');
F=[F zeros(20,1)];
theta2=subspacea(F,G,'Afunc');
G=[G zeros(20,2)];
theta3=subspacea(F,G,'Afunc');
diff=norm(theta1-theta2)+norm(theta1-theta3)
%Adjustment makes tol consistent with experimental results
tol=max(size(QF,2),size(QG,2))*eps^.936
if diff>tol
  status_error
  numb_error=numb_error+1;
  diff
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 18, make sure output theta is a column vector'
F=rand(20,1);
G=rand(20,7);
theta1=subspacea(F,G);
theta2=subspacea(G,F);
if size(theta1,2)>1 | size(theta2,2)>1
  status_error
  numb_error=numb_error+1;
else
  status_pass
  numb_pass=numb_pass+1;
end

A=eye(size(F,1));
theta1=subspacea(F,G,A);
theta2=subspacea(G,F,A);
if size(theta1,2)>1 | size(theta2,2)>1
  status_error
  numb_error=numb_error+1;
else
  status_pass
  numb_pass=numb_pass+1;
end

A=eye(size(F,1));
theta1=subspacea(F,G,'Afunc');
theta2=subspacea(G,F,'Afunc');
if size(theta1,2)>1 | size(theta2,2)>1
  status_error
  numb_error=numb_error+1;
else
  status_pass
  numb_pass=numb_pass+1;
end

'Test 19, test to detect higher precision arithmetic'
%Determine if we are a processor with higher precision arithmetic
precision_test=[1 eps/1024 -1]*[1 1 1]'; 
high_precision=(precision_test>0)
halfsize=10;
F=gallery('vander',[1:2*halfsize]);
F=F(:,halfsize+1:2*halfsize);
G=eye(2*halfsize);
G=G(:,1:halfsize);

A=(10^(-14))*eye(2*halfsize)+hilb(2*halfsize);
[theta,U,V] = subspacea(G,F,A);
error1=norm(V'*(A*V) - eye(size(V,2)))+...
norm(U'*(A*U) - eye(size(U,2)))+...
norm((diag(cos(theta))) - U'*(A*V))
%Adjustment makes tol consistent with experimental results
if (high_precision & error1>eps^.3996) |...
   (~high_precision & error1>eps^.2334)
  status_warning
  numb_warning=numb_warning+1;
  error1
else
  status_pass
  numb_pass=numb_pass+1;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numb_pass
numb_warning
numb_error

diary off

fid = fopen('testspa.txt','a');

fclose(fid);
return