
% for inA and pdB; attain upper boundary

clc
clear

density = 1e-2;
repeats = 10;
% 0 for easy case, 1 for hard case 1, 2 for hard case 2

% options for eigifp
opts.useprecon = 'NO';
opts.disp = 0;
opts.maxit = 10000;

% options for galahad
control.out = 6;
control.unitm = 0;
control.print_level = -1;
fprintf(' Generating\n')

fid1 = fopen('pdBindAe.txt','w');
fid2 = fopen('pdBindAh1.txt','w');
fid3 = fopen('pdBindAh2.txt','w');

for n = 10000*[1 1.5 2]
  
  for cases = 0:2
    
    iter1 = 0; iter2 = 0; iter3 = 0;
    time1 = 0; time2 = 0; time3 = 0;
    fval1 = 0; fval2 = 0; fval3 = 0;
    feas1 = 0; feas2 = 0; feas3 = 0;  
    
    for repeat = 1:repeats
      
      tst = tic;
      B = sprandsym(n,density,0.1,2); % positive definite with cond(B) approximately 1/0.1
      A = sprandsym(n,density); 
      
      opts.initialvec = sum(A)';
      [lambda0,va] = eigifp(A,B,1,opts);
      fprintf(' Done generation, time %g\n', toc(tst))
      
      if cases == 0  
    
        lambda = lambda0 - 5*rand(1) - 5; % avoid hard case.
        x = randn(n,1)/10; 
        a = (A - lambda*B)*x;
        s = x'*B*x;
        ub = s*1.2;
        lb = s*0.8;  
      
      elseif cases == 1
        
        lambda = lambda0;
        x = randn(n,1)/10;
        a = (A - lambda*B)*x;
        s = x'*B*x;
        ub = s*0.6; % make this small to likely get hard case 1
        lb = ub*0.6;
    
      elseif cases == 2 % one dimensional nullspace
      
        lambda = lambda0;
        x = randn(n,1)/10;
        a = (A - lambda*B)*x;
        s = x'*B*x;
        ub = s*1.2;
        lb = s*1.1;
    
      end
      fprintf(' RW \n')
      [y,iter,fval,feas,dval,k_slope,time,lambda] = gtrsalgo_ul(A,a,B,lb,ub,1e-13);
      feas = y'*B*y - ub;
      iter1 = iter1 + iter;
      time1 = time1 + time;
      fval1 = fval1 + fval;
      feas1 = feas1 + abs(feas);
    
      fprintf(' GLTR \n')
      tstart = tic;
      [x,obj,inform] = galahad_gltr(A,-a,0,sqrt(ub),control,B);
      time = toc(tstart); fval = 2*obj;  
      feas = x'*B*x - ub;
      iter = inform.iter;
      iter2 = iter2 + iter;
      time2 = time2 + time;
      fval2 = fval2 + fval;
      feas2 = feas2 + abs(feas);
    
      fprintf(' Newton \n')
      [u,lambda,fval,feas,iter,time] = Newton_gtrs(A,a,B,ub,1e-12,10);
      iter3 = iter3 + iter;
      time3 = time3 + time;
      fval3 = fval3 + fval;
      feas3 = feas3 + abs(feas);
      
    end
  
    if cases == 0; 
      fprintf(fid1,'%5d & %3.0f/%6.2f/%10.9e/%2.1e &', n, iter1/repeats, time1/repeats, fval1/repeats, feas1/repeats);
      fprintf(fid1,'%3.0f/%6.2f/%10.9e/%2.1e &', iter2/repeats, time2/repeats, fval2/repeats, feas2/repeats);
      fprintf(fid1,'%3.0f/%6.2f/%10.9e/%2.1e\n', iter3/repeats, time3/repeats, fval3/repeats, feas3/repeats);
    elseif cases == 1;
      fprintf(fid2,'%5d & %3.0f/%6.2f/%10.9e/%2.1e &', n, iter1/repeats, time1/repeats, fval1/repeats, feas1/repeats);
      fprintf(fid2,'%3.0f/%6.2f/%10.9e/%2.1e &', iter2/repeats, time2/repeats, fval2/repeats, feas2/repeats);
      fprintf(fid2,'%3.0f/%6.2f/%10.9e/%2.1e\n', iter3/repeats, time3/repeats, fval3/repeats, feas3/repeats);
    elseif cases == 2;
      fprintf(fid3,'%5d & %3.0f/%6.2f/%10.9e/%2.1e &', n, iter1/repeats, time1/repeats, fval1/repeats, feas1/repeats);
      fprintf(fid3,'%3.0f/%6.2f/%10.9e/%2.1e &', iter2/repeats, time2/repeats, fval2/repeats, feas2/repeats);
      fprintf(fid3,'%3.0f/%6.2f/%10.9e/%2.1e\n', iter3/repeats, time3/repeats, fval3/repeats, feas3/repeats);
    end
  
  end

end

fclose(fid1);
fclose(fid2);
fclose(fid3);