Sources Esope | Cast3M

dvdiag
C DVDIAG    SOURCE    KICH      19/09/26    21:15:05     10311          
      subroutine dvdiag(a,d)
c
c     ***************
c     diagonalisation
c     ***************
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension sig(6),a(3,3),d(3),e(3)
      n  = 3
      np = 3
        do 18 i=n,2,-1
          l=i-1
          h=0.
          scale=0.
          if(l.gt.1)then
            do 11 k=1,l
              scale=scale+abs(a(i,k))
 11         continue
            if(scale.eq.0.)then
              e(i)=a(i,l)
            else
              do 12 k=1,l
                a(i,k)=a(i,k)/scale
                h=h+a(i,k)**2
 12           continue
              f=a(i,l)
              g=-sign(sqrt(h),f)
              e(i)=scale*g
              h=h-f*g
              a(i,l)=f-g
              f=0.
              do 15 j=1,l
                a(j,i)=a(i,j)/h
                g=0.
                do 13 k=1,j
                  g=g+a(j,k)*a(i,k)
 13             continue
                if(l.gt.j)then
                  do 14 k=j+1,l
                    g=g+a(k,j)*a(i,k)
 14               continue
                endif
                e(j)=g/h
                f=f+e(j)*a(i,j)
 15           continue
              hh=f/(h+h)
              do 17 j=1,l
                f=a(i,j)
                g=e(j)-hh*f
                e(j)=g
                do 16 k=1,j
                  a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
 16             continue
 17           continue
            endif
          else
            e(i)=a(i,l)
          endif
          d(i)=h
 18     continue
      d(1)=0.
      e(1)=0.
      do 23 i=1,n
        l=i-1
        if(d(i).ne.0.)then
          do 21 j=1,l
            g=0.
            do 19 k=1,l
              g=g+a(i,k)*a(k,j)
 19         continue
            do 20 k=1,l
              a(k,j)=a(k,j)-g*a(k,i)
 20         continue
 21       continue
        endif
        d(i)=a(i,i)
        a(i,i)=1.
        if(l.ge.1)then
          do 22 j=1,l
            a(i,j)=0.
            a(j,i)=0.
 22       continue
        endif
 23   continue
        do 911 i=2,n
          e(i-1)=e(i)
 911     continue
        e(n)=0.
        do 915 l=1,n
          iter=0
 91        do 912 m=l,n-1
            dd=abs(d(m))+abs(d(m+1))
            if (abs(e(m))+dd.eq.dd) go to 92
 912       continue
          m=n
 92        if(m.ne.l)then
            if(iter.eq.30) then
              write(6,*) 'too many iterations'
              call erreur(5)
              return
            endif
            iter=iter+1
            g=(d(l+1)-d(l))/(2.*e(l))
            r=sqrt(g**2+1.)
            g=d(m)-d(l)+e(l)/(g+sign(r,g))
            s=1.
            c=1.
            p=0.
            do 914 i=m-1,l,-1
              f=s*e(i)
              b=c*e(i)
              if(abs(f).ge.abs(g))then
                c=g/f
                r=sqrt(c**2+1.)
                e(i+1)=f*r
                s=1./r
                c=c*s
              else
                s=f/g
                r=sqrt(s**2+1.)
                e(i+1)=g*r
                c=1./r
                s=s*c
              endif
              g=d(i+1)-p
              r=(d(i)-g)*s+2.*c*b
              p=s*r
              d(i+1)=g+p
              g=c*r-b
              do 913 k=1,n
                f=a(k,i+1)
                a(k,i+1)=s*a(k,i)+c*f
                a(k,i)=c*a(k,i)-s*f
 913           continue
 914         continue
            d(l)=d(l)-p
            e(l)=g
            e(m)=0.
            go to 91
          endif
 915     continue
      return
      end