PRACE NPT

      Subroutine cg( n1, n2, n3, m, mp, a, al, b, x, gamma, maxit,

     &               conv, tol, prec )

! ---------------------------------------------------------------------

! --- cg does an iterative solve of ax = b.

!     Real(l_) a(n1*n2*n3,0:3) : Main diagonal and the 3 upper

!                                diagonals of the matrix.

!     Real(l_) al(n1*n2*n3,1:3): The 3 lower diagonals. Identical to

!                                the upper diagonals but used for

!                                efficiency reasons.

!     Real(l_) b(n1*n2*n3)     : The right hand side.

!     Real(l_) x(n1*n2*n3)     : On input the initial estimate of the

!                                solution, on output the solution

!                                (hopefully).

!     Real(l_) gamma(m+1)      : Polynomial coefficients used in the

!                                preconditioning.

!     Integer  maxit           : On input the maximum number of

!                                iterations allowed.

!     Real(l_) conv(maxit)     : Residuals of the sequence of iter's.

!     Real(l_) tol             : Tolerance used as a stop criterion.

!     External prec            : Name of subroutine that implements the

!                                preconditioner.

! ---------------------------------------------------------------------

      Use         numerics

      Use         floptime

      Use         mpi_module

      Implicit    None

      Integer  :: n1, n2, n3, m, mp

      Real(l_) :: a(m,0:3), al(m,1:3), b(n1*n2*n3), x(n1*n2*n3)

      Real(l_) :: gamma(mp+1)

      Integer  :: maxit

      Real(l_) :: conv(maxit)

      Real(l_) :: tol

      External    prec

      Integer  :: i, it, its, j, ntot

      Real(l_) :: ap(n1*n2*n3), p(n1*n2*n3), r(n1*n2*n3)

      Real(l_) :: alpha, beta, nr0, nrm, nnrm, pap, tol2

      Real(l_) :: dotpr

      External    dotpr

! ---------------------------------------------------------------------

      ntot = n1*n2*n3

      tol2 = tol*tol

      Call sym7mxv( n1, n2, n3, m, a, al, x, r )

      r(lb:gub) = b(lb:gub) - r(lb:gub)

! ---------------------------------------------------------------------

! --- Precondition of initial residual.

      Call prec( n1, n2, n3, m, mp, a, al, nr0, r, p, gamma )

      nrm = nr0

! ---------------------------------------------------------------------

! --- Iterate at most maxit times.  

      its = maxit

      Do i = 1, its

         Call sym7mxv( n1, n2, n3, m, a, al, p, ap )

         pap = dotpr( ntot, p, ap )

         alpha = nrm/pap

         x(lb:gub) = x(lb:gub) + alpha*p(lb:gub)

         r(lb:gub) = r(lb:gub) - alpha*ap(lb:gub)

         Call prec( n1, n2, n3, m, mp, a, al, nnrm, r, ap, gamma )

         conv(i) = Sqrt( nnrm )

         If ( nnrm < nr0*tol2 ) Then

            maxit = i

            Go To 1000

         End If

         beta = nnrm/nrm

         nrm = nnrm

         p(lb:gub) = ap(lb:gub) + beta*p(lb:gub)

      End Do

!     Print *, 'CG: no convergence in ', maxit, ' iterations.'

 1000 flops = flops + maxit*( 6*m + 4 ) + m + 1

! ---------------------------------------------------------------------

      End Subroutine cg