矩阵的特殊结构很容易通过基于 Givens 旋转的定制 这导致了一种复杂度的方法。下面是一个 Octave/Matlab 代码，它为大小为的循环双对角矩阵。$QR$$O(n)$$Mx=rhs$$M$$n \ge 2$

% Create a cyclic bidiagonal testmatrix
a=[2 -3 -1 5 3]';
b=[7 2 -4 -1 6]';
rhs=[-1 -3 2 6 2]';
n=5;

% Full matrix solution
M=diag(a)+diag(b(1:n-1),1);  M(n,1)=b(n);
x_full_matrix=M\rhs;
res_nrn_full_matrix=norm(M*x_full_matrix-rhs)

% QR factorization using Givens rotations (see Golub and Van Loan), without 
% forming the full Matrix M. M is transformed to an upper triangular 
% matrix R (which is also not formed as a full matrix)
% a and b are overwritten. d stores fill-in entries of R
c=zeros(n-1,1);
s=zeros(n-1,1);
d=zeros(n-2,1);

z=b(n);
for i=1:n-1,
  % Compute Givens rotation
  if z==0, c(i)=1; s(i)=0;
  else
    if abs(z)>abs(a(i)), tau=-a(i)/z; s(i)=1/sqrt(1+tau*tau); c(i)=s(i)*tau;
    else                 tau=-z/a(i); c(i)=1/sqrt(1+tau*tau); s(i)=c(i)*tau;
    end
  end
  % Apply Givens rotation
  a(i)=c(i)*a(i)-s(i)*z;
  a_n=a(n);
  if i==n-1, break; end
  z   =s(i)*b(i);
  b(i)=c(i)*b(i);
  a(n)=c(i)*a_n;
  d(i)=-s(i)*a_n;
end
a(n)  =s(n-1)*b(n-1)+c(n-1)*a_n;
b(n-1)=c(n-1)*b(n-1)-s(n-1)*a_n;

% Apply Givens rotations to rhs
x=rhs;
for i=1:n-1,
  x_i=x(i);
  x(i)=c(i)*x_i-s(i)*x(n);
  x(n)=s(i)*x_i+c(i)*x(n);
end

% Triangular solve
x(n)=a(n)\x(n);
x(n-1)=a(n-1)\(x(n-1)-b(n-1)*x(n));
for i=n-2:-1:1,
  x(i)=a(i)\(x(i)-d(i)*x(n)-b(i)*x(i+1));
end

% Check the answer
res_nrm_cycl_bidiag_QR=norm(M*x-rhs)
diff=norm(x_full_matrix-x)

显然，如果计算效率很重要，应该将此 Octave/Matlab 代码翻译成 Fortran 或 C 等。即使对于小的，例如，这种方法的这种实现也可能比一般的稀疏分解或密集的 LAPACK/BLAS分解快得多。$n$$n=100$$LU$$LU$

“小/中”有多小？N=10？N=10000？N=1000000？您必须解决多少个这样的方程组（一个？几百个？几百万个？）

如果 N 很小并且您必须解决的这些系统的总数很小，那么没有必要做任何比使用密集 LU 分解ala LAPACK/BLAS 更复杂的事情。

如果 N 真的很大并且你必须解决的这些问题的总数很小，那么使用稀疏的 LU 分解例程来完成这项工作。

只有当 N 真的很大并且你有很多要解决的问题时，尝试为问题开发自定义解决方案才有任何意义。