tfqmr.f

网络带宽测试工具

      Subroutine tfqmr( 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.
!     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, im, im0, iter
      Real(l_) :: alpha, beta, c, eta, eta0, kappa, rho, rho0, rhstp,
     &            sigma, tau, tau0, theta, theta0
      Real(l_) :: d(n), g(n), h(n), p(n), r(n), rt(n), v(n), w(n),
     &            y(n), y0(n), z(n)
      Real(l_) :: dotpr, nrm2
      External    dotpr, nrm2
! ----------------------------------------------------------------------
      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 )
!$omp parallel do
      Do il = 1, n
         w(il) = r(il)
         y(il) = r(il)
      End Do
      Call spmxv( n, nel, indx, rowp, matvals, y, z )
      Call prec( n, nel, m, indx, rowp, matvals, z, g, gamma )
!$omp parallel do
      Do il = 1, n
         v(il) = g(il)
         d(il) = 0.0_l_
      End Do
      tau = Sqrt( nrm2( n, r ) )
      theta = 0.0_l_
      eta   = 0.0_l_
!$omp parallel do
      Do il = 1, n
         rt(il) = r(il)
      End Do
!
! --- We can use 'nrm2' instead of dotpr because rt = r.
!
      rho = nrm2( n, r )
      rhstp = tol*Sqrt( nrm2( n, b ) )
      im0 = 1
      Do i = 1, maxit
         iter = i
         sigma = dotpr( n, rt, v )
         If ( sigma == 0.0_l_ ) Then
            Print *, 'Stop: sigma = 0'; Stop
         End If
         alpha = rho/sigma
!$omp parallel do
         Do il = 1, n
            y0(il) = y(il)
            y(il)  = y(il) - alpha*v(il)
         End Do
         Call spmxv( n, nel, indx, rowp, matvals, y, z )
         Call prec( n, nel, m, indx, rowp, matvals, z, h, gamma )
         Do im = im0, im0 + 1
!$omp parallel do
            Do il = 1, n
               w(il) = w(il) - alpha*g(il)
            End Do
            theta0 = theta
            tau0   = tau
            If ( tau0 == 0.0_l_ ) Then
               Print *, 'Stop: tau0 = 0'; Stop
            End If
            theta = Sqrt( nrm2( n, w ) )/tau
            c = 1.0_l_/Sqrt( 1.0_l_ + theta*theta )
            tau  = tau0*theta*c
            eta0 = eta
            eta  = c*c*alpha
            If ( alpha == 0.0_l_ ) Then
               Print *, 'Stop: alpha = 0'; Stop
            End If
!$omp parallel do
            Do il = 1, n
               d(il) = y0(il) + (theta0*theta0*eta0/alpha)*d(il)
               x(il) = x(il) + eta*d(il)
            End Do
            exnrm = Sqrt( nrm2( n, r ) )
            kappa = Sqrt( Real( im + 1, l_ ) )*tau
            If ( kappa < tol ) Then
               Call spmxv( n, nel, indx, rowp, matvals, x, p )
!$omp parallel do
               Do il = 1, n
                  z(il) = b(il) - p(il)
               End Do
               Call prec( n, nel, m, indx, rowp, matvals, z, r, gamma )
               exnrm = Sqrt( nrm2( n, r ) )
               If ( exnrm < rhstp ) Go To 10  ! <--- Convergence.
               flops = flops + n + 9
            End If
!$omp parallel do
            Do il = 1, n
               y0(il) = y(il)
               g(il)  = h(il)
            End Do
         End Do                               ! <--- im-loop
         rho0 = rho
         rho  = dotpr( n, rt, w )
         If ( rho0 == 0.0_l_ ) Then
            Print *, 'Stop: rho0 = 0'; Stop
         End If
         beta = rho/rho0
!$omp parallel do
         Do il = 1, n
            y(il) = w(il) + beta*y0(il)
         End Do
         Call spmxv( n, nel, indx, rowp, matvals, y, z )
         Call prec( n, nel, m, indx, rowp, matvals, z, g, gamma )
!$omp parallel do
         Do il = 1, n
            v(il) = g(il) + beta*( h(il) + beta*v(il) )
         End Do
         im0 = im0 + 2
      End Do                                 ! <--- 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
   10 flops = flops + n + 10 + iter*( 22*n + 79 ) ! <--- # of flops.
! ----------------------------------------------------------------------
      End Subroutine tfqmr