view release on metacpan or  search on metacpan

reference/qralg/qralg.f  view on Meta::CPAN

************************************************************************
*
*  Solve the real coefficient monic algebraic equation by the QR-method.
*
************************************************************************
      subroutine qr_algeq_solver(n,c,macheps,
     &                 zr,zi,detil,
     &                 a,cnt,rtn_code)
      implicit none
      integer n
      double precision c(0:n-1)
      double precision macheps
*
      double precision zr(n),zi(n)
      double precision detil
      double precision a(n,n)
      integer cnt(n)
      integer rtn_code
Comments:
*     n       : degree of the monic polynomial.
*     c       : non-principal coefficients of the polynomial.
*               i-th degree coefficients is stored in c(i).
*     macheps : double precision machine epsilon.
*     zr,zi   : Re and Im part of output roots.
*     detil   : The accuracy hint.
*     a       : work matrix.
*     cnt     : work area for counting the qr-iterations.
*     rtn_code: return code from hqr_eigen_hessenberg.
*end comments;
*================
      integer i
      integer iter
      double precision afnorm
*
      if(n.le.0)then
        print *, 'qr_alg_solver: wrong arguments. (n<=0)'
        rtn_code=3
        return
      endif
*
* Build the companion matrix A.
      call build_companion(n,a,c)
*
* Balancing the A itself.
      call balance_companion(n,a)
*
* Compute the Frobenius norm of the balanced companion matrix A.
      call frobenius_norm_companion(n,a,afnorm)
*
* QR iterations from A.
      call hqr_eigen_hessenberg(n,a,macheps,  zr,zi,cnt,rtn_code)
      if(rtn_code.ne.0)then
        print*, 'in calling hqr_eigen_hessenberg abnormal condition'
        print*, 'rtn_code=',rtn_code
        if(rtn_code.eq.1) print *, 'matrix is completely zero.'
        if(rtn_code.eq.2) print *, 'QR iteration does not converged.'
        if(rtn_code.gt.3) print *, 'arguments violate the condition.'
        return
      endif
*
* Count the total QR iteration.
      iter=0
      do i=1,n
        if(cnt(i).gt.0)iter=iter+cnt(i)
      enddo 
      print *, 'iteration ',iter,' times.'
*
* Calculate the accuracy hint.
      detil=macheps*n*iter*afnorm
      print 600, 'O(del-e-tilda)=',detil
  600 format(1x,a,1pe8.2)
*
      end
************************************************************************
*
*  build the Companion Matrix of the polynomial.
*
************************************************************************
      subroutine build_companion(n,a,c)
      implicit none
      integer n
      double precision a(n,n)
      double precision c(0:n-1)
*
      integer i,j
*
      do j=1,n
        do i=1,n
          a(i,j)=0.d0
        enddo
      enddo
*
      do i=2,n
        a(i,i-1)=1.d0
      enddo
*
      do i=1,n
        a(i,n)=-c(i-1)
      enddo
*
      end
************************************************************************
*
*     Blancing the unsymmetric matrix A.
*
*  This fortran code is based on the Algol code "balance" from paper:
*   "Balancing a Matrixfor Calculation of Eigenvalues and Eigenvectors"
*   by B.N.Parlett and C.Reinsch, Numer. Math. 13, 293-304(1969).
*
*  Note: The only non-zero elements of the companion matrix are touched.

reference/qralg/qralg.f  view on Meta::CPAN

                a(j,i)=a(j,i)*f
              enddo
            endif
C end specific code.
          endif
* Comment: 
*     the j joops may be done by exponent modification
*     in machine Language;
  100   enddo
      if(noconv) goto 1
      end
************************************************************************
*
*  Calculate the Frobenius norm of the companion-like matrix.
*  Note: The only non-zero elements of the companion matrix are touched.
*
************************************************************************
      subroutine frobenius_norm_companion(n,a,afnorm)
      implicit none
      integer n
      double precision a(n,n)
      double precision afnorm
*
      integer i,j
      double precision sum
*
C     sum=0.d0
C     do j=1,n
C       do i=1,n
C         sum=sum+a(i,j)**2
C       enddo
C     enddo
*
      sum=0.d0
      do j=1,n-1
        sum=sum+a(j+1,j)**2
      enddo
      do i=1,n
        sum=sum+a(i,n)**2
      enddo
*
      afnorm=sqrt(sum)
      end
************************************************************************
*
*     Eigenvalue Computation by the Householder QR method 
*     for the Real Hessenberg matrix.
* This fortran code is based on the Algol code "hqr" from the paper:
*       "The QR Algorithm for Real Hessenberg Matrices"
*       by R.S.Martin, G.Peters and J.H.Wilkinson,
*       Numer. Math. 14, 219-231(1970).
*
************************************************************************
      subroutine hqr_eigen_hessenberg(n0,h,macheps,  wr,wi,cnt,rtn_code)
      implicit none
*Comment: Finds the eigenvalues of a real upper Hessenberg matrix,
*         H, stored in the array h(1:n0,1:n0), and stores the real
*         parts in the array wr(1:n0) and the imaginary parts in the
*         array wi(1:n0). macheps is the relative machine precision.
*         The procedure fails if any eigenvalue takes more than 
*         30 iterations.
*
      integer n0
      double precision h(n0,n0)
      double precision macheps
*
      double precision wr(n0),wi(n0)
      integer cnt(n0)
      integer rtn_code
*===
      integer i,j,k,l,m,na,its
      double precision p,q,r,s,t,w,x,y,z
      logical notlast
      integer n
*===
* Note: n is changing in this subroutine.
      n=n0
*
      rtn_code=0
      t=0.d0
*label nextw:
    1 continue
      if(n.eq.0)goto 9
      its=0
      na=n-1
*Comment: look for single small sub-diagonal element;
*label nextit:
    2 continue
      do l=n,2,-1
       if(abs(h(l,l-1)).le.macheps*(abs(h(l-1,l-1))+abs(h(l,l))))then
         go to 3
       endif
      enddo
      l=1
*label cont1:
    3 continue
      x=h(n,n)
      if(l.eq.n)goto 6
      y=h(na,na)
      w=h(n,na)*h(na,n)
      if(l.eq.na)goto 7
      if(its.eq.30)then
        rtn_code=1
        return
      endif
      if(its.eq.10.or.its.eq.20)then
* Comment: form exceptional shift;
        t=t+x
        do i=1,n
          h(i,i)=h(i,i)-x
        enddo
        s=abs(h(n,na))+abs(h(na,n-2))
        y=0.75*s
        x=y
        w=-0.4375*s**2
      endif
      its=its+1
* Comment: look for two consecutive small sub-diagonal 
*          elements;
      do m=n-2,l,-1
        z=h(m,m)