c---------------------------------------------------------------------
c Заполнение трехмерной матрици нулями
c---------------------------------------------------------------------
c @param double precission :: z(n1 ,n2 ,n3) - некоторая матрица
c---------------------------------------------------------------------	 
      subroutine zero3(z,n1,n2,n3)
c---------------------------------------------------------------------
!DVM$ INHERIT z
c---------------------------------------------------------------------

      implicit none      

      integer n1, n2, n3
      double precision z(n1,n2,n3)
      integer i1, i2, i3
      
!!DVM$ INTERVAL 3   
!DVM$ REGION  
!DVM$ PARALLEL (i3,i2,i1) ON z(i1,i2,i3), PRIVATE(i3,i2,i1)
!DVM$&, cuda_block(32,6)    
      do  i3=1,n3      
         do  i2=1,n2
            do  i1=1,n1
               z(i1,i2,i3)=0.0D0
            enddo
         enddo
      enddo
!DVM$ END REGION      
!!DVM$ END INTERVAL
      
      return
      end


c----- end of program ------------------------------------------------

c -- DEBUG --------------------------------------------------------      

      subroutine printMatrix(a,n1,n2,n3)
c---------------------------------------------------------------------
!DVM$ INHERIT a
      integer n1,n2,n3,i1,i2,i3
      double precision a(n1,n2,n3), z(n2)
      integer m1, m2, m3

      write(*,*) 'MATRIX ------------------'
      write(*,*) a
      write(*,*) ' '
!      m1 = min(n1,18)
!      m2 = min(n2,14)
!      m3 = min(n3,18)

!      write(*,*)'  '
!      do  i3=1,m3
!         do  i1=1,m1
!            do i2=1,m2
!               z(i2) = a(i1,i2,i3)
!            enddo
!            write(*,6)(z(i2),i2=1,m2)
!         enddo
!         write(*,*)' - - - - - - - '
!      enddo
!      write(*,*)'  '
! 6    format(15f6.3)

      return 
      end 
      
      subroutine printMatrixNN(a,n1,n2,n3)
c---------------------------------------------------------------------
!DVM$ INHERIT a
      integer n1, n2, n3
      double precision a(n1,n2,n3), z

      do  i3=1,n3
         do  i2=1,n2
            do  i1=1,n1
               z = a(i1,i2,i3)
               if (z.ne.0) then
                write(*,*) '(',i1,',',i2,',',i3,')=',z
               endif  
            enddo
         enddo
      enddo
 
      end  
      

c -- EXTERNAL --------------------------------------------------------      
      
      subroutine timer_clear(n)

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

      implicit none
      integer n
      
      double precision start(64), elapsed(64)
      common /tt/ start, elapsed

      elapsed(n) = 0.0
      return
      end

      subroutine timer_start(n)

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

      implicit none
      external         elapsed_time
      double precision elapsed_time
      integer n
      double precision start(64), elapsed(64)
      common /tt/ start, elapsed

      start(n) = elapsed_time()
      return
      end
     
      subroutine timer_stop(n)

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

      implicit none
      external         elapsed_time
      double precision elapsed_time
      integer n
      double precision start(64), elapsed(64)
      common /tt/ start, elapsed
      double precision t, now

      now = elapsed_time()
      t = now - start(n)
      elapsed(n) = elapsed(n) + t
      return
      end
      
      
      double precision function timer_read(n)

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

      implicit none
      integer n
      double precision start(64), elapsed(64)
      common /tt/ start, elapsed
    
      timer_read = elapsed(n)
      return
      end

      double precision function elapsed_time()

c---------------------------------------------------------------------
c---------------------------------------------------------------------
      implicit none

        double precision t
        double precision dvtime

        data t/0.d0/ 
c This function must measure wall clock time, not CPU time. 
c Since there is no portable timer in Fortran (77)
c we call a routine compiled in C (though the C source may have
c to be tweaked). 
      t = dvtime()
c The following is not ok for "official" results because it reports
c CPU time not wall clock time. It may be useful for developing/testing
c on timeshared Crays, though. 
c     call second(t)

      elapsed_time = t

      return
      end

      subroutine print_results(name, class, n1, n2, n3, niter, 
     >               t, mops, optype, verified, npbversion) 
c               ,compiletime, cs1, cs2, cs3, cs4, cs5, cs6, cs7)
      
      implicit none
      character*2 name
      character*1 class
      integer n1, n2, n3, niter, j
      double precision t, mops
      character optype*24, size*13
      logical verified
      character*5 npbversion
c     >          , compiletime, cs1, cs2, cs3, cs4, cs5, cs6, cs7

         write (*, 2) name 
 2       format(//, ' ', A2, ' Benchmark Completed.')

         write (*, 3) Class
 3       format(' Class           = ', 12x, a12)

c   If this is not a grid-based problem (EP, FT, CG), then
c   we only print n1, which contains some measure of the
c   problem size. In that case, n2 and n3 are both zero.
c   Otherwise, we print the grid size n1xn2xn3

         if ((n2 .eq. 0) .and. (n3 .eq. 0)) then
            if (name(1:2) .eq. 'EP') then
               write(size, '(f12.0)' ) 2.d0**n1
               do j =13,1,-1
                  if (size(j:j) .eq. '.') size(j:j) = ' '
               end do
               write (*,42) size
 42            format(' Size            = ',12x, a14)
            else
               write (*,44) n1
 44            format(' Size            = ',12x, i12)
            endif
         else
            write (*, 4) n1,n2,n3
 4          format(' Size            =  ',12x, i3,'x',i3,'x',i3)
         endif

         write (*, 5) niter
 5       format(' Iterations      = ', 12x, i12)
         
         write (*, 6) t
 6       format(' Time in seconds = ',12x, f12.2)
         
         write (*,9) mops
 9       format(' Mop/s total     = ',12x, f12.2)

         write(*, 11) optype
 11      format(' Operation type  = ', a24)

         if (verified) then 
            write(*,12) '  SUCCESSFUL'
         else
            write(*,12) 'UNSUCCESSFUL'
         endif
 12      format(' Verification    = ', 12x, a)

         write(*,13) npbversion
 13      format(' Version         = ', 12x, a12)

c         write(*,14) compiletime
c 14      format(' Compile date    = ', 12x, a12)


c         write (*,121) cs1
c 121     format(/, ' Compile options:', /, 
c     >          '    F77          = ', A)

c         write (*,122) cs2
c 122     format('    FLINK        = ', A)

c         write (*,123) cs3
c 123     format('    F_LIB        = ', A)
c
c         write (*,124) cs4
c 124     format('    F_INC        = ', A)
c
c         write (*,125) cs5
c 125     format('    FFLAGS       = ', A)
c
c         write (*,126) cs6
c 126     format('    FLINKFLAGS   = ', A)
c
c         write(*, 127) cs7
c 127     format('    RAND         = ', A)
        
         write (*,130)
 130     format(//' Please send the results of this run to:'//
     >            ' NPB Development Team '/
     >            ' Internet: npb@nas.nasa.gov'/
     >            ' '/
     >            ' If email is not available, send this to:'//
     >            ' MS T27A-1'/
     >            ' NASA Ames Research Center'/
     >            ' Moffett Field, CA  94035-1000'//
     >            ' Fax: 415-604-3957'//)


         return
         end

   
      double precision function randlc (x, a)

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

c---------------------------------------------------------------------
c
c   This routine returns a uniform pseudorandom double precision number in the
c   range (0, 1) by using the linear congruential generator
c
c   x_{k+1} = a x_k  (mod 2^46)
c
c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
c   before repeating.  The argument A is the same as 'a' in the above formula,
c   and X is the same as x_0.  A and X must be odd double precision integers
c   in the range (1, 2^46).  The returned value RANDLC is normalized to be
c   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
c   the new seed x_1, so that subsequent calls to RANDLC using the same
c   arguments will generate a continuous sequence.
c
c   This routine should produce the same results on any computer with at least
c   48 mantissa bits in double precision floating point data.  On 64 bit
c   systems, double precision should be disabled.
c
c   David H. Bailey     October 26, 1990
c
c---------------------------------------------------------------------

      implicit none

      double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
      parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
     >  t46 = t23 ** 2)

c---------------------------------------------------------------------
c   Break A into two parts such that A = 2^23 * A1 + A2.
c---------------------------------------------------------------------

      t1 = r23 * a
      a1 = int (t1)
      a2 = a - t23 * a1

c---------------------------------------------------------------------
c   Break X into two parts such that X = 2^23 * X1 + X2, compute
c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
c   X = 2^23 * Z + A2 * X2  (mod 2^46).
c---------------------------------------------------------------------
      t1 = r23 * x
      x1 = int (t1)
      x2 = x - t23 * x1
      t1 = a1 * x2 + a2 * x1
      t2 = int (r23 * t1)
      z = t1 - t23 * t2
      t3 = t23 * z + a2 * x2
      t4 = int (r46 * t3)
      x = t3 - t46 * t4
      randlc = r46 * x

      return
      end

      subroutine vranlc (n, x, a, y)

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

c---------------------------------------------------------------------
c
c   This routine generates N uniform pseudorandom double precision numbers in
c   the range (0, 1) by using the linear congruential generator
c
c   x_{k+1} = a x_k  (mod 2^46)
c
c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
c   before repeating.  The argument A is the same as 'a' in the above formula,
c   and X is the same as x_0.  A and X must be odd double precision integers
c   in the range (1, 2^46).  The N results are placed in Y and are normalized
c   to be between 0 and 1.  X is updated to contain the new seed, so that
c   subsequent calls to VRANLC using the same arguments will generate a
c   continuous sequence.  If N is zero, only initialization is performed, and
c   the variables X, A and Y are ignored.
c
c   This routine is the standard version designed for scalar or RISC systems.
c   However, it should produce the same results on any single processor
c   computer with at least 48 mantissa bits in double precision floating point
c   data.  On 64 bit systems, double precision should be disabled.
c
c---------------------------------------------------------------------

      implicit none

      integer i,n
      double precision y,r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
      dimension y(*)
      parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
     >  t46 = t23 ** 2)


c---------------------------------------------------------------------
c   Break A into two parts such that A = 2^23 * A1 + A2.
c---------------------------------------------------------------------
      t1 = r23 * a
      a1 = int (t1)
      a2 = a - t23 * a1

c---------------------------------------------------------------------
c   Generate N results.   This loop is not vectorizable.
c---------------------------------------------------------------------
      do i = 1, n

c---------------------------------------------------------------------
c   Break X into two parts such that X = 2^23 * X1 + X2, compute
c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
c   X = 2^23 * Z + A2 * X2  (mod 2^46).
c---------------------------------------------------------------------
        t1 = r23 * x
        x1 = int (t1)
        x2 = x - t23 * x1
        t1 = a1 * x2 + a2 * x1
        t2 = int (r23 * t1)
        z = t1 - t23 * t2
        t3 = t23 * z + a2 * x2
        t4 = int (r46 * t3)
        x = t3 - t46 * t4
        y(i) = r46 * x
      enddo
      return
      end
c---------------------------------------------------------------------      