mod2cr.f

网络带宽测试工具

      Program mod2cr
! ----------------------------------------------------------------------
! **********************************************************************
! *** This program is part of the Euroben Benchmark                  ***
! ***                                                                ***
! *** Copyright: EuroBen Group p/o                                   ***
! ***            Utrecht University, Computational Physics Dept.     ***
! ***            P.O. Box 80.000                                     ***
! ***            3508 TA Utrecht                                     ***
! ***            The Netherlands                                     ***
! ***                                                                ***
! *** Author of this program: Aad van der Steen                      ***
! *** Date                    Spring 2005                            ***
! **********************************************************************
! ----------------------------------------------------------------------
! --- Version 1.0 (Shared-memory parallel, OpenMP).

! --- Purpose of program mod2cr
! -----------------------------
! --- Solver for sparse linear systems with iterative methods and
!     two types of preconditioners. The systems are not actually
!     solved but rather a predefined number of iterations is performed
!     to assess the speed of the solver-preconditioner combination.
!     Two types of systems are considered:
!     1. - Systems stemming from 3-D finite difference problems
!          resulting in 7-banded matrices. We consider only the 
!          symmetric type in this program 'mod2cr'.
!     2. - Finite element type, irregularly filled matrices.
!          They are stored in CRS format. These systems are addressed
!          in ANOTHER program: 'mod2ci'. 
!
! --- In this Program 'mod2cr' we address the regular, banded type
!     of problem. For the irregular CRS matrix-based type we have
!     program 'mod2ci'.
!
!     Solvers used:
!     -  For the symmetric banded systems in 'mod2cr': CG with ILU(0)
!        and polynomial preconditioning.
!     -  For irregular systems in 'mod2ci': RGMRES and TFQMR with
!        polynomial preconditioning.
! ----------------------------------------------------------------------
      Use                       numerics
      Use                       floptime
      Implicit                  None

      Real(l_), Allocatable  :: areg(:,:), aorg(:,:), b(:), borg(:),
     &                          conv(:), fac(:), x(:), q(:)
      Real(l_)               :: corr, maxdf, mindf, res
      Real(l_)               :: time, mflops
      Integer                :: maxit = 50
      Integer, Parameter     :: m = 1
      Real(l_)               :: gamma(m+1)
      Integer                :: nprocs, omp_get_num_threads
      Real(l_), Parameter    :: tol = 1.0e-12_l_
      Real(l_), Parameter    :: micro = 1.0e-6_l_, two = 2.0_l_,
     &                          perc = 1.0e2_l_
      Integer                :: n1, n2, n3, ntot, nrep
      Integer                :: i
      Logical                :: ok
      External               :: sym7icg0, sym7pol
! ----------------------------------------------------------------------
      Call state( 'mod2cr  ' )
      Open( 1, File = 'mod2cr.in' )
!$omp parallel shared(nprocs)
      nprocs = omp_get_num_threads()
!$omp end parallel
      Print 1000, nprocs
! ----------------------------------------------------------------------
! --- Do 7-band case with ILU(0) preconditioning.
 
      Print 1005
      Call pcoefs( m, gamma )
   10 Read( 1, *, End = 20 ) n1, n2, n3, nrep
         ntot = n1*n2*n3
         Allocate( areg(ntot,0:3), aorg(ntot,0:3), b(ntot), borg(ntot),
     &             conv(maxit), fac(ntot), x(ntot), q(ntot) ) 
         Call sym7gen( n1, n2, n3, aorg, borg )
         flops = 0
         q = 1.0_l_
         time = wclock()
         Do i = 1, nrep
            Call zero( ntot, x )
            areg = aorg
            b    = borg
            Call sym7pre( n1, n2, n3, areg, b, fac )
            Call cg( n1, n2, n3, m, areg, b, x, q, gamma, maxit, conv,
     &               tol, sym7icg0 )
            x = x*fac
            flops = flops + ntot
         End Do
         time   = wclock() - time
         corr = wclock()
         Do i = 1, nrep
            Call zero( ntot, x )
            areg = aorg
            b    = borg
         End Do
         time = time - wclock() + corr
         mflops = micro*Real( flops, l_ )/time
         time   = time/Real( nrep, l_ )
         res = conv(maxit)
         ok = .TRUE.
! ----------------------------------------------------------------------
! --- Routine 'checkreg' is only used for debugging but has no function
!     in the finished program.
!        Call checkreg( n1, n2, n3, areg, b, x, mindf, maxdf )

         Deallocate ( areg, aorg, b, borg, conv, fac, x, q )
         Print 1010, n1, n2, n3, time, mflops, res
      Go To 10
   20 Print 1020
      Print 1025
! ----------------------------------------------------------------------
! --- Do 7-band case with von Neumann polynomial preconditioning.

      Rewind 1
      Print 1030, nprocs
      Print 1005
      Call pcoefs( m, gamma )
   30 Read ( 1, *, End = 40 ) n1, n2, n3, nrep
         Allocate( areg(ntot,0:3), aorg(ntot,0:3), b(ntot), borg(ntot),
     &             conv(maxit), fac(ntot), x(ntot), q(ntot) )
         Call sym7gen( n1, n2, n3, aorg, borg )
         flops = 0
         q = 1.0_l_
         time = wclock()
         Do i = 1, nrep
            Call zero( ntot, x )
            areg = aorg
            b    = borg
            Call sym7pre( n1, n2, n3, areg, b, fac )
            Call cg( n1, n2, n3, m, areg, b, x, q, gamma, maxit, conv,
     &               tol, sym7pol )
            x = x*fac
            flops = flops + ntot
         End Do
         time = wclock() - time
         corr = wclock()
         Do i = 1, nrep
            Call zero( ntot, x )
            areg = aorg
            b    = borg
         End Do
         time = time - wclock() + corr
         mflops = micro*Real( flops, l_ )/time
         time   = time/Real( nrep, l_ )
         res = conv(maxit)
         ok = .TRUE.
! ----------------------------------------------------------------------
! --- Routine 'checkreg' is only used for debugging but has no function
!     in the finished program.
!        Call checkreg( n1, n2, n3, areg, b, x, mindf, maxdf )

         Deallocate ( aorg, areg, b, borg, conv, fac, x, q )
         Print 1010, n1, n2, n3, time, mflops, res
      Go To 30
   40 Print 1020
! ----------------------------------------------------------------------
 1000 Format( 'Program mod2cr: Sparse iterative solver test'/
     &        'Symmetric, 7-band: CG with ILU(0) preconditioner.'/
     &        'No. of proc.s = ', i4 )
 1005 Format( '-------------------------------------------------------',
     &        '-----'/
     &        '  n1 |  n2 |  n3 |   Time(s)   |   Mflop/s   |   Residu',
     &        'e   |'/
     &        '-------------------------------------------------------',
     &        '----|' )
 1010 Format( i5, '|', i5, '|', i5, '|', g13.5, '|', g13.5, '|', g13.5,
     &        '|' )
 1020 Format( '-------------------------------------------------------', 
     &        '-----' )
 1025 Format ( / )
 1030 Format( '-------------------------------------------------------',
     &        '-----'/
     &        'Program mod2cr: Sparse iterative solver test'/
     &        'Symmetric, 7-band: CG with polynomial preconditioner.'/
     &        'No. of proc.s = ', i4 )
! ----------------------------------------------------------------------
      End Program mod2cr