Trefethen 和 Bau 的书Numerical Linear Algebra在第 10 章中有 Householder QR 算法，它是考虑一般矩形矩阵编写的。它也在Golub 和 van Loan的矩阵计算中。这两本书或多或少都以 Matlab 代码的形式给出了算法，因此您必须在两者之间进行一些翻译。

您可以在LAPACK中找到 QR 的 Fortran 实现，以及许多其他算法。

MINPACK 使用带旋转的非方形 QR 分解，即它不是真正的 QR 分解（如 LAPACK 中的分解），而是更广泛的更稳健的 QRP 分解（其中 P 是表示旋转的置换矩阵）。

我有这些子程序，其中一个对 nxm 矩阵 A 进行 qr 分解，并将户主矩阵的上三角 R 和 w 存储在 A 中，同时将 R 的对角线保存在 d 向量中。然后我尝试形成 Q ^ Tb 这样我就可以用反向替换解决 Rx=Q^tb。我用它用最小二乘法进行数据拟合（所以使用反向替换我应该采用系数），但它给了我错误的结果。如果您对我做错了什么有任何想法，那将有很大帮助，因为我看不出问题出在哪里！谢谢您的帮助！

子程序 qrdecomposition(a,c,d) 隐式 无

实数，维度（n，m） :: a 实数，维度（m） :: c 实数，维度（m） :: d 逻辑 :: sing !integer :: i,j,k,n,m 实数 :: scale1 ,sigma,sum1,tau sing=.false。do k=1,m scale1=0 do i=k,n scale1=max(scale1,abs(a(i,k))) end do if(scale1.eq.0.)then sing=.true。c(k)=0。d(k)=0。否则做 i=k,n a(i,k)=a(i,k)/scale1 end 做 sum1=0。do i=k,n sum1=sum1+a(i,k)*a(i,k) end do sigma=sign(sqrt(sum1),a(k,k)) a(k,k)=a( k,k)+sigma c(k)=sigma*a(k,k) d(k)=-scale1*sigma do j=k+1,m sum1=0。do i=k,n sum1=sum1+a(i,k)*a(i,j) end do tau=sum1/c(k) do i=k,n a(i,j)=a(i, j)-tau*a(i,k) end do end do end if end do d(m)=a(m,m) if(d(m).eq.0.)sing=.true。返回结束子程序 qrdecomposition

 subroutine qrsolution(a,c,d,b)
  implicit none
  !integer :: n,m
 real, dimension(n,m) :: a 
  real, dimension(n) :: b
  real, dimension(m) :: c
  real, dimension(m) :: d

  !integer :: i,j
  real :: sum1,tau
  do  j=1,m-1
    sum1=0.
   do  i=j,n
      sum1=sum1+a(i,j)*b(i)

end do tau=sum1/c(j) do i=j,m b(i)=b(i)-tau*a(i,j) end do end do

做 i=m,1,-1 sum1=0。do j=i+1,m sum1=sum1+a(i,j)*b(j) end do b(i)=(b(i)-sum1)/d(i) end do return

结束子程序