% This one use infinite norms for all orthogonality checks and also
% calculates uu' and u'u

clc;
clear all;
close all;

% matrix dimensions
n=100:100:1500;
% Sparse matrix generation
r=10:5:20; %in percentage the number of non zero elemements

% initializing
t1=zeros(length(r),length(n));
t2=t1;
t3=t1;
e1=t1;
e2=t1;
e3=t1;
orth1=t1;
orth2=t1;
orth3=t1;
orth4=t1;
orth5=t1;
orth6=t1;
orth1_1=t1;
orth2_1=t1;
orth3_1=t1;
orth4_1=t1;
orth5_1=t1;
orth6_1=t1;
esval1=t1;
esval2=t1;

% loop starts from here
for i=1:length(r)
    for j=1:length(n)
        % number of non zero elements
        m=n(j)*n(j)*r(i)/100;
        % randomizing indices
        x=randi(n(j),1,m);
        y=randi(n(j),1,m);
        % generating elements
        z=randn(1,m);
        % Matrix generation
        A=sparse(x,y,z,n(j),n(j));
        A=full(A);
        
        % Decompositions
        
        % built-in SVD decomposition
        tic
        [u1,d1,v1]=svd(A);
        t1(i,j) = toc;
        [x,y]=size(d1);
        
        % SVDecom
        tic
        [u2,d2,v2]=SVDecom(A);
        t2(i,j) = toc;
        
        % mp_mysvd
        tic
        [u3,d3,v3]=mp_mysvd(A);
        t3(i,j) = toc;
        d31=d3(1:x,1:y);
        
        % errors in decompositions
        e1(i,j)=norm(u1*d1*v1'-A,inf);
        e2(i,j)=norm(u2*d2*v2'-A,inf);
        e3(i,j)=norm(u3*d31*v3'-A,inf);
        
        % errors in singular values
        esval1(i,j)=norm(d2-d1, inf);
        esval2(i,j)=norm(d31-d1, inf);
        
        % orthogonality checks
        orth1(i,j)=norm(u1*u1'-eye(n(j)),inf);
        orth2(i,j)=norm(v1*v1'-eye(n(j)),inf);
        orth3(i,j)=norm(u2*u2'-eye(n(j)),inf);
        orth4(i,j)=norm(v2*v2'-eye(n(j)),inf);
        orth5(i,j)=norm(u3*u3'-eye(n(j)),inf);
        orth6(i,j)=norm(v3*v3'-eye(n(j)),inf);
        orth1_1(i,j)=norm(u1'*u1-eye(n(j)),inf);
        orth2_1(i,j)=norm(v1'*v1-eye(n(j)),inf);
        orth3_1(i,j)=norm(u2'*u2-eye(n(j)),inf);
        orth4_1(i,j)=norm(v2'*v2-eye(n(j)),inf);
        orth5_1(i,j)=norm(u3'*u3-eye(n(j)),inf);
        orth6_1(i,j)=norm(v3'*v3-eye(n(j)),inf);
        
    end
end

save('t1.mat','t1');
save('t2.mat','t2');
save('t3.mat','t3');
save('e1.mat','e1');
save('e2.mat','e2');
save('e3.mat','e3');
save('orth1.mat','orth1');
save('orth2.mat','orth2');
save('orth3.mat','orth3');
save('orth4.mat','orth4');
save('orth5.mat','orth5');
save('orth6.mat','orth6');
save('orth1_1.mat','orth1_1');
save('orth2_1.mat','orth2_1');
save('orth3_1.mat','orth3_1');
save('orth4_1.mat','orth4_1');
save('orth5_1.mat','orth5_1');
save('orth6_1.mat','orth6_1');
save('esval1.mat','esval1');
save('esval2.mat','esval2');