PRACE NPT

      Subroutine rgmres( n, nl, nel, mp, indx, rowp, matvals, q, x, b,

     &                   gamma, maxit, tol, exnrm, prec )

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

! --- rgmres does an iterative solve on a local part of Ax = b, where A

!     is in CRS format:

!     Integer  indx(nel),

!     Integer  rowp(n+1), and

!     Real(l_) matvals(nel). 

! --- The vectors b, x, and q are completely in memory but only a local

!     of x is updated (and combined later).

!     Real(l_) b(n)       : The righthand side.

!     Real(l_) x(n)       : On input the initial guess of the solution

!                           On convergence 'x' contains the solution.

!     Real(l_) q(n)       : Contains a (left) preconditioning vector.

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

!                           preconditioning.

!     Integer  maxit      : The maximum number of iterations allowed.

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

!     Real(l_) exnrm      : Contains the norm of the residual on exit.

!     External prec       : Name of subroutine performing the

!                           preconditioning.

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

      Use                    numerics

      Use                    floptime

      Use                    mpi_module

      Implicit               None

      Integer             :: n, nl, nel, mp, maxit

      Integer             :: indx(nel), rowp(nl+1)

      Real(l_)            :: matvals(nel)

      Real(l_)            :: q(n), x(n), b(n)

      Real(l_)            :: gamma(mp+1), tol, exnrm

      External               prec

      Integer             :: i, iter, j, jc, k, k0, k1

      Logical             :: conv

      Integer, Parameter  :: ib = 10, ir = 50

      Real(l_)            :: g(ir+2), rho(ir), rm(ir+1,ir+1), v(ib*n)

      Real(l_)            :: alpha, beta, eta, rhstp, tau1, tau2

      Real(l_)            :: r(n), w(n), z(n)

      Real(l_)            :: dotpr, nrm2

      External               dotpr, nrm2

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

      conv = .FALSE.

      rhstp = Sqrt( nrm2( nl, b(lb) ) )

      Call spmxv( nl, nel, indx, rowp, matvals, x, w(lb) )

      z(lb:gub) = b(lb:gub) - w(lb:gub)

      Call prec( nl, nel, mp, indx, rowp, matvals, z, r, gamma )

      beta = Sqrt( nrm2( nl, r(lb) ) )

      Do i = 1, maxit

         iter = i

         g(1) = beta

         g(2) = beta

         If ( beta == 0.0_l_ ) Then

            Print *, 'Stop: beta = 0'; Stop

         End If

         v(lb:gub) = r(lb:gub)/beta

         k0 = 0

         Do j = 1, ib

            jc = j

            Call spmxv( nl, nel, indx, rowp, matvals, v(k0+1:k0+n),

     &                  w(lb) )

            Call prec( nl, nel, mp, indx, rowp, matvals, w, z, gamma )

            Call MPI_Allgatherv( z(lb), nl, rtyp, z, sizes, offset,

     &                           rtyp, comm, ierr )

            k1 = 0

            Do k = 1, j

               rm(k,j) = dotpr( n, v(k1+1), z )

               k1 = k1 + n

            End Do

            k1 = 0

            w = 0.0_l_

            Do k = 1, j

               w(lb:gub) = w(lb:gub) + rm(k,j)*v(k1+lb:k1+gub)

               k1 = k1 + n

               flops = flops + 2*nl

            End Do

            w(lb:gub) = z(lb:gub) - w(lb:gub)

            rm(j+1,j) = Sqrt( nrm2( nl, w(lb) ) )

            If ( rm(j+1,j) == 0.0_l_ ) Then

               Print *, ' Stop: r(j+1,j) = 0'; Stop

            End If

            k0 = k0 + n

            w(lb:gub) = w(lb:gub)/rm(j+1,j)

            v(k0+lb:k0+gub) = w(lb:gub)

            Do k = 1, j - 1

               Call rotf( rho(k), alpha, eta )

               tau1 = rm(k,j)

               tau2 = rm(k+1,j)

               rm(k,j)   = alpha*tau1 - eta*tau2

               rm(k+1,j) = eta*tau1   + alpha*tau2

               flops = flops + 4

            End Do

            Call givens( rm(j,j), rm(j+1,j), alpha, eta )

            tau1 = rm(j,j)

            tau2 = rm(j+1,j)

            rm(j,j)   = alpha*tau1 - eta*tau2

            rm(j+1,j) = eta*tau1   + alpha*tau2

            flops = flops + 4

            Call rotb( rho(j), alpha, eta )

            g(j)   = g(j)*alpha

            g(j+1) = g(j+1)*eta

            g(j+2) = g(j+1)

!

! --- Test for convergence.

!

            exnrm = Abs( g(j+1) )

            conv = exnrm < rhstp

            If ( conv ) Exit

         End Do                                      ! <--- End j-loop

         Call lsqslv( jc, ir+2, rm, g )      

         k1 = 0

         Do k = 1, jc

            x(lb:gub) = x(lb:gub) + g(k)*v(k1+lb:k1+gub)

            k1 = k1 + n

            flops = flops + 2*nl

         End Do

         If ( conv ) Return                          ! <--- Convergence

         Call spmxv( nl, nel, indx, rowp, matvals, x, w(lb) )

         z(lb:gub) = b(lb:gub) - w(lb:gub)

         Call prec( nl, nel, mp, indx, rowp, matvals, z, r, gamma )

         beta = Sqrt( nrm2( nl, r(lb) ) )

      End Do                                         ! <--- End i-loop

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

! --- Normally we would end up here and with no convergence issue a

!     warning. Because we only want to see the residual value and

!     the speed for benchmarking purposes, we comment out the following

!     lines:

!

!     If ( iter >= maxit ) Then

!          iter = maxit

!          Print *, 'No convergence in ', maxit, ' iterations;'

!          Print *, 'Norm of residual =', exnrm

!          Stop

!     End If

      flops = flops + nl + 9 + iter*( 2*nl + 9 ) ! <--- # of rest flops.

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

      End Subroutine rgmres