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