C PDIRE3    SOURCE    CHAT      05/01/13    02:11:26     5004
c---------------------------------------------------------------------
c
      SUBROUTINE PDIRE3 (S, XL, V)
c
c=====================================================================
c                                                                    =
c   This routine computes the eigenvectors.                          =
c                                                                    =
c    Input :     xl  (3)   eigenvalues                               =
c                s   (6)   original matrix                           =
c    Output:     v (3,3)   eigenvectors (normalised)                 =
c                                                                    =
c    Note:   s  = (Sxx, Syy, Szz, Sxy, Sxz, Syz)                     =
c            xl = (S11, S22, S33)                                    =
c                                                                    =
c=====================================================================
      IMPLICIT INTEGER(I-N)
      real*8    s (6), xl (3), v (3,3)
c
      real*8    toler, comp, compi, comp2, a1, a2, a3, dot, xn
c
      parameter (toler = 1.0 d-07)
c
      comp2 = xl (1)**2 + xl (2)**2 + xl (3)**2
      compi = 1.0 d0 / sqrt(comp2)
c
      a1 = abs ((xl(2)-xl(3))*compi)
      a2 = abs ((xl(1)-xl(3))*compi)
      a3 = abs ((xl(1)-xl(2))*compi)
c=====----------------------------------------------------------------
c  case-a: three equal eigenvalues                                   |
c=====----------------------------------------------------------------
      if (a1 .lt. toler .and. a2 .lt. toler)  then
          v (1,1) = 1.0 d0
          v (2,2) = 1.0 d0
          v (3,3) = 1.0 d0
          v (1,2) = 0.0 d0
          v (1,3) = 0.0 d0
          v (2,1) = 0.0 d0
          v (2,3) = 0.0 d0
          v (3,1) = 0.0 d0
          v (3,2) = 0.0 d0
          return
      end if
c=====----------------------------------------------------------------
c  case-b: two equal eigenvalues                                     |
c=====----------------------------------------------------------------
      if (a1 .lt. toler)  then
          call pvecto (s, xl(1), v(1,1), comp2)
          call vecbas (v(1,1), v(1,2), v(1,3))
          return
      end if
      if (a2 .lt. toler)  then
          call pvecto (s, xl(2), v(1,2), comp2)
          call vecbas (v(1,2), v(1,3), v(1,1))
          return
      end if
      if (a3 .lt. toler)  then
          call pvecto (s, xl (3), v(1,3), comp2)
          call vecbas (v(1,3), v(1,1), v(1,2))
          return
      end if
c=====----------------------------------------------------------------
c  case-c: three different eigenvalues                               |
c=====----------------------------------------------------------------
      call pvecto (s, xl(1), v(1,1), comp2)
      call pvecto (s, xl(2), v(1,2), comp2)
c
      dot = v (1,1) * v (1,2) + v (2,1) * v (2,2) + v (3,1) * v (3,2)
c
      v (1,2) = v (1,2) - dot * v (1,1)
      v (2,2) = v (2,2) - dot * v (2,1)
      v (3,2) = v (3,2) - dot * v (3,1)
c
      xn = 1.0 d0 / SQRT (v(1,2)**2+v(2,2)**2+v(3,2)**2)
c
      v (1,2) = v (1,2) * xn
      v (2,2) = v (2,2) * xn
      v (3,2) = v (3,2) * xn
      v (1,3) = v (2,1) * v (3,2) - v (3,1) * v (2,2)
      v (2,3) = v (3,1) * v (1,2) - v (1,1) * v (3,2)
      v (3,3) = v (1,1) * v (2,2) - v (2,1) * v (1,2)
c
      return
      end