 
! *** generated by SAPFOR with version 1756 and build date: Mar 23 2021 12:41:48 
 
!-------------------------------------------------------------------------!   
!                                                                         !   
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3         !   
!                                                                         !   
!                                   E P                                   !   
!                                                                         !   
!-------------------------------------------------------------------------!   
!                                                                         !   
!    This benchmark is part of the NAS Parallel Benchmark 3.3 suite.      !   
!    It is described in NAS Technical Reports 95-020 and 02-007           !   
!                                                                         !   
!    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                                        !   
!                                                                         !   
!-------------------------------------------------------------------------!   
!---------------------------------------------------------------------   
!        
! Authors: P. O. Frederickson    
!          D. H. Bailey   
!          A. C. Woo   
!          R. F. Van der Wijngaart   
!---------------------------------------------------------------------   
!---------------------------------------------------------------------   
      program embar 
 
!---------------------------------------------------------------------   
!        
!   This is the MPI version of the APP Benchmark 1,   
!   the "embarassingly parallel" benchmark.   
!        
!        
!   M is the Log_2 of the number of complex pairs of uniform (0, 1) random   
!   numbers.  MK is the Log_2 of the size of each batch of uniform random   
!   numbers.  MK can be set for convenience on a given system, since it does   
!   not affect the results.   
!---------------------------------------------------------------------  
!---------------------------------------------------------------------  
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
!--------------------------------------------------------------------- 
      implicit none 
 
      include 'mpinpb.h'   
      include 'npbparams.h'   
      double precision ::  mops,epsilon,a,s,t1,t2,t3,t4,x,x1,x2,q,sx,sy, 
     &tm,an,tt,gc,dum(3),timer_read 
      double precision ::  sx_verify_value,sy_verify_value,sx_err,sy_err 
      integer ::  mk,mm,nn,nk,nq,np,ierr,node,no_nodes,i,ik,kk,l,k,nit,i 
     &errcode,no_large_nodes,np_add,k_offset,j 
      logical ::  verified,timers_enabled 
      external randlc,timer_read 
      double precision ::  randlc,qq 
      character*15 ::  size 
      integer ::  fstatus 
      integer ::  t_total,t_gpairs,t_randn,t_rcomm,t_last 
      parameter (t_total = 1,t_gpairs = 2,t_randn = 3,t_rcomm = 4,t_last 
     & = 4) 
      double precision ::  tsum(t_last + 2),t1m(t_last + 2),tming(t_last 
     & + 2),tmaxg(t_last + 2) 
      character ::  t_recs(t_last + 2)*8 
      parameter (mk = 16,mm = m - mk,nn = 2** mm,nk = 2** mk,nq = 10,eps 
     &ilon = 1.d-8,a = 1220703125.d0,s = 271828183.d0) 
      common /storage/x(2 * nk),q(0:nq - 1),qq(10000) 
      data             dum /1.d0, 1.d0, 1.d0/   
      data t_recs/'total', 'gpairs', 'randn', 'rcomm',            ' totc 
     &omp', ' totcomm'/   
      double precision ::  r23_0,r23_1,r46_0,r46_1,t23_0,t23_1,t46_0,t46 
     &_1 
      parameter (r23_0 = 0.5d0** 23,t23_0 = 2.d0** 23,t23_1 = 2.d0** 23, 
     &r23_1 = 0.5d0** 23,r46_0 = r23_0** 2,t46_0 = t23_0** 2,t46_1 = t23 
     &_1** 2,r46_1 = r23_1** 2) 
      double precision ::  a1_0,a1_2,a2_0,a2_2,randlc_0,t1_0,t1_2,t2_0,t 
     &2_2,t3_0,t3_2,t4_0,t4_2,x1_0,x1_2,x2_0,x2_2,z_0,z_2,arg0 
      integer ::  arg_0,i_0 
	  intrinsic int
      call mpi_init(ierr) 
      call mpi_comm_rank(mpi_comm_world,node,ierr) 
      call mpi_comm_size(mpi_comm_world,no_nodes,ierr) 
      root = 0 
      if (.not.(convertdouble)) then 
         dp_type = mpi_double_precision 
      else   
         dp_type = mpi_real 
      endif   
      if (node .eq. root) then 
 
!   Because the size of the problem is too large to store in a 32-bit   
!   integer for some classes, we put it into a string (for printing).   
!   Have to strip off the decimal point put in there by the floating   
!   point print statement (internal file)   
         write (unit = *,fmt = 1000)  
         write (unit = size,fmt = '(f15.0)') 2.d0** (m + 1) 
         j = 15 
         if (size(j:j) .eq. '.')  j = j - 1 
         write (unit = *,fmt = 1001) size(1:j) 
         write (unit = *,fmt = 1003) no_nodes 
1000              format(/,' NAS Parallel Benchmarks 3.3 -- EP Benchmark 
     &',/)   
1001                  format(' Number of random numbers generated: ', a1 
     &5)   
1003                  format(' Number of active processes:         ', 2x 
     &, i13, /)  
         open (unit = 2,file = 'timer.flag',status = 'old',iostat = fsta 
     &tus) 
         timers_enabled = .FALSE. 
         if (fstatus .eq. 0) then 
            timers_enabled = .TRUE. 
            close (unit = 2) 
         endif   
      endif   
      call mpi_bcast(timers_enabled,1,mpi_logical,root,mpi_comm_world,ie 
     &rr) 
      verified = .FALSE. 
 
!   Compute the number of "batches" of random number pairs generated    
!   per processor. Adjust if the number of processors does not evenly    
!   divide the total number   
      np = nn / no_nodes 
      no_large_nodes = mod (nn,no_nodes) 
      if (node .lt. no_large_nodes) then 
         np_add = 1 
      else   
         np_add = 0 
      endif   
      np = np + np_add 
      if (np .eq. 0) then 
         write (unit = 6,fmt = 1) no_nodes,nn 
1                    format ('Too many nodes:',2i6)   
         ierrcode = 1 
         call mpi_abort(mpi_comm_world,ierrcode,ierr) 
         stop 
      endif   
 
!   Call the random number generator functions and initialize   
!   the x-array to reduce the effects of paging on the timings.   
!   Also, call all mathematical functions that are used. Make   
!   sure these initializations cannot be eliminated as dead code.   
!DVM$ GET_ACTUAL (dum) 
      call vranlc(0,dum(1),dum(2),dum(3)) 
!DVM$ ACTUAL (dum) 
      dum(1) = randlc (dum(2),dum(3)) 
!DVM$ ACTUAL (dum) 
!DVM$ REGION  
!DVM$ PARALLEL (i), PRIVATE (i),TIE (x(i)) 
      do  i = 1,2 * nk 
         x(i) = (-(1.d99)) 
      enddo   
!DVM$ END REGION 
      mops = log (sqrt (abs (max (1.d0,1.d0)))) 
 
!---------------------------------------------------------------------   
!      Synchronize before placing time stamp   
!---------------------------------------------------------------------   
      do  i = 1,t_last 
         call timer_clear(i) 
      enddo   
      call mpi_barrier(mpi_comm_world,ierr) 
      call timer_start(1) 
      t1 = a 
      call vranlc(0,t1,a,x) 
 
!   Compute AN = A ^ (2 * NK) (mod 2^46).   
      t1 = a 
      do  i = 1,mk + 1 
         t2 = randlc (t1,t1) 
      enddo   
      an = t1 
      tt = s 
      gc = 0.d0 
      sx = 0.d0 
      sy = 0.d0 
      do  i = 0,nq - 1 
         q(i) = 0.d0 
      enddo   
 
!   Each instance of this loop may be performed independently. We compute   
!   the k offsets separately to take into account the fact that some nodes   
!   have more numbers to generate than others   
      if (np_add .eq. 1) then 
         k_offset = node * np - 1 
      else   
         k_offset = no_large_nodes * (np + 1) + (node - no_large_nodes)  
     &* np - 1 
      endif   
!DVM$ REGION	  
!DVM$ PARALLEL(k),private( kk, t1, t2, i, ik, t1_0,a1_0,a2_0,x1_0,x2_0,
!DVM$& t2_0, z_0, t3_0, t4_0, randlc_0, t3, arg_0, t1_2, a1_2, a2_2, 
!DVM$& i_0, x1_2, x2_2, t2_2, z_2, t3_2, t4_2,x1,x2,arg0,t4,l),
!DVM$& reduction (sum(q),sum(sy),sum(sx))
      do  k = 1,np 
         kk = k_offset + k 
         t1 = s 
         t2 = an 
 
!        Find starting seed t1 for this kk.   
         do  i = 1,100 
            ik = kk / 2 
            if (2 * ik .ne. kk) then 
               t1_0 = r23_0 * t2 
               a1_0 = int (t1_0) 
               a2_0 = t2 - t23_0 * a1_0 
               t1_0 = r23_0 * t1 
               x1_0 = int (t1_0) 
               x2_0 = t1 - t23_0 * x1_0 
               t1_0 = a1_0 * x2_0 + a2_0 * x1_0 
               t2_0 = int (r23_0 * t1_0) 
               z_0 = t1_0 - t23_0 * t2_0 
               t3_0 = t23_0 * z_0 + a2_0 * x2_0 
               t4_0 = int (r46_0 * t3_0) 
               t1 = t3_0 - t46_0 * t4_0 
               randlc_0 = r46_0 * t1 
               t3 = randlc_0 
            endif   
            if (ik .eq. 0)  goto 130 
            t1_0 = r23_0 * t2 
            a1_0 = int (t1_0) 
            a2_0 = t2 - t23_0 * a1_0 
            t1_0 = r23_0 * t2 
            x1_0 = int (t1_0) 
            x2_0 = t2 - t23_0 * x1_0 
            t1_0 = a1_0 * x2_0 + a2_0 * x1_0 
            t2_0 = int (r23_0 * t1_0) 
            z_0 = t1_0 - t23_0 * t2_0 
            t3_0 = t23_0 * z_0 + a2_0 * x2_0 
            t4_0 = int (r46_0 * t3_0) 
            t2 = t3_0 - t46_0 * t4_0 
            randlc_0 = r46_0 * t2 
            t3 = randlc_0 
            kk = ik 
         enddo   
 
!        Compute uniform pseudorandom numbers.   
130      continue 
         arg_0 = 2 * nk 
         t1_2 = r23_1 * a 
         a1_2 = int (t1_2) 
         a2_2 = a - t23_1 * a1_2 
         do  i_0 = 1,nk 
            t1_2 = r23_1 * t1 
            x1_2 = int (t1_2) 
            x2_2 = t1 - t23_1 * x1_2 
            t1_2 = a1_2 * x2_2 + a2_2 * x1_2 
            t2_2 = int (r23_1 * t1_2) 
            z_2 = t1_2 - t23_1 * t2_2 
            t3_2 = t23_1 * z_2 + a2_2 * x2_2 
            t4_2 = int (r46_1 * t3_2) 
            t1 = t3_2 - t46_1 * t4_2 
            x1 = r46_1 * t1 
            t1_2 = r23_1 * t1 
            x1_2 = int (t1_2) 
            x2_2 = t1 - t23_1 * x1_2 
            t1_2 = a1_2 * x2_2 + a2_2 * x1_2 
            t2_2 = int (r23_1 * t1_2) 
            z_2 = t1_2 - t23_1 * t2_2 
            t3_2 = t23_1 * z_2 + a2_2 * x2_2 
            t4_2 = int (r46_1 * t3_2) 
            t1 = t3_2 - t46_1 * t4_2 
            x2 = r46_1 * t1 
 
!            x1 = 2.d0 * x(2 * i_0 - 1) - 1.d0 
!            x2 = 2.d0 * x(2 * i_0) - 1.d0 
            x1 = 2.d0 * x1 - 1.d0 
            x2 = 2.d0 * x2 - 1.d0 
            arg0 = x1** 2 + x2** 2 
            if (arg0 .le. 1.d0) then 
               t2 = sqrt ((-(2.d0)) * log (arg0) / arg0) 
               t3 = x1 * t2 
               t4 = x2 * t2 
               l = max (abs (t3),abs (t4)) 
               q(l) = q(l) + 1.d0 
               sx = sx + t3 
               sy = sy + t4 
            endif   
         enddo   
 
!         if (timers_enabled)  call timer_stop(t_gpairs)   
      enddo   
!DVM$ END REGION	  
      if (timers_enabled)  call timer_start(t_rcomm) 
!DVM$ GET_ACTUAL (x) 
      call mpi_allreduce(sx,x,1,dp_type,mpi_sum,mpi_comm_world,ierr) 
!DVM$ ACTUAL (x) 
      sx = x(1) 
!DVM$ GET_ACTUAL (x) 
      call mpi_allreduce(sy,x,1,dp_type,mpi_sum,mpi_comm_world,ierr) 
!DVM$ ACTUAL (x) 
      sy = x(1) 
!DVM$ GET_ACTUAL (x) 
      call mpi_allreduce(q,x,nq,dp_type,mpi_sum,mpi_comm_world,ierr) 
!DVM$ ACTUAL (x) 
      if (timers_enabled)  call timer_stop(t_rcomm) 
!DVM$ REGION  
!DVM$ PARALLEL (i), PRIVATE (i),TIE (x(i)) 
      do  i = 1,nq 
         q(i - 1) = x(i) 
      enddo   
!DVM$ END REGION 
      do  i = 0,nq - 1 
         gc = gc + q(i) 
      enddo   
      call timer_stop(1) 
      tm = timer_read (1) 
!DVM$ GET_ACTUAL (x) 
      call mpi_allreduce(tm,x,1,dp_type,mpi_max,mpi_comm_world,ierr) 
!DVM$ ACTUAL (x) 
      tm = x(1) 
      if (node .eq. root) then 
         nit = 0 
         verified = .TRUE. 
         if (m .eq. 24) then 
            sx_verify_value = (-(3.247834652034740d+3)) 
            sy_verify_value = (-(6.958407078382297d+3)) 
         else if (m .eq. 25) then   
            sx_verify_value = (-(2.863319731645753d+3)) 
            sy_verify_value = (-(6.320053679109499d+3)) 
         else if (m .eq. 28) then   
            sx_verify_value = (-(4.295875165629892d+3)) 
            sy_verify_value = (-(1.580732573678431d+4)) 
         else if (m .eq. 30) then   
            sx_verify_value = 4.033815542441498d+4 
            sy_verify_value = (-(2.660669192809235d+4)) 
         else if (m .eq. 32) then   
            sx_verify_value = 4.764367927995374d+4 
            sy_verify_value = (-(8.084072988043731d+4)) 
         else if (m .eq. 36) then   
            sx_verify_value = 1.982481200946593d+5 
            sy_verify_value = (-(1.020596636361769d+5)) 
         else if (m .eq. 40) then   
            sx_verify_value = (-(5.319717441530d+05)) 
            sy_verify_value = (-(3.688834557731d+05)) 
         else   
            verified = .FALSE. 
         endif   
         if (verified) then 
            sx_err = abs ((sx - sx_verify_value) / sx_verify_value) 
            sy_err = abs ((sy - sy_verify_value) / sy_verify_value) 
            verified = sx_err .le. epsilon .and. sy_err .le. epsilon 
         endif   
         mops = 2.d0** (m + 1) / tm / 1000000.d0 
         write (unit = 6,fmt = 11) tm,m,gc,sx,sy,(i,q(i), i = 0,nq - 1) 
11                   format ('EP Benchmark Results:'//'CPU Time =',f10.4 
     &/'N = 2^',           i5/'No. Gaussian Pairs =',f15.0/'Sums = ',1p, 
     &2d25.15/           'Counts:'/(i3,0p,f15.0))   
         call print_results('EP',class,32 + 1,0,0,nit,npm,no_nodes,tm,mo 
     &ps,'Random numbers generated',verified,npbversion,compiletime,cs1, 
     &cs2,cs3,cs4,cs5,cs6,cs7) 
      endif   
      if (.not.(timers_enabled))  goto 999 
!DVM$ GET_ACTUAL (t1m) 
      do  i = 1,t_last 
         t1m(i) = timer_read (i) 
      enddo   
!DVM$ ACTUAL (t1m) 
      t1m(t_last + 2) = t1m(t_rcomm) 
!DVM$ ACTUAL (t1m(t_last + 2)) 
      t1m(t_last + 1) = t1m(t_total) - t1m(t_last + 2) 
!DVM$ ACTUAL (t1m(t_last + 1)) 
      call mpi_reduce(t1m,tsum,4 + 2,dp_type,mpi_sum,0,mpi_comm_world,ie 
     &rr) 
!DVM$ GET_ACTUAL (t1m) 
      call mpi_reduce(t1m,tming,4 + 2,dp_type,mpi_min,0,mpi_comm_world,i 
     &err) 
!DVM$ GET_ACTUAL (t1m) 
      call mpi_reduce(t1m,tmaxg,4 + 2,dp_type,mpi_max,0,mpi_comm_world,i 
     &err) 
      if (node .eq. 0) then 
         write (unit = *,fmt = 800) no_nodes 
!DVM$    GET_ACTUAL (t_recs,tmaxg,tming,tsum) 
         do  i = 1,t_last + 2 
            tsum(i) = tsum(i) / no_nodes 
            write (unit = *,fmt = 810) i,t_recs(i),tming(i),tmaxg(i),tsu 
     &m(i) 
         enddo   
!DVM$    ACTUAL (tsum) 
      endif   
800   format(' nprocs =', i6, 11x, 'minimum', 5x, 'maximum',        5x,  
     &'average')   
810   format(' timer ', i2, '(', A8, ') :', 3(2x,f10.4))   
999   continue 
      call mpi_finalize(ierr) 
      end 
 
