
% for pdA and pdB; it looks at instances that attain lower 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('pdBpdAetwo.txt','w');
fid2 = fopen('pdBpdAh1two.txt','w');
fid3 = fopen('pdBpdAh2two.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;
    
    time0 = 0;
    
    repeat = 1;
    while repeat <= repeats
      
      tst = tic;
      B = sprandsym(n,density,0.1,2); % positive definite with cond(B) approximately 1/0.1
      A = sprandsym(n,density,0.1,2) + 10*B; 
      
      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; % used to divide by 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 % meanwhile 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,side] = gtrsalgo_ul(A,a,B,lb,ub,1e-13);
      feas = y'*B*y - lb;
      iter1 = iter1 + iter;
      time1 = time1 + time;
      fval1 = fval1 + fval;
      feas1 = feas1 + abs(feas);
    
      if abs(side - lb)/abs(lb) <= 100*eps
        fprintf(' Newton \n')
        [u,~,fval,feas,iter,time] = Newton_gtrs(A,a,B,lb,1e-12,10);
        iter3 = iter3 + iter;
        time3 = time3 + time;
        fval3 = fval3 + fval;
        feas3 = feas3 + abs(feas);
        repeat = repeat + 1;
      else
        fprintf('\n\n REGENERATE !!!\n\n')
        iter1 = iter1 - iter;
        time1 = time1 - time;
        fval1 = fval1 - fval;
        feas1 = feas1 - abs(feas);
      end
      
      if abs(side - lb)/abs(lb) <= 100*eps
        opts.useprecon = 'NO';
        opts.disp = 0;
        opts.initialvec = sum(A)';
        opts.tol = 1e-8;
        fprintf(' GLTR \n')
        
        tstart = tic;
          tstart2 = tic;
          lambda0 = eigifp(A,B,1,opts);
          timepre = toc(tstart2);
        lambda1 = lambda0 + 1;
        Atmp = A - lambda1*B; % makes the matrix indefinite
        [x,obj,inform] = galahad_gltr(Atmp,-a,0,sqrt(lb),control,B);
        time = toc(tstart); fval = 2*obj + lambda1*lb;  iter = inform.iter;
        feas = x'*B*x - lb;
        
        iter2 = iter2 + iter;
        time2 = time2 + time;  time0 = time0 + timepre;
        fval2 = fval2 + fval;
        feas2 = feas2 + abs(feas);
      end
      
    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 (%6.2f)/%10.9e/%2.1e &', iter2/repeats, time2/repeats, time0/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 (%6.2f)/%10.9e/%2.1e &', iter2/repeats, time2/repeats, time0/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 (%6.2f)/%10.9e/%2.1e &', iter2/repeats, time2/repeats, time0/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);