PRACE NPT

      Subroutine rgmres( n, nel, m, indx, rowp, matvals, q, x, b, gamma,

     &                  maxit, tol, exnrm, prec )

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

! --- tfqmr does an iterative solve of Ax = b, where A is in CRS format:

!     Integer  indx(nel),

!     Integer  rowp(n+1), and

!     Real(l_) matvals(nel).

!     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

      Implicit               None

      Integer             :: n, nel, m, maxit

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

      Real(l_)            :: matvals(nel)

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

      Real(l_)            :: gamma(m+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( n, b ) )

      Call spmxv( n, nel, indx, rowp, matvals, x, w )

      z = b - w

      Call prec( n, nel, m, indx, rowp, matvals, z, r, gamma )

      beta = Sqrt( nrm2( n, r ) )

      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(1:n) = r/beta

         k0 = 0

         Do j = 1, ib

            jc = j

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

            Call prec( n, nel, m, indx, rowp, matvals, w, z, gamma )

            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 = w + rm(k,j)*v(k1+1:k1+n)

               k1 = k1 + n

               flops = flops + 2*n

            End Do

            w = z - w

            rm(j+1,j) = Sqrt( nrm2( n, w ) )

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

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

            End If

            k0 = k0 + n

            v(k0+1:k0+n) = w/rm(j+1,j)

            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 k-loop

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

         k1 = 0

         Do k = 1, jc

            x = x + g(k)*v(k1+1:k1+n)

            k1 = k1 + n

            flops = flops + 2*n

         End Do

         If ( conv ) Return                          ! <--- Convergence

         Call spmxv( n, nel, indx, rowp, matvals, x, w )

         z = b - w

         Call prec( n, nel, m, indx, rowp, matvals, z, r, gamma )

         beta = Sqrt( nrm2( n, r ) )

      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 + n + 9 + iter*( 2*n + 9 ) ! <--- # of rest flops.

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

      End Subroutine rgmres