用Givens旋转进行矩阵QR分解

不多废话，直接贴代码

function [A_, T_, T] = my_qr_givens( A )
%利用givens旋转进行qr分解
%输出
%A_ 每次变换后的A矩阵
%T_ 对应于A_的变换矩阵
    A_ = sym([]);
    T_ = sym([]);
    A = sym(A);
    n = size(A,1);
    T = sym(eye(n));
    sum = 1;
    for j = 1:n-1 %从第一列到第n-1列，全部变换到第j列第j个元素上
        for i = j+1:n %从第j+1行到第n行
            if A(i,j) ~= 0
                r = (A(j,j)*A(j,j)'+ A(i,j)*A(i,j)')^(1/2);
                %r = simplify(sym(real(A(j,j))^2 + imag(A(j,j))^2 + real(A(i,j))^2 + imag(A(i,j))^2)^(1/2));
                T_(:,:,sum) = eye(n);
                T_(j,j,sum) = A(j,j)/r;
                T_(i,i,sum) =  A(j,j)/r;
                T_(j,i,sum) = A(i,j)/r;
                T_(i,j,sum) = -A(i,j)/r;
                A_(:,:,sum) = T_(:,:,sum)*A;
                T = simplify(T_(:,:,sum)*T);
                A = A_(:,:,sum);
                sum = sum + 1;
            end
        end
    end
end