
      SUBROUTINE DVEC (SIGN, USITES, U, DU)
C     -------------------------------------

      INCLUDE 'qcd2.inc'

      INTEGER  SIGN, USITES
      DOUBLE PRECISION  U( VSIZE ), DU( VSIZE ), V( VSIZE ),
     *      USHIFT( VSIZE )

      INTEGER NINIT
      PARAMETER ( NINIT = NSITES * NCOLOR**(NCOLOR-1) )
      DOUBLE PRECISION     GAUGE ( NINIT, 4, 0:1 )
      COMMON / LINKS / GAUGE

      INTEGER     I, TDIR,  XDIR,  YDIR,  ZDIR,  EVEN,  ODD,
     *            PSENSE,  NSENSE,  PLUS,  MINUS,  RE,  IM,  VSITES
      PARAMETER ( TDIR = 1,    XDIR = 2,    YDIR = 3,    ZDIR = 4,
     *            EVEN = 0,     ODD = 1,   PSENSE= 0,   NSENSE= 1,
     *            PLUS = 0,   MINUS = 1,      RE = 1,      IM = 2 )

C     -----------------------------------------------------------------------
C     |       Performs the [matrix] * {vector} product :                    |
C     |                                                                     |
C     |              {Du} = +/- [D] {u}                                     |
C     |                                                                     |
C     |  [D]         is the Dirac matrix                                    |
C     |  {u}, {Du}   are fermion vectors (defined only at even or odd sites)|
C     |  sign      = 'plus' or 'minus' (the sign on the rhs of above eqn)   |
C     |  usites    = 'odd'  or 'even'  (the parity of the sites at which    |
C     |                                                   {u} is defined).  |
C     |                                                                     |
C     |  Explanation                                                        |
C     |  -----------                                                        |
C     |      [D] only has non-zero elements between sites at opposite ends  |
C     |  of the same link.                                                  |
C     |      Consider a general site 'i'.  Let '(i + mu)' represent the site|
C     |  at the end of the link emerging from site 'i' in the POSITIVE  'mu'|
C     |  direction,  and '(i - mu)' represent the site at the end of the    |
C     |  NEGATIVE 'mu' link.  Then, {Du} at site 'i' is the sum of 8 terms  |
C     |  (one for each link from 'i'):                                      |
C     |                                                                     |
C     |                                                                     |
C     |  {Du}i  =  SUM  D(i,i+mu) * u(i+mu)   +   SUM  D(i,i-mu) * u(i-mu)  |
C     |         +ve dirns                      -ve dirns                    |
C     |             mu                             mu                       |
C     |                                                                     |
C     |                                                                     |
C     |       Each term in the sum is evaluated (for all sites 'i') by a    |
C     |  single call to DVECDIR.   DVECDIR contains a loop over sites,      |
C     |  which may be vectorised.                                           |
C     |---------------------------------------------------------------------|
C     |  Negative directions                                                |
C     |  -------------------                                                |
C     |                                                                     |
C     |       v (i) = D (i,j) * u (j),   where site j = (i-mu).             |
C     |                                                                     |
C     |                    u(j)             v(i)                            |
C     |                      <---- D(i,j) -----                             |
C     |                      j                i                             |
C     |                                                                     |
C     |       In this case :                                                |
C     |               D(i,j) = - herm. conj. [ gauge(j,i) ],                |
C     |  so the relevant gauge matrix is associated with the site 'j' at    |
C     |  the end of the link (where u (j) is defined).                      |
C     |       To vectorise over sites, the elements of {u}, {v} and         |
C     |  [gauge] must be aligned (i.e. corresponding elements stored at     |
C     |  corresponding positions).  Hence:                                  |
C     |                                                                     |
C     |  (1) Call DVECDIR to evaluate :                                     |
C     |                                                                     |
C     |         'vshifted (j)'  =  - herm.conj. [gauge(j,i)] * u(j)         |
C     |                                                                     |
C     |    for all 'usites' j (in DVECDIR).  (N.B. By 'usites' we mean      |
C     |    the sites, even or odd, at which vector {u} is defined).         |
C     |      Now the result for site 'i',  v (i),  is held in the vector    |
C     |    {vshifted} at a position corresponding to site 'j'.              |
C     |                                                                     |
C     |  (2) Shift all the 'vshifted (j)' in the positive direction         |
C     |    (along the link) to the corresponding sites 'i':                 |
C     |    vshifted (j)  -->  v (i)    (using SHIFTVEC).                    |
C     |                                                                     |
C     |  (3) Accumulate {v} in {Du}.                                        |
C     |                                                                     |
C     |---------------------------------------------------------------------|
C     |  Positive directions                                                |
C     |  -------------------                                                |
C     |                                                                     |
C     |       v (i) = D (i,j) * u (j),   where site j = (i+mu).             |
C     |                                                                     |
C     |                    v(i)             u(j)                            |
C     |                      ----- D(i,j) ---->                             |
C     |                      i                j                             |
C     |                                                                     |
C     |       In this case,  D(i,j) = gauge(i,j),  so the appropriate       |
C     |   gauge matrix is associated with site 'i'  (at the start of the    |
C     |   link, the same as v (i)).                                         |
C     |       To vectorise over sites, the elements of {u}, {v} and         |
C     |  [gauge] must be aligned (i.e. corresponding elements stored at     |
C     |  corresponding positions).  Hence:                                  |
C     |                                                                     |
C     |  (1) Shift all the u (j) in the negative direction (along the link) |
C     |    to the corresponding sites 'i':    u (j) --> ushifted (i)        |
C     |    (using SHIFTVEC).                                                |
C     |                                                                     |
C     |  (2) Perform                                                        |
C     |                 v(i)  =  gauge(i,j) * ushifted(i)                   |
C     |                                                                     |
C     |    for all 'vsites' i (in DVECDIR).  (N.B. 'vsites' means the       |
C     |    sites, even or odd, at which vector {v} is defined).             |
C     |      The result v (i)  for site 'i' is now held at the correct      |
C     |    position in vector {v}, so no further shift is necessary.        |
C     |                                                                     |
C     |  (3) Accumulate {v} in {Du}.                                        |
C     |                                                                     |
C     |---------------------------------------------------------------------|
C     |  Communications                                                     |
C     |  --------------                                                     |
C     |       When a vector is shifted in one of the dimensions that are    |
C     |  distributed over several processors (viz. x, y or z), the 'edge    |
C     |  sites in that direction have to be shifted to and from the         |
C     |  neighbour processors in that direction.  Eg, for a +x shift of     |
C     |  vector u (t,x,y,z..), the 'edge', u (t,segx,y,z..) for all t,y,z.. |
C     |  has to be moved TO the +x neighbour, and the 'edge' u (t,1,y,z..)  |
C     |  is input FROM the -x neighbour.  These edge communications are     |
C     |  INITIATED by calls to  SEND,  but they do not have to be           |
C     |  COMPLETED until the corresponding RECEIVE.                         |
C     |  This allows edge communications to be overlapped with internal     |
C     |  shifts.                                                            |
C     -----------------------------------------------------------------------



C     -------------------------------
C     |   Positive directed links   |
C     -------------------------------
      VSITES = 1 - USITES

C------ + t direction :

      CALL  SHIFTV ( U, USITES, TDIR, NSENSE, USHIFT )
      CALL  DVDIR ( PSENSE, GAUGE (1, TDIR, VSITES), USHIFT, DU )

C------ - t direction :

      CALL  DVDIR ( NSENSE, GAUGE (1, TDIR, USITES), U, USHIFT )
      CALL  SHIFTV ( USHIFT, USITES, TDIR, PSENSE, V )
      DO 20 I = 1, VSIZE
 20     DU(I) = DU(I) + V(I)

C------ + x direction

      CALL  SHIFTV ( U, USITES, XDIR, NSENSE, USHIFT )
C     CALL  blshft ( u, usites, xdir, nsense, ushift )

      CALL  DVDIR ( PSENSE, GAUGE (1, XDIR, VSITES), USHIFT, V )
      DO 120 I = 1, VSIZE
 120    DU(I) = DU(I) + V(I)

C------ + y direction

      CALL  SHIFTV ( U, USITES, YDIR, NSENSE, USHIFT )
C     CALL  blshft ( u, usites, ydir, nsense, ushift )

      CALL  DVDIR ( PSENSE, GAUGE (1, YDIR, VSITES), USHIFT, V )
      DO 220 I = 1, VSIZE
 220    DU(I) = DU(I) + V(I)

C------ + z direction

C     CALL  shiftv ( u, usites, zdir, nsense, ushift )
      CALL  BLSHFT ( U, USITES, ZDIR, NSENSE, USHIFT )

      CALL  DVDIR ( PSENSE, GAUGE (1, ZDIR, VSITES), USHIFT, V )
      DO 320 I = 1, VSIZE
 320    DU(I) = DU(I) + V(I)


C     -------------------------------
C     |   Negative directed links   |
C     -------------------------------

C------ - x direction :

      CALL  DVDIR ( NSENSE, GAUGE (1, XDIR, USITES), U, USHIFT )

      CALL  SHIFTV ( USHIFT, USITES, XDIR, PSENSE, V )
C     CALL  blshft ( ushift, usites, xdir, psense, v )

      DO 420 I = 1, VSIZE
 420    DU(I) = DU(I) + V(I)

C------ - y direction :

      CALL  DVDIR ( NSENSE, GAUGE (1, YDIR, USITES), U, USHIFT )

      CALL  SHIFTV ( USHIFT, USITES, YDIR, PSENSE, V )
C     CALL  blshft ( ushift, usites, ydir, psense, v )

      DO 520 I = 1, VSIZE
 520    DU(I) = DU(I) + V(I)

C------ - z direction :

      CALL  DVDIR ( NSENSE, GAUGE (1, ZDIR, USITES), U, USHIFT )

C     CALL  shiftv ( ushift, usites, zdir, psense, v )
      CALL  BLSHFT ( USHIFT, USITES, ZDIR, PSENSE, V )

      DO 620 I = 1, VSIZE
 620    DU(I) = DU(I) + V(I)


C------ Result has now been accumulated in 'Du'.  All that remains is to 
C------ scale 'Du' by  +/-(1/2).

      IF  (SIGN  .EQ.  PLUS)  THEN
        DO 700 I =1, VSIZE
 700      DU(I) =   0.5D0 * DU(I)
      ELSE 
        DO 720 I =1, VSIZE
 720      DU(I) =  -0.5D0 * DU(I)
      ENDIF

      END
