c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole Publ., 1990 c c Section 5.3 c c Example of QR-factorization using Householder transformations c c c file: ex1s53.f c parameter (n=3,ia=4,m=4) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(i,1),i=1,m) /63.,42.,0.,126./ data (A(i,2),i=1,m) /41.,60.,-28.,82./ data (A(i,3),i=1,m) /-88.,51.,56.,-71./ c print * print *,' QR-factorization using Householder transformations' print *,' Section 5.3, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c call QRfac(n,ia,m,A,V,U,Q,W) c print *,' Matrix Q' call prtmtx(m,ia,m,Q) call mult(n,ia,m,Q,A) print *,' Check A = QR' call prtmtx(n,ia,m,A) c stop end c subroutine QRfac(n,ia,m,A,V,U,Q,W) c c QR-factorization c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c call setoI(m,ia,m,Q) do 2 k=1,n mm = m - k + 1 nn = n - k + 1 call setoI(m,ia,m,W) call UtimesA(nn,ia,mm,A(k,k),V(k),U(k,k),W(k,k)) print *,' Matrix A' call prtmtx(n,ia,m,A) print *,' Matrix W' call prtmtx(m,ia,m,W) call mult(m,ia,m,W,Q) print *,' Matrix Q' call prtmtx(m,ia,m,Q) 2 continue call trans(m,ia,m,Q) c return end c subroutine setoI(n,ia,m,Q) c c set Q to I c dimension Q(ia,n) c do 3 i=1,m do 2 j=1,n Q(i,j) = 0.0 2 continue Q(i,i) = 1.0 3 continue c return end c subroutine prtmtx(n,ia,m,A) c c print array A c dimension A(ia,n) c do 2 i=1,m print 3,(A(i,j),j=1,n) 2 continue print * c return 3 format(2x,4(e13.6,2x)) end c subroutine UtimesA(n,ia,m,A,V,U,W) c c compute UA c dimension A(ia,n),U(ia,m),V(m),W(ia,m) c call findV(n,ia,m,A,V) print *,' V vector' call prtmtx(1,ia,m,V) call formU(m,ia,V,U) print *,' Matrix U' call prtmtx(m,ia,m,U) call mult(n,ia,m,U,A) call formW(m,ia,U,W) c return end c subroutine findV(n,ia,m,A,V) c c determine vector V c dimension A(ia,n),V(m) c beta = - sqrt(prod(m,A(1,1),A(1,1))) V(1) = A(1,1) - beta dem = V(1)*V(1) + prod(m-1,A(2,1),A(2,1)) alpha = sqrt(2.0/dem) V(1) = alpha*(A(1,1) - beta) do 2 i=2,m V(i) = alpha*A(i,1) 2 continue c return end c function prod(n,x,y) dimension x(n),y(n) c sum = 0.0 do 2 i=1,n sum = sum + x(i)*y(i) 2 continue prod = sum c return end c subroutine formU(m,ia,V,U) c c compute the unitary factor U c dimension U(ia,m),V(m) c do 3 i=1,m do 2 j=1,m U(i,j) = - V(i)*V(j) 2 continue U(i,i) = 1.0 + U(i,i) 3 continue c return end c subroutine formW(m,ia,U,W) c c store U in W c dimension U(ia,m),W(ia,m) c do 2 i=1,m do 2 j=1,m W(i,j)=U(i,j) 2 continue c return end c subroutine trans(n,ia,m,A) c c compute transpose of the matrix A c dimension A(ia,n),C(4,3) c do 2 i=1,m do 2 j=1,n C(i,j) = A(j,i) 2 continue c do 3 i=1,m do 3 j=1,n A(i,j) = C(i,j) 3 continue c return end c subroutine mult(n,ia,m,B,A) c c compute the matrix-matrix product c dimension B(ia,m),A(ia,n),C(4,3) c do 3 j=1,n do 3 i=1,m sum = 0.0 do 2 k=1,m sum = sum + B(i,k)*A(k,j) 2 continue C(i,j) = sum 3 continue c do 4 i=1,m do 4 j=1,n A(i,j) = C(i,j) 4 continue c return end