
      SUBROUTINE INITPD( PSI, PSIDOT, RANDOM)  
C     ---------------------------------------

      INTEGER VSIZE
      COMMON /VECLEN/ VSIZE

      DOUBLE PRECISION  PSI(VSIZE), PSIDOT(VSIZE), RANDOM(VSIZE)

C     --------------------------------------------------------------------
C     |  Initialise the fermion vector  {psidot}  at t = 0 .             |
C     |  The correct distribution is obtained by solving :               |
C     |                                                                  |
C     |               [-D + m] {psidot}  =  {z1 + i.z2}                  |
C     |                                                                  |
C     |  where z1, z2 are normally distributed with std.dev. SQRT (2).   |
C     |  [D] is the Dirac matrix, which is a function of the gauge       |
C     |  fields (see the introduction).                                  |
C     |                                                                  |
C     |    This is solved for {psidot} using Conjugate Gradient (CG)     |
C     |  iteration.  However, CG can only solve for hermitian matrices,  |
C     |  so need to multiply both sides by  hc.[-D + m] = [ D + m],      |
C     |  and solve :                                                     |
C     |                                                                  |
C     |           [-D*D + m*m] {psidot}  =  [D + m] {z1 + i.z2}          |
C     |                                                                  |
C     |  N.B. !! Need the gauge fields U to be defined already !!        |
C     --------------------------------------------------------------------

      INTEGER     TDIR,  XDIR,  YDIR,  ZDIR,  EVEN,  ODD,
     *            PSENSE,  NSENSE,  PLUS,  MINUS,  RE,  IM
      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 )

      INTEGER  NSWEEP
      INTEGER  MAXIT, SIGFIG
      DOUBLE PRECISION     H,  HALFH,  MQUARK,  MSQ
      COMMON / PARAMS / H,HALFH,MQUARK,MSQ,NSWEEP,MAXIT,SIGFIG

      INTEGER SEED, INNODE 
      DOUBLE PRECISION N2
      LOGICAL  STORED   
      COMMON / RANDOM / N2, SEED, INNODE, STORED

C     ---------------------------------------------------------------------
C     |(i)                                                                |
C     |                 {random} = [D] {z1 + i.z2}                        |
C     |                                                                   |
C     |  {z1 + i.z2}  is put in  {psidot}  for temporary storage.  It is  |
C     |  defined at ODD sites, so  {random}  is defined at EVEN sites.    |
C     ---------------------------------------------------------------------

      CALL RANDV (PSIDOT)                   
      CALL DVEC (PLUS, ODD, PSIDOT, RANDOM)
C  @@@@  Communication Structure  @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C  @                                                                  @
      IF( INNODE .EQ. 0) PRINT *, 'Initial Dvec finished'
C  @                                                                  @
C  @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
                             
C     ---------------------------------------------------------------------
C     |(ii)                                                               |
C     |                 {random} = {random} + m {z1 + i.z2}               |
C     |                                                                   |
C     |  where {z1 + i.z2} is defined on EVEN sites this time and so is   |
C     |  a different random vector.  (Again it is put in  {psidot}  for   |
C     |  temporary storage.                                               |
C     |    Therefore,                                                     |
C     |                 {random} = [D + m] {z1 + i.z2}                    |
C     |                                                                   |
C     |  defined at even sites.                                           |
C     ---------------------------------------------------------------------

      CALL RANDV (PSIDOT)
      CALL ADDSCV (RANDOM, RANDOM, MQUARK, PSIDOT)

C     ----------------------------------------
C     |(iii)                                 |
C     |  Solve    [F] {psidot} = {random}    |
C     ----------------------------------------

      CALL SOLVE (PSIDOT, RANDOM)

      END
