!-------------------------------------------------------------------------!
!                                                                         !
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3         !
!                                                                         !
!                       D V M     V E R S I O N                           !
!                                                                         !
!                                   M G                                   !
!                                                                         !
!-------------------------------------------------------------------------!
!                                                                         !
!    This benchmark is an OpenMP version of the NPB MG code.              !
!    It is described in NAS Technical Report 99-011.                      !
!                                                                         !
!    Permission to use, copy, distribute and modify this software         !
!    for any purpose with or without fee is hereby granted.  We           !
!    request, however, that all derived work reference the NAS            !
!    Parallel Benchmarks 3.3. This software is provided "as is"           !
!    without express or implied warranty.                                 !
!                                                                         !
!    Information on NPB 3.3, including the technical report, the          !
!    original specifications, source code, results and information        !
!    on how to submit new results, is available at:                       !
!                                                                         !
!           http://www.nas.nasa.gov/Software/NPB/                         !
!                                                                         !
!    Send comments or suggestions to  npb@nas.nasa.gov                    !
!                                                                         !
!          NAS Parallel Benchmarks Group                                  !
!          NASA Ames Research Center                                      !
!          Mail Stop: T27A-1                                              !
!          Moffett Field, CA   94035-1000                                 !
!                                                                         !
!          E-mail:  npb@nas.nasa.gov                                      !
!          Fax:     (650) 604-3957                                        !
!                                                                         !
!-------------------------------------------------------------------------!


c---------------------------------------------------------------------
c
c Authors: 
c      Original:
c          E. Barszcz
c          P. Frederickson
c          A. Woo
c          M. Yarrow
c          H. Jin
c      DVM/DVMH vesion:  
c          A. Shubert
c      Optimized for DVM/DVMH:
c          A. Kolganov
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c Программа решающая уравнение Пуасона многосеточным методом (V-цикл)
c---------------------------------------------------------------------
      program mg	  
c---------------------------------------------------------------------
      implicit none

      include 'globals.h'
      include 'dvmvars.h'   
      
c---------------------------------------------------------------------------c
c k is the current level. It is passed down through subroutine args
c and is NOT global. it is the current iteration
c---------------------------------------------------------------------------c

      integer k, it, pdim, pi
      
      external timer_read
      double precision t, tinit, mflops, timer_read

c---------------------------------------------------------------------------c
c These arrays are in common because they are quite large
c and probably shouldn't be allocated on the stack. They
c are always passed as subroutine args. 
c---------------------------------------------------------------------------c      

      double precision rnm2, rnmu, epsilon
      integer n1, n2, n3, nit
      double precision nn, verify_value, err, t_1,t_2
      logical verified

      integer i, fstatus
      character t_names(t_last)*8
      double precision tmax, elapsed_time
	  
      do i = T_init, T_last
         call timer_clear(i)
      end do

c---------------------------------------------------------------------
c Read in and broadcast input data
c---------------------------------------------------------------------

      open(unit=7,file='timer.flag', status='old', iostat=fstatus)
      if (fstatus .eq. 0) then
         timeron = .true.
         t_names(t_init) = 'init'
         t_names(t_bench) = 'benchmk'
         t_names(t_mg3P) = 'mg3P'
         t_names(t_psinv) = 'psinv'
         t_names(t_resid) = 'resid'
         t_names(t_rprj3) = 'rprj3'
         t_names(t_interp) = 'interp'
         t_names(t_norm2) = 'norm2'
         t_names(t_comm3) = 'comm3'
         close(7)
      else
         timeron = .false.
      endif

      write (*, 1000) 

      open(unit=7,file="mg.input", status="old", iostat=fstatus)
      if (fstatus .eq. 0) then
         write(*,50) 
 50      format(' Reading from input file mg.input')
         read(7,*) lt
         read(7,*) nx(lt), ny(lt), nz(lt)
         read(7,*) nit
         read(7,*) (debug_vec(i),i=0,7)
      else
         write(*,51) 
 51      format(' No input file. Using compiled defaults ')
         lt = lt_default
         nit = nit_default
         nx(lt) = nx_default
         ny(lt) = ny_default
         nz(lt) = nz_default
         do i = 0,7
            debug_vec(i) = debug_default
         end do
      endif


      if ( (nx(lt) .ne. ny(lt)) .or. (nx(lt) .ne. nz(lt)) ) then
         Class = 'U' 
      else if( nx(lt) .eq. 32 .and. nit .eq. 4 ) then
         Class = 'S'
      else if( nx(lt) .eq. 128 .and. nit .eq. 4 ) then
         Class = 'W'
      else if( nx(lt) .eq. 256 .and. nit .eq. 4 ) then  
         Class = 'A'
      else if( nx(lt) .eq. 256 .and. nit .eq. 20 ) then
         Class = 'B'
      else if( nx(lt) .eq. 512 .and. nit .eq. 20 ) then  
         Class = 'C'
      else if( nx(lt) .eq. 1024 .and. nit .eq. 50 ) then  
         Class = 'D'
      else if( nx(lt) .eq. 2048 .and. nit .eq. 50 ) then  
         Class = 'E'
      else
         Class = 'U'
      endif

c---------------------------------------------------------------------
c  Use these for debug info:
c---------------------------------------------------------------------
c     debug_vec(0) = 1 !=> report all norms
c     debug_vec(1) = 1 !=> some setup information
c     debug_vec(1) = 2 !=> more setup information
c     debug_vec(2) = k => at level k or below, show result of resid
c     debug_vec(3) = k => at level k or below, show result of psinv
c     debug_vec(4) = k => at level k or below, show result of rprj
c     debug_vec(5) = k => at level k or below, show result of interp
c     debug_vec(6) = 1 => (unused)
c     debug_vec(7) = 1 => (unused)
c---------------------------------------------------------------------
      a(0) = -8.0D0/3.0D0 
      a(1) =  0.0D0 
      a(2) =  1.0D0/6.0D0 
      a(3) =  1.0D0/12.0D0
      
      if(Class .eq. 'A' .or. Class .eq. 'S'.or. Class .eq.'W') then
c---------------------------------------------------------------------
c     Coefficients for the S(a) smoother
c---------------------------------------------------------------------
         c(0) =  -3.0D0/8.0D0
         c(1) =  +1.0D0/32.0D0
         c(2) =  -1.0D0/64.0D0
         c(3) =   0.0D0
      else
c---------------------------------------------------------------------
c     Coefficients for the S(b) smoother
c---------------------------------------------------------------------
         c(0) =  -3.0D0/17.0D0
         c(1) =  +1.0D0/33.0D0
         c(2) =  -1.0D0/61.0D0
         c(3) =   0.0D0
      endif
      lb = 1
      k  = lt
      
c**********************************************************************
c**********************************************************************
c********************* START HERE *************************************
c**********************************************************************
c**********************************************************************
       
      call setup(n1,n2,n3,k)
      call setupDVM(0)

      call timer_start(T_init)

      call zero3(u(pu),n1,n2,n3)
      call zran3(v(pv),n1,n2,n3,nx(lt),ny(lt),k,class)    

      write (*, 1001) nx(lt),ny(lt),nz(lt), Class
      write (*, 1002) nit
      write (*, *)

      pdim = PROCESSORS_RANK()
      write (*, 310) pdim
      do pi=1, pdim
          write (*, 311) pi, PROCESSORS_SIZE(pi)
      enddo           
      write (*, *)  ' '
      
 310  format(' Processors grid rank: ', i4)
 311  format(' Grid dimension [', i4, '] size: ', i4)
 
 1000 format(//,' NAS Parallel Benchmarks (NPB3.3-DVMH)',
     >          ' - MG Benchmark', /)
 1001 format(' Size: ', i4, 'x', i4, 'x', i4, '  (class ', A, ')' )
 1002 format(' Iterations:                  ', i5)

      call resid(u(pu),v(pv),r(pr),n1,n2,n3,a,k)
      
c---------------------------------------------------------------------
c     One iteration for startup
c---------------------------------------------------------------------
      call mg3P(n1,n2,n3,k)      
      call resid(u(pu),v(pv),r(pr),n1,n2,n3,a,k)
      call norm2u3(r(pr),n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))         
      call setup(n1,n2,n3,k)            
      
      call zero3(u(pu),n1,n2,n3)
      call zran3(v(pv),n1,n2,n3,nx(lt),ny(lt),k,class)
      
      call timer_stop(T_init)
      tinit = timer_read(T_init)

      write( *,'(A,F15.3,A/)' ) 
     >     ' Initialization time: ',tinit, ' seconds'

      do i = T_bench, T_last
         call timer_clear(i)
      end do

      call timer_start(T_bench)      
      call resid(u(pu),v(pv),r(pr),n1,n2,n3,a,k)      
      call norm2u3(r(pr),n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))      
   
! временный вызов
!      call timer_start(T_bench)  
c*****************************************************************
c********************MAIN LOOP ***********************************
c*****************************************************************

      do  it = 1, nit
         if (it.eq.1 .or. it.eq.nit .or. mod(it,5).eq.0) then
            write(*,80) it
   80       format('  iter ',i3)
         endif
         call mg3P(n1,n2,n3,k)
         call resid(u(pu),v(pv),r(pr),n1,n2,n3,a,k)
      enddo

      call norm2u3(r(pr),n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))   
      
      call timer_stop(T_bench)

      t = timer_read(T_bench)

      verified = .FALSE.
      verify_value = 0.0

      write(*,100)
 100  format(/' Benchmark completed ')

      epsilon = 1.d-8
      if (Class .ne. 'U') then
         if(Class.eq.'S') then
            verify_value = 0.5307707005734d-04
         elseif(Class.eq.'W') then
            verify_value = 0.6467329375339d-05
         elseif(Class.eq.'A') then
            verify_value = 0.2433365309069d-05
         elseif(Class.eq.'B') then
            verify_value = 0.1800564401355d-05
         elseif(Class.eq.'C') then
            verify_value = 0.5706732285740d-06
         elseif(Class.eq.'D') then
            verify_value = 0.1583275060440d-09
         elseif(Class.eq.'E') then
            verify_value = 0.5630442584711d-10
         endif

         err = abs( rnm2 - verify_value ) / verify_value
         if( err .le. epsilon ) then
            verified = .TRUE.
            write(*, 200)
            write(*, 201) rnm2
            write(*, 202) err
 200        format(' VERIFICATION SUCCESSFUL ')
 201        format(' L2 Norm is ', E20.13)
 202        format(' Error is   ', E20.13)
         else
            verified = .FALSE.
            write(*, 300) 
            write(*, 301) rnm2
            write(*, 302) verify_value
 300        format(' VERIFICATION FAILED')
 301        format(' L2 Norm is             ', E20.13)
 302        format(' The correct L2 Norm is ', E20.13)
         endif
      else
         verified = .FALSE.
         write (*, 400)
         write (*, 401)
         write (*, 201) rnm2
 400     format(' Problem size unknown')
 401     format(' NO VERIFICATION PERFORMED')
      endif

      nn = 1.0d0*nx(lt)*ny(lt)*nz(lt)

      if( t .ne. 0. ) then
         mflops = 58.*nit*nn*1.0D-6 /t
      else
         mflops = 0.0
      endif

      call print_results('MG', class, nx(lt), ny(lt), nz(lt), 
     >                   nit, t,
     >                   mflops, '          floating point', 
     >                   verified, npbversion, compiletime,
     >                   cs1, cs2, cs3, cs4, cs5, cs6, cs7)


 600  format( i4, 2e19.12)

c---------------------------------------------------------------------
c      More timers
c---------------------------------------------------------------------
      if (.not.timeron) goto 999

      tmax = timer_read(t_bench)
      if (tmax .eq. 0.0) tmax = 1.0

      write(*,800)
 800  format('  SECTION   Time (secs)')
      do i=t_bench, t_last
         t = timer_read(i)
         if (i.eq.t_resid2) then
            t = timer_read(T_resid) - t
            write(*,820) 'mg-resid', t, t*100./tmax
         else
            write(*,810) t_names(i), t, t*100./tmax
         endif
 810     format(2x,a8,':',f9.3,'  (',f6.2,'%)')
 820     format('    --> ',a8,':',f9.3,'  (',f6.2,'%)')
      end do

 999  continue

      end
