%%MGS.m ======== Implement the modified Gram-Schmidt method for A=QR
%%Usage [1]:  [Q R]=mgs(A)      for a n x n matrix A=Q*R
%%      [2]:  [A R]=mgs(A, k)   for step k of MGS method (Sequential!)
%%      [3]:     A =mgs(A, k)   for step k of MGS method (Sequential!)
%%e.g. _____________________________________________________________________
%  n1=5;n2=3; rand('seed',0); A = rand(n1,n2);  A_orig = A;  
%        [A R] = mgs(A,1),  [A R] = mgs(A,2),  [Q R] = mgs(A,3),
%        A=A_orig,      R = Q'*A %----------Short Q Version --------qr(A,0)
%  extra = n2+1:n1;Af=[];  Af(:,1:n2)=A;  %_____________ <Method a> ________
%                    Af(extra,extra)=eye(n1-n2); 
%        [Q1 R1] = mgs(Af)       % Do QR for Af=[A I_e] effectively
%  disp(['Only the top block in R1 to match R i.e. 1:' num2str(n2)])
%  extra = n2+1:n1;Af=[];  Af(:,1:n2)=Q;  %_____________ <Method b> ________
%                    Af(extra,extra)=eye(n1-n2); 
%        [Q2 R2] = mgs(Af)        % Do QR for Af=[Q I_e] effectively
%  disp(['Only the top block in R2 to match R i.e. 1:' num2str(n2)])
%  extra = n2+1:n1;Ae=zeros(n1,n1-n2);    %_____________ <Method c> ________
%      Ae=-Q*Q(extra,:)';
%                    Ae(extra,:)=eye(n1-n2) + Ae(extra,:); 
%        [Q3 R3] = mgs(Ae,1:n1-n2)        % Do QR for Af=[A I_e] effectively
%         Q_full = [Q Q3] %---------------- Get FULL Q_{n1 x n1} -------qr(A)
function [Q, R] = Modified_GS(A,step)
if nargin<1, help mgs; return, end
   [m n]= size(A);
if nargin<2, step=1:n; end  % Complete steps
if step > n, disp('Error in mgs.m:  Cannot have so many steps!'),
           Size_Input_Matrix = size(A), step, return,end
 R=zeros(m,n); 
for j=step;  U=zeros(1,n); % U a ROW of R ===========================
fprintf('mgs.m ------------ step  j = %d \n', j)

r_j = norm(A(:,j));
if r_j < 1.0E-10 & j<n
   disp('Matrix A (input) may be singular!'), break;
else
   q_found = A(:,j) / r_j;   A(:,j) = q_found;
end 

for k=j:n % by row (modified version)
    U(k) = A(:,k)'*q_found;
end
    U(j)=r_j;  %% Row position j re-computed
for k=j+1:n
    A(:,k) = A(:,k) - q_found*U(k); % Modification of remaining a_j's
end

%if j+1>n
%fprintf('2-norm of a_%d = r_%d = %f (See changes to the last col a_%d) \n  [Q^T * A = R  or A = Q*R] \n',...
%         j,j,r_j,  n)
%else
%fprintf('2-norm of a_%d = r_%d = %f (See changes in a_%d to a_%d)\n',...
%         j,j,r_j,  j+1,n)
%end
    R(j,:)=U;  
end %j =========================== =========================== ==========
    Q = A;
