This is the mail archive of the gcc-prs@gcc.gnu.org mailing list for the GCC project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

optimization/7563: g77 3.1.1 internal error: Segmentation Fault at -O2 -funroll-loops -m64


>Number:         7563
>Category:       optimization
>Synopsis:       g77 3.1.1 internal error: Segmentation Fault at -O2 -funroll-loops -m64
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    unassigned
>State:          open
>Class:          ice-on-legal-code
>Submitter-Id:   net
>Arrival-Date:   Fri Aug 09 16:56:00 PDT 2002
>Closed-Date:
>Last-Modified:
>Originator:     David G Hough
>Release:        3.1.1
>Organization:
Sun Microsystems
>Environment:
System: SunOS dgh 5.8 Generic_108528-15 sun4u sparc SUNW,Ultra-60
Architecture: sun4

host: sparc-sun-solaris2.8
build: sparc-sun-solaris2.8
target: sparc-sun-solaris2.8
configured with: /net/validator/export/set/validator18/hwvalid/compdir/gcc-3.1.1/configure --prefix=/net/validator/export/set/validator18/hwvalid/compdir/g311 --enable-languages=c,f77
>Description:
g77 -c qcd.f -O2 -funroll-loops -m64
qcd.f: In program `qcd2':
qcd.f:290: internal error: Segmentation Fault

>How-To-Repeat:

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                THIS SOFTWARE COPYRIGHTED BY CALTECH, 1983, 1986, 1987
C
C                PURE GAUGE SU(3) MONTE CARLO
C                THIS CODE HAS A LONG HISTORY:
C                    EARLY VERSIONS OF THIS DEVELOPED FOR THE MARK I
C                    HYPERCUBE BY STEVE OTTO AND PAUL STOLORZ (19
C                    JON FLOWER AND STEVE OTTO EXTENDED AND IMPROVED THE
C                    CODE FOR THE MARK II MACHINES. (1984-86)
C                    THE CURRENT INCARNATION IS A "PORT" TO THE NCUBE
C                    DONE BY CLIVE BAILLIE, ED FELTEN, STEVE OTTO,
C                    AND DAVID WALKER  (1987)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PURE GAUGE SU(3) THEORY IN 3+1 DIMENSIONS
C               LINK VARIABLE UPDATING DONE VIA CABIBBO-MARINARI HEAT BA
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE
C     AND DAVID WALKER, SUMMER 1987
C
      PROGRAM qcd2
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      PARAMETER (mbar=18*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
      PARAMETER (lismax=(l1*l2*l3*l4)/2+1)
      COMMON /colcom/ lisred(lismax) , lisblk(lismax) , nred , nblack
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      DOUBLEPRECISION ubar(mbar) , sbar(mbar)
      COMMON /comubr/ ubar , sbar
C
      INTEGER numstg
      COMMON /comfil/ numstg
C
      INTEGER i , j , k , numoff , avgnum , flg
      INTEGER r , t , numbar , dir , rmin , rmax , tmin(13) , tmax(13)
      INTEGER r1max , r2max , r1 , r2 , r3 , tmx , tmn
      CHARACTER c
      INTEGER sweeps , pinval
C
      DOUBLEPRECISION fc(2)
c                      isvalid, unit 8 changed to unit 0, etc. for /valid
      CHARACTER*7 isvalid
C
       write(6,*) ' BEGIN perfect.qcd2.F,1.12,99/03/12 '
C

        call validflush
	call startvalidtime

C
      CALL prnset(42)
C
c      WRITE(8,201)
c201   FORMAT(1X,'OUTPUT FOR PERFECT CLUB BENCHMARK: QCD'/
c     *       1X,'$Revision: 1.12 $ $Author: lpointer $'/)
c      WRITE(8,202)
c202   FORMAT(1X,'VALIDATION PARAMETERS:'/)
C
      WRITE (6,*) 'PURE GAUGE QCD USING PSEUDO HEAT BATH'
      WRITE (6,*) '-------------------------------------'
      WRITE (6,*) ' '
      WRITE (6,*) 'ENTER COMMAND (? FOR HELP)'
      WRITE (6,*) ' '
  100 CONTINUE
      READ (5,99008) c
C
      IF ( c.EQ.'R' ) THEN
         CALL setsed()
      ELSEIF ( c.EQ.'D' ) THEN
         CALL dmpsed()
      ELSEIF ( c.EQ.'S' ) THEN
         CALL setpar()
      ELSEIF ( c.EQ.'P' ) THEN
         CALL projec()
      ELSEIF ( c.EQ.'U' ) THEN
         WRITE (6,*) 
     &              'ENTER NUMBER OF TEMPORAL LINK HITS TO AVERAGE OVER'
         READ (5,*) avgnum
         WRITE (6,*) avgnum
C       CLOCK = SECNDS(0)
         CALL linkbr(avgnum,3,ubar)
C       WRITE(6,101) SECNDS(0)-CLOCK
         WRITE (6,*) 'UBARS MADE.'
      ELSEIF ( c.EQ.'M' ) THEN
         WRITE (6,*) 'ENTER DIR (0<->X,1<->Y,2<->Z), RMIN, RMAX'
         READ (5,*) dir , rmin , rmax
         WRITE (6,*) dir , rmin , rmax
         DO r = rmin , rmax
            WRITE (6,*) 'FOR R = ' , r , ' ENTER TMIN, TMAX '
            READ (5,*) tmin(r+1) , tmax(r+1)
            WRITE (6,*) tmin(r+1) , tmax(r+1)
         ENDDO
C
         DO r = 1 , 13
            DO t = 1 , 13
               which(r,t) = 0
            ENDDO
         ENDDO
C
         DO r = rmin , rmax
            DO t = 1 , 13
               IF ( (t.GE.tmin(r+1)) .AND. (t.LE.tmax(r+1)) )
     &               which(r+1,t+1) = 1
            ENDDO
         ENDDO
C
         WRITE (6,*) 'ENTER NUMBER OF SPATIAL LINK HITS TO AVERAGE OVER'
         READ (5,*) avgnum
         WRITE (6,*) avgnum
         WRITE (6,*) 'MEASURING WILSON (ON-AXIS) LOOPS IN DIRECTION ' , 
     &               dir
C       CLOCK = SECNDS(0)
         CALL measur(dir,rmin,rmax,avgnum)
C       WRITE(6,102) SECNDS(0)-CLOCK
         fc(1) = numstg
         DO r = rmin , rmax
            DO t = tmin(r+1) , tmax(r+1)
               fc(2) = wloop(r+1,t+1)
               fc(2) = fc(2)/fc(1)
               WRITE (6,1) r , t , fc(2)
 1               format ('W' , i3 , 'X' , i3 , ':' , e19.9)
c            WRITE(8,*) 'W',R,'X',T,':',FC(2)
C
            ENDDO
         ENDDO
C
         ival = 0
         IF ( abs((fc(2)-.633259)/.633259).GE.1E-2 ) ival = 1
C
      ELSEIF ( c.EQ.'I' ) THEN
         numoff = 0
      ELSEIF ( c.EQ.'A' ) THEN
         WRITE (6,*) 'ENTER R1,R2,R3,TMN,TMX'
         READ (5,*) r1 , r2 , r3 , tmn , tmx
         WRITE (6,*) r1 , r2 , r3 , tmn , tmx
         WRITE (6,*) 'MEASURING BARYON LOOPS ...'
C       CLOCK = SECNDS(0)
         DO t = tmn , tmx
            CALL qqqmea(r1,r2,r3,t)
         ENDDO
C       WRITE(6,103) SECNDS(0)-CLOCK
C


	write(6,11) fc(2),0.633259d0
 11	format(' fc(2) ',f10.6,' vs ', f10.6)

         IF ( ival.EQ.0 ) THEN
c           WRITE(8,203)
99001       FORMAT (//1X,'RESULTS FOR THIS RUN ARE:   VALID')
            isvalid = 'VALID'

	write(6,*) 'qcd2=LG results are VALID'


         ELSE
c           WRITE(8,204)
            isvalid = 'INVALID'
99002       FORMAT (//1X,'RESULTS FOR THIS RUN ARE: INVALID')

	write(6,*) 'qcd2=LG results are INVALID'
	call validexitbad

         ENDIF
C
      ELSEIF ( c.EQ.'E' ) THEN
         WRITE (6,*) 'ENTER DIR (0<->X,1<->Y,2<->Z)'
         READ (5,*) dir
         WRITE (6,*) dir
         WRITE (6,*) 'ENTER NUMBER OF SPATIAL LINK HITS TO AVERAGE OVER'
         READ (5,*) numbar
         WRITE (6,*) numbar
         DO i = 1 , 13
            DO j = 1 , 13
               DO k = 1 , 6
                  res(i,j,k) = 0.0
               ENDDO
            ENDDO
         ENDDO
         WRITE (6,*) 'ENTER R1MAX AND R2MAX :'
         READ (5,*) r1max , r2max
         WRITE (6,*) r1max , r2max
         WRITE (6,*) 'NOW ENTER THE QUADRANT TABLE OF TMINS ;'
         WRITE (6,*) r1max , 'ROWS &' , r2max , 'COLS (0 ==> OMIT)'
         READ (5,*) ((which(i,j),j=1,r2max),i=1,r1max)
         WRITE (6,*) 'MEASURING LOOPS IN DIRECTION ' , dir
C
C       CLOCK = SECNDS(0)
         CALL rotmea(dir,numbar,r1max,r2max)
C       WRITE(6,105) SECNDS(0)-CLOCK
         DO i = 1 , r1max
            DO j = 1 , r2max
               IF ( which(i,j).NE.0 ) THEN
                  tmx = min(12,which(i,j)+5)
                  DO t = which(i,j) , tmx - 1
                     WRITE (6,2) i + 1 , j + 1 , 
     &                           t , res(i,j,t-which(i,j)+1)
     &                           /real(numstg)
 2                     format ('R (' , i3 , ',' , i3 , ') X' , 
     &                           i3 , ':' , e19.9)
c                WRITE(8,*) 'R (',I+1,',',J+1,') X',T,':',
c     .            RES(I,J,T-WHICH(I,J)+1)/DOUBLEPRECISION(NUMSTG)
                  ENDDO
               ENDIF
            ENDDO
         ENDDO
C
         IF ( abs((res(1,1,1-which(1,1)+1)/real(numstg)-.434777)
     &        /.434777).GE.1E-2 ) ival = 1
         IF ( abs((res(1,1,2-which(1,1)+1)/real(numstg)-.299319)
     &        /.299319).GE.1E-2 ) ival = 1
         IF ( abs((res(1,1,3-which(1,1)+1)/real(numstg)-.205511)
     &        /.205511).GE.1E-2 ) ival = 1
         IF ( abs((res(1,1,4-which(1,1)+1)/real(numstg)-.136702)
     &        /.136702).GE.1E-2 ) ival = 1
         IF ( abs((res(1,1,5-which(1,1)+1)/real(numstg)-9.14717E-02)
     &        /9.14717E-02).GE.1E-2 ) ival = 1
C
      ELSEIF ( c.EQ.'C' ) THEN
         WRITE (6,*) '0 <-> COLD START, 1 <-> RESTART'
         READ (5,*) flg
         WRITE (6,*) flg
         CALL init(flg)
         WRITE (6,*) 'INITIALIZED SUCCESSFULLY.'
      ELSEIF ( c.EQ.'B' ) THEN
         CALL lattot()
         WRITE (6,*) 'BACKED UP SUCCESSFULLY'
      ELSEIF ( c.EQ.'F' ) THEN
         WRITE (6,*) 'INPUT # OF SWEEPS AND PROJECT INTERVAL'
         READ (5,*) sweeps , pinval
         WRITE (6,*) sweeps , pinval
         WRITE (6,*) 'RUNNING...'
C       CLOCK = SECNDS(0)
         CALL update(sweeps,pinval)
C       WRITE(6,104) SECNDS(0)-CLOCK
         WRITE (6,*) 'LATTICE UPDATED'
      ELSEIF ( c.EQ.'?' ) THEN
         WRITE (6,*) 'COMMANDS (AND THEIR PARAMETERS) ARE:'
         WRITE (6,*) '===================================='
         WRITE (6,*) 'C :             INITIALIZE (COLD OR READ)'
         WRITE (6,*) 'D :             UPDATES (SWEEPS, PINVAL)'
         WRITE (6,*) 'B :             BACK UP LATTICE (FILENAME)'
         WRITE (6,*) 'S :             SET PARAMETERS'
         WRITE (6,*) 'R :             RESET RANDOM NUMBER SEEDS'
         WRITE (6,*) 'D :             DUMP CURRENT SEED TO FILE'
         WRITE (6,*) 'P :             PROJECT MATRICES ONTO SU(3)'
         WRITE (6,*) 'M :             MEASURE ON AXIS LOOPS'
         WRITE (6,*) 'E :             MEASURE OFF AXIS LOOP'
         WRITE (6,*) 'A :             MEASURE BARYON LOOPS'
         WRITE (6,*) 'U :             MAKE U BARS'
         WRITE (6,*) 'Q :             QUIT'
      ELSEIF ( c.EQ.'Q' ) THEN
         WRITE (6,*) 'FINISHED.'
C

	call printvalidtime(1.0d0,"perfect.qcd2.F,1.12,99/03/12")

C
c      WRITE(8,*)
         OPEN (UNIT=1,FILE='SEED',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
         CLOSE (1,STATUS='delete')
         STOP
      ELSE
         WRITE (6,*) 'UNKNOWN COMMAND (? FOR HELP)'
      ENDIF
      GOTO 100
C
99003 FORMAT (/5X,'Time for LINKBR = ',F11.4,' seconds')
99004 FORMAT (/5X,'Time for MEASUR = ',F11.4,' seconds')
99005 FORMAT (/5X,'Time for QQQMEA = ',F11.4,' seconds')
99006 FORMAT (/5X,'Time for UPDATE  = ',F11.4,' seconds')
99007 FORMAT (/5X,'Time for ROTMEA = ',F11.4,' seconds')
99008 FORMAT (A1)
      END
**==cfill.f  
 
 
 
 
C
      SUBROUTINE cfill(a,b)
      DOUBLEPRECISION a(*) , b(*)
C
      DO i = 1 , 12
         a(i) = b(i)
      ENDDO
      a(13) = a(3)*a(11) - a(4)*a(12) - a(5)*a(9) + a(6)*a(10)
      a(14) = -a(4)*a(11) - a(3)*a(12) + a(6)*a(9) + a(5)*a(10)
      a(15) = a(5)*a(7) - a(6)*a(8) - a(1)*a(11) + a(2)*a(12)
      a(16) = -a(6)*a(7) - a(5)*a(8) + a(2)*a(11) + a(1)*a(12)
      a(17) = a(1)*a(9) - a(2)*a(10) - a(3)*a(7) + a(4)*a(8)
      a(18) = -a(2)*a(9) - a(1)*a(10) + a(4)*a(7) + a(3)*a(8)
C
      RETURN
      END
**==choos.f  
 
C     HEAT BATH - USES CABIBBO-MARINARI ALGORITHM INSTEAD OF
C                 METROPOLIS
C                 3 SUBGROUP VERSION
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE
C     AND DAVID WALKER, SUMMER 1987
C
      SUBROUTINE choos(sst,pu)
      DOUBLEPRECISION sst(*)
      DOUBLEPRECISION pu(*)
C
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER dims
      PARAMETER (dims=4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION r0 , r1 , r2 , r3 , r 
      DOUBLEPRECISION a0 , a1 , a2 , a3 , xx , efac 
      DOUBLEPRECISION rfac , norm
      DOUBLEPRECISION part(18) , plaq(18) , mat(18) 
      DOUBLEPRECISION dis , b0 , b1 , b2 , b3 
      DOUBLEPRECISION temp(18)
      INTEGER i , j
C
      DO i = 1 , 18
         part(i) = 0.0
         mat(i) = pu(i)
      ENDDO
C
      DO j = 1 , 6
         DO i = 1 , 18
            part(i) = part(i) + sst(i+18*(j-1))
         ENDDO
      ENDDO
C
      DO i = 1 , 3
         CALL mult(mat,part,plaq)
         IF ( i.EQ.1 ) THEN
            r0 = 0.5*(plaq(1)+plaq(9))
            r1 = 0.5*(plaq(2)-plaq(10))
            r2 = 0.5*(plaq(3)-plaq(7))
            r3 = 0.5*(plaq(4)+plaq(8))
         ELSEIF ( i.EQ.2 ) THEN
            r0 = 0.5*(plaq(9)+plaq(17))
            r1 = 0.5*(plaq(10)-plaq(18))
            r2 = 0.5*(plaq(11)-plaq(15))
            r3 = 0.5*(plaq(12)+plaq(16))
         ELSE
            r0 = 0.5*(plaq(1)+plaq(17))
            r1 = 0.5*(plaq(2)-plaq(18))
            r2 = 0.5*(plaq(5)-plaq(13))
            r3 = 0.5*(plaq(6)+plaq(14))
         ENDIF
         r = r0*r0 + r1*r1 + r2*r2 + r3*r3
         r = sqrt(r)
         r0 = r0/r
         r1 = r1/r
         r2 = r2/r
         r3 = r3/r
         xx = beta*r*2.0/3.0
         efac = exp(-2.0*xx)
         rfac = 1.0 - efac
C
   50    CONTINUE
         a0 = 1.0 + log(efac+rfac*pranf(0))/xx
         xtemp = pranf(0)
         IF ( (xtemp*xtemp).GT.(1.0-a0*a0) ) GOTO 50
         norm = sqrt(1.0-a0*a0)
C
  100    CONTINUE
         a1 = 2.0*pranf(0) - 1.0
         a2 = 2.0*pranf(0) - 1.0
         a3 = 2.0*pranf(0) - 1.0
         dis = a1*a1 + a2*a2 + a3*a3
         IF ( dis.GT.1.0 ) GOTO 100
C
         norm = norm/sqrt(dis)
         a1 = a1*norm
         a2 = a2*norm
         a3 = a3*norm
         b0 = a0*r0 + a1*r1 + a2*r2 + a3*r3
         b1 = -a0*r1 + a1*r0 + a2*r3 - a3*r2
         b2 = a2*r0 + a3*r1 - a0*r2 - a1*r3
         b3 = -a2*r1 + a3*r0 - a0*r3 + a1*r2
         IF ( i.EQ.1 ) THEN
            temp(1) = b0*mat(1) - b1*mat(2) - b2*mat(3) - b3*mat(4)
            temp(2) = b1*mat(1) + b0*mat(2) + b3*mat(3) - b2*mat(4)
            temp(3) = b2*mat(1) - b3*mat(2) + b0*mat(3) + b1*mat(4)
            temp(4) = b3*mat(1) + b2*mat(2) - b1*mat(3) + b0*mat(4)
            temp(7) = b0*mat(7) - b1*mat(8) - b2*mat(9) - b3*mat(10)
            temp(8) = b1*mat(7) + b0*mat(8) + b3*mat(9) - b2*mat(10)
            temp(9) = b2*mat(7) - b3*mat(8) + b0*mat(9) + b1*mat(10)
            temp(10) = b3*mat(7) + b2*mat(8) - b1*mat(9) + b0*mat(10)
            temp(13) = b0*mat(13) - b1*mat(14) - b2*mat(15) - b3*mat(16)
            temp(14) = b1*mat(13) + b0*mat(14) + b3*mat(15) - b2*mat(16)
            temp(15) = b2*mat(13) - b3*mat(14) + b0*mat(15) + b1*mat(16)
            temp(16) = b3*mat(13) + b2*mat(14) - b1*mat(15) + b0*mat(16)
            mat(1) = temp(1)
            mat(2) = temp(2)
            mat(3) = temp(3)
            mat(4) = temp(4)
            mat(7) = temp(7)
            mat(8) = temp(8)
            mat(9) = temp(9)
            mat(10) = temp(10)
            mat(13) = temp(13)
            mat(14) = temp(14)
            mat(15) = temp(15)
            mat(16) = temp(16)
         ELSEIF ( i.EQ.2 ) THEN
            temp(3) = b0*mat(3) - b1*mat(4) - b2*mat(5) - b3*mat(6)
            temp(4) = b1*mat(3) + b0*mat(4) + b3*mat(5) - b2*mat(6)
            temp(5) = b2*mat(3) - b3*mat(4) + b0*mat(5) + b1*mat(6)
            temp(6) = b3*mat(3) + b2*mat(4) - b1*mat(5) + b0*mat(6)
            temp(9) = b0*mat(9) - b1*mat(10) - b2*mat(11) - b3*mat(12)
            temp(10) = b1*mat(9) + b0*mat(10) + b3*mat(11) - b2*mat(12)
            temp(11) = b2*mat(9) - b3*mat(10) + b0*mat(11) + b1*mat(12)
            temp(12) = b3*mat(9) + b2*mat(10) - b1*mat(11) + b0*mat(12)
            temp(15) = b0*mat(15) - b1*mat(16) - b2*mat(17) - b3*mat(18)
            temp(16) = b1*mat(15) + b0*mat(16) + b3*mat(17) - b2*mat(18)
            temp(17) = b2*mat(15) - b3*mat(16) + b0*mat(17) + b1*mat(18)
            temp(18) = b3*mat(15) + b2*mat(16) - b1*mat(17) + b0*mat(18)
            mat(3) = temp(3)
            mat(4) = temp(4)
            mat(5) = temp(5)
            mat(6) = temp(6)
            mat(9) = temp(9)
            mat(10) = temp(10)
            mat(11) = temp(11)
            mat(12) = temp(12)
            mat(15) = temp(15)
            mat(16) = temp(16)
            mat(17) = temp(17)
            mat(18) = temp(18)
         ELSE
            temp(1) = b0*mat(1) - b1*mat(2) - b2*mat(5) - b3*mat(6)
            temp(2) = b1*mat(1) + b0*mat(2) + b3*mat(5) - b2*mat(6)
            temp(5) = b2*mat(1) - b3*mat(2) + b0*mat(5) + b1*mat(6)
            temp(6) = b3*mat(1) + b2*mat(2) - b1*mat(5) + b0*mat(6)
            temp(7) = b0*mat(7) - b1*mat(8) - b2*mat(11) - b3*mat(12)
            temp(8) = b1*mat(7) + b0*mat(8) + b3*mat(11) - b2*mat(12)
            temp(11) = b2*mat(7) - b3*mat(8) + b0*mat(11) + b1*mat(12)
            temp(12) = b3*mat(7) + b2*mat(8) - b1*mat(11) + b0*mat(12)
            mat(1) = temp(1)
            mat(2) = temp(2)
            mat(5) = temp(5)
            mat(6) = temp(6)
            mat(7) = temp(7)
            mat(8) = temp(8)
            mat(11) = temp(11)
            mat(12) = temp(12)
         ENDIF
      ENDDO
C
      CALL cfill(pu,mat)
C
      RETURN
      END
**==cpymat.f  
 
      SUBROUTINE cpymat(outmat,inmat,n)
      DOUBLEPRECISION inmat(*) , outmat(*)
      INTEGER n
C
      DO i = 1 , n
         outmat(i) = inmat(i)
      ENDDO
C
      RETURN
      END
**==icpymt.f  
 
      SUBROUTINE icpymt(outmat,inmat,n)
      INTEGER inmat(*) , outmat(*) , n
 
C     ADDED BY LPOINTER 4-3-89
 
      DO i = 1 , n
         outmat(i) = inmat(i)
      ENDDO
 
      RETURN
      END
**==dmpsed.f  
 
      SUBROUTINE dmpsed()
C
      OPEN (UNIT=1,FILE='SEED',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
      REWIND (1)
      ino = igetsd(0)
      WRITE (1,99001) ino
99001 FORMAT (I12)
      WRITE (6,*) 'O.K. ...... DONE IT !'
      CLOSE (1)
      END
**==getpos.f  
C
CCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE DMPSED CCCCCCCCCCCCCCCCCCCC
C
 
 
 
 
      SUBROUTINE getpos(latt,coord,index)
      INTEGER coord(4) , latt(4) , index
C
      j = index
      DO i = 1 , 4
         coord(i) = mod(j,latt(i))
         j = j/latt(i)
      ENDDO
C
      RETURN
      END
**==hit.f  
 
      SUBROUTINE hit(coord,origin,dir)
      INTEGER coord(4) , dir , origin
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mbar=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION  beta , legs(mbar)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      INTEGER othend , site , cnt
      INTEGER pu
      DOUBLEPRECISION st(18,6)
      INTEGER ptemp(24)
C
      INTEGER pstr(24,4)
      SAVE pstr
C
      DATA pstr/2 , 5 , 6 , 14 , 6 , 5 , 2 , 14 , 3 , 5 , 7 , 14 , 7 , 
     &     5 , 3 , 14 , 4 , 5 , 8 , 14 , 8 , 5 , 4 , 14 , 1 , 6 , 5 , 
     &     14 , 5 , 6 , 1 , 14 , 3 , 6 , 7 , 14 , 7 , 6 , 3 , 14 , 4 , 
     &     6 , 8 , 14 , 8 , 6 , 4 , 14 , 1 , 7 , 5 , 14 , 5 , 7 , 1 , 
     &     14 , 2 , 7 , 6 , 14 , 6 , 7 , 2 , 14 , 4 , 7 , 8 , 14 , 8 , 
     &     7 , 4 , 14 , 1 , 8 , 5 , 14 , 5 , 8 , 1 , 14 , 2 , 8 , 6 , 
     &     14 , 6 , 8 , 2 , 14 , 3 , 8 , 7 , 14 , 7 , 8 , 3 , 14/
C
C      CALL CPYMAT(PTEMP,PSTR(1,DIR+1),24) LPOINTER 4-3-89
      CALL icpymt(ptemp,pstr(1,dir+1),24)
C
      IF ( coord(dir+1).EQ.latt1(dir+1) ) THEN
         othend = 0
      ELSE
         othend = coord(dir+1) + 1
      ENDIF
C
      coord(dir+1) = othend
      site = coord(1)*mov(1) + coord(2)*mov(2) + coord(3)*mov(3)
     &        + coord(4)*mov(4)
      cnt = 1
      DO j = 1 , 6
         coord(dir+1) = othend
         CALL syslop(coord,site,ptemp(cnt),st(1,j),nchar)
         cnt = nchar + cnt
      ENDDO
C
      pu = origin + rot(dir+1)
      CALL choos(st(1,1),u1(pu+1))
C
      RETURN
      END
**==init.f  
 
 
 
      SUBROUTINE init(flag)
C
      INTEGER numstg
      COMMON /comfil/ numstg
C
      INTEGER dims
      PARAMETER (dims=4)
      INTEGER umax
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3,umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      DOUBLEPRECISION facfi , facif
      COMMON /comfac/ facfi , facif
C
      INTEGER flag , bignum , movx , buff , rem , bufnum
      INTEGER pu
      CHARACTER*5 bname
C
      movx = 72
      buff = 256
      bignum = 4*latt(1)*latt(2)*latt(3)*latt(4)
C
      IF ( flag.EQ.1 ) THEN
         WRITE (6,*) 'BACKUP FILE NAME ? : '
         READ (5,99001) bname
         OPEN (UNIT=1,FILE=bname,ACCESS='SEQUENTIAL',STATUS='UNKNOWN',
     &         FORM='UNFORMATTED')
         REWIND (1)
         bignum = movx*numstg
         bufnum = bignum/buff
         rem = mod(bignum,buff)
         WRITE (6,*) 'BACKUP FILE NAME = ' , bname , 'BUFNUM = ' , 
     &               bufnum , 'REM = ' , rem
         pu = 1
         DO i = 0 , bufnum - 1
            READ (1) (u1(j),j=pu,pu+buff)
            pu = pu + buff
            WRITE (6,*) 'B'
         ENDDO
         WRITE (6,*) 'E'
         CLOSE (1)
      ELSE
         pu = 1
         DO i = 0 , bignum - 1
            u1(pu) = 1.0
            u1(pu+1) = 0.0
            u1(pu+2) = 0.0
            u1(pu+3) = 0.0
            u1(pu+4) = 0.0
            u1(pu+5) = 0.0
            u1(pu+6) = 0.0
            u1(pu+7) = 0.0
            u1(pu+8) = 1.0
            u1(pu+9) = 0.0
            u1(pu+10) = 0.0
            u1(pu+11) = 0.0
            u1(pu+12) = 0.0
            u1(pu+13) = 0.0
            u1(pu+14) = 0.0
            u1(pu+15) = 0.0
            u1(pu+16) = 1.0
            u1(pu+17) = 0.0
            pu = pu + rot(2)
         ENDDO
      ENDIF
C
      RETURN
99001 FORMAT (A)
      END
**==lattot.f  
 
C
C     THE ROUTINES USED TO INITIALIZE AND DUMP THE LATTICE.
C
C     CONVERTED FORM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, SUMMER 1987
C
      SUBROUTINE lattot()
C
      INTEGER umax
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (umax=72*l1*l2*l3*l4)
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      INTEGER numstg
      COMMON /comfil/ numstg
C
      INTEGER bignum
      INTEGER bufnum , rem
      INTEGER pu
      INTEGER buff , movx
      CHARACTER*5 bname
C
      WRITE (6,*) 'ENTER 5-LETTER BACKUP FILE NAME ? : '
      READ (5,99002) bname
      WRITE (6,*) bname
      OPEN (UNIT=1,FILE=bname,ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
C     OPEN(UNIT=1,FILE = BNAME,ACCESS = 'SEQUENTIAL',STATUS='UNKNOWN',
C    .FORM='UNFORMATTED')
      REWIND (1)
      movx = 72
      buff = 256
      bignum = movx*numstg
      bufnum = bignum/buff
      rem = mod(bignum,buff)
      WRITE (6,*) 'BACKUP FILE NAME = ' , bname , 'BUFNUM = ' , bufnum , 
     &            'REM = '
      pu = 1
      DO i = 0 , bufnum - 1
C       WRITE (1) (U1(J),J=PU,PU+BUFF)
         WRITE (1,99001) (u1(j),j=pu,pu+buff)
99001    FORMAT (9F8.3)
         pu = pu + buff
         WRITE (6,*) 'B'
      ENDDO
      WRITE (6,*) 'E'
      CLOSE (1)
      RETURN
99002 FORMAT (A)
      END
**==legmak.f  
 
C
C     THIS ROUTINE MAKES BARRED LEGS ONLY IN THE ARRAY LEGS
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE legmak(t,tt)
      INTEGER tt , t
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      INTEGER strb(26)
C
      INTEGER site , coord(4) , plegs
C
      DO i = 1 , t
         strb(i) = 18
      ENDDO
C
      DO i = t + 1 , 2*t + 1
         strb(i) = 30
      ENDDO
C
      strb(2*t+1) = 14
C
      coord(4) = tt
      site = tt*mov(4)
      plegs = 1
      DO i = 0 , latt(3) - 1
         coord(3) = i
         DO j = 0 , latt(2) - 1
            coord(2) = j
            DO k = 0 , latt(1) - 1
               coord(1) = k
               CALL syslop(coord,site,strb,legs(plegs),nchar)
               plegs = plegs + 18
               site = site + mov(1)
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
**==linkbr.f  
 
C     LINKBR MULTI-HITS AND AVERAGES U MATRICES IN ANY DIRECTION
C     AND PUTS THE RESULTS IN THE ARRAY POINTED TO BY "WHERE".
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE linkbr(n,dir,where)
      DOUBLEPRECISION where(*)
      INTEGER n , dir
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      INTEGER coord(4) , othend , site , endsit , pt , cnt
      INTEGER pu
      INTEGER pstr(24)
      DOUBLEPRECISION avmat(18) , mat2(18) , st(18,6)
C
      INTEGER env(24,4)
      SAVE env
C
      DATA env/2 , 5 , 6 , 14 , 6 , 5 , 2 , 14 , 3 , 5 , 7 , 14 , 7 , 
     &     5 , 3 , 14 , 4 , 5 , 8 , 14 , 8 , 5 , 4 , 14 , 1 , 6 , 5 , 
     &     14 , 5 , 6 , 1 , 14 , 3 , 6 , 7 , 14 , 7 , 6 , 3 , 14 , 4 , 
     &     6 , 8 , 14 , 8 , 6 , 4 , 14 , 1 , 7 , 5 , 14 , 5 , 7 , 1 , 
     &     14 , 2 , 7 , 6 , 14 , 6 , 7 , 2 , 14 , 4 , 7 , 8 , 14 , 8 , 
     &     7 , 4 , 14 , 1 , 8 , 5 , 14 , 5 , 8 , 1 , 14 , 2 , 8 , 6 , 
     &     14 , 6 , 8 , 2 , 14 , 3 , 8 , 7 , 14 , 7 , 8 , 3 , 14/
C
      pt = 1
      DO icord4 = 0 , latt(4) - 1
         coord(4) = icord4
         DO icord3 = 0 , latt(3) - 1
            coord(3) = icord3
            DO icord2 = 0 , latt(2) - 1
               coord(2) = icord2
               DO icord1 = 0 , latt(1) - 1
                  coord(1) = icord1
                  cnt = 1
                  site = icord1*mov(1) + icord2*mov(2) + icord3*mov(3)
     &                    + icord4*mov(4)
                  pu = site + rot(dir+1)
                  CALL cpymat(avmat,u1(pu+1),18)
                  IF ( n.GT.1 ) THEN
                     DO j = 1 , 18
                        mat2(j) = u1(pu+j)
                     ENDDO
                     IF ( coord(dir+1).EQ.latt1(dir+1) ) THEN
                        othend = 0
                        endsit = site - latt1(dir+1)*mov(dir+1)
                     ELSE
                        othend = coord(dir+1) + 1
                        endsit = site + mov(dir+1)
                     ENDIF
C                CALL CPYMAT(PSTR,ENV(1,DIR+1),24) LPOINTER 4-3-89
                     CALL icpymt(pstr,env(1,dir+1),24)
C
                     DO j = 1 , 6
                        coord(dir+1) = othend
                        CALL syslop(coord,endsit,pstr(cnt),st(1,j),
     &                              nchar)
                        cnt = cnt + nchar
                     ENDDO
                     DO j = 1 , n - 1
                        CALL choos(st(1,1),mat2)
                        DO i = 1 , 18
                           avmat(i) = avmat(i) + mat2(i)
                        ENDDO
                     ENDDO
                     DO i = 1 , 18
                        avmat(i) = avmat(i)/n
                     ENDDO
                  ENDIF
                  DO k = 1 , 18
                     where(pt+k-1) = avmat(k)
                  ENDDO
                  pt = pt + 18
               ENDDO
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
**==lpstrm.f  
 
 
 
 
C
C     S. W. OTTO AND J. W. FLOWER, SUMMER 1985.
C
C     CONSTRUCTS THE STRINGS FOR MEASURING PLANAR RxT WILSON LOOPS.
C     IF R=1 THEN INSTEAD OF FETCHING A BARRED LEG FOR THE T DIRECTION W
C     CONSTRUCT IT EXPLICITLY USING THE REGULAR LINKS.
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE lpstrm(pstrt,pstru,pstrb,r,t,dir)
      INTEGER pstrt(60) , pstru(60) , pstrb(60)
      INTEGER r , t , dir
C
      INTEGER rr , tt , cnt , temp
C
C  FIRST GET THE TOP OF THE LOOP BARRING AS MANY LINKS AS POSSIBLE
C
      cnt = 1
      DO tt = 0 , t - 1
         pstrt(cnt) = 26
         cnt = cnt + 1
      ENDDO
C
      temp = 23 + dir
      DO rr = 0 , r - 1
         pstrt(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      temp = 5 + dir
      pstrt(cnt) = temp
      cnt = cnt + 1
C
      temp = 19 + dir
      DO rr = 1 , r - 2
         pstrt(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      IF ( r.GE.2 ) THEN
         temp = 5 + dir
         pstrt(cnt) = temp
         cnt = cnt + 1
      ENDIF
C
      DO tt = 0 , t - 1
         pstrt(cnt) = 30
         cnt = cnt + 1
      ENDDO
      pstrt(cnt) = 14
C
C  NOW CONSTRUCT THE TIMELIKE LEG. IF R=1 MAKE EXPLICITLY
C
      cnt = 1
      temp = 23 + dir
      DO rr = 0 , r - 1
         pstru(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      IF ( r.GT.1 ) THEN
         pstru(cnt) = 36
         cnt = cnt + 1
      ELSE
         DO tt = 0 , t - 1
            pstru(cnt) = 4
            cnt = cnt + 1
         ENDDO
         DO tt = 0 , t - 1
            pstru(cnt) = 30
            cnt = cnt + 1
         ENDDO
      ENDIF
C
      temp = 27 + dir
      DO rr = 0 , r - 1
         pstru(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      pstru(cnt) = 14
      cnt = cnt + 1
C
C  FINALLY MAKE THE BOTTOM OF THE LOOP
C
      cnt = 1
      temp = 1 + dir
      pstrb(cnt) = temp
C
      cnt = cnt + 1
      temp = 15 + dir
      DO rr = 1 , r - 2
         pstrb(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      IF ( r.GE.2 ) THEN
         temp = 1 + dir
         pstrb(cnt) = temp
         cnt = cnt + 1
      ENDIF
C
      temp = 27 + dir
      DO rr = 0 , r - 1
         pstrb(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      pstrb(cnt) = 14
C
      RETURN
      END
**==make.f  
 
 
 
 
C
C     STEVE OTTO (AND JON FLOWER) JULY 1984
C
C     SETS UP INITIAL MOVEMENT TABLES.
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, SUMMER 1987
C
      SUBROUTINE make
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DO i = 1 , 4
         rot(i) = 18*(i-1)
      ENDDO
C
      mov(1) = 72
      mov(2) = mov(1)*latt(1)
      mov(3) = mov(2)*latt(2)
      mov(4) = mov(3)*latt(3)
C
      DO i = 1 , 4
         latt1(i) = latt(i) - 1
      ENDDO
C
      RETURN
      END
**==measur.f  
 
C
C     JON FLOWER, STEVE OTTO   MARCH 1985.  MOD: JUNE 85
C
C     MEASURES WILSON LOOPS UP TO RxT BARRING EVERYTHING IN SIGHT....
C     ALSO MAKES LEGS.
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE measur(dir,rmin,rmax,numhit)
      INTEGER dir , rmin , rmax , numhit
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      PARAMETER (mbar=18*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION ubar(mbar) , sbar(mbar)
      COMMON /comubr/ ubar , sbar
C
      INTEGER r , nchar , coord(4) , tflag , site , t , plegs
      INTEGER strt(60) , stru(60) , strb(60)
      DOUBLEPRECISION up(18) , down(18) , top(18) 
      DOUBLEPRECISION  bottom(18) , stemp(18,2)
C
      DO r = 0 , 12
         DO t = 0 , 12
            wloop(r+1,t+1) = 0.0
         ENDDO
      ENDDO
C
      CALL linkbr(numhit,dir,sbar)
C
C     ONE TIME SLICE AT A TIME
C
      DO icord4 = 0 , latt(4) - 1
         coord(4) = icord4
         DO t = 1 , 12
            tflag = 0
            DO r = rmin , rmax
               IF ( which(r+1,t+1).EQ.1 ) tflag = 1
            ENDDO
C
C    TFLAG=1 MEANS WE NEED LEGS OF THIS T
C
            IF ( tflag.EQ.1 ) CALL legmak(t,icord4)
C
            DO r = rmin , rmax
               IF ( (r.LE.t) .AND. (which(r+1,t+1).EQ.1) ) THEN
                  CALL lpstrm(strt,stru,strb,r,t,dir)
                  site = icord4*mov(4)
                  plegs = 1
                  DO icord3 = 0 , latt(3) - 1
                     coord(3) = icord3
                     DO icord2 = 0 , latt(2) - 1
                        coord(2) = icord2
                        DO icord1 = 0 , latt(1) - 1
                           coord(1) = icord1
C
C                    FETCH AND DAGGER THE DOWN LEG
C
                           CALL udag(legs(plegs),down)
C
C                    FIND THE TOP
C
                           CALL syslop(coord,site,strt,top,nchar)
C
C                    FIND THE UP
C
                           CALL syslop(coord,site,stru,up,nchar)
C
C                    FIND THE BOTTOM
C
                           CALL syslop(coord,site,strb,bottom,nchar)
C
C                    PUT THEM ALL TOGETHER
C
                           CALL mult(bottom,up,stemp(1,1))
                           CALL mult(stemp(1,1),top,stemp(1,2))
                           CALL mult(stemp(1,2),down,stemp(1,1))
                           wloop(r+1,t+1) = wloop(r+1,t+1)
     &                            + (stemp(1,1)+stemp(9,1)+stemp(17,1))
     &                           /3.0
                           plegs = plegs + 18
                           site = site + mov(1)
                        ENDDO
                     ENDDO
                  ENDDO
               ENDIF
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
**==mult.f  
 
      SUBROUTINE mult(a,b,c)
      DOUBLEPRECISION a(18) , b(18) , c(18)
C
      c(1) = a(1)*b(1) - a(2)*b(2) + a(7)*b(3) - a(8)*b(4) + a(13)*b(5)
     &        - a(14)*b(6)
      c(2) = a(1)*b(2) + a(2)*b(1) + a(7)*b(4) + a(8)*b(3) + a(13)*b(6)
     &        + a(14)*b(5)
      c(3) = a(3)*b(1) - a(4)*b(2) + a(9)*b(3) - a(10)*b(4) + a(15)*b(5)
     &        - a(16)*b(6)
      c(4) = a(3)*b(2) + a(4)*b(1) + a(9)*b(4) + a(10)*b(3) + a(15)*b(6)
     &        + a(16)*b(5)
      c(5) = a(5)*b(1) - a(6)*b(2) + a(11)*b(3) - a(12)*b(4) + a(17)
     &       *b(5) - a(18)*b(6)
      c(6) = a(5)*b(2) + a(6)*b(1) + a(11)*b(4) + a(12)*b(3) + a(17)
     &       *b(6) + a(18)*b(5)
      c(7) = a(1)*b(7) - a(2)*b(8) + a(7)*b(9) - a(8)*b(10) + a(13)
     &       *b(11) - a(14)*b(12)
      c(8) = a(1)*b(8) + a(2)*b(7) + a(7)*b(10) + a(8)*b(9) + a(13)
     &       *b(12) + a(14)*b(11)
      c(9) = a(3)*b(7) - a(4)*b(8) + a(9)*b(9) - a(10)*b(10) + a(15)
     &       *b(11) - a(16)*b(12)
      c(10) = a(3)*b(8) + a(4)*b(7) + a(9)*b(10) + a(10)*b(9) + a(15)
     &        *b(12) + a(16)*b(11)
      c(11) = a(5)*b(7) - a(6)*b(8) + a(11)*b(9) - a(12)*b(10) + a(17)
     &        *b(11) - a(18)*b(12)
      c(12) = a(5)*b(8) + a(6)*b(7) + a(11)*b(10) + a(12)*b(9) + a(17)
     &        *b(12) + a(18)*b(11)
      c(13) = a(1)*b(13) - a(2)*b(14) + a(7)*b(15) - a(8)*b(16) + a(13)
     &        *b(17) - a(14)*b(18)
      c(14) = a(1)*b(14) + a(2)*b(13) + a(7)*b(16) + a(8)*b(15) + a(13)
     &        *b(18) + a(14)*b(17)
      c(15) = a(3)*b(13) - a(4)*b(14) + a(9)*b(15) - a(10)*b(16) + a(15)
     &        *b(17) - a(16)*b(18)
      c(16) = a(3)*b(14) + a(4)*b(13) + a(9)*b(16) + a(10)*b(15) + a(15)
     &        *b(18) + a(16)*b(17)
      c(17) = a(5)*b(13) - a(6)*b(14) + a(11)*b(15) - a(12)*b(16)
     &         + a(17)*b(17) - a(18)*b(18)
      c(18) = a(5)*b(14) + a(6)*b(13) + a(11)*b(16) + a(12)*b(15)
     &         + a(17)*b(18) + a(18)*b(17)
C
      RETURN
      END
**==prnse2.f  
C
      SUBROUTINE prnse2(iseed)
      PARAMETER (mult=1103515245,iadd=12345)
      INTEGER aaa , bbb , randx
      COMMON /compra/ aaa , bbb , randx
      SAVE /compra/
      aaa = mult
      bbb = iadd
      randx = iseed
      RETURN
      END
**==observ.f  
C
C
C    J. W. FLOWER   AUG 1985
C
C     CONSTRUCTS THE THREE QUARK OBSERVABLE :
C
C                            (1)  (2)  (3)
C       EPSILON   EPSILON    U    U    U
C              IJK       PQR  IP   JQ   KR
C
C    THIS CODE IS CERTAINLY NON-OPTIMAL BUT SACRIFICES
C            SPEED FOR CODING SIMPLICITY
C
C    CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C    AND DAVID WALKER, 1987
C
      SUBROUTINE observ(tot,u1,u2,u3)
      DOUBLEPRECISION tot(2) , u1(2,9) , u2(2,9) 
      DOUBLEPRECISION u3(2,9)
C
      INTEGER p , q , r , epsilo(3,3,3)
      DOUBLEPRECISION fac
C
      DO i = 1 , 3
         DO j = 1 , 3
            DO k = 1 , 3
               epsilo(i,j,k) = 0
            ENDDO
         ENDDO
      ENDDO
C
      epsilo(1,2,3) = 1
      epsilo(2,3,1) = 1
      epsilo(3,1,2) = 1
      epsilo(2,1,3) = -1
      epsilo(1,3,2) = -1
      epsilo(3,2,1) = -1
C
      DO i = 0 , 2
         DO p = 0 , 2
            DO j = 0 , 2
               DO q = 0 , 2
                  DO k = 0 , 2
                     IF ( epsilo(i+1,j+1,k+1).NE.0 ) THEN
                        DO r = 0 , 2
                           IF ( epsilo(p+1,q+1,r+1).NE.0 ) THEN
                              fac = epsilo(i+1,j+1,k+1)
     &                              *epsilo(p+1,q+1,r+1)
                              tot(1) = tot(1) + fac*u1(1,3*i+p+1)
     &                                 *u2(1,3*j+q+1)*u3(1,3*k+r+1)
                              tot(1) = tot(1) - fac*u1(2,3*i+p+1)
     &                                 *u2(2,3*j+q+1)*u3(1,3*k+r+1)
                              tot(1) = tot(1) - fac*u1(1,3*i+p+1)
     &                                 *u2(2,3*j+q+1)*u3(2,3*k+r+1)
                              tot(1) = tot(1) - fac*u1(2,3*i+p+1)
     &                                 *u2(1,3*j+q+1)*u3(2,3*k+r+1)
                              tot(2) = tot(2) + fac*u1(1,3*i+p+1)
     &                                 *u2(1,3*j+q+1)*u3(2,3*k+r+1)
                              tot(2) = tot(2) + fac*u1(1,3*i+p+1)
     &                                 *u2(2,3*j+q+1)*u3(1,3*k+r+1)
                              tot(2) = tot(2) + fac*u1(2,3*i+p+1)
     &                                 *u2(1,3*j+q+1)*u3(1,3*k+r+1)
                              tot(2) = tot(2) - fac*u1(2,3*i+p+1)
     &                                 *u2(2,3*j+q+1)*u3(2,3*k+r+1)
                           ENDIF
                        ENDDO
                     ENDIF
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
**==prnset.f  
C
      SUBROUTINE prnset(iseed)
      COMMON /ran   / irand0 , irand1 , iab0 , iab1 , icb0 , icb1
      SAVE /ran/
C
      PARAMETER (mult0=20077,mult1=16838,iadd0=12345,iadd1=0)
      PARAMETER (ipow16=2**16)
C
      iseed0 = mod(iseed,ipow16)
      iseed1 = iseed/ipow16
      iab0 = mult0
      iab1 = mult1
      icb0 = iadd0
      icb1 = iadd1
      CALL lmult(j0,j1,iseed0,iseed1,iab0,iab1)
      CALL ladd(irand0,irand1,j0,j1,icb0,icb1)
      RETURN
      END
**==pranf.f  
C
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE PRNSET CCCCCCCCCCCCCCCCCCCCCCC C
C
      FUNCTION pranf(idummy)
      COMMON /ran   / irand0 , irand1 , iab0 , iab1 , icb0 , icb1
      SAVE /ran/
C
      DOUBLE PRECISION randb , divfac , pow16
      PARAMETER (divfac=2.147483648D9,pow16=6.5536D4)
C
      randb = (dble(irand0)+pow16*dble(irand1))/divfac
      pranf = real(randb)
      CALL lmult(j0,j1,irand0,irand1,iab0,iab1)
      CALL ladd(irand0,irand1,j0,j1,icb0,icb1)
C
      RETURN
      END
**==ladd.f  
C
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE PRANF CCCCCCCCCCCCCCCCCCCCCCC C
C
      SUBROUTINE ladd(i0,i1,j0,j1,k0,k1)
C
      PARAMETER (ipow16=2**16,ipow15=2**15)
C
      i0 = mod(k0+j0,ipow16)
      i1 = mod((k0+j0)/ipow16+k1+j1,ipow15)
C
      RETURN
      END
**==lmult.f  
C
CCCCCCCCCCCCCCCCCCCCC END OF ROUTINE LADD CCCCCCCCCCCCCCCCCCCCCCC C
C
      SUBROUTINE lmult(i0,i1,j0,j1,k0,k1)
C
      PARAMETER (ipow30=2**30,ipow16=2**16,ipow15=2**15)
C
      k = k0*j0
      IF ( k.LT.0 ) k = (k+ipow30) + ipow30
      i0 = mod(k,ipow16)
      l = k0*j1 + k1*j0
      IF ( l.LT.0 ) l = (l+ipow30) + ipow30
      k = k/ipow16 + l
      IF ( k.LT.0 ) k = (k+ipow30) + ipow30
      i1 = mod(k,ipow15)
C
      RETURN
      END
**==igetsd.f  
C
CCCCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE LMULT CCCCCCCCCCCCCCCCCCCCCCC
C
      FUNCTION igetsd(idummy)
C
      COMMON /ran   / irand0 , irand1 , iab0 , iab1 , icb0 , icb1
      SAVE /ran/
C
      igetsd = irand0 + irand1*2**16
C
      RETURN
      END
**==projec.f  
C
CCCCCCCCCCCCCCCCCCCCCCCCCCC  END OF ROUTINE IGETSD  CCCCCCCCCCCCCCCCCCCC
C
 
C
C     STEVE OTTO  AND JON FLOWER, JULY 1984.
C
C     DUE TO THE STRANGE SHAPE OF THE SU(3) GROUP MANIFOLD (WHATEVER IT
C     THE LINK MATRICES TEND TO DRIFT OFF AND MUST BE PROJECTED BACK
C     PERIODICALLY. THE SITUATION IS OBVIOUSLY NOT HELPED BY THE LIMITED
C     ACCURACY OF THE FIXED POINT FORMAT.                             */
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE projec
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      DOUBLEPRECISION x(12) , z(6) , fac
      INTEGER site
      INTEGER pu
C
      numlin = 4*latt(1)*latt(2)*latt(3)*latt(4)
C
      DO i = 0 , numlin - 1
         site = rot(2)*i
         pu = site
         DO j = 1 , 12
            x(j) = u1(pu+j)
         ENDDO
         z(1) = x(3)*x(11) - x(4)*x(12) - x(5)*x(9) + x(6)*x(10)
         z(2) = -x(4)*x(11) - x(3)*x(12) + x(6)*x(9) + x(5)*x(10)
         z(3) = x(5)*x(7) - x(6)*x(8) - x(1)*x(11) + x(2)*x(12)
         z(4) = -x(6)*x(7) - x(5)*x(8) + x(2)*x(11) + x(1)*x(12)
         z(5) = x(1)*x(9) - x(2)*x(10) - x(3)*x(7) + x(4)*x(8)
         z(6) = -x(2)*x(9) - x(1)*x(10) + x(4)*x(7) + x(3)*x(8)
         fac = 0.0
         DO j = 1 , 6
            fac = fac + x(j)*x(j)
         ENDDO
         fac = 1.0/sqrt(fac)
         DO j = 1 , 6
            x(j) = x(j)*fac
         ENDDO
         x(7) = z(3)*x(5) - z(4)*x(6) - z(5)*x(3) + z(6)*x(4)
         x(8) = -z(3)*x(6) - z(4)*x(5) + z(5)*x(4) + z(6)*x(3)
         x(9) = z(5)*x(1) - z(6)*x(2) - z(1)*x(5) + z(2)*x(6)
         x(10) = -z(5)*x(2) - z(6)*x(1) + z(1)*x(6) + z(2)*x(5)
         x(11) = z(1)*x(3) - z(2)*x(4) - z(3)*x(1) + z(4)*x(2)
         x(12) = -z(1)*x(4) - z(2)*x(3) + z(3)*x(2) + z(4)*x(1)
         fac = 0.0
         DO j = 1 , 6
            fac = fac + x(6+j)*x(6+j)
         ENDDO
         fac = 1.0/sqrt(fac)
         DO j = 1 , 6
            x(6+j) = x(6+j)*fac
         ENDDO
         CALL cfill(u1(pu+1),x)
      ENDDO
C
      RETURN
      END
**==qqqlpc.f  
 
C
C     J. W. FLOWER & STEVE OTTO, SEPTEMBER 1984.
C
C     USED TO ACTUALLY CALCULATE THE THREE QUARK OBSERVABLES.       */
C
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE qqqlpc(r1,r2,r3,t,lpstr,tot)
      INTEGER r1 , r2 , r3 , t
      INTEGER lpstr(512)
      DOUBLEPRECISION tot(2)
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      INTEGER coord(4) , odir , site , cnt
      INTEGER pstr(512)
      DOUBLEPRECISION u1(18) , u2(18) , u3(18) 
      DOUBLEPRECISION bits(18,3) , temp(18)
C
      site = t*mov(4)
C
      coord(4) = t
      DO icord3 = 0 , latt(3) - 1
         coord(3) = icord3
         DO icord2 = 0 , latt(2) - 1
            coord(2) = icord2
            DO icord1 = 0 , latt(1) - 1
               coord(1) = icord1
C
C     FIRST GET THE MATRIX CORRESPONDING TO DIRECTION R1.
C
C            CALL CPYMAT(PSTR,LPSTR,512) LPOINTER 4-3-89
               CALL icpymt(pstr,lpstr,512)
               cnt = 1
               DO i = 1 , 3
                  CALL syslop(coord,site,pstr(cnt),bits(1,i),nchar)
                  cnt = nchar + cnt
               ENDDO
               IF ( r1.GT.0 ) THEN
                  CALL mult(bits(1,1),bits(1,2),temp)
                  CALL mult(temp,bits(1,3),u1)
               ELSE
                  DO i = 1 , 18
                     u1(i) = bits(i,2)
                  ENDDO
               ENDIF
C
C      NEXT GET THE MATRIX CORRESPONDING TO DIRECTION R2.
C
               DO i = 1 , 3
                  CALL syslop(coord,site,pstr(cnt),bits(1,i),nchar)
                  cnt = cnt + nchar
               ENDDO
               IF ( r2.GT.0 ) THEN
                  CALL mult(bits(1,1),bits(1,2),temp)
                  CALL mult(temp,bits(1,3),u2)
               ELSE
                  DO i = 1 , 18
                     u2(i) = bits(i,2)
                  ENDDO
               ENDIF
C
C       NOW LOOP OVER THE TWO (CURRENTLY) DIRECTIONS FOR THE R3 LEG
C       CALCULATING THE OBSERVABLES AS WE GO.
C
               DO odir = 1 , 2
                  DO i = 1 , 3
                     CALL syslop(coord,site,pstr(cnt),bits(1,i),nchar)
                     cnt = cnt + nchar
                  ENDDO
                  IF ( r3.GT.0 ) THEN
                     CALL mult(bits(1,1),bits(1,2),temp)
                     CALL mult(temp,bits(1,3),u3)
                  ELSE
                     DO i = 1 , 18
                        u3(i) = bits(i,2)
                     ENDDO
                  ENDIF
C
C        CALCULATE OBSERVABLES
C
                  CALL observ(tot,u1,u2,u3)
C
               ENDDO
C
C        MOVE TO NEXT SITE IN LATTICE.
C
               site = site + mov(1)
            ENDDO
         ENDDO
      ENDDO
C
      RETURN
      END
**==qqqlps.f  
 
C
C     J .W .FLOWER  AUG 1985
C
C
C     FROM ABOVE THE THREE QUARK ORIENTATION IS SPECIFIED AS FOLLOWS :
C
C
C                                   * ^
C                                   |
C                                   | R3
C                                   |
C                     * ------------|-------------------- *
C                     <---- R1 ---->|<-------- R2 -------->
C
C     THE LIMB WITH R1 & R2 IS TAKEN TO LIE IN DIRECTION "DIR" AND THE
C     OTHER LIMB IS TAKEN IN THE TWO TRANSVERSE DIRECTIONS SUCCESSIVELY.
C     NOTE : SINCE WE HAVE ACCESS TO ONLY BARRED 'LEGS' WE NEED TO MAKE
C     SOME TIMELIKE PIECES ON THE FLY.
C                                                                     */
C     CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE qqqlps(r1,r2,r3,t,dir,str)
      INTEGER r1 , r2 , r3 , t , dir
      INTEGER str(*)
      INTEGER ctr , odir , temp
C
C  DO R1 LIMB - GO BACK IN "DIR" FIRST.
C
      ctr = 1
      temp = 5 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      temp = 23 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
      temp = 27 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      IF ( r1.NE.0 ) THEN
         str(ctr) = 36
         ctr = ctr + 1
      ELSE
         DO i = 0 , t - 1
            str(ctr) = 4
            ctr = ctr + 1
         ENDDO
         DO i = 0 , t - 1
            str(ctr) = 30
            ctr = ctr + 1
         ENDDO
      ENDIF
C
      temp = 23 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
      temp = 27 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      DO i = 0 , t - 1
         str(ctr) = 26
         ctr = ctr + 1
      ENDDO
C
      temp = 1 + dir
      DO i = 0 , r1 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      DO i = 0 , t - 1
         str(ctr) = 30
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
C   DO R2 LIMB - GO FORWARD IN "DIR" FIRST.
C
      temp = 1 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      temp = 27 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
      temp = 23 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      IF ( r2.NE.0 ) THEN
         str(ctr) = 36
         ctr = ctr + 1
      ELSE
         DO i = 0 , t - 1
            str(ctr) = 4
            ctr = ctr + 1
         ENDDO
         DO i = 0 , t - 1
            str(ctr) = 30
            ctr = ctr + 1
         ENDDO
      ENDIF
C
      temp = 27 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
      temp = 23 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      DO i = 0 , t - 1
         str(ctr) = 26
         ctr = ctr + 1
      ENDDO
C
      temp = 5 + dir
      DO i = 0 , r2 - 1
         str(ctr) = temp
         ctr = ctr + 1
      ENDDO
C
      DO i = 0 , t - 1
         str(ctr) = 30
         ctr = ctr + 1
      ENDDO
C
      str(ctr) = 14
      ctr = ctr + 1
C
C  NOW GET OTHER LEG ... LOOP OVER TRANSVERSE DIRECTIONS
C   .... AND ALSO GO FORWARD AND (EVENTUALLY) BACK TO R3.
C
      DO ntrans = 1 , 2
         odir = mod((dir+ntrans),3)
         temp = 1 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         temp = 27 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         str(ctr) = 14
         ctr = ctr + 1
C
         temp = 23 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         IF ( r3.NE.0 ) THEN
            str(ctr) = 36
            ctr = ctr + 1
         ELSE
            DO i = 0 , t - 1
               str(ctr) = 4
               ctr = ctr + 1
            ENDDO
            DO i = 0 , t - 1
               str(ctr) = 30
               ctr = ctr + 1
            ENDDO
         ENDIF
         temp = 27 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         str(ctr) = 14
         ctr = ctr + 1
C
         temp = 23 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         DO i = 0 , t - 1
            str(ctr) = 26
            ctr = ctr + 1
         ENDDO
         temp = 5 + odir
         DO i = 0 , r3 - 1
            str(ctr) = temp
            ctr = ctr + 1
         ENDDO
         DO i = 0 , t - 1
            str(ctr) = 30
            ctr = ctr + 1
         ENDDO
         str(ctr) = 14
C
         ctr = ctr + 1
      ENDDO
C
      RETURN
      END
**==qqqmea.f  
 
 
 
 
C
C   J. W. FLOWER & STEVE OTTO  SEPTEMBER 1984
C
C   MEASURES THREE QUARK OBSERVABLES USING UBARS.
C
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C   AND DAVID WALKER, 1987
C
      SUBROUTINE qqqmea(r1,r2,r3,tt)
      INTEGER r1 , r2 , r3 , tt
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      INTEGER numstg
      COMMON /comfil/ numstg
C
      INTEGER dir , t
CMWB 6.22.88      CHARACTER*512 LPSTR
      INTEGER lpstr(512)
      DOUBLEPRECISION tot(2) , vals(2,15)
C
      vals(1,tt) = 0.0
      vals(2,tt) = 0.0
C
      DO t = 0 , latt(4) - 1
         CALL legmak(tt,t)
         tot(1) = 0.0
         tot(2) = 0.0
         DO dir = 0 , 2
            CALL qqqlps(r1,r2,r3,tt,dir,lpstr)
            CALL qqqlpc(r1,r2,r3,t,lpstr,tot)
         ENDDO
         vals(1,tt) = vals(1,tt) + tot(1)
         vals(2,tt) = vals(2,tt) + tot(2)
      ENDDO
C
      vals(1,tt) = vals(1,tt)/real(numstg*3*2)
      vals(2,tt) = vals(2,tt)/real(numstg*3*2)
      WRITE (6,1) r1 , r2 , r3 , tt , vals(1,tt) , vals(2,tt)
 1      format   ('Q (' , 3i3, ')' , '  X ' , i3 , ' : ' , 
     &            2e19.9)
C
      RETURN
      END
**==rotlps.f  
 
C   MAKE THE STRINGS USED TO MEASURE OFF-AXIS WILSON LOOPS.
C   STRING COMES IN THREE PARTS.
C    (i)  OUTWARD MOVES AND FAKE BACK
C    (ii) FETCH LEG AND FAKE BACK.
C    (iii GET TOP PART AND FAKE BACK.
C   BARS ALL T LINKS EXCEPT IF R1+R2=1.
C
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C   AND DAVID WALKER, 1987
C
      SUBROUTINE rotlps(pstr,r1,r2,t,dir1,dir2)
      INTEGER pstr(256)
      INTEGER r1 , r2 , t , dir1 , dir2
C
      INTEGER r , dir , cnt , temp
C
C   GO OUT IN SPATIAL DIRECTION.
C
      num1 = 0
      num2 = 0
      cnt = 1
  100 CONTINUE
      IF ( (num1.LT.r1) .AND. (num2.LT.r2) ) THEN
         temp = 1 + dir2
         pstr(cnt) = temp
         cnt = cnt + 1
         num2 = num2 + 1
         IF ( (num1+num2).LT.(r1+r2-1) ) THEN
            temp = 15 + dir1
         ELSE
            temp = 1 + dir1
         ENDIF
         pstr(cnt) = temp
         cnt = cnt + 1
         num1 = num1 + 1
         GOTO 100
      ENDIF
C
      IF ( num1.LT.r1 ) THEN
         dir = dir1
      ELSE
         dir = dir2
      ENDIF
C
      DO r = (num1+num2) , (r1+r2-1)
         IF ( (dir.EQ.dir1) .AND. (r.LT.r1+r2-1) .AND. (r.GT.0) ) THEN
            temp = 15 + dir
         ELSE
            temp = 1 + dir
         ENDIF
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
C   NOW FAKE BACK... THIS PART IS EASY !
C
      temp = 27 + dir1
      DO i = 0 , r1 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      temp = 27 + dir2
      DO i = 0 , r2 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      pstr(cnt) = 14
      cnt = cnt + 1
C
C   NOW GO FETCH THE LEG. IF R1+R2 = 1 THEN WE HAVE TO DO
C   THE LEG ON THE FLY
C
      temp = 23 + dir1
      DO i = 0 , r1 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      temp = 23 + dir2
      DO i = 0 , r2 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      IF ( (r1+r2).NE.1 ) THEN
         pstr(cnt) = 36
         cnt = cnt + 1
      ELSE
         DO i = 0 , t - 1
            pstr(cnt) = 4
            cnt = cnt + 1
         ENDDO
         DO i = 0 , t - 1
            pstr(cnt) = 30
            cnt = cnt + 1
         ENDDO
      ENDIF
      temp = 27 + dir1
      DO i = 0 , r1 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      temp = 27 + dir2
      DO i = 0 , r2 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      pstr(cnt) = 14
      cnt = cnt + 1
C
C   NOW FETCH THE TOP PART OF THE LOOP.
C
      temp = 23 + dir1
      DO i = 0 , r1 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      temp = 23 + dir2
      DO i = 0 , r2 - 1
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
      DO i = 0 , t - 1
         pstr(cnt) = 26
         cnt = cnt + 1
      ENDDO
C
      num1 = 0
      num2 = 0
  200 CONTINUE
      IF ( (num1.LT.r1) .AND. (num2.LT.r2) ) THEN
         temp = 5 + dir2
         pstr(cnt) = temp
         cnt = cnt + 1
         num2 = num2 + 1
         IF ( (num1+num2).LT.(r1+r2-1) ) THEN
            temp = 19 + dir1
         ELSE
            temp = 5 + dir1
         ENDIF
         pstr(cnt) = temp
         cnt = cnt + 1
         num1 = num1 + 1
         GOTO 200
      ENDIF
C
      IF ( num1.LT.r1 ) THEN
         dir = dir1
      ELSE
         dir = dir2
      ENDIF
      DO r = (num1+num2) , (r1+r2-1)
         IF ( (dir.EQ.dir1) .AND. (r.LT.(r1+r2-1)) .AND. (r.GT.0) ) THEN
            temp = 19 + dir
         ELSE
            temp = 5 + dir
         ENDIF
         pstr(cnt) = temp
         cnt = cnt + 1
      ENDDO
C
      DO i = 0 , t - 1
         pstr(cnt) = 30
         cnt = cnt + 1
      ENDDO
C
      pstr(cnt) = 14
C
      RETURN
      END
**==rotmea.f  
 
 
 
 
C
C   J. W. FLOWER  OCT 1985
C
C   THIS CODE MEASURES THE OFF AXIS WILSON LOOPS.
C
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C   AND DAVID WALKER, 1987
C
      SUBROUTINE rotmea(dir,numhit,r1max,r2max)
      INTEGER dir , numhit , r1max , r2max
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      PARAMETER (mbar=18*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION ubar(mbar) , sbar(mbar)
      COMMON /comubr/ ubar , sbar
C
      INTEGER tmax , t , plegs , cnt
      INTEGER coord(4) , nchar , tflag , site , tmp
      DOUBLEPRECISION bot(18) , top(18) , up(18) 
      DOUBLEPRECISION  down(18) , stemp(18,2)
      DOUBLEPRECISION temp
      INTEGER pstr(256) , lpstr(256)
C
      CALL linkbr(numhit,dir,sbar)
C
      DO icord4 = 0 , latt(4) - 1
         coord(4) = icord4
         DO t = 1 , 12
            tflag = 0
            DO i = 1 , r1max
               DO j = 1 , r2max
                  IF ( which(i,j).NE.0 ) THEN
                     tmax = min(12,which(i,j)+5)
                     IF ( (which(i,j).LE.t) .AND. (t.LE.tmax) )
     &                     tflag = 1
                  ENDIF
               ENDDO
            ENDDO
            IF ( tflag.EQ.1 ) THEN
               CALL legmak(t,icord4)
C
C    NOW LOOP THROUGH THE REQUIRED LOOPS.
C
               DO i = 1 , r1max
                  DO j = 1 , r2max
                     IF ( which(i,j).NE.0 ) THEN
                        wloop(i,j) = 0.
                        tmax = min(12,which(i,j)+5)
                        IF ( (t.GE.which(i,j)) .AND. (t.LE.tmax) ) THEN
                           CALL rotlps(lpstr,i,j,t,dir,mod((dir+1),3))
                           site = coord(4)*mov(4)
                           plegs = 1
                           DO icord3 = 0 , latt(3) - 1
                              coord(3) = icord3
                              DO icord2 = 0 , latt(2) - 1
                                 coord(2) = icord2
                                 DO icord1 = 0 , latt(1) - 1
                                    coord(1) = icord1
                                    CALL udag(legs(plegs),down)
C                        CALL CPYMAT(PSTR,LPSTR,256) LPOINTER 4-3-89
                                    CALL icpymt(pstr,lpstr,256)
                                    cnt = 1
                                    CALL syslop(coord,site,pstr(cnt),
     &                                    bot,nchar)
                                    cnt = cnt + nchar
                                    CALL syslop(coord,site,pstr(cnt),up,
     &                                    nchar)
                                    cnt = cnt + nchar
                                    CALL syslop(coord,site,pstr(cnt),
     &                                    top,nchar)
                                    CALL mult(bot,up,stemp(1,1))
                                    CALL mult(stemp(1,1),top,stemp(1,2))
                                    CALL mult(stemp(1,2),down,stemp(1,1)
     &                                    )
                                    wloop(i,j) = wloop(i,j)
     &                                     + (stemp(1,1)+stemp(9,1)
     &                                    +stemp(17,1))/3.0
                                    plegs = plegs + 18
                                    site = site + mov(1)
                                 ENDDO
                              ENDDO
                           ENDDO
                        ENDIF
                     ENDIF
                  ENDDO
               ENDDO
C
               DO i = 1 , r1max
                  DO j = 1 , r2max
                     IF ( which(i,j).NE.0 ) THEN
                        temp = wloop(i,j)
                        tmp = t - which(i,j) + 1
                        res(i,j,tmp) = res(i,j,tmp) + temp
                     ENDIF
                  ENDDO
               ENDDO
            ENDIF
         ENDDO
      ENDDO
C
      RETURN
      END
**==setcol.f  
 
      SUBROUTINE setcol(latt)
      INTEGER latt(4)
C
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (lismax=(l1*l2*l3*l4)/2+1)
      COMMON /colcom/ lisred(lismax) , lisblk(lismax) , nred , nblack
C
      nred = 0
      nblack = 0
      iptr = 0
      index4 = 0
      DO i = 1 , latt(4)
         index3 = index4
         DO j = 1 , latt(3)
            index2 = index3
            DO k = 1 , latt(2)
               index1 = index2
               DO l = 1 , latt(1)
                  IF ( index1.EQ.0 ) THEN
                     nred = nred + 1
                     lisred(nred) = iptr
                  ELSE
                     nblack = nblack + 1
                     lisblk(nblack) = iptr
                  ENDIF
                  iptr = iptr + 72
                  index1 = 1 - index1
               ENDDO
               index2 = 1 - index2
            ENDDO
            index3 = 1 - index3
         ENDDO
         index4 = 1 - index4
      ENDDO
C
      RETURN
      END
**==setpar.f  
 
C
C   STEVE OTTO AND JON FLOWER   JULY 1984.
C
C   RECEIVES  PARAMETER INPUT FOR THE SU3 ROUTINES.
C
C   CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C   AND DAVID WALKER, 1987
C
      SUBROUTINE setpar
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      INTEGER numstg
      COMMON /comfil/ numstg
C
  100 CONTINUE
      WRITE (6,*) 'ENTER X,Y,Z,T LATTICE SIZE ?'
      READ (5,*) latt
      WRITE (6,*) latt
      num = latt(1)*latt(2)*latt(3)*latt(4)
C
      IF ( num*72.GT.umax ) THEN
         WRITE (6,*) 'NOT ENOUGH ROOM FOR GAUGE FIELD'
         WRITE (6,*) '*******************************'
         WRITE (6,*) 'TRY AGAIN'
         GOTO 100
      ENDIF
C
      numsit = latt(1)*latt(2)*latt(3)*latt(4)
      numstg = latt(1)*latt(2)*latt(3)*latt(4)
      WRITE (6,*) 'ENTER VALUE OF BETA '
      READ (5,*) beta
      WRITE (6,*) beta
C
      CALL make()
C
      CALL setcol(latt)
C
      WRITE (6,*) 'PARAMETERS NOW SET.'
      WRITE (6,*) 'NOW NEED TO INITIALIZE..........'
C
      RETURN
      END
**==setsed.f  
 
      SUBROUTINE setsed()
C
      INTEGER seed , mode
C
      WRITE (6,*) 'SELECT MODE OF SEED SETTING: '
      WRITE (6,*) '0  <->  BY HAND'
      WRITE (6,*) '1  <->  READ IN FROM STANDARD FILE, SEED'
      READ (5,*) mode
      IF ( mode.EQ.0 ) THEN
         WRITE (6,*) 'INPUT SEED'
         READ (5,*) seed
         CALL prnset(seed)
      ELSE
         OPEN (UNIT=1,FILE='SEED',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
         REWIND (1)
         READ (1,99001) seed
99001    FORMAT (I12)
         CALL prnset(seed)
         CLOSE (1)
      ENDIF
      WRITE (6,*) 'SEED SET'
      RETURN
      END
**==syslop.f  
C
CCCCCCCCCCCCCCCCCCCCCCCCC END OF ROUTINE SETSED CCCCCCCCCCCCCCCCCCCC
C
 
 
 
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                     C
C          THIS SOFTWARE COPYRIGHTED BY CALTECH, 1983, 1986, 1987     C
C                                                                     C
C          PURE GAUGE SU(3) MONTE CARLO                               C
C                                                                     C
C          THIS CODE HAS A LONG HISTORY:                              C
C              EARLY VERSIONS OF THIS DEVELOPED FOR THE MARK I        C
C              HYPERCUBE BY STEVE OTTO AND PAUL STOLORZ (1981-83)     C
C              JON FLOWER AND STEVE OTTO EXTENDED AND IMPROVED THE    C
C              CODE FOR THE MARK II MACHINES. (1984-86)               C
C              THE CURRENT INCARNATION IS A "PORT" TO THE NCUBE       C
C              DONE BY CLIVE BAILLIE, ED FELTEN, STEVE OTTO,          C
C              AND DAVID WALKER  (1987)                               C
C                                                                     C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       SYSLOP SYSTEMATICALLY CALCULATES WILSON LOOPS
C       FOR SU3 THEORY IN 3+1 DIM WITH THE INT STRINGS
C       AS INPUT, ALSO KEEPS TRACK OF THE LENGTH OF THE STRING.
C
C       CODING OF LATTICE PATHS:
C                1, 2, 3, 4 +X,+Y,+Z,+T
C                5, 6, 7, 8 -X,-Y,-Z,-T
C               15,16,17,18  X,+Y,+Z,+T USING BARS.
C               19,20,21,22 -X,-Y,-Z,-T USING BARS.
C               23,24,25,26 +X,+Y,+Z,+T FAKE STEPS
C               27,28,29,30 -X,-Y,-Z,-T FAKE STEPS
C               36          RETURNS LEG
C
C       CONVERTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C       AND DAVID WALKER, 1987
C
      SUBROUTINE syslop(coord,site,ptr,ploop,nn)
      INTEGER coord(4) , nn , site , sites
      INTEGER ptr(*)
      DOUBLEPRECISION ploop(18)
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      PARAMETER (mbar=18*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      DOUBLEPRECISION ubar(mbar) , sbar(mbar)
      COMMON /comubr/ ubar , sbar
C
      DOUBLEPRECISION ftemp(18,2) , filmat(18) 
      DOUBLEPRECISION extra(18)
      INTEGER pu , pbar
      INTEGER plegs , tindex , setflg
C
      sites = site
      index = 0
      setflg = 0
      nn = 0
C
  100 CONTINUE
      nn = nn + 1
      ic = ptr(nn)
      IF ( ic.EQ.14 ) THEN
C
C
         DO i = 1 , 18
            ploop(i) = ftemp(i,index+1)
         ENDDO
C
         RETURN
      ELSE
C
         IF ( ic.LE.4 ) THEN
            ind = ic
            pu = sites + rot(ind)
            IF ( setflg.EQ.0 ) THEN
               CALL cpymat(ftemp(1,index+1),u1(pu+1),18)
               setflg = 1
            ELSE
               CALL cpymat(filmat,u1(pu+1),18)
               tindex = 1 - index
               CALL mult(ftemp(1,index+1),filmat,ftemp(1,tindex+1))
               index = tindex
            ENDIF
            coord(ind) = mod(coord(ind)+1,latt1(ind)+1)
            IF ( coord(ind).EQ.0 ) THEN
               sites = sites - mov(ind)*latt1(ind)
            ELSE
               sites = sites + mov(ind)
            ENDIF
C
         ELSEIF ( ic.LE.8 ) THEN
            ind = ic - 4
            IF ( coord(ind).EQ.0 ) THEN
               coord(ind) = latt1(ind)
               sites = sites + mov(ind)*latt1(ind)
            ELSE
               coord(ind) = coord(ind) - 1
               sites = sites - mov(ind)
            ENDIF
            pu = sites + rot(ind)
            IF ( setflg.EQ.0 ) THEN
               CALL cpymat(filmat,u1(pu+1),18)
               CALL udag(filmat,ftemp(1,index+1))
               setflg = 1
            ELSE
               CALL cpymat(extra,u1(pu+1),18)
               CALL udag(extra,filmat)
               tindex = 1 - index
               CALL mult(ftemp(1,index+1),filmat,ftemp(1,tindex+1))
               index = tindex
            ENDIF
C
         ELSEIF ( ic.EQ.36 ) THEN
            lsite = sites - coord(4)*mov(4)
            plegs = (lsite/72)*18
            DO i = 1 , 18
               ftemp(i,index+1) = legs(plegs+i)
            ENDDO
            setflg = 1
C
         ELSEIF ( ic.LE.18 ) THEN
            ind = ic - 14
            pbar = (sites/72)*18
            IF ( ind.EQ.4 ) THEN
               CALL cpymat(filmat,ubar(pbar+1),18)
            ELSE
               CALL cpymat(filmat,sbar(pbar+1),18)
            ENDIF
            IF ( setflg.EQ.0 ) THEN
               DO j = 1 , 18
                  ftemp(j,index+1) = filmat(j)
               ENDDO
               setflg = 1
            ELSE
               tindex = 1 - index
               CALL mult(ftemp(1,index+1),filmat,ftemp(1,tindex+1))
               index = tindex
            ENDIF
            coord(ind) = mod(coord(ind)+1,latt1(ind)+1)
            IF ( coord(ind).EQ.0 ) THEN
               sites = sites - mov(ind)*latt1(ind)
            ELSE
               sites = sites + mov(ind)
            ENDIF
C
         ELSEIF ( ic.LE.22 ) THEN
            ind = ic - 18
            IF ( coord(ind).EQ.0 ) THEN
               coord(ind) = latt1(ind)
               sites = sites + mov(ind)*latt1(ind)
            ELSE
               coord(ind) = coord(ind) - 1
               sites = sites - mov(ind)
            ENDIF
            pbar = (sites/72)*18
            IF ( ind.EQ.4 ) THEN
               CALL cpymat(extra,ubar(pbar+1),18)
            ELSE
               CALL cpymat(extra,sbar(pbar+1),18)
            ENDIF
            IF ( setflg.EQ.0 ) THEN
               CALL udag(extra,ftemp(1,index+1))
               setflg = 1
            ELSE
               CALL udag(extra,filmat)
               tindex = 1 - index
               CALL mult(ftemp(1,index+1),filmat,ftemp(1,tindex+1))
               index = tindex
            ENDIF
C
         ELSEIF ( ic.LE.26 ) THEN
            ind = ic - 22
            coord(ind) = mod(coord(ind)+1,latt1(ind)+1)
            IF ( coord(ind).EQ.0 ) THEN
               sites = sites - mov(ind)*latt1(ind)
            ELSE
               sites = sites + mov(ind)
            ENDIF
C
         ELSEIF ( ic.LE.30 ) THEN
            ind = ic - 26
            IF ( coord(ind).EQ.0 ) THEN
               coord(ind) = latt1(ind)
               sites = sites + mov(ind)*latt1(ind)
            ELSE
               coord(ind) = coord(ind) - 1
               sites = sites - mov(ind)
            ENDIF
C
         ENDIF
         GOTO 100
      ENDIF
      END
**==udag.f  
 
C
C   UDAG EVALUATES THE COMPLEX CONJUGATE OF THE INPUT 3x3 MATRIX
C
      SUBROUTINE udag(pin,pout)
      DOUBLEPRECISION pin(18) , pout(18)
C
      pout(1) = pin(1)
      pout(7) = pin(3)
      pout(13) = pin(5)
      pout(2) = -pin(2)
      pout(8) = -pin(4)
      pout(14) = -pin(6)
      pout(3) = pin(7)
      pout(9) = pin(9)
      pout(15) = pin(11)
      pout(4) = -pin(8)
      pout(10) = -pin(10)
      pout(16) = -pin(12)
      pout(5) = pin(13)
      pout(11) = pin(15)
      pout(17) = pin(17)
      pout(6) = -pin(14)
      pout(12) = -pin(16)
      pout(18) = -pin(18)
C
      RETURN
      END
**==update.f  
 
 
 
C
C     J. W. FLOWER AND STEVE OTTO   JULY 1984.
C
C     UPDATES THE LATICE "SWEEPS" TIMES WITH PROJECTIONS EVERY "PINVAL"
C
C     CONVERRTED FROM C TO FORTRAN BY KEN BARISH, CLIVE BAILLIE,
C     AND DAVID WALKER, 1987
C
      SUBROUTINE update(sweeps,pinval)
      INTEGER sweeps , pinval
C
      INTEGER dims
      PARAMETER (dims=4)
      PARAMETER (l1=12,l2=12,l3=12,l4=12)
      PARAMETER (mleg=18*l1*l2*l3)
      INTEGER umax
      PARAMETER (umax=72*l1*l2*l3*l4)
      INTEGER latt(dims) , latt1(dims) , mov(dims) , rot(dims) , 
     &        which(13,13)
      DOUBLEPRECISION wloop(13,13) , res(13,13,6) 
      DOUBLEPRECISION beta , legs(mleg)
      COMMON /comsys/ beta , wloop , res , legs , latt , latt1 , mov , 
     &                rot , which
      PARAMETER (lismax=(l1*l2*l3*l4)/2+1)
      COMMON /colcom/ lisred(lismax) , lisblk(lismax) , nred , nblack
C
C
      DOUBLEPRECISION u1(umax)
      COMMON /comu1 / u1
C
      INTEGER coord(4) , dir , origin
C
      DO loop = 0 , sweeps - 1
C
C  FIRST DO RED POINTS
C
         DO i = 1 , nred
            origin = lisred(i)
            index = origin/72
            CALL getpos(latt,coord,index)
            DO dir = 0 , 3
               CALL hit(coord,origin,dir)
            ENDDO
         ENDDO
C
C  NOW DO BLACK POINTS
C
         DO i = 1 , nblack
            origin = lisblk(i)
            index = origin/72
            CALL getpos(latt,coord,index)
            DO dir = 0 , 3
               CALL hit(coord,origin,dir)
            ENDDO
         ENDDO
C
         IF ( mod(loop,pinval).EQ.(pinval-1) ) CALL projec()
      ENDDO
C
      RETURN
      END

>Fix:
>Release-Note:
>Audit-Trail:
>Unformatted:


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]