jacobi.F90

Go to the documentation of this file.

!
! $Id: jacobi.F90 1907 2013-11-26 13:10:46Z lguez $
!
      SUBROUTINE jacobi(A,N,NP,D,V,NROT)
      implicit none
! Arguments:
      integer,intent(in) :: N
      integer,intent(in) :: NP
      integer,intent(out) :: NROT
      real,intent(inout) :: A(np,np)
      real,intent(out) :: D(np)
      real,intent(out) :: V(np,np)

! local variables:
      integer :: IP,IQ,I,J
      real :: SM,TRESH,G,H,T,THETA,C,S,TAU
      real :: B(n)
      real :: Z(n)
      
      DO ip=1,n
        DO iq=1,n
          v(ip,iq)=0.
        ENDDO
        v(ip,ip)=1.
      ENDDO
      DO ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
      ENDDO
      nrot=0
      DO i=1,50 ! 50? I suspect this should be NP
                !     but convergence is fast enough anyway
        sm=0.
        DO ip=1,n-1
          DO iq=ip+1,n
            sm=sm+abs(a(ip,iq))
          ENDDO
        ENDDO
        IF(sm.EQ.0.)RETURN
        IF(i.LT.4)THEN
          tresh=0.2*sm/n**2
        ELSE
          tresh=0.
        ENDIF
        DO ip=1,n-1
          DO iq=ip+1,n
            g=100.*abs(a(ip,iq))
            IF((i.GT.4).AND.(abs(d(ip))+g.EQ.abs(d(ip))) &
               .AND.(abs(d(iq))+g.EQ.abs(d(iq))))THEN
              a(ip,iq)=0.
            ELSE IF(abs(a(ip,iq)).GT.tresh)THEN
              h=d(iq)-d(ip)
              IF(abs(h)+g.EQ.abs(h))THEN
                t=a(ip,iq)/h
              ELSE
                theta=0.5*h/a(ip,iq)
                t=1./(abs(theta)+sqrt(1.+theta**2))
                IF(theta.LT.0.)t=-t
              ENDIF
              c=1./sqrt(1+t**2)
              s=t*c
              tau=s/(1.+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.
              DO j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
             ENDDO
              DO j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
              ENDDO
              DO j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
              ENDDO
              DO j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
              ENDDO
              nrot=nrot+1
            ENDIF
          ENDDO
        ENDDO
        DO ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.
        ENDDO
      ENDDO ! of DO I=1,50
      stop 'Jacobi: 50 iterations should never happen'
      RETURN
      END SUBROUTINE jacobi