pyeig4.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyeig4.f

 
C*********************************************************************
 
C...PYEIG4
C...Finds eigenvalues and eigenvectors to a 4 * 4 matrix.
C...Specific application: mixing in neutralino sector.
 
      SUBROUTINE pyeig4(A,W,Z)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
 
C...Arrays: in call and local.
      dimension a(4,4),w(4),z(4,4),x(4),d(4,4),e(4)
 
C...Coefficients of fourth-degree equation from matrix.
C...x**4 + b3 * x**3 + b2 * x**2 + b1 * x + b0 = 0.
      b3=-(a(1,1)+a(2,2)+a(3,3)+a(4,4))
      b2=0d0
      DO 110 i=1,3
        DO 100 j=i+1,4
          b2=b2+a(i,i)*a(j,j)-a(i,j)*a(j,i)
  100   CONTINUE
  110 CONTINUE
      b1=0d0
      b0=0d0
      DO 120 i=1,4
        i1=mod(i,4)+1
        i2=mod(i+1,4)+1
        i3=mod(i+2,4)+1
        b1=b1+a(i,i)*(-a(i1,i1)*a(i2,i2)+a(i1,i2)*a(i2,i1)+
     &  a(i1,i3)*a(i3,i1)+a(i2,i3)*a(i3,i2))-
     &  a(i,i1)*a(i1,i2)*a(i2,i)-a(i,i2)*a(i2,i1)*a(i1,i)
        b0=b0+(-1d0)**(i+1)*a(1,i)*(
     &  a(2,i1)*(a(3,i2)*a(4,i3)-a(3,i3)*a(4,i2))+
     &  a(2,i2)*(a(3,i3)*a(4,i1)-a(3,i1)*a(4,i3))+
     &  a(2,i3)*(a(3,i1)*a(4,i2)-a(3,i2)*a(4,i1)))
  120 CONTINUE
 
C...Coefficients of third-degree equation needed for
C...separation into two second-degree equations.
C...u**3 + c2 * u**2 + c1 * u + c0 = 0.
      c2=-b2
      c1=b1*b3-4d0*b0
      c0=-b1**2-b0*b3**2+4d0*b0*b2
      cq=c1/3d0-c2**2/9d0
      cr=c1*c2/6d0-c0/2d0-c2**3/27d0
      cqr=cq**3+cr**2
 
C...Cases with one or three real roots.
      IF(cqr.GE.0d0) THEN
        s1=(cr+sqrt(cqr))**(1d0/3d0)
        s2=(cr-sqrt(cqr))**(1d0/3d0)
        u=s1+s2-c2/3d0
      ELSE
        sabs=sqrt(-cq)
        the=acos(cr/sabs**3)/3d0
        sre=sabs*cos(the)
        u=2d0*sre-c2/3d0
      ENDIF
 
C...Find and solve two second-degree equations.
      p1=b3/2d0-sqrt(b3**2/4d0+u-b2)
      p2=b3/2d0+sqrt(b3**2/4d0+u-b2)
      q1=u/2d0+sqrt(u**2/4d0-b0)
      q2=u/2d0-sqrt(u**2/4d0-b0)
      IF(abs(p1*q1+p2*q2-b1).LT.abs(p1*q2+p2*q1-b1)) THEN
        qsav=q1
        q1=q2
        q2=qsav
      ENDIF
      x(1)=-p1/2d0+sqrt(p1**2/4d0-q1)
      x(2)=-p1/2d0-sqrt(p1**2/4d0-q1)
      x(3)=-p2/2d0+sqrt(p2**2/4d0-q2)
      x(4)=-p2/2d0-sqrt(p2**2/4d0-q2)
 
C...Order eigenvalues in asceding mass.
      w(1)=x(1)
      DO 150 i1=2,4
        DO 130 i2=i1-1,1,-1
          IF(abs(x(i1)).GE.abs(w(i2))) goto 140
          w(i2+1)=w(i2)
  130   CONTINUE
  140   w(i2+1)=x(i1)
  150 CONTINUE
 
C...Find equation system for eigenvectors.
      DO 250 i=1,4
        DO 170 j1=1,4
          d(j1,j1)=a(j1,j1)-w(i)
          DO 160 j2=j1+1,4
            d(j1,j2)=a(j1,j2)
            d(j2,j1)=a(j2,j1)
  160     CONTINUE
  170   CONTINUE
 
C...Find largest element in matrix.
        damax=0d0
        DO 190 j1=1,4
          DO 180 j2=1,4
            IF(abs(d(j1,j2)).LE.damax) goto 180
            ja=j1
            jb=j2
            damax=abs(d(j1,j2))
  180     CONTINUE
  190   CONTINUE
 
C...Subtract others by multiple of row selected above.
        damax=0d0
        DO 210 j3=ja+1,ja+3
          j1=j3-4*((j3-1)/4)
          rl=d(j1,jb)/d(ja,jb)
          DO 200 j2=1,4
            d(j1,j2)=d(j1,j2)-rl*d(ja,j2)
            IF(abs(d(j1,j2)).LE.damax) goto 200
            jc=j1
            jd=j2
            damax=abs(d(j1,j2))
  200     CONTINUE
  210   CONTINUE
 
C...Do one more subtraction of a row.
        damax=0d0
        DO 230 j3=jc+1,jc+3
          j1=j3-4*((j3-1)/4)
          IF(j1.EQ.ja) goto 230
          rl=d(j1,jd)/d(jc,jd)
          DO 220 j2=1,4
            IF(j2.EQ.jb) goto 220
            d(j1,j2)=d(j1,j2)-rl*d(jc,j2)
            IF(abs(d(j1,j2)).LE.damax) goto 220
            je=j1
            damax=abs(d(j1,j2))
  220     CONTINUE
  230   CONTINUE
 
C...Construct unnormalized eigenvector.
        jf1=jd+1-4*(jd/4)
        jf2=jd+2-4*((jd+1)/4)
        IF(jf1.EQ.jb) jf1=jd+3-4*((jd+2)/4)
        IF(jf2.EQ.jb) jf2=jd+3-4*((jd+2)/4)
        e(jf1)=-d(je,jf2)
        e(jf2)=d(je,jf1)
        e(jd)=-(d(jc,jf1)*e(jf1)+d(jc,jf2)*e(jf2))/d(jc,jd)
        e(jb)=-(d(ja,jf1)*e(jf1)+d(ja,jf2)*e(jf2)+d(ja,jd)*e(jd))/
     &  d(ja,jb)
 
C...Normalize and fill in final array.
        ea=sqrt(e(1)**2+e(2)**2+e(3)**2+e(4)**2)
        sgn=(-1d0)**int(pyr(0)+0.5d0)
        DO 240 j=1,4
          z(i,j)=sgn*e(j)/ea
  240   CONTINUE
  250 CONTINUE
 
      RETURN
      END