
       PROGRAM QCD2
C   ==================================================================
C   ===                                                            ===
C   ===           GENESIS Distributed Memory Benchmarks            ===
C   ===                                                            ===
C   ===            Conjugate Gradient Kernel from QCD              ===
C   ===          Single Node Standard Fortran77 version            ===
C   ===                                                            ===
C   ===     Original author:        John Merlin                    ===
C   ===     Mofified by    :        Ivan Wolton                    ===
C   ===     Department of Electronics and Computer Science         ===
C   ===               University of Southampton                    ===
C   ===               Southampton SO9 5NH, U.K.                    ===
C   ===     fax.:+44-703-593045   e-mail:icw@uk.ac.soton.ecs       ===
C   ===                                                            ===
C   ===          Last update: April 1991; Release: 2.0             ===
C   ===                                                            ===
C   ==================================================================
      INCLUDE 'qcd2.inc'

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

      DOUBLE PRECISION  PSI( VSIZE ), PSIDOT( VSIZE )
      COMMON / SITES / PSI, PSIDOT

      INTEGER CVSIZE
      COMMON /VECLEN/ CVSIZE

      DOUBLE PRECISION   RHS(VSIZE), TEMP(VSIZE), DVWORK(VSIZE)
      INTEGER  ISWEEP, NLATT
      LOGICAL LASTSW
      DATA  LASTSW /.FALSE./

      CVSIZE = VSIZE
      NLATT = NSITES

C     --------------------------------------      
C     |  Input and initialise fields :     |
C     --------------------------------------      
      CALL SETPAR
      CALL HEAD
      CALL INIT( PSI, PSIDOT, TEMP, NCOLOR)
      CALL OUTPUT( PSI, PSIDOT, TEMP, DVWORK, NLATT, LASTSW)

      DO 1088   ISWEEP = 1, NSWEEP
C       ---------------------------------------------------------------
C       |  Update the vector {psi} :                                  |
C       |                                                             |
C       |     {psi (1/2)} = {psi (-1/2)} + h {psidot (0)}             |
C       |                                                             |
C       |  (The numbers in brackets represent the timestep).          |
C       ---------------------------------------------------------------
        CALL ADDSCV ( PSI, PSI, H, PSIDOT )


C       ---------------------------------------------------------------
C       |  To update the vector {psidot}, solve                       |
C       |                                                             |
C       |     [F] {psidot (1)} = [F] {psidot (0)} - h * psi (1/2)}    |
C       |                                                             |
C       |  (using the CG algorithm)                                   |
C       ---------------------------------------------------------------
        CALL FVEC   ( PSIDOT, RHS, DVWORK)
        CALL ADDSCV ( RHS, RHS, - H, PSI )
        CALL SOLVE  ( PSIDOT, RHS)

        IF( ISWEEP .EQ. NSWEEP ) LASTSW = .TRUE.
        CALL OUTPUT( PSI, PSIDOT, TEMP, DVWORK, NLATT, LASTSW)

 1088 CONTINUE
      STOP
      END
