% 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');