www.pudn.com > NumericalComputingwithMatlabCode.zip > qrsteps.m, change:2003-10-28,size:1458b

```function [A,b] = qrsteps(A,b)
%QRSTEPS  Orthogonal-triangular decomposition.
%   Demonstrates M-file version of built-in QR function.
%   R = QRSTEPS(A) is the upper trapezoidal matrix R that
%   results from the orthogonal transformation, R = Q'*A.
%   With no output argument, QRSTEPS(A) shows the steps in
%   the computation of R.  Press <enter> key after each step.
%   [R,bout]  = QRSTEPS(A,b) also applies the transformation to b.

[m,n] = size(A);
if nargin < 2, b = zeros(m,0); end

% Householder reduction to triangular form.

if nargout == 0
clc
A
if ~isempty(b), b, end
pause
end

for k = 1:min(m-1,n)
% Introduce zeros below the diagonal in the k-th column.
% Use Householder transformation, I - rho*u*u'.
i = k:m;
u = A(i,k);
sigma = norm(u);

% Skip transformation if column is already zero.
if sigma ~= 0

if u(1) ~= 0, sigma = sign(u(1))*sigma; end
u(1) = u(1) + sigma;
rho = 1/(conj(sigma)*u(1));

% Update the k-th column.
A(i,k) = 0;
A(k,k) = -sigma;

% Apply the transformation to remaining columns of A.
j = k+1:n;
v = rho*(u'*A(i,j));
A(i,j) = A(i,j) - u*v;

% Apply the transformation to b.
if ~isempty(b)
tau = rho*(u'*b(i));
b(i) = b(i) - tau*u;
end
end
if nargout == 0
clc
A
if ~isempty(b), b, end
pause
end
end
```