function xLS = ls_rotation(A,y)
%  xLS = ls_rotation(A,y)
%
%  Use rotation matrices to compute 
%    xLS = min_x ||A*x - y||_2

%  Error checking.

  [my,ny] = size(y);
  if ny ~= 1
    fprintf(' *** Error! y must be a column vector.\n');
    xLS = [];  %  Return empty vector.
    return
  end
  [mA,nA] = size(A);
  if my ~= mA
    fprintf(' *** Error! Row dimension of A and y must be the same.\n');
    xLS = [];
    return
  end
  if nA>mA
  fprintf(' *** Error! A cannot have more columns than rows.\n');
    xLS = [];
    return
  end
   
%  Reduce A to upper triangular form using rotation matrices.

  for j = 1:nA

    %  Reduce column j of A below the diagonal to zeros.

    for i = j+1:mA

      %  Build vector v = [A(j,j)]
      %                   [A(i,j)]

      row_indices = [j,i];
      v = A(row_indices,j);

      %  Build rotation matrix Q which reduces v to form (r,0).

      [c,s,Q] = rotation_mat(v);

      %  Apply rotation matrix to rows j and i of A and to entries 
      %  j and i of y.

      A(row_indices,:) = Q'*A(row_indices,:);
      y(row_indices) = Q'*y(row_indices);

    end  %  end loop on i
  end  %  end loop on j

%  Solve R*xLS = y1, where R is the upper triangle of A and y1 consists 
%  of the first nA entries of y. Note that the entries of y have been 
%  modified in the reduction process.

  R = triu(A);
  R = R(1:nA,1:nA);  %  If A isn't square, need to extract nA x nA block.
  y1 = y(1:nA);  %  Extract 1st nA entries from y.
  xLS = R\y1;    %  Solve R*x = y1 to get xLS.