Numérotation des lignes :

C EIG35     SOURCE    CHAT      05/01/12    23:28:43     5004c**************************************************************************      SUBROUTINE EIG35(V,D,ROT)c**************************************************************************      IMPLICIT INTEGER(I-N)      integer rot, its, i,j,k      real*8  g,h,aij,sm,thresh,t,c,s,tau,v(3,3),d(3),a(3),b(3),z(3)      a(1) = v(1,2)      a(2) = v(2,3)      a(3) = v(3,1)      do i = 1,3        d(i) = v(i,i)        b(i) = d(i)        z(i) = 0.0d0        do j = 1,3          v(i,j) = 0.0d0        end do        v(i,i) = 1.d0      end do      rot = 0      do its = 1,50c.... set convergence test and threshold        sm = abs(a(1)) + abs(a(2)) + abs(a(3))        if (sm.eq.0.0d0) return        if(its.lt.4) then          thresh = 0.011d0*sm        else          thresh = 0.0d0        end ifc.... perform sweeps for rotations        do i = 1,3          j = mod(i,3) + 1          k = mod(j,3) + 1          aij  = a(i)          g    = 100.d0*abs(aij)          if(abs(d(i))+g.ne.abs(d(i)) .or.     1       abs(d(j))+g.ne.abs(d(j))) then            if(abs(aij).gt.thresh) then              a(i) = 0.0d0              h    = d(j) - d(i)              if(abs(h)+g.eq.abs(h)) then                t = aij/h              else                t = sign(2.d0,h/aij)/(abs(h/aij)+sqrt(4.d0+(h/aij)**2))              endifc.... set rotation parameters              c    = 1.d0/sqrt(1.d0+t*t)              s    = t*c              tau  = s/(1.d0+c)c.... rotate diagonal terms              h    = t*aij              z(i) = z(i) - h              z(j) = z(j) + h              d(i) = d(i) - h              d(j) = d(j) + hc.... rotate off-diagonal terms              h    = a(j)              g    = a(k)              a(j) = h + s*(g - h*tau)              a(k) = g - s*(h + g*tau)c.... rotate eigenvectors              do k = 1,3                g      = v(k,i)                h      = v(k,j)                v(k,i) = g - s*(h + g*tau)                v(k,j) = h + s*(g - h*tau)              end do              rot = rot + 1            end if          else            a(i) = 0.0d0          end if        end doc.... update the diagonal terms        do i = 1,3          b(i) = b(i) + z(i)          d(i) = b(i)          z(i) = 0.0d0        end do      end do      end

© Cast3M 2003 - Tous droits réservés.
Mentions légales