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.5 c c Example of shifted QR-factorization c c c file: ex3s55.f c parameter (ia=4,m=4,n=4,itmax=4) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(1,j),j=1,m) /1.,2.,3.,4./ data (A(2,j),j=1,m) /4.,5.,6.,7./ data (A(3,j),j=1,m) /2.,1.,5.,0./ data (A(4,j),j=1,m) /4.,2.,1.,0./ c print * print *,' Shifted QR-factorization example' print *,' Section 5.5, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c call hess(n,ia,A,V,U) print *,' After Hess' c do 2 nn=n,2,-1 print *,' Matrix A' call prtmtx(nn,ia,nn,A) call submtx(nn,ia,nn,A,itmax,Q,U,W,V) print *,' Eigenvalue=',A(nn,nn) 2 continue print *,' Eigenvalue=',A(1,1) stop end c subroutine submtx(n,ia,m,A,itmax,Q,U,W,V) c c shifted QR on deflated matrix c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c do 2 k=1,itmax print *,' k =',k z = A(n,n) call shiftA(n,ia,m,A,z) print *,' Shifted A-zI' call prtmtx(n,ia,m,A) call QRfac(n,ia,m,A,V,U,Q,W) c 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) c call copy(n,ia,m,A,W) call mult(n,ia,m,Q,W) print *,' Check A = QR + zI' call unshiftA(n,ia,m,W,z) call prtmtx(n,ia,m,W) c c A = R*Q + zI c call mult(n,ia,m,A,Q) call copy(n,ia,m,Q,A) call unshiftA(n,ia,m,A,z) print *,' Matrix A(k+1)' call prtmtx(n,ia,m,A) 2 continue c return end c subroutine shiftA(n,ia,m,A,z) c c form shifted matrix c dimension A(ia,n) c do 2 i=1,n a(i,i) = a(i,i) - z 2 continue c return end c subroutine unshiftA(n,ia,m,A,z) c c form unshifted matrix c dimension A(ia,n) c do 2 i=1,n a(i,i) = a(i,i) + z 2 continue c return end c subroutine hess(n,ia,A,V,U) c c compute Hessenberg matrix of A c dimension A(ia,n),U(ia,n),V(n) c do 2 k=1,n-2 nn = n - k call findU(nn,ia,A(k+1,k),V,U) call mult(nn+1,ia,nn,U,A(k+1,k)) print *,' Matrix UA' call prtmtx(n,ia,n,A) call trans(nn,ia,nn,U) call mult(nn,ia,n,A(1,k+1),U) call copy(nn,ia,n,U,A(1,k+1)) print *,' Matrix UAU*' call prtmtx(n,ia,n,A) 2 continue c return end c subroutine QRfac(n,ia,m,A,V,U,Q,W) c c compute 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-1 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 U times A, store into A 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 findU(n,ia,A,V,U) c c find matrix U c dimension A(ia,n),U(ia,n),V(n) c call findV(n,ia,n,A,V) print *,' V vector' call prtmtx(1,ia,n,V) call formU(n,ia,V,U) print *,' Matrix U' call prtmtx(n,ia,n,U) 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 print *,' beta=',beta,' alpha=',alpha 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 form matrix 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 form matrix 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(10,10) c do 2 i=1,m do 2 j=1,n C(i,j) = A(j,i) 2 continue 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(10,10) 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 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 copy A into 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 Q & R (positive diagonals in 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