c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole., 1990 c c Section 5.5 c c Example of QR-factorization to reduce given matrix to upper c Hessenberg form c c c file: ex2s55.f c parameter (n=4,ia=4,m=4,itmax=10) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(i,1),i=1,m) /1.,-6.,0.,0./ data (A(i,2),i=1,m) /-5.,8.5556,-5.9182,0./ data (A(i,3),i=1,m) /-1.5395,0.10848,-2.1689,-0.14276/ data (A(i,4),i=1,m) /-1.2767,3.6986,-1.1428,3.6133/ c print * print *,' QR-factorization example' print *,' Section 5.5, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c do 2 k=1,itmax call QRfac(n,ia,m,A,V,U,Q,W) call scale(n,ia,m,Q,A) print *,' Matrix Q' call prtmtx(m,ia,m,Q) print *,' Matrix R' call prtmtx(n,ia,m,A) call copy(m,ia,m,Q,W) call mult(n,ia,m,A,W) c call mult(n,ia,m,Q,A) print *,' Check A = QR' call prtmtx(n,ia,m,A) c call copy(n,ia,m,W,A) print *,' Matrix A(k+1)' call prtmtx(n,ia,m,A) 2 continue 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,3 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) 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,4) 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,4) 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 c subroutine copy(n,ia,m,A,B) c c store A in B c dimension A(ia,n),B(ia,n) c do 2 i=1,m do 2 j=1,n B(i,j) = A(i,j) 2 continue c return end c subroutine scale(n,ia,m,Q,R) c c scale the matrix Q & R c dimension Q(ia,m),R(ia,n),D(4) c do 2 k=1,m D(k) = sign(1.0,R(k,k)) 2 continue do 3 j=1,n do 3 i=1,m Q(i,j) = Q(i,j)*D(i) 3 continue do 4 i=1,m do 4 j=1,n R(i,j) = D(j)*R(i,j) 4 continue c return end