rgmres.f

网络带宽测试工具

      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, il, 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 )
!$omp parallel do
      Do il = 1, n
         z(il) = b(il) - w(il)
      End Do
      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
!$omp parallel do
         Do il = 1, n
            v(il) = r(il)/beta
         End Do
         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
!$omp parallel do
            Do il = 1, n
               w(il) = 0.0_l_
            End Do
            Do k = 1, j
!$omp parallel do
               Do il = 1, n
                  w(il) = w(il) + rm(k,j)*v(k1+il)
               End Do
               k1 = k1 + n
               flops = flops + 2*n
            End Do
!$omp parallel do
            Do il = 1, n
               w(il) = z(il) - w(il)
            End Do
            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
!$omp parallel do
            Do il = 1, n
               w(il) = w(il)/rm(j+1,j)
               v(k0+il) = w(il)
            End Do
            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
!$omp parallel do
            Do il = 1, n
               x(il) = x(il) + g(k)*v(k1+il)
            End Do
            k1 = k1 + n
            flops = flops + 2*n
         End Do
         If ( conv ) Return                          ! <--- Convergence
         Call spmxv( n, nel, indx, rowp, matvals, x, w )
!$omp parallel do
         Do il = 1, n
            z(il) = b(il) - w(il)
         End Do
         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