2
$\begingroup$

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!

Improve this question
$\endgroup$

3 Answers 3

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Just to emphasize this answer, the QR routines in lapack do handle nonsquare inputs out of the box. The relevant routines for e.g. double numbers are dgeqrf (for factoring A=QR) and dormqr (for applying Q*C to some other matrix C). You should be able to backsolve by R using dtrsm in the "tall-skinny case", I'd have to consult the documentation myself for what to do about the trapezoidal R that arises in the "short-fat" case. $\endgroup$
    – rchilton1980
    Apr 25, 2014 at 15:01
0
$\begingroup$

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).

Improve this answer
$\endgroup$
0
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.