!-------------------------------------------------------------------------!
!                                                                         !
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3.1       !
!                                                                         !
!                       D V M H     V E R S I O N                         !
!                                                                         !
!                                   S P                                   !
!                                                                         !
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!

c---------------------------------------------------------------------
c
c Authors:
c      Original:
c          R. Van der Wijngaart
c          W. Saphir
c          H. Jin
c      Optimize for DVMH:
c          Kolganov A.S.
c---------------------------------------------------------------------

c---------------------------------------------------------------------
       program SP
c---------------------------------------------------------------------

       include  'header.h'
      
       integer          i, niter, step, fstatus, n3
       external         timer_read
       double precision mflops, t, tmax, timer_read, trecs(t_last)
       logical          verified
       character        class
       character        t_names(t_last)*8

c---------------------------------------------------------------------
c      Read input file (if it exists), else take
c      defaults from parameters
c---------------------------------------------------------------------
          
       open (unit=2,file='timer.flag',status='old', iostat=fstatus)
       if (fstatus .eq. 0) then
         timeron = .true.
         t_names(t_total) = 'total'
         t_names(t_rhsx) = 'rhsx'
         t_names(t_rhsy) = 'rhsy'
         t_names(t_rhsz) = 'rhsz'
         t_names(t_rhs) = 'rhs'
         t_names(t_xsolve) = 'xsolve'
         t_names(t_ysolve) = 'ysolve'
         t_names(t_zsolve) = 'zsolve'
         t_names(t_rdis1) = 'redist1'
         t_names(t_rdis2) = 'redist2'
         t_names(t_tzetar) = 'tzetar'
         t_names(t_ninvr) = 'ninvr'
         t_names(t_pinvr) = 'pinvr'
         t_names(t_txinvr) = 'txinvr'
         t_names(t_add) = 'add'
         close(2)
       else
         timeron = .false.
       endif

       write(*, 1000)
       open (unit=2,file='inputsp.data',status='old', iostat=fstatus)

       if (fstatus .eq. 0) then
         write(*,233) 
 233     format(' Reading from input file inputsp.data')
         read (2,*) niter
         read (2,*) dt
         read (2,*) grid_points(1), grid_points(2), grid_points(3)
         close(2)
       else
         write(*,234) 
         niter = niter_default
         dt    = dt_default
         grid_points(1) = problem_size
         grid_points(2) = problem_size
         grid_points(3) = problem_size
       endif
 234   format(' No input file inputsp.data. Using compiled defaults')
       open (unit = 2,file = 'inputStage',status = 'old',iostat = fstat
     &us)	
      if (fstatus .eq. 0) then
         read (unit = 2,fmt = *) stage_n
         close (unit = 2)
      else  
         stage_n = 0
      endif 
      write(*,*) 'stage = ', stage_n
	  
       write(*, 1001) problem_size, problem_size, problem_size
       write(*, 1002) niter, dt
       write(*, *)

 1000  format(//, ' NAS Parallel Benchmarks (NPB3.3.1-DVMH)',
     >            ' - SP Benchmark', /)
 1001  format(' Size: ', i4, 'x', i4, 'x', i4)
 1002  format(' Iterations: ', i4, '    dt:  ', F11.7)
 1003  format(' Number of available threads: ', i5)

       if ( (problem_size .gt. IMAX) .or.
     >      (problem_size .gt. JMAX) .or.
     >      (problem_size .gt. KMAX) ) then
             print *, (grid_points(i),i=1,3)
             print *,' Problem size too big for compiled array sizes'
             goto 999
       endif
       nx2 = problem_size - 2
       ny2 = problem_size - 2
       nz2 = problem_size - 2

       call set_constants
       call exact_rhs

       call initialize
       call adi_first
       call adi_first
       call initialize	   
	   
       do i = 1, t_last
          call timer_clear(i)
       end do
       call timer_start(1)
!DVM$ BARRIER
       do  step = 1, niter

          if (mod(step, 20) .eq. 0 .or. step .eq. 1) then
             write(*, 200) step
 200         format(' Time step ', i4)
          endif

          call adi

       end do
       call timer_stop(1)
       tmax = timer_read(1)

       call verify(niter, class, verified)

       if( tmax .ne. 0. ) then
          n3 = problem_size*problem_size*problem_size
          t = (problem_size+problem_size+problem_size)/3.0
          mflops = (881.174 * float( n3 )
     >             -4683.91 * t**2
     >             +11484.5 * t
     >             -19272.4) * float( niter ) / (tmax*1000000.0d0)
       else
          mflops = 0.0
       endif

      call print_results('SP', class, problem_size, 
     >     problem_size, problem_size, niter, 
     >     tmax, mflops, '          floating point', 
     >     verified, npbversion,compiletime, cs1, cs2, cs3, cs4, cs5, 
     >     cs6, '(none)')

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

       do i=1, t_last
          trecs(i) = timer_read(i)
       end do
       if (tmax .eq. 0.0) tmax = 1.0

       write(*,800)
 800   format('  SECTION   Time (secs)')

       do i=1, t_last
          write(*,810) t_names(i), trecs(i), trecs(i)*100./tmax
          if (i.eq.t_rhs) then
             t = trecs(t_rhsx) + trecs(t_rhsy) + trecs(t_rhsz)
             write(*,820) 'sub-rhs', t, t*100./tmax
             t = trecs(t_rhs) - t
             write(*,820) 'rest-rhs', t, t*100./tmax
          elseif (i.eq.t_zsolve) then
             t = trecs(t_zsolve) - trecs(t_rdis1) - trecs(t_rdis2)
             write(*,820) 'sub-zsol', t, t*100./tmax
          elseif (i.eq.t_rdis2) then
             t = trecs(t_rdis1) + trecs(t_rdis2)
             write(*,820) 'redist', 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
	   
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine  adi_first

c---------------------------------------------------------------------
c---------------------------------------------------------------------

       call compute_rhs(1)
       call x_solve
       call y_solve
       call z_solve

       return
       end
	   
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine  adi

c---------------------------------------------------------------------
c---------------------------------------------------------------------
!DVM$ interval 1
       call compute_rhs(1)
!DVM$ end interval	   
!DVM$ interval 12	   
       call x_solve
!DVM$ end interval	   
!DVM$ interval 13	   
       call y_solve
!DVM$ end interval	   	   
!DVM$ interval 14	   
       call z_solve
!DVM$ end interval	   

       return
       end	   