
      SUBROUTINE SOLVE (X, B)
C     -----------------------

      INCLUDE 'qcd2.inc'

      DOUBLE PRECISION  X(VSIZE), B(VSIZE)

C     -------------------------------------------------------------------
C     |  Uses the CG algorithm to solve :                               |
C     |                        [F] {x} = {b}                            |
C     |  for {x}, given {b}.                                            |
C     |                                                                 |
C     |  {}    denotes a fermion vector (defined on even sites)         |
C     |  [F] = (- [D][D] + msq),  where [D] is the Dirac matrix.        |
C     |                                                                 |
C     |  At start, {x} has the initial trial solution.                  |
C     -------------------------------------------------------------------

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

      INTEGER  FLOPS, NITER
      DOUBLE PRECISION   RSQ,  RSQTOL,  SCPROD, TSTART, TEND, TIME
      COMMON / RESULT /  TIME, FLOPS, NITER

      DOUBLE PRECISION   P( VSIZE ), R( VSIZE )
      DOUBLE PRECISION   FV(VSIZE), DV(VSIZE)
 
C     ---------------------------------------------------
C     |  Initialise {p}, {r}, rsq, rsqtol  and  niter   |
C     |                                                 |
C     |     {p}  =  {r}  =  {b} - {Fx}                  |
C     |     rsq  =  {r}.{r}                             |
C     ---------------------------------------------------

      CALL FVEC ( X, R, DV )
      CALL SUBVEC ( R, B, R )
      CALL COPVEC ( P, R )
      RSQTOL  =  SCPROD (B, B)  * (1.0D-02 ** (SIGFIG))
      RSQ     =  SCPROD (R, R)
      NITER   =  0

C     ------------------------------------------
C     |  Do CG iterations until {x} converged  |
C     ------------------------------------------

      CALL TIMER (TSTART)

 1    CONTINUE
        NITER  =  NITER + 1
        CALL CGITER ( X, P, R, FV, DV, RSQ )
      IF  ((RSQ .GT. RSQTOL)  .AND.  (NITER .LT. MAXIT)) GO TO 1

      CALL TIMER (TEND)

C------ 'time'   = average time (in secs) per CG iteration.

      TIME  =  (TEND - TSTART) / DBLE (NITER)
      IF (TIME .LE. 0.0) TIME   =  - 1.0D20

      IF (NCOLOR .EQ. 2) THEN
         FLOPS = 280 * NSITES
      ELSE IF (NCOLOR .EQ. 3) THEN
         FLOPS = 612 * NSITES
      ENDIF

      FLOPS = FLOPS * PROCX * PROCY * PROCZ

      END
