% sparse matrix factorization
clear all;

rand('seed',12);
randn('seed',12);

fid_1 = fopen('sparse_matrix_factorization.txt','w');

repeats = 20; % repeats > 1

tol = 1e-9;
maxiter = 10000;
sigma = 1e-4;

%% begin loop for n and r
for n = 100:100:400
    for r = 0.6:0.1:0.9
        s = n*r;
        
        progress_1 = []; progress_2 = [];
        fopen(fid_1);
       
        for repeat = 1: repeats
            % preprocessing
            tstart_pre = tic;
            A_0 = randn(n);
            A_0 = sparsematrix(A_0,s);
            A = A_0*A_0';
                       
            [U,S] = eig(A);
            V = U*diag(sqrt(diag(S)))*U';
        
            
            Q_0 = eye(n);
            
            time_eig = 0;
            if n > 2000
                clear opts
                opts.issym = 1;
                tstart = tic;
                L = eigs(V*V',1,'LM',opts);
                time_eig = time_eig + toc(tstart);
            else
                tstart = tic;
                L = norm(V*V');
                time_eig = time_eig + toc(tstart);
            end
            
            fprintf('\n Lipschitz constant L = %g, time_eig = %g\n', L, time_eig)
           
            L_0 = L + sigma;
            
            time0 = toc(tstart_pre);
            fprintf('******** Done preprocessing ******\n')
            
            fprintf('******** NPG_BB for SMF ********\n')
            
            clear opts
            
            opts.maxiter = maxiter;
            opts.freq = 10;
            opts.disp = 0;
            opts.M = 5;
            opts.tol = tol;
            tstart = tic;
            [ Q_1, fnew_1, iter_1, ~] = NPG_BB(V, Q_0, s,opts);
            time_1 = toc(tstart);
           
            progress_1 = [progress_1; sqrt(fnew_1*2) iter_1 time_1];
            
            fprintf('******** NPG_constant for SMF ********\n')
            
            clear opts
            
            opts.maxiter = maxiter;
            opts.freq = 1;
            opts.disp = 0;
            opts.tol = tol;
            tstart = tic;
            [ Q_2, fnew_2, iter_2, ~] = NPG_constant(V, Q_0, L_0, s,opts);
            time_2 = toc(tstart);
            progress_2 = [progress_2; sqrt(fnew_2*2) iter_2 time_2];
            
        end
        fprintf(fid_1, '%5d & %5.1f &%1.0e & %5.0f & %5.1f &%1.0e & %5.0f & %5.1f\n', n, r, mean(progress_1), mean(progress_2));
    end
  
end
fclose(fid_1);

