function x = lsnormal(A,b)
% LSNORMAL   Solve linear least squares by normal equations.
% Input: 
%   A     coefficient matrix (m by n, m>n)
%   b     right-hand side (m by 1)
% Output:
%   x     minimizer of || b-Ax ||

B = A'*A;  c = A'*b;
R = chol(B);
z = forwardsub(R',c);                   % solve R'z=c
x = backsub(R,z);                       % solve Rx=z