pyldcm.f

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

 
C*********************************************************************
 
C...PYLDCM
C...Auxiliary to PYSIGH, for technicolor corrections to QCD 2 -> 2
C...processes.
 
      SUBROUTINE pyldcm(A,N,NP,INDX,D)
      IMPLICIT NONE
      INTEGER n,np,indx(n)
      REAL*8 d,tiny
      COMPLEX*16 a(np,np)
      parameter(tiny=1.0d-20)
      INTEGER i,imax,j,k
      REAL*8 aamax,vv(6),dum
      COMPLEX*16 sum,dumc
 
      d=1d0
      DO 110 i=1,n
        aamax=0d0
        DO 100 j=1,n
          IF (abs(a(i,j)).GT.aamax) aamax=abs(a(i,j))
  100   CONTINUE
        IF (aamax.EQ.0d0) CALL pyerrm(28,'(PYLDCM:) singular matrix')
        vv(i)=1d0/aamax
  110 CONTINUE
      DO 180 j=1,n
        DO 130 i=1,j-1
          sum=a(i,j)
          DO 120 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
  120     CONTINUE
          a(i,j)=sum
  130   CONTINUE
        aamax=0d0
        DO 150 i=j,n
          sum=a(i,j)
          DO 140 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
  140     CONTINUE
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          IF (dum.GE.aamax) THEN
            imax=i
            aamax=dum
          ENDIF
  150   CONTINUE
        IF (j.NE.imax)THEN
          DO 160 k=1,n
            dumc=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dumc
  160     CONTINUE
          d=-d
          vv(imax)=vv(j)
        ENDIF
        indx(j)=imax
        IF(abs(a(j,j)).EQ.0d0) a(j,j)=dcmplx(tiny,0d0)
        IF(j.NE.n)THEN
          DO 170 i=j+1,n
            a(i,j)=a(i,j)/a(j,j)
  170     CONTINUE
        ENDIF
  180 CONTINUE
 
      RETURN
      END