pythiaeRHIC/pythia/pyeig4.f

0001
0002 C*********************************************************************
0003
0004 C...PYEIG4
0005 C...Finds eigenvalues and eigenvectors to a 4 * 4 matrix.
0006 C...Specific application: mixing in neutralino sector.
0007
0008       SUBROUTINE PYEIG4(A,W,Z)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014
0015 C...Arrays: in call and local.
0016       DIMENSION A(4,4),W(4),Z(4,4),X(4),D(4,4),E(4)
0017
0018 C...Coefficients of fourth-degree equation from matrix.
0019 C...x**4 + b3 * x**3 + b2 * x**2 + b1 * x + b0 = 0.
0020       B3=-(A(1,1)+A(2,2)+A(3,3)+A(4,4))
0021       B2=0D0
0022       DO 110 I=1,3
0023         DO 100 J=I+1,4
0024           B2=B2+A(I,I)*A(J,J)-A(I,J)*A(J,I)
0025   100   CONTINUE
0026   110 CONTINUE
0027       B1=0D0
0028       B0=0D0
0029       DO 120 I=1,4
0030         I1=MOD(I,4)+1
0031         I2=MOD(I+1,4)+1
0032         I3=MOD(I+2,4)+1
0033         B1=B1+A(I,I)*(-A(I1,I1)*A(I2,I2)+A(I1,I2)*A(I2,I1)+
0034      &  A(I1,I3)*A(I3,I1)+A(I2,I3)*A(I3,I2))-
0035      &  A(I,I1)*A(I1,I2)*A(I2,I)-A(I,I2)*A(I2,I1)*A(I1,I)
0036         B0=B0+(-1D0)**(I+1)*A(1,I)*(
0037      &  A(2,I1)*(A(3,I2)*A(4,I3)-A(3,I3)*A(4,I2))+
0038      &  A(2,I2)*(A(3,I3)*A(4,I1)-A(3,I1)*A(4,I3))+
0039      &  A(2,I3)*(A(3,I1)*A(4,I2)-A(3,I2)*A(4,I1)))
0040   120 CONTINUE
0041
0042 C...Coefficients of third-degree equation needed for
0043 C...separation into two second-degree equations.
0044 C...u**3 + c2 * u**2 + c1 * u + c0 = 0.
0045       C2=-B2
0046       C1=B1*B3-4D0*B0
0047       C0=-B1**2-B0*B3**2+4D0*B0*B2
0048       CQ=C1/3D0-C2**2/9D0
0049       CR=C1*C2/6D0-C0/2D0-C2**3/27D0
0050       CQR=CQ**3+CR**2
0051
0052 C...Cases with one or three real roots.
0053       IF(CQR.GE.0D0) THEN
0054         S1=(CR+SQRT(CQR))**(1D0/3D0)
0055         S2=(CR-SQRT(CQR))**(1D0/3D0)
0056         U=S1+S2-C2/3D0
0057       ELSE
0058         SABS=SQRT(-CQ)
0059         THE=ACOS(CR/SABS**3)/3D0
0060         SRE=SABS*COS(THE)
0061         U=2D0*SRE-C2/3D0
0062       ENDIF
0063
0064 C...Find and solve two second-degree equations.
0065       P1=B3/2D0-SQRT(B3**2/4D0+U-B2)
0066       P2=B3/2D0+SQRT(B3**2/4D0+U-B2)
0067       Q1=U/2D0+SQRT(U**2/4D0-B0)
0068       Q2=U/2D0-SQRT(U**2/4D0-B0)
0069       IF(ABS(P1*Q1+P2*Q2-B1).LT.ABS(P1*Q2+P2*Q1-B1)) THEN
0070         QSAV=Q1
0071         Q1=Q2
0072         Q2=QSAV
0073       ENDIF
0074       X(1)=-P1/2D0+SQRT(P1**2/4D0-Q1)
0075       X(2)=-P1/2D0-SQRT(P1**2/4D0-Q1)
0076       X(3)=-P2/2D0+SQRT(P2**2/4D0-Q2)
0077       X(4)=-P2/2D0-SQRT(P2**2/4D0-Q2)
0078
0079 C...Order eigenvalues in asceding mass.
0080       W(1)=X(1)
0081       DO 150 I1=2,4
0082         DO 130 I2=I1-1,1,-1
0083           IF(ABS(X(I1)).GE.ABS(W(I2))) GOTO 140
0084           W(I2+1)=W(I2)
0085   130   CONTINUE
0086   140   W(I2+1)=X(I1)
0087   150 CONTINUE
0088
0089 C...Find equation system for eigenvectors.
0090       DO 250 I=1,4
0091         DO 170 J1=1,4
0092           D(J1,J1)=A(J1,J1)-W(I)
0093           DO 160 J2=J1+1,4
0094             D(J1,J2)=A(J1,J2)
0095             D(J2,J1)=A(J2,J1)
0096   160     CONTINUE
0097   170   CONTINUE
0098
0099 C...Find largest element in matrix.
0100         DAMAX=0D0
0101         DO 190 J1=1,4
0102           DO 180 J2=1,4
0103             IF(ABS(D(J1,J2)).LE.DAMAX) GOTO 180
0104             JA=J1
0105             JB=J2
0106             DAMAX=ABS(D(J1,J2))
0107   180     CONTINUE
0108   190   CONTINUE
0109
0110 C...Subtract others by multiple of row selected above.
0111         DAMAX=0D0
0112         DO 210 J3=JA+1,JA+3
0113           J1=J3-4*((J3-1)/4)
0114           RL=D(J1,JB)/D(JA,JB)
0115           DO 200 J2=1,4
0116             D(J1,J2)=D(J1,J2)-RL*D(JA,J2)
0117             IF(ABS(D(J1,J2)).LE.DAMAX) GOTO 200
0118             JC=J1
0119             JD=J2
0120             DAMAX=ABS(D(J1,J2))
0121   200     CONTINUE
0122   210   CONTINUE
0123
0124 C...Do one more subtraction of a row.
0125         DAMAX=0D0
0126         DO 230 J3=JC+1,JC+3
0127           J1=J3-4*((J3-1)/4)
0128           IF(J1.EQ.JA) GOTO 230
0129           RL=D(J1,JD)/D(JC,JD)
0130           DO 220 J2=1,4
0131             IF(J2.EQ.JB) GOTO 220
0132             D(J1,J2)=D(J1,J2)-RL*D(JC,J2)
0133             IF(ABS(D(J1,J2)).LE.DAMAX) GOTO 220
0134             JE=J1
0135             DAMAX=ABS(D(J1,J2))
0136   220     CONTINUE
0137   230   CONTINUE
0138
0139 C...Construct unnormalized eigenvector.
0140         JF1=JD+1-4*(JD/4)
0141         JF2=JD+2-4*((JD+1)/4)
0142         IF(JF1.EQ.JB) JF1=JD+3-4*((JD+2)/4)
0143         IF(JF2.EQ.JB) JF2=JD+3-4*((JD+2)/4)
0144         E(JF1)=-D(JE,JF2)
0145         E(JF2)=D(JE,JF1)
0146         E(JD)=-(D(JC,JF1)*E(JF1)+D(JC,JF2)*E(JF2))/D(JC,JD)
0147         E(JB)=-(D(JA,JF1)*E(JF1)+D(JA,JF2)*E(JF2)+D(JA,JD)*E(JD))/
0148      &  D(JA,JB)
0149
0150 C...Normalize and fill in final array.
0151         EA=SQRT(E(1)**2+E(2)**2+E(3)**2+E(4)**2)
0152         SGN=(-1D0)**INT(PYR(0)+0.5D0)
0153         DO 240 J=1,4
0154           Z(I,J)=SGN*E(J)/EA
0155   240   CONTINUE
0156   250 CONTINUE
0157
0158       RETURN
0159       END