hhlqr.m

来自「％Jacobi迭代法求方程组等程序」· M 代码 · 共 33 行

function [Q,R]=hhlqr(A)
%用Householder变换求解长方阵A的正交三角分解，即QR分解
%其中Q为正交矩阵,R为与A同型的上三角矩阵，使得A=Q*R
[m,n]=size(A);
if  m<n
    error('the rows of the input matrix must be not smaller than its columns!')
end
for j=1:n                                           %compute the inner product of every column of matrix A
    s(j)=sum(A(:,j).*A(:,j));
end
H=eye(m);
for k=1:n
    sigma(k)=sign(A(k,k))*sqrt(s(k));
    beta(k)=sigma(k)*(sigma(k)-A(k,k));
    A(k,k)=A(k,k)-sigma(k);   
    for j=k+1:n
        gamma=sum(A(k:m,k).*A(k:m,j));
        mu=gamma/beta(k);
        A(k:m,j)=A(k:m,j)-mu*A(k:m,k);
        s(j)=s(j)-A(k,j)^2;        
    end
    for j=1:m
        gamma=sum(A(k:m,k).*H(k:m,j));
        mu=gamma/beta(k);
        H(k:m,j)=H(k:m,j)-mu*A(k:m,k);    
    end
    A(k,k)=sigma(k);
    A(k+1:m,k)=0;
    k
    A
end
Q=H';
R=A;