% Define matrices A and B
%A = sprandn(4, 3,0.5);
A = complex(randn(4,3),randn(4,3));
%B = sprandn(4, 1,0.5);
B = complex(randn(4,1),randn(4,1));
display(A);
display(B);

[m, n] = size(A);
x=lsqr(A,B); 
disp(x);

% Set initial guess for the solution vector and tolerance
n = size(A, 2); % Define the size of the solution vector prior
tol = 1e-4;
max_iters = 50;
s = zeros(n, 1); % Initialize the solution vector
rk = B; % Initialize the residual vector

% Perform the two-loop iteration

 for k = 1: max_iters
   
   i=1;
   while i <= n
       % Solve for the update vector sk using the least-squares method
        ai = A(:, i);
        delta_sk = (ai'*rk) / (norm(ai)^2);
        %disp(delta_sk);
        s(i) = s(i)+ delta_sk; % Update the update vector
       
        rk = rk-(ai*delta_sk);% Update the residual vector rk
        % disp(rk);
        c = norm(rk);
       
        i= i+1;
   end
  % Check for convergence
        if c < tol
            break;
        end
  plot(c,k);
 end 
 
disp("Solution vector:");
disp(s);

disp("Norm of the residual:");
disp(c);