% Givens旋转变换的QL/R分解脚本（直接运行版）
% 功能：无需函数调用，通过Givens旋转变换实现QL分解（X=1）或QR分解（X=2），输出变换矩阵Q、分解后矩阵A及迭代次数k
% 原理：根据X选择上三角（QL）或下三角（QR）分解，通过Givens旋转逐步消去非三角元素，直到收敛
% 运行方式：直接运行此m脚本，无需命令行操作，结果自动输出
clear; clc;  % 清空工作区变量与命令行窗口，避免历史数据干扰


%% 1. 输入计算所需参数（可根据需求替换）
% 待分解的方阵A（示例为3阶矩阵，可替换为任意方阵，确保迭代收敛）
A = [5 2 1; 
     2 6 3; 
     1 3 7];  

% 分解模式选择：X=1为QL分解（基于上三角）；X=2为QR分解（基于下三角）
X = 1;  % 可修改为2切换至QR分解

% 输出输入信息，便于后续对比
fprintf('待分解的矩阵A：\n');
disp(A);
if X == 1
    fprintf('\n分解模式：QL分解（基于Givens旋转变换，目标化为下三角矩阵）\n');
elseif X == 2
    fprintf('\n分解模式：QR分解（基于Givens旋转变换，目标化为上三角矩阵）\n');
else
    error('分解模式X必须为1（QL）或2（QR）！');  % 检查输入合法性
end
[m,n] = size(A);
if m ~= n
    error('输入矩阵A必须为方阵！');  % QL/R分解仅适用于方阵
end


%% 2. Givens旋转变换QL/R分解核心逻辑（整合原R_QLR函数核心代码）
C2 = 1;          % 收敛判断指标（初始设为1，满足循环条件）
Q = eye(n);      % 初始化变换矩阵Q为n阶单位矩阵（累积Givens旋转）
k = 0;           % 迭代次数计数器

% 根据X选择分解模式（QL或QR），设置对应的矩阵提取函数和参数
switch X
    case 1  % QL分解：基于上三角部分，目标将A化为下三角矩阵
        fun = @triu;       % 提取上三角矩阵（含对角线）
        option = 1;        % triu的第二个参数，1表示提取主对角线及以上
        opbreak = 1;       % 停止条件的列索引（当最大元素在第1列时停止）
    case 2  % QR分解：基于下三角部分，目标将A化为上三角矩阵
        fun = @tril;       % 提取下三角矩阵（含对角线）
        option = -1;       % tril的第二个参数，-1表示提取主对角线及以下
        opbreak = n;       % 停止条件的列索引（当最大元素在第n列时停止）
end

% 迭代循环：当最大非三角元素的绝对值小于1e-10时收敛
while abs(C2) > 1e-10
    % 步骤1：提取当前模式下的三角部分（上三角或下三角）
    A1 = fun(A, option);
    % 步骤2：寻找三角部分中绝对值最大的元素位置（i,j）
    [C1, I1] = max(abs(A1));  % C1为每列最大绝对值，I1为对应行索引
    [C2, j] = max(C1);        % C2为全局最大绝对值，j为对应列索引
    i = I1(j);                % i为全局最大元素的行索引
    % 步骤3：判断是否达到停止条件（最大元素在目标列时停止）
    if j == opbreak
        break;
    end
    % 步骤4：构造Givens旋转矩阵R（用于消去A(i,j)元素）
    d = A(i,j) / A(j,j);               % 计算旋转参数d
    c = 1 / sqrt(1 + d^2);             % 旋转矩阵的余弦分量
    s = -d * c;                        % 旋转矩阵的正弦分量
    R = eye(n);                        % 初始化n阶单位矩阵
    R(i,i) = c; R(i,j) = s;            % 填充旋转矩阵的(i,i)和(i,j)位置
    R(j,i) = -s; R(j,j) = c;           % 填充旋转矩阵的(j,i)和(j,j)位置
    % 步骤5：更新变换矩阵Q和待分解矩阵A
    Q = Q * R';    % Q累积旋转（R'为R的转置，因Givens矩阵正交，R'=R逆）
    A = R * A;     % 左乘R，逐步将A化为目标三角矩阵
    k = k + 1;     % 迭代次数加1
end


%% 3. 输出计算结果
fprintf('\n==================== Givens旋转变换分解结果 ====================\n');
fprintf('1. 迭代次数 k = %d\n', k);  % 输出总迭代次数

fprintf('\n2. 正交变换矩阵 Q（满足Q''=Q逆）：\n');
disp(Q);

if X == 1
    fprintf('\n3. QL分解后的下三角矩阵 A（满足A = QL，L为下三角）：\n');
else
    fprintf('\n3. QR分解后的上三角矩阵 A（满足A = QR，R为上三角）：\n');
end
disp(A);


%% 4. 结果验证（验证分解正确性）
fprintf('\n==================== 结果验证 ====================\n');
% 验证1：Q是否为正交矩阵（Q'*Q应接近单位矩阵）
ortho_check = Q' * Q;
fprintf('Q的正交性验证（Q''*Q 应接近单位矩阵）：\n');
disp(ortho_check);

% 验证2：分解正确性（QL或QR应接近原始矩阵）
if X == 1  % QL分解：原始A ≈ Q * 分解后的A（因A=QL，L为分解后的下三角）
    recover_A = Q * A;
    fprintf('\nQL分解验证（原始A ≈ Q * 分解后A）：\n');
else  % QR分解：原始A ≈ Q * 分解后的A（因A=QR，R为分解后的上三角）
    recover_A = Q * A;
    fprintf('\nQR分解验证（原始A ≈ Q * 分解后A）：\n');
end
disp(recover_A);

% 量化误差（Frobenius范数）
error = norm(recover_A - A, 'fro');
fprintf('\n分解误差（||恢复矩阵 - 原始矩阵||_F）：%.6e\n', error);
fprintf('说明：误差为浮点计算误差，数值≤1e-10时表明分解正确。\n');