cg.f

网络带宽测试工具

      Subroutine cg( n1, n2, n3, m, a, b, x, q, gamma, maxit, conv,
     &               tol, prec )
! ---------------------------------------------------------------------
! --- cg does an iterative solve of ax = b.
!     Real(l_) a(n1*n2*n3,0:3) : 7-diagonal symmetric matrix of which
!                                only the 4 upper diagonals are stored.
!     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_) q(n1*n2*n3)     : Contains (left) preconditioning vector
!     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
      Implicit    None

      Integer  :: n1, n2, n3, m
      Real(l_) :: a(n1*n2*n3,0:3), b(n1*n2*n3), x(n1*n2*n3), q(n1*n2*n3)
      Real(l_) :: gamma(m+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, a, x, r )
!$omp parallel do
      Do j = 1, ntot
         r(j) = b(j) - r(j)
      End Do
! ---------------------------------------------------------------------
! --- Precondition of initial residual.

      Call prec( n1, n2, n3, m, a, nr0, q, r, p, gamma )
      nrm = nr0
! ---------------------------------------------------------------------
! --- Iterate at most maxit times.  

      its = maxit
      Do i = 1, its
         Call sym7mxv( n1, n2, n3, a, p, ap )
         pap = dotpr( ntot, p, ap )
         alpha = nrm/pap

!$omp parallel do
         Do j = 1, ntot
            x(j) = x(j) + alpha*p(j)
            r(j) = r(j) - alpha*ap(j)
         End Do

         Call prec( n1, n2, n3, m, a, nnrm, q, r, ap, gamma )
         conv(i) = Sqrt( nnrm )
!     PRINT *, 'i = ', i, '; norm = ', conv(i)
         If ( nnrm < nr0*tol2 ) Then
            maxit = i
            Go To 1000
         End If
         beta = nnrm/nrm
         nrm = nnrm
!$omp parallel do
         Do j = 1, ntot
            p(j) = ap(j) + beta*p(j)
         End Do
      End Do
!     Print *, 'CG: no convergence in ', maxit, ' iterations.'
 1000 flops = flops + maxit*( 8*ntot + 4 ) + ntot + 1
! ---------------------------------------------------------------------
      End Subroutine cg