fortran code-algorithm for qr decomposition of non-square matrix

Question

I am trying to implement QR factorization of a non-square matrix in FORTRAN. I have the algorithm for a square matrix but not for a non-square. I use Housholder matrices. Do you know where I could find the whole algorithm or code for the non-square case so that I can be sure I am doing it the right way (i am new in programming!)? I have to solve (backsubstitute) for the upper trangular R and the orthogonal Q at the last part but the decomposition must be done for the non-square case. Thank you!

Answer 1

Trefethen and Bau's book Numerical Linear Algebra has the Householder QR algorithm in chapter 10, and it's written considering general rectangular matrices. It's also in Matrix Computations by Golub and van Loan. Both books give algorithms more or less in the form of Matlab code, so you'll have to do some translation between the two.

You can find a Fortran implementation of QR in LAPACK, as well as numerous other algorithms.

Answer 2

MINPACK uses a non-square QR factorization with pivoting, i.e. it's not a real QR decomposition (like the one in LAPACK), but a generalized more robust QRP decomposition (where P is the permutation matrix representing the pivoting).

Answer 3

I have these subroutines the one does the qr factorization for the nxm matrix A and it stores the upper triangular R and the w of the householder matrices in A while it saves the diagonals of R in d vector.then i try to form the Q^Tb so that i can solve the Rx=Q^tb with backsubstitution.I use it to do a data fitting with least square method(so with the backsubstitution i am supposed to take the coefficients) but it gives me wrong results. If you have any ideas of what i am doing wrong it would be of great help because I cannot see where the problem is!Thank you for your help!

subroutine qrdecomposition(a,c,d) implicit none

real, dimension(n,m) :: a real, dimension(m) :: c real, dimension(m) :: d logical :: sing !integer :: i,j,k,n,m real :: 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. else do i=k,n a(i,k)=a(i,k)/scale1 end do 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. return end subroutine 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

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

end subroutine qrsolution