SP Report

SP
Exec Time (%)	Codelet Name	Error (%)
99.2	sp_MAIN___150	95.0
1.2	pinvr_pinvr__19	1.7
3.1	rhs_compute_rhs__21	10.9
3.4	rhs_compute_rhs__351	0.4
6.4	rhs_compute_rhs__273	0.9
1.2	ninvr_ninvr__19	5.7
2.0	add_add__18	5.6
2.5	txinvr_txinvr__21	7.6
2.9	tzetar_tzetar__22	4.6
21.1	z_solve_z_solve__27	6.8
1.2	rhs_compute_rhs__387	5.2
8.3	rhs_compute_rhs__64	0.4
21.0	y_solve_y_solve__27	0.8
13.4	x_solve_x_solve__26	0.4
2.2	rhs_compute_rhs__50	7.4
8.7	rhs_compute_rhs__166	6.3

Exec Time : 1.012064e+11 cycles

Matching : 94.5%

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

       subroutine compute_rhs

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

       include 'header.h'

       integer i, j, k, m
       double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
     >                  wijk, wp1, wm1


c---------------------------------------------------------------------
c      compute the reciprocal of density, and the kinetic energy, 
c      and the speed of sound. 
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                rho_inv = 1.0d0/u(1,i,j,k)
                rho_i(i,j,k) = rho_inv
                us(i,j,k) = u(2,i,j,k) * rho_inv
                vs(i,j,k) = u(3,i,j,k) * rho_inv
                ws(i,j,k) = u(4,i,j,k) * rho_inv
                square(i,j,k)     = 0.5d0* (
     >                        u(2,i,j,k)*u(2,i,j,k) + 
     >                        u(3,i,j,k)*u(3,i,j,k) +
     >                        u(4,i,j,k)*u(4,i,j,k) ) * rho_inv
                qs(i,j,k) = square(i,j,k) * rho_inv
c---------------------------------------------------------------------
c               (don't need speed and ainx until the lhs computation)
c---------------------------------------------------------------------
                aux = c1c2*rho_inv* (u(5,i,j,k) - square(i,j,k))
                speed(i,j,k) = dsqrt(aux)
c                ainv(i,j,k)  = 1.0d0/speed(i,j,k)
             end do
          end do
       end do

c---------------------------------------------------------------------
c copy the exact forcing term to the right hand side;  because 
c this forcing term is known, we can store it on the whole grid
c including the boundary                   
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                do    m = 1, 5
                   rhs(m,i,j,k) = forcing(m,i,j,k)
                end do
             end do
          end do
       end do

c---------------------------------------------------------------------
c      compute xi-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsx)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                uijk = us(i,j,k)
                up1  = us(i+1,j,k)
                um1  = us(i-1,j,k)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dx1tx1 * 
     >                    (u(1,i+1,j,k) - 2.0d0*u(1,i,j,k) + 
     >                     u(1,i-1,j,k)) -
     >                    tx2 * (u(2,i+1,j,k) - u(2,i-1,j,k))

                rhs(2,i,j,k) = rhs(2,i,j,k) + dx2tx1 * 
     >                    (u(2,i+1,j,k) - 2.0d0*u(2,i,j,k) + 
     >                     u(2,i-1,j,k)) +
     >                    xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
     >                    tx2 * (u(2,i+1,j,k)*up1 - 
     >                           u(2,i-1,j,k)*um1 +
     >                           (u(5,i+1,j,k)- square(i+1,j,k)-
     >                            u(5,i-1,j,k)+ square(i-1,j,k))*
     >                            c2)

                rhs(3,i,j,k) = rhs(3,i,j,k) + dx3tx1 * 
     >                    (u(3,i+1,j,k) - 2.0d0*u(3,i,j,k) +
     >                     u(3,i-1,j,k)) +
     >                    xxcon2 * (vs(i+1,j,k) - 2.0d0*vs(i,j,k) +
     >                              vs(i-1,j,k)) -
     >                    tx2 * (u(3,i+1,j,k)*up1 - 
     >                           u(3,i-1,j,k)*um1)

                rhs(4,i,j,k) = rhs(4,i,j,k) + dx4tx1 * 
     >                    (u(4,i+1,j,k) - 2.0d0*u(4,i,j,k) +
     >                     u(4,i-1,j,k)) +
     >                    xxcon2 * (ws(i+1,j,k) - 2.0d0*ws(i,j,k) +
     >                              ws(i-1,j,k)) -
     >                    tx2 * (u(4,i+1,j,k)*up1 - 
     >                           u(4,i-1,j,k)*um1)

                rhs(5,i,j,k) = rhs(5,i,j,k) + dx5tx1 * 
     >                    (u(5,i+1,j,k) - 2.0d0*u(5,i,j,k) +
     >                     u(5,i-1,j,k)) +
     >                    xxcon3 * (qs(i+1,j,k) - 2.0d0*qs(i,j,k) +
     >                              qs(i-1,j,k)) +
     >                    xxcon4 * (up1*up1 -       2.0d0*uijk*uijk + 
     >                              um1*um1) +
     >                    xxcon5 * (u(5,i+1,j,k)*rho_i(i+1,j,k) - 
     >                              2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                              u(5,i-1,j,k)*rho_i(i-1,j,k)) -
     >                    tx2 * ( (c1*u(5,i+1,j,k) - 
     >                             c2*square(i+1,j,k))*up1 -
     >                            (c1*u(5,i-1,j,k) - 
     >                             c2*square(i-1,j,k))*um1 )
             end do

c---------------------------------------------------------------------
c      add fourth order xi-direction dissipation               
c---------------------------------------------------------------------

             i = 1
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) +
     >                            u(m,i+2,j,k))
             end do

             i = 2
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i-1,j,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i+1,j,k) + u(m,i+2,j,k))
             end do

             do  i = 3, nx2-2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) + 
     >                         u(m,i+2,j,k) )
                end do
             end do

             i = nx2-1
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) )
             end do

             i = nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.d0*u(m,i-1,j,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsx)

c---------------------------------------------------------------------
c      compute eta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsy)
       do     k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                vijk = vs(i,j,k)
                vp1  = vs(i,j+1,k)
                vm1  = vs(i,j-1,k)
                rhs(1,i,j,k) = rhs(1,i,j,k) + dy1ty1 * 
     >                   (u(1,i,j+1,k) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j-1,k)) -
     >                   ty2 * (u(3,i,j+1,k) - u(3,i,j-1,k))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dy2ty1 * 
     >                   (u(2,i,j+1,k) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j-1,k)) +
     >                   yycon2 * (us(i,j+1,k) - 2.0d0*us(i,j,k) + 
     >                             us(i,j-1,k)) -
     >                   ty2 * (u(2,i,j+1,k)*vp1 - 
     >                          u(2,i,j-1,k)*vm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dy3ty1 * 
     >                   (u(3,i,j+1,k) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j-1,k)) +
     >                   yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
     >                   ty2 * (u(3,i,j+1,k)*vp1 - 
     >                          u(3,i,j-1,k)*vm1 +
     >                          (u(5,i,j+1,k) - square(i,j+1,k) - 
     >                           u(5,i,j-1,k) + square(i,j-1,k))
     >                          *c2)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dy4ty1 * 
     >                   (u(4,i,j+1,k) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j-1,k)) +
     >                   yycon2 * (ws(i,j+1,k) - 2.0d0*ws(i,j,k) + 
     >                             ws(i,j-1,k)) -
     >                   ty2 * (u(4,i,j+1,k)*vp1 - 
     >                          u(4,i,j-1,k)*vm1)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dy5ty1 * 
     >                   (u(5,i,j+1,k) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j-1,k)) +
     >                   yycon3 * (qs(i,j+1,k) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j-1,k)) +
     >                   yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk + 
     >                             vm1*vm1) +
     >                   yycon5 * (u(5,i,j+1,k)*rho_i(i,j+1,k) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j-1,k)*rho_i(i,j-1,k)) -
     >                   ty2 * ((c1*u(5,i,j+1,k) - 
     >                           c2*square(i,j+1,k)) * vp1 -
     >                          (c1*u(5,i,j-1,k) - 
     >                           c2*square(i,j-1,k)) * vm1)
             end do
          end do

c---------------------------------------------------------------------
c      add fourth order eta-direction dissipation         
c---------------------------------------------------------------------

          j = 1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) +
     >                            u(m,i,j+2,k))
             end do
          end do

          j = 2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j-1,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j+1,k) + u(m,i,j+2,k))
             end do
          end do

          do    j = 3, ny2-2
             do  i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) + 
     >                         u(m,i,j+2,k) )
                end do
             end do
          end do
 
          j = ny2-1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) )
             end do
          end do

          j = ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.d0*u(m,i,j-1,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsy)

c---------------------------------------------------------------------
c      compute zeta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsz)
       do    k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                wijk = ws(i,j,k)
                wp1  = ws(i,j,k+1)
                wm1  = ws(i,j,k-1)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dz1tz1 * 
     >                   (u(1,i,j,k+1) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j,k-1)) -
     >                   tz2 * (u(4,i,j,k+1) - u(4,i,j,k-1))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dz2tz1 * 
     >                   (u(2,i,j,k+1) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j,k-1)) +
     >                   zzcon2 * (us(i,j,k+1) - 2.0d0*us(i,j,k) + 
     >                             us(i,j,k-1)) -
     >                   tz2 * (u(2,i,j,k+1)*wp1 - 
     >                          u(2,i,j,k-1)*wm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dz3tz1 * 
     >                   (u(3,i,j,k+1) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j,k-1)) +
     >                   zzcon2 * (vs(i,j,k+1) - 2.0d0*vs(i,j,k) + 
     >                             vs(i,j,k-1)) -
     >                   tz2 * (u(3,i,j,k+1)*wp1 - 
     >                          u(3,i,j,k-1)*wm1)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dz4tz1 * 
     >                   (u(4,i,j,k+1) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j,k-1)) +
     >                   zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
     >                   tz2 * (u(4,i,j,k+1)*wp1 - 
     >                          u(4,i,j,k-1)*wm1 +
     >                          (u(5,i,j,k+1) - square(i,j,k+1) - 
     >                           u(5,i,j,k-1) + square(i,j,k-1))
     >                          *c2)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dz5tz1 * 
     >                   (u(5,i,j,k+1) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j,k-1)) +
     >                   zzcon3 * (qs(i,j,k+1) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j,k-1)) +
     >                   zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk + 
     >                             wm1*wm1) +
     >                   zzcon5 * (u(5,i,j,k+1)*rho_i(i,j,k+1) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j,k-1)*rho_i(i,j,k-1)) -
     >                   tz2 * ( (c1*u(5,i,j,k+1) - 
     >                            c2*square(i,j,k+1))*wp1 -
     >                           (c1*u(5,i,j,k-1) - 
     >                            c2*square(i,j,k-1))*wm1)
             end do
          end do
       end do

c---------------------------------------------------------------------
c      add fourth order zeta-direction dissipation                
c---------------------------------------------------------------------

       k = 1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) +
     >                            u(m,i,j,k+2))
             end do
          end do
       end do

       k = 2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j,k-1) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j,k+1) + u(m,i,j,k+2))
             end do
          end do
       end do

       do     k = 3, nz2-2
          do     j = 1, ny2
             do     i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) + 
     >                         u(m,i,j,k+2) )
                end do
             end do
          end do
       end do
 
       k = nz2-1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) )
             end do
          end do
       end do

       k = nz2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.d0*u(m,i,j,k-1) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsz)

       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                do    m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) * dt
                end do
             end do
          end do
       end do
    
       return
       end

invivorhs_compute_rhs__50

Real time: 2.248011e+09 -- Predicted time: 2.079801e+09

Callcount: 402

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	5.979548e+06	1.018011e+07	41.2
160	### 2	401.264282868583	5.168218e+06	5.582352e+06	7.4

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

       subroutine compute_rhs

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

       include 'header.h'

       integer i, j, k, m
       double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
     >                  wijk, wp1, wm1


c---------------------------------------------------------------------
c      compute the reciprocal of density, and the kinetic energy, 
c      and the speed of sound. 
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                rho_inv = 1.0d0/u(1,i,j,k)
                rho_i(i,j,k) = rho_inv
                us(i,j,k) = u(2,i,j,k) * rho_inv
                vs(i,j,k) = u(3,i,j,k) * rho_inv
                ws(i,j,k) = u(4,i,j,k) * rho_inv
                square(i,j,k)     = 0.5d0* (
     >                        u(2,i,j,k)*u(2,i,j,k) + 
     >                        u(3,i,j,k)*u(3,i,j,k) +
     >                        u(4,i,j,k)*u(4,i,j,k) ) * rho_inv
                qs(i,j,k) = square(i,j,k) * rho_inv
c---------------------------------------------------------------------
c               (don't need speed and ainx until the lhs computation)
c---------------------------------------------------------------------
                aux = c1c2*rho_inv* (u(5,i,j,k) - square(i,j,k))
                speed(i,j,k) = dsqrt(aux)
c                ainv(i,j,k)  = 1.0d0/speed(i,j,k)
             end do
          end do
       end do

c---------------------------------------------------------------------
c copy the exact forcing term to the right hand side;  because 
c this forcing term is known, we can store it on the whole grid
c including the boundary                   
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                do    m = 1, 5
                   rhs(m,i,j,k) = forcing(m,i,j,k)
                end do
             end do
          end do
       end do

c---------------------------------------------------------------------
c      compute xi-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsx)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                uijk = us(i,j,k)
                up1  = us(i+1,j,k)
                um1  = us(i-1,j,k)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dx1tx1 * 
     >                    (u(1,i+1,j,k) - 2.0d0*u(1,i,j,k) + 
     >                     u(1,i-1,j,k)) -
     >                    tx2 * (u(2,i+1,j,k) - u(2,i-1,j,k))

                rhs(2,i,j,k) = rhs(2,i,j,k) + dx2tx1 * 
     >                    (u(2,i+1,j,k) - 2.0d0*u(2,i,j,k) + 
     >                     u(2,i-1,j,k)) +
     >                    xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
     >                    tx2 * (u(2,i+1,j,k)*up1 - 
     >                           u(2,i-1,j,k)*um1 +
     >                           (u(5,i+1,j,k)- square(i+1,j,k)-
     >                            u(5,i-1,j,k)+ square(i-1,j,k))*
     >                            c2)

                rhs(3,i,j,k) = rhs(3,i,j,k) + dx3tx1 * 
     >                    (u(3,i+1,j,k) - 2.0d0*u(3,i,j,k) +
     >                     u(3,i-1,j,k)) +
     >                    xxcon2 * (vs(i+1,j,k) - 2.0d0*vs(i,j,k) +
     >                              vs(i-1,j,k)) -
     >                    tx2 * (u(3,i+1,j,k)*up1 - 
     >                           u(3,i-1,j,k)*um1)

                rhs(4,i,j,k) = rhs(4,i,j,k) + dx4tx1 * 
     >                    (u(4,i+1,j,k) - 2.0d0*u(4,i,j,k) +
     >                     u(4,i-1,j,k)) +
     >                    xxcon2 * (ws(i+1,j,k) - 2.0d0*ws(i,j,k) +
     >                              ws(i-1,j,k)) -
     >                    tx2 * (u(4,i+1,j,k)*up1 - 
     >                           u(4,i-1,j,k)*um1)

                rhs(5,i,j,k) = rhs(5,i,j,k) + dx5tx1 * 
     >                    (u(5,i+1,j,k) - 2.0d0*u(5,i,j,k) +
     >                     u(5,i-1,j,k)) +
     >                    xxcon3 * (qs(i+1,j,k) - 2.0d0*qs(i,j,k) +
     >                              qs(i-1,j,k)) +
     >                    xxcon4 * (up1*up1 -       2.0d0*uijk*uijk + 
     >                              um1*um1) +
     >                    xxcon5 * (u(5,i+1,j,k)*rho_i(i+1,j,k) - 
     >                              2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                              u(5,i-1,j,k)*rho_i(i-1,j,k)) -
     >                    tx2 * ( (c1*u(5,i+1,j,k) - 
     >                             c2*square(i+1,j,k))*up1 -
     >                            (c1*u(5,i-1,j,k) - 
     >                             c2*square(i-1,j,k))*um1 )
             end do

c---------------------------------------------------------------------
c      add fourth order xi-direction dissipation               
c---------------------------------------------------------------------

             i = 1
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) +
     >                            u(m,i+2,j,k))
             end do

             i = 2
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i-1,j,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i+1,j,k) + u(m,i+2,j,k))
             end do

             do  i = 3, nx2-2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) + 
     >                         u(m,i+2,j,k) )
                end do
             end do

             i = nx2-1
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) )
             end do

             i = nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.d0*u(m,i-1,j,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsx)

c---------------------------------------------------------------------
c      compute eta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsy)
       do     k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                vijk = vs(i,j,k)
                vp1  = vs(i,j+1,k)
                vm1  = vs(i,j-1,k)
                rhs(1,i,j,k) = rhs(1,i,j,k) + dy1ty1 * 
     >                   (u(1,i,j+1,k) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j-1,k)) -
     >                   ty2 * (u(3,i,j+1,k) - u(3,i,j-1,k))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dy2ty1 * 
     >                   (u(2,i,j+1,k) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j-1,k)) +
     >                   yycon2 * (us(i,j+1,k) - 2.0d0*us(i,j,k) + 
     >                             us(i,j-1,k)) -
     >                   ty2 * (u(2,i,j+1,k)*vp1 - 
     >                          u(2,i,j-1,k)*vm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dy3ty1 * 
     >                   (u(3,i,j+1,k) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j-1,k)) +
     >                   yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
     >                   ty2 * (u(3,i,j+1,k)*vp1 - 
     >                          u(3,i,j-1,k)*vm1 +
     >                          (u(5,i,j+1,k) - square(i,j+1,k) - 
     >                           u(5,i,j-1,k) + square(i,j-1,k))
     >                          *c2)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dy4ty1 * 
     >                   (u(4,i,j+1,k) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j-1,k)) +
     >                   yycon2 * (ws(i,j+1,k) - 2.0d0*ws(i,j,k) + 
     >                             ws(i,j-1,k)) -
     >                   ty2 * (u(4,i,j+1,k)*vp1 - 
     >                          u(4,i,j-1,k)*vm1)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dy5ty1 * 
     >                   (u(5,i,j+1,k) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j-1,k)) +
     >                   yycon3 * (qs(i,j+1,k) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j-1,k)) +
     >                   yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk + 
     >                             vm1*vm1) +
     >                   yycon5 * (u(5,i,j+1,k)*rho_i(i,j+1,k) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j-1,k)*rho_i(i,j-1,k)) -
     >                   ty2 * ((c1*u(5,i,j+1,k) - 
     >                           c2*square(i,j+1,k)) * vp1 -
     >                          (c1*u(5,i,j-1,k) - 
     >                           c2*square(i,j-1,k)) * vm1)
             end do
          end do

c---------------------------------------------------------------------
c      add fourth order eta-direction dissipation         
c---------------------------------------------------------------------

          j = 1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) +
     >                            u(m,i,j+2,k))
             end do
          end do

          j = 2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j-1,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j+1,k) + u(m,i,j+2,k))
             end do
          end do

          do    j = 3, ny2-2
             do  i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) + 
     >                         u(m,i,j+2,k) )
                end do
             end do
          end do
 
          j = ny2-1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) )
             end do
          end do

          j = ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.d0*u(m,i,j-1,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsy)

c---------------------------------------------------------------------
c      compute zeta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsz)
       do    k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                wijk = ws(i,j,k)
                wp1  = ws(i,j,k+1)
                wm1  = ws(i,j,k-1)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dz1tz1 * 
     >                   (u(1,i,j,k+1) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j,k-1)) -
     >                   tz2 * (u(4,i,j,k+1) - u(4,i,j,k-1))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dz2tz1 * 
     >                   (u(2,i,j,k+1) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j,k-1)) +
     >                   zzcon2 * (us(i,j,k+1) - 2.0d0*us(i,j,k) + 
     >                             us(i,j,k-1)) -
     >                   tz2 * (u(2,i,j,k+1)*wp1 - 
     >                          u(2,i,j,k-1)*wm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dz3tz1 * 
     >                   (u(3,i,j,k+1) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j,k-1)) +
     >                   zzcon2 * (vs(i,j,k+1) - 2.0d0*vs(i,j,k) + 
     >                             vs(i,j,k-1)) -
     >                   tz2 * (u(3,i,j,k+1)*wp1 - 
     >                          u(3,i,j,k-1)*wm1)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dz4tz1 * 
     >                   (u(4,i,j,k+1) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j,k-1)) +
     >                   zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
     >                   tz2 * (u(4,i,j,k+1)*wp1 - 
     >                          u(4,i,j,k-1)*wm1 +
     >                          (u(5,i,j,k+1) - square(i,j,k+1) - 
     >                           u(5,i,j,k-1) + square(i,j,k-1))
     >                          *c2)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dz5tz1 * 
     >                   (u(5,i,j,k+1) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j,k-1)) +
     >                   zzcon3 * (qs(i,j,k+1) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j,k-1)) +
     >                   zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk + 
     >                             wm1*wm1) +
     >                   zzcon5 * (u(5,i,j,k+1)*rho_i(i,j,k+1) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j,k-1)*rho_i(i,j,k-1)) -
     >                   tz2 * ( (c1*u(5,i,j,k+1) - 
     >                            c2*square(i,j,k+1))*wp1 -
     >                           (c1*u(5,i,j,k-1) - 
     >                            c2*square(i,j,k-1))*wm1)
             end do
          end do
       end do

c---------------------------------------------------------------------
c      add fourth order zeta-direction dissipation                
c---------------------------------------------------------------------

       k = 1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) +
     >                            u(m,i,j,k+2))
             end do
          end do
       end do

       k = 2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j,k-1) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j,k+1) + u(m,i,j,k+2))
             end do
          end do
       end do

       do     k = 3, nz2-2
          do     j = 1, ny2
             do     i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) + 
     >                         u(m,i,j,k+2) )
                end do
             end do
          end do
       end do
 
       k = nz2-1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) )
             end do
          end do
       end do

       k = nz2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.d0*u(m,i,j,k-1) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsz)

       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                do    m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) * dt
                end do
             end do
          end do
       end do
    
       return
       end

invivorhs_compute_rhs__166

Real time: 8.819109e+09 -- Predicted time: 9.420095e+09

Callcount: 402

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
121	### 1	401.700290809342	2.345056e+07	2.194015e+07	6.4

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

       subroutine z_solve

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

c---------------------------------------------------------------------
c this function performs the solution of the approximate factorization
c step in the z-direction for all five matrix components
c simultaneously. The Thomas algorithm is employed to solve the
c systems for the z-lines. Boundary conditions are non-periodic
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k, k1, k2, m
       double precision ru1, fac1, fac2, rtmp(5,0:KMAX)


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

       if (timeron) call timer_start(t_zsolve)
       do   j = 1, ny2
          do   i = 1, nx2

c---------------------------------------------------------------------
c Computes the left hand side for the three z-factors   
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c first fill the lhs for the u-eigenvalue                          
c---------------------------------------------------------------------

             do   k = 0, nz2+1
                ru1 = c3c4*rho_i(i,j,k)
                cv(k) = ws(i,j,k)
                rhos(k) = dmax1(dz4 + con43 * ru1,
     >                          dz5 + c1c5 * ru1,
     >                          dzmax + ru1,
     >                          dz1)
             end do

             call lhsinit(grid_points(3)-1)
             do   k =  1, nz2
                lhs(1,k) =  0.0d0
                lhs(2,k) = -dttz2 * cv(k-1) - dttz1 * rhos(k-1)
                lhs(3,k) =  1.0 + c2dttz1 * rhos(k)
                lhs(4,k) =  dttz2 * cv(k+1) - dttz1 * rhos(k+1)
                lhs(5,k) =  0.0d0
             end do

c---------------------------------------------------------------------
c      add fourth order dissipation                                  
c---------------------------------------------------------------------

             k = 1
             lhs(3,k) = lhs(3,k) + comz5
             lhs(4,k) = lhs(4,k) - comz4
             lhs(5,k) = lhs(5,k) + comz1

             k = 2
             lhs(2,k) = lhs(2,k) - comz4
             lhs(3,k) = lhs(3,k) + comz6
             lhs(4,k) = lhs(4,k) - comz4
             lhs(5,k) = lhs(5,k) + comz1

             do    k = 3, nz2-2
                lhs(1,k) = lhs(1,k) + comz1
                lhs(2,k) = lhs(2,k) - comz4
                lhs(3,k) = lhs(3,k) + comz6
                lhs(4,k) = lhs(4,k) - comz4
                lhs(5,k) = lhs(5,k) + comz1
             end do

             k = nz2-1
             lhs(1,k) = lhs(1,k) + comz1
             lhs(2,k) = lhs(2,k) - comz4
             lhs(3,k) = lhs(3,k) + comz6
             lhs(4,k) = lhs(4,k) - comz4

             k = nz2
             lhs(1,k) = lhs(1,k) + comz1
             lhs(2,k) = lhs(2,k) - comz4
             lhs(3,k) = lhs(3,k) + comz5


c---------------------------------------------------------------------
c      subsequently, fill the other factors (u+c), (u-c) 
c---------------------------------------------------------------------
             do    k = 1, nz2
                lhsp(1,k) = lhs(1,k)
                lhsp(2,k) = lhs(2,k) - 
     >                            dttz2 * speed(i,j,k-1)
                lhsp(3,k) = lhs(3,k)
                lhsp(4,k) = lhs(4,k) + 
     >                            dttz2 * speed(i,j,k+1)
                lhsp(5,k) = lhs(5,k)
                lhsm(1,k) = lhs(1,k)
                lhsm(2,k) = lhs(2,k) + 
     >                            dttz2 * speed(i,j,k-1)
                lhsm(3,k) = lhs(3,k)
                lhsm(4,k) = lhs(4,k) - 
     >                            dttz2 * speed(i,j,k+1)
                lhsm(5,k) = lhs(5,k)
             end do


c---------------------------------------------------------------------
c Load a row of K data
c---------------------------------------------------------------------

             do   k = 0, nz2+1
	        rtmp(1,k) = rhs(1,i,j,k)
	        rtmp(2,k) = rhs(2,i,j,k)
	        rtmp(3,k) = rhs(3,i,j,k)
	        rtmp(4,k) = rhs(4,i,j,k)
	        rtmp(5,k) = rhs(5,i,j,k)
	     end do

c---------------------------------------------------------------------
c                          FORWARD ELIMINATION  
c---------------------------------------------------------------------

             do    k = 0, grid_points(3)-3
                k1 = k  + 1
                k2 = k  + 2
                fac1      = 1.d0/lhs(3,k)
                lhs(4,k)  = fac1*lhs(4,k)
                lhs(5,k)  = fac1*lhs(5,k)
                do    m = 1, 3
                   rtmp(m,k) = fac1*rtmp(m,k)
                end do
                lhs(3,k1) = lhs(3,k1) -
     >                         lhs(2,k1)*lhs(4,k)
                lhs(4,k1) = lhs(4,k1) -
     >                         lhs(2,k1)*lhs(5,k)
                do    m = 1, 3
                   rtmp(m,k1) = rtmp(m,k1) -
     >                         lhs(2,k1)*rtmp(m,k)
                end do
                lhs(2,k2) = lhs(2,k2) -
     >                         lhs(1,k2)*lhs(4,k)
                lhs(3,k2) = lhs(3,k2) -
     >                         lhs(1,k2)*lhs(5,k)
                do    m = 1, 3
                   rtmp(m,k2) = rtmp(m,k2) -
     >                         lhs(1,k2)*rtmp(m,k)
                end do
             end do

c---------------------------------------------------------------------
c      The last two rows in this grid block are a bit different, 
c      since they do not have two more rows available for the
c      elimination of off-diagonal entries
c---------------------------------------------------------------------
             k  = grid_points(3)-2
             k1 = grid_points(3)-1
             fac1      = 1.d0/lhs(3,k)
             lhs(4,k)  = fac1*lhs(4,k)
             lhs(5,k)  = fac1*lhs(5,k)
             do    m = 1, 3
                rtmp(m,k) = fac1*rtmp(m,k)
             end do
             lhs(3,k1) = lhs(3,k1) -
     >                      lhs(2,k1)*lhs(4,k)
             lhs(4,k1) = lhs(4,k1) -
     >                      lhs(2,k1)*lhs(5,k)
             do    m = 1, 3
                rtmp(m,k1) = rtmp(m,k1) -
     >                      lhs(2,k1)*rtmp(m,k)
             end do
c---------------------------------------------------------------------
c               scale the last row immediately
c---------------------------------------------------------------------
             fac2      = 1.d0/lhs(3,k1)
             do    m = 1, 3
                rtmp(m,k1) = fac2*rtmp(m,k1)
             end do

c---------------------------------------------------------------------
c      do the u+c and the u-c factors               
c---------------------------------------------------------------------
             do    k = 0, grid_points(3)-3
             	k1 = k  + 1
             	k2 = k  + 2
	     	m = 4
             	fac1	   = 1.d0/lhsp(3,k)
             	lhsp(4,k)  = fac1*lhsp(4,k)
             	lhsp(5,k)  = fac1*lhsp(5,k)
             	rtmp(m,k)  = fac1*rtmp(m,k)
             	lhsp(3,k1) = lhsp(3,k1) -
     >       		    lhsp(2,k1)*lhsp(4,k)
             	lhsp(4,k1) = lhsp(4,k1) -
     >       		    lhsp(2,k1)*lhsp(5,k)
             	rtmp(m,k1) = rtmp(m,k1) -
     >       		    lhsp(2,k1)*rtmp(m,k)
             	lhsp(2,k2) = lhsp(2,k2) -
     >       		    lhsp(1,k2)*lhsp(4,k)
             	lhsp(3,k2) = lhsp(3,k2) -
     >       		    lhsp(1,k2)*lhsp(5,k)
             	rtmp(m,k2) = rtmp(m,k2) -
     >       		    lhsp(1,k2)*rtmp(m,k)
	     	m = 5
             	fac1	   = 1.d0/lhsm(3,k)
             	lhsm(4,k)  = fac1*lhsm(4,k)
             	lhsm(5,k)  = fac1*lhsm(5,k)
             	rtmp(m,k)  = fac1*rtmp(m,k)
             	lhsm(3,k1) = lhsm(3,k1) -
     >       		    lhsm(2,k1)*lhsm(4,k)
             	lhsm(4,k1) = lhsm(4,k1) -
     >       		    lhsm(2,k1)*lhsm(5,k)
             	rtmp(m,k1) = rtmp(m,k1) -
     >       		    lhsm(2,k1)*rtmp(m,k)
             	lhsm(2,k2) = lhsm(2,k2) -
     >       		    lhsm(1,k2)*lhsm(4,k)
             	lhsm(3,k2) = lhsm(3,k2) -
     >       		    lhsm(1,k2)*lhsm(5,k)
             	rtmp(m,k2) = rtmp(m,k2) -
     >       		    lhsm(1,k2)*rtmp(m,k)
             end do

c---------------------------------------------------------------------
c         And again the last two rows separately
c---------------------------------------------------------------------
             k  = grid_points(3)-2
             k1 = grid_points(3)-1
	     m = 4
             fac1	= 1.d0/lhsp(3,k)
             lhsp(4,k)  = fac1*lhsp(4,k)
             lhsp(5,k)  = fac1*lhsp(5,k)
             rtmp(m,k)  = fac1*rtmp(m,k)
             lhsp(3,k1) = lhsp(3,k1) -
     >       		 lhsp(2,k1)*lhsp(4,k)
             lhsp(4,k1) = lhsp(4,k1) -
     >       		 lhsp(2,k1)*lhsp(5,k)
             rtmp(m,k1) = rtmp(m,k1) -
     >       		 lhsp(2,k1)*rtmp(m,k)
	     m = 5
             fac1	= 1.d0/lhsm(3,k)
             lhsm(4,k)  = fac1*lhsm(4,k)
             lhsm(5,k)  = fac1*lhsm(5,k)
             rtmp(m,k)  = fac1*rtmp(m,k)
             lhsm(3,k1) = lhsm(3,k1) -
     >       		 lhsm(2,k1)*lhsm(4,k)
             lhsm(4,k1) = lhsm(4,k1) -
     >       		 lhsm(2,k1)*lhsm(5,k)
             rtmp(m,k1) = rtmp(m,k1) -
     >       		 lhsm(2,k1)*rtmp(m,k)
c---------------------------------------------------------------------
c               Scale the last row immediately (some of this is overkill
c               if this is the last cell)
c---------------------------------------------------------------------
             rtmp(4,k1) = rtmp(4,k1)/lhsp(3,k1)
             rtmp(5,k1) = rtmp(5,k1)/lhsm(3,k1)


c---------------------------------------------------------------------
c                         BACKSUBSTITUTION 
c---------------------------------------------------------------------

             k  = grid_points(3)-2
             k1 = grid_points(3)-1
             do   m = 1, 3
                rtmp(m,k) = rtmp(m,k) -
     >                             lhs(4,k)*rtmp(m,k1)
             end do

             rtmp(4,k) = rtmp(4,k) -
     >                             lhsp(4,k)*rtmp(4,k1)
             rtmp(5,k) = rtmp(5,k) -
     >                             lhsm(4,k)*rtmp(5,k1)

c---------------------------------------------------------------------
c      Whether or not this is the last processor, we always have
c      to complete the back-substitution 
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c      The first three factors
c---------------------------------------------------------------------
             do   k = grid_points(3)-3, 0, -1
                k1 = k  + 1
                k2 = k  + 2
                do   m = 1, 3
                   rtmp(m,k) = rtmp(m,k) - 
     >                          lhs(4,k)*rtmp(m,k1) -
     >                          lhs(5,k)*rtmp(m,k2)
                end do

c---------------------------------------------------------------------
c      And the remaining two
c---------------------------------------------------------------------
                rtmp(4,k) = rtmp(4,k) - 
     >                          lhsp(4,k)*rtmp(4,k1) -
     >                          lhsp(5,k)*rtmp(4,k2)
                rtmp(5,k) = rtmp(5,k) - 
     >                          lhsm(4,k)*rtmp(5,k1) -
     >                          lhsm(5,k)*rtmp(5,k2)
             end do

c---------------------------------------------------------------------
c      Store result
c---------------------------------------------------------------------
             do   k = 0, nz2+1
	        rhs(1,i,j,k) = rtmp(1,k)
	        rhs(2,i,j,k) = rtmp(2,k)
	        rhs(3,i,j,k) = rtmp(3,k)
	        rhs(4,i,j,k) = rtmp(4,k)
	        rhs(5,i,j,k) = rtmp(5,k)
	     end do

          end do
       end do
       if (timeron) call timer_stop(t_zsolve)

       if (timeron) call timer_start(t_tzetar)
       call tzetar
       if (timeron) call timer_stop(t_tzetar)

       return
       end

invivoz_solve_z_solve__27

Real time: 2.145082e+10 -- Predicted time: 2.303774e+10

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
206	### 1	400.870219834059	5.746932e+07	5.340264e+07	7.0

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

       subroutine  add

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

c---------------------------------------------------------------------
c addition of update to the vector u
c---------------------------------------------------------------------

       include 'header.h'

       integer i,j,k,m

       do k = 1, nz2
          do j = 1, ny2
             do i = 1, nx2
                do m = 1, 5
                   u(m,i,j,k) = u(m,i,j,k) + rhs(m,i,j,k)
                end do
	     end do
	  end do
       end do

       return
       end

invivoadd_add__18

Real time: 2.046508e+09 -- Predicted time: 1.931876e+09

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
204	### 1	394.090049716406	4.840534e+06	4.912664e+06	1.4
366	### 2	4.99815140053304	4.855808e+06	5.069784e+06	4.2

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

       subroutine y_solve

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

c---------------------------------------------------------------------
c this function performs the solution of the approximate factorization
c step in the y-direction for all five matrix components
c simultaneously. The Thomas algorithm is employed to solve the
c systems for the y-lines. Boundary conditions are non-periodic
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k, j1, j2, m
       double precision ru1, fac1, fac2


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

       if (timeron) call timer_start(t_ysolve)
       do  k = 1, grid_points(3)-2
          do  i = 1, grid_points(1)-2

c---------------------------------------------------------------------
c Computes the left hand side for the three y-factors   
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c      first fill the lhs for the u-eigenvalue         
c---------------------------------------------------------------------

             do  j = 0, grid_points(2)-1
                ru1 = c3c4*rho_i(i,j,k)
                cv(j) = vs(i,j,k)
                rhoq(j) = dmax1( dy3 + con43 * ru1,
     >                           dy5 + c1c5*ru1,
     >                           dymax + ru1,
     >                           dy1)
             end do
            
             call lhsinit(grid_points(2)-1)
             do  j = 1, grid_points(2)-2
                lhs(1,j) =  0.0d0
                lhs(2,j) = -dtty2 * cv(j-1) - dtty1 * rhoq(j-1)
                lhs(3,j) =  1.0 + c2dtty1 * rhoq(j)
                lhs(4,j) =  dtty2 * cv(j+1) - dtty1 * rhoq(j+1)
                lhs(5,j) =  0.0d0
             end do

c---------------------------------------------------------------------
c      add fourth order dissipation                             
c---------------------------------------------------------------------

             j = 1

             lhs(3,j) = lhs(3,j) + comz5
             lhs(4,j) = lhs(4,j) - comz4
             lhs(5,j) = lhs(5,j) + comz1
       
             lhs(2,j+1) = lhs(2,j+1) - comz4
             lhs(3,j+1) = lhs(3,j+1) + comz6
             lhs(4,j+1) = lhs(4,j+1) - comz4
             lhs(5,j+1) = lhs(5,j+1) + comz1

             do   j=3, grid_points(2)-4

                lhs(1,j) = lhs(1,j) + comz1
                lhs(2,j) = lhs(2,j) - comz4
                lhs(3,j) = lhs(3,j) + comz6
                lhs(4,j) = lhs(4,j) - comz4
                lhs(5,j) = lhs(5,j) + comz1
             end do

             j = grid_points(2)-3
             lhs(1,j) = lhs(1,j) + comz1
             lhs(2,j) = lhs(2,j) - comz4
             lhs(3,j) = lhs(3,j) + comz6
             lhs(4,j) = lhs(4,j) - comz4

             lhs(1,j+1) = lhs(1,j+1) + comz1
             lhs(2,j+1) = lhs(2,j+1) - comz4
             lhs(3,j+1) = lhs(3,j+1) + comz5

c---------------------------------------------------------------------
c      subsequently, do the other two factors                    
c---------------------------------------------------------------------
             do    j = 1, grid_points(2)-2
                lhsp(1,j) = lhs(1,j)
                lhsp(2,j) = lhs(2,j) - 
     >                            dtty2 * speed(i,j-1,k)
                lhsp(3,j) = lhs(3,j)
                lhsp(4,j) = lhs(4,j) + 
     >                            dtty2 * speed(i,j+1,k)
                lhsp(5,j) = lhs(5,j)
                lhsm(1,j) = lhs(1,j)
                lhsm(2,j) = lhs(2,j) + 
     >                            dtty2 * speed(i,j-1,k)
                lhsm(3,j) = lhs(3,j)
                lhsm(4,j) = lhs(4,j) - 
     >                            dtty2 * speed(i,j+1,k)
                lhsm(5,j) = lhs(5,j)
             end do


c---------------------------------------------------------------------
c                          FORWARD ELIMINATION  
c---------------------------------------------------------------------

             do    j = 0, grid_points(2)-3
                j1 = j  + 1
                j2 = j  + 2
                fac1      = 1.d0/lhs(3,j)
                lhs(4,j)  = fac1*lhs(4,j)
                lhs(5,j)  = fac1*lhs(5,j)
                do    m = 1, 3
                   rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
                end do
                lhs(3,j1) = lhs(3,j1) -
     >                         lhs(2,j1)*lhs(4,j)
                lhs(4,j1) = lhs(4,j1) -
     >                         lhs(2,j1)*lhs(5,j)
                do    m = 1, 3
                   rhs(m,i,j1,k) = rhs(m,i,j1,k) -
     >                         lhs(2,j1)*rhs(m,i,j,k)
                end do
                lhs(2,j2) = lhs(2,j2) -
     >                         lhs(1,j2)*lhs(4,j)
                lhs(3,j2) = lhs(3,j2) -
     >                         lhs(1,j2)*lhs(5,j)
                do    m = 1, 3
                   rhs(m,i,j2,k) = rhs(m,i,j2,k) -
     >                         lhs(1,j2)*rhs(m,i,j,k)
                end do
             end do

c---------------------------------------------------------------------
c      The last two rows in this grid block are a bit different, 
c      since they do not have two more rows available for the
c      elimination of off-diagonal entries
c---------------------------------------------------------------------

             j  = grid_points(2)-2
             j1 = grid_points(2)-1
             fac1      = 1.d0/lhs(3,j)
             lhs(4,j)  = fac1*lhs(4,j)
             lhs(5,j)  = fac1*lhs(5,j)
             do    m = 1, 3
                rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             end do
             lhs(3,j1) = lhs(3,j1) -
     >                      lhs(2,j1)*lhs(4,j)
             lhs(4,j1) = lhs(4,j1) -
     >                      lhs(2,j1)*lhs(5,j)
             do    m = 1, 3
                rhs(m,i,j1,k) = rhs(m,i,j1,k) -
     >                      lhs(2,j1)*rhs(m,i,j,k)
             end do
c---------------------------------------------------------------------
c            scale the last row immediately 
c---------------------------------------------------------------------
             fac2      = 1.d0/lhs(3,j1)
             do    m = 1, 3
                rhs(m,i,j1,k) = fac2*rhs(m,i,j1,k)
             end do

c---------------------------------------------------------------------
c      do the u+c and the u-c factors                 
c---------------------------------------------------------------------
             do    j = 0, grid_points(2)-3
             	j1 = j  + 1
             	j2 = j  + 2
	     	m = 4
             	fac1	   = 1.d0/lhsp(3,j)
             	lhsp(4,j)  = fac1*lhsp(4,j)
             	lhsp(5,j)  = fac1*lhsp(5,j)
             	rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             	lhsp(3,j1) = lhsp(3,j1) -
     >       		    lhsp(2,j1)*lhsp(4,j)
             	lhsp(4,j1) = lhsp(4,j1) -
     >       		    lhsp(2,j1)*lhsp(5,j)
             	rhs(m,i,j1,k) = rhs(m,i,j1,k) -
     >       		    lhsp(2,j1)*rhs(m,i,j,k)
             	lhsp(2,j2) = lhsp(2,j2) -
     >       		    lhsp(1,j2)*lhsp(4,j)
             	lhsp(3,j2) = lhsp(3,j2) -
     >       		    lhsp(1,j2)*lhsp(5,j)
             	rhs(m,i,j2,k) = rhs(m,i,j2,k) -
     >       		    lhsp(1,j2)*rhs(m,i,j,k)
	     	m = 5
             	fac1	   = 1.d0/lhsm(3,j)
             	lhsm(4,j)  = fac1*lhsm(4,j)
             	lhsm(5,j)  = fac1*lhsm(5,j)
             	rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             	lhsm(3,j1) = lhsm(3,j1) -
     >       		    lhsm(2,j1)*lhsm(4,j)
             	lhsm(4,j1) = lhsm(4,j1) -
     >       		    lhsm(2,j1)*lhsm(5,j)
             	rhs(m,i,j1,k) = rhs(m,i,j1,k) -
     >       		    lhsm(2,j1)*rhs(m,i,j,k)
             	lhsm(2,j2) = lhsm(2,j2) -
     >       		    lhsm(1,j2)*lhsm(4,j)
             	lhsm(3,j2) = lhsm(3,j2) -
     >       		    lhsm(1,j2)*lhsm(5,j)
             	rhs(m,i,j2,k) = rhs(m,i,j2,k) -
     >       		    lhsm(1,j2)*rhs(m,i,j,k)
             end do

c---------------------------------------------------------------------
c         And again the last two rows separately
c---------------------------------------------------------------------
             j  = grid_points(2)-2
             j1 = grid_points(2)-1
	     m = 4
             fac1	= 1.d0/lhsp(3,j)
             lhsp(4,j)  = fac1*lhsp(4,j)
             lhsp(5,j)  = fac1*lhsp(5,j)
             rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             lhsp(3,j1) = lhsp(3,j1) -
     >       		 lhsp(2,j1)*lhsp(4,j)
             lhsp(4,j1) = lhsp(4,j1) -
     >       		 lhsp(2,j1)*lhsp(5,j)
             rhs(m,i,j1,k)   = rhs(m,i,j1,k) -
     >       		 lhsp(2,j1)*rhs(m,i,j,k)
	     m = 5
             fac1	= 1.d0/lhsm(3,j)
             lhsm(4,j)  = fac1*lhsm(4,j)
             lhsm(5,j)  = fac1*lhsm(5,j)
             rhs(m,i,j,k) = fac1*rhs(m,i,j,k)
             lhsm(3,j1) = lhsm(3,j1) -
     >       		 lhsm(2,j1)*lhsm(4,j)
             lhsm(4,j1) = lhsm(4,j1) -
     >       		 lhsm(2,j1)*lhsm(5,j)
             rhs(m,i,j1,k)   = rhs(m,i,j1,k) -
     >       		 lhsm(2,j1)*rhs(m,i,j,k)
c---------------------------------------------------------------------
c               Scale the last row immediately 
c---------------------------------------------------------------------
             rhs(4,i,j1,k)   = rhs(4,i,j1,k)/lhsp(3,j1)
             rhs(5,i,j1,k)   = rhs(5,i,j1,k)/lhsm(3,j1)


c---------------------------------------------------------------------
c                         BACKSUBSTITUTION 
c---------------------------------------------------------------------

             j  = grid_points(2)-2
             j1 = grid_points(2)-1
             do   m = 1, 3
                rhs(m,i,j,k) = rhs(m,i,j,k) -
     >                           lhs(4,j)*rhs(m,i,j1,k)
             end do

             rhs(4,i,j,k) = rhs(4,i,j,k) -
     >                           lhsp(4,j)*rhs(4,i,j1,k)
             rhs(5,i,j,k) = rhs(5,i,j,k) -
     >                           lhsm(4,j)*rhs(5,i,j1,k)

c---------------------------------------------------------------------
c      The first three factors
c---------------------------------------------------------------------
             do   j = grid_points(2)-3, 0, -1
                j1 = j  + 1
                j2 = j  + 2
                do   m = 1, 3
                   rhs(m,i,j,k) = rhs(m,i,j,k) - 
     >                          lhs(4,j)*rhs(m,i,j1,k) -
     >                          lhs(5,j)*rhs(m,i,j2,k)
                end do

c---------------------------------------------------------------------
c      And the remaining two
c---------------------------------------------------------------------
                rhs(4,i,j,k) = rhs(4,i,j,k) - 
     >                          lhsp(4,j)*rhs(4,i,j1,k) -
     >                          lhsp(5,j)*rhs(4,i,j2,k)
                rhs(5,i,j,k) = rhs(5,i,j,k) - 
     >                          lhsm(4,j)*rhs(5,i,j1,k) -
     >                          lhsm(5,j)*rhs(5,i,j2,k)
             end do
          end do
       end do
       if (timeron) call timer_stop(t_ysolve)

       if (timeron) call timer_start(t_pinvr)
       call pinvr
       if (timeron) call timer_stop(t_pinvr)

       return
       end

invivoy_solve_y_solve__27

Real time: 2.131032e+10 -- Predicted time: 2.150313e+10

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
270	### 1	401.242968082245	5.359130e+07	5.332230e+07	0.5

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

       subroutine  ninvr

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

c---------------------------------------------------------------------
c   block-diagonal matrix-vector multiplication              
c---------------------------------------------------------------------

       include 'header.h'

       integer  i, j, k
       double precision r1, r2, r3, r4, r5, t1, t2

       do k = 1, nz2
          do j = 1, ny2
             do i = 1, nx2

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)
               
                t1 = bt * r3
                t2 = 0.5d0 * ( r4 + r5 )

                rhs(1,i,j,k) = -r2
                rhs(2,i,j,k) =  r1
                rhs(3,i,j,k) = bt * ( r4 - r5 )
                rhs(4,i,j,k) = -t1 + t2
                rhs(5,i,j,k) =  t1 + t2
             enddo    
          enddo
       enddo

       return
       end

invivoninvr_ninvr__19

Real time: 1.303880e+09 -- Predicted time: 1.229155e+09

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
4	### 1	7.92277729696289	3.069106e+06	3.711292e+06	17.3
355	### 2	393.699167132387	3.060303e+06	3.242292e+06	5.6

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

       subroutine compute_rhs

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

       include 'header.h'

       integer i, j, k, m
       double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
     >                  wijk, wp1, wm1


c---------------------------------------------------------------------
c      compute the reciprocal of density, and the kinetic energy, 
c      and the speed of sound. 
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                rho_inv = 1.0d0/u(1,i,j,k)
                rho_i(i,j,k) = rho_inv
                us(i,j,k) = u(2,i,j,k) * rho_inv
                vs(i,j,k) = u(3,i,j,k) * rho_inv
                ws(i,j,k) = u(4,i,j,k) * rho_inv
                square(i,j,k)     = 0.5d0* (
     >                        u(2,i,j,k)*u(2,i,j,k) + 
     >                        u(3,i,j,k)*u(3,i,j,k) +
     >                        u(4,i,j,k)*u(4,i,j,k) ) * rho_inv
                qs(i,j,k) = square(i,j,k) * rho_inv
c---------------------------------------------------------------------
c               (don't need speed and ainx until the lhs computation)
c---------------------------------------------------------------------
                aux = c1c2*rho_inv* (u(5,i,j,k) - square(i,j,k))
                speed(i,j,k) = dsqrt(aux)
c                ainv(i,j,k)  = 1.0d0/speed(i,j,k)
             end do
          end do
       end do

c---------------------------------------------------------------------
c copy the exact forcing term to the right hand side;  because 
c this forcing term is known, we can store it on the whole grid
c including the boundary                   
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                do    m = 1, 5
                   rhs(m,i,j,k) = forcing(m,i,j,k)
                end do
             end do
          end do
       end do

c---------------------------------------------------------------------
c      compute xi-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsx)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                uijk = us(i,j,k)
                up1  = us(i+1,j,k)
                um1  = us(i-1,j,k)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dx1tx1 * 
     >                    (u(1,i+1,j,k) - 2.0d0*u(1,i,j,k) + 
     >                     u(1,i-1,j,k)) -
     >                    tx2 * (u(2,i+1,j,k) - u(2,i-1,j,k))

                rhs(2,i,j,k) = rhs(2,i,j,k) + dx2tx1 * 
     >                    (u(2,i+1,j,k) - 2.0d0*u(2,i,j,k) + 
     >                     u(2,i-1,j,k)) +
     >                    xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
     >                    tx2 * (u(2,i+1,j,k)*up1 - 
     >                           u(2,i-1,j,k)*um1 +
     >                           (u(5,i+1,j,k)- square(i+1,j,k)-
     >                            u(5,i-1,j,k)+ square(i-1,j,k))*
     >                            c2)

                rhs(3,i,j,k) = rhs(3,i,j,k) + dx3tx1 * 
     >                    (u(3,i+1,j,k) - 2.0d0*u(3,i,j,k) +
     >                     u(3,i-1,j,k)) +
     >                    xxcon2 * (vs(i+1,j,k) - 2.0d0*vs(i,j,k) +
     >                              vs(i-1,j,k)) -
     >                    tx2 * (u(3,i+1,j,k)*up1 - 
     >                           u(3,i-1,j,k)*um1)

                rhs(4,i,j,k) = rhs(4,i,j,k) + dx4tx1 * 
     >                    (u(4,i+1,j,k) - 2.0d0*u(4,i,j,k) +
     >                     u(4,i-1,j,k)) +
     >                    xxcon2 * (ws(i+1,j,k) - 2.0d0*ws(i,j,k) +
     >                              ws(i-1,j,k)) -
     >                    tx2 * (u(4,i+1,j,k)*up1 - 
     >                           u(4,i-1,j,k)*um1)

                rhs(5,i,j,k) = rhs(5,i,j,k) + dx5tx1 * 
     >                    (u(5,i+1,j,k) - 2.0d0*u(5,i,j,k) +
     >                     u(5,i-1,j,k)) +
     >                    xxcon3 * (qs(i+1,j,k) - 2.0d0*qs(i,j,k) +
     >                              qs(i-1,j,k)) +
     >                    xxcon4 * (up1*up1 -       2.0d0*uijk*uijk + 
     >                              um1*um1) +
     >                    xxcon5 * (u(5,i+1,j,k)*rho_i(i+1,j,k) - 
     >                              2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                              u(5,i-1,j,k)*rho_i(i-1,j,k)) -
     >                    tx2 * ( (c1*u(5,i+1,j,k) - 
     >                             c2*square(i+1,j,k))*up1 -
     >                            (c1*u(5,i-1,j,k) - 
     >                             c2*square(i-1,j,k))*um1 )
             end do

c---------------------------------------------------------------------
c      add fourth order xi-direction dissipation               
c---------------------------------------------------------------------

             i = 1
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) +
     >                            u(m,i+2,j,k))
             end do

             i = 2
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i-1,j,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i+1,j,k) + u(m,i+2,j,k))
             end do

             do  i = 3, nx2-2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) + 
     >                         u(m,i+2,j,k) )
                end do
             end do

             i = nx2-1
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) )
             end do

             i = nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.d0*u(m,i-1,j,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsx)

c---------------------------------------------------------------------
c      compute eta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsy)
       do     k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                vijk = vs(i,j,k)
                vp1  = vs(i,j+1,k)
                vm1  = vs(i,j-1,k)
                rhs(1,i,j,k) = rhs(1,i,j,k) + dy1ty1 * 
     >                   (u(1,i,j+1,k) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j-1,k)) -
     >                   ty2 * (u(3,i,j+1,k) - u(3,i,j-1,k))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dy2ty1 * 
     >                   (u(2,i,j+1,k) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j-1,k)) +
     >                   yycon2 * (us(i,j+1,k) - 2.0d0*us(i,j,k) + 
     >                             us(i,j-1,k)) -
     >                   ty2 * (u(2,i,j+1,k)*vp1 - 
     >                          u(2,i,j-1,k)*vm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dy3ty1 * 
     >                   (u(3,i,j+1,k) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j-1,k)) +
     >                   yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
     >                   ty2 * (u(3,i,j+1,k)*vp1 - 
     >                          u(3,i,j-1,k)*vm1 +
     >                          (u(5,i,j+1,k) - square(i,j+1,k) - 
     >                           u(5,i,j-1,k) + square(i,j-1,k))
     >                          *c2)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dy4ty1 * 
     >                   (u(4,i,j+1,k) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j-1,k)) +
     >                   yycon2 * (ws(i,j+1,k) - 2.0d0*ws(i,j,k) + 
     >                             ws(i,j-1,k)) -
     >                   ty2 * (u(4,i,j+1,k)*vp1 - 
     >                          u(4,i,j-1,k)*vm1)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dy5ty1 * 
     >                   (u(5,i,j+1,k) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j-1,k)) +
     >                   yycon3 * (qs(i,j+1,k) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j-1,k)) +
     >                   yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk + 
     >                             vm1*vm1) +
     >                   yycon5 * (u(5,i,j+1,k)*rho_i(i,j+1,k) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j-1,k)*rho_i(i,j-1,k)) -
     >                   ty2 * ((c1*u(5,i,j+1,k) - 
     >                           c2*square(i,j+1,k)) * vp1 -
     >                          (c1*u(5,i,j-1,k) - 
     >                           c2*square(i,j-1,k)) * vm1)
             end do
          end do

c---------------------------------------------------------------------
c      add fourth order eta-direction dissipation         
c---------------------------------------------------------------------

          j = 1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) +
     >                            u(m,i,j+2,k))
             end do
          end do

          j = 2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j-1,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j+1,k) + u(m,i,j+2,k))
             end do
          end do

          do    j = 3, ny2-2
             do  i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) + 
     >                         u(m,i,j+2,k) )
                end do
             end do
          end do
 
          j = ny2-1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) )
             end do
          end do

          j = ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.d0*u(m,i,j-1,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsy)

c---------------------------------------------------------------------
c      compute zeta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsz)
       do    k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                wijk = ws(i,j,k)
                wp1  = ws(i,j,k+1)
                wm1  = ws(i,j,k-1)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dz1tz1 * 
     >                   (u(1,i,j,k+1) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j,k-1)) -
     >                   tz2 * (u(4,i,j,k+1) - u(4,i,j,k-1))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dz2tz1 * 
     >                   (u(2,i,j,k+1) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j,k-1)) +
     >                   zzcon2 * (us(i,j,k+1) - 2.0d0*us(i,j,k) + 
     >                             us(i,j,k-1)) -
     >                   tz2 * (u(2,i,j,k+1)*wp1 - 
     >                          u(2,i,j,k-1)*wm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dz3tz1 * 
     >                   (u(3,i,j,k+1) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j,k-1)) +
     >                   zzcon2 * (vs(i,j,k+1) - 2.0d0*vs(i,j,k) + 
     >                             vs(i,j,k-1)) -
     >                   tz2 * (u(3,i,j,k+1)*wp1 - 
     >                          u(3,i,j,k-1)*wm1)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dz4tz1 * 
     >                   (u(4,i,j,k+1) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j,k-1)) +
     >                   zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
     >                   tz2 * (u(4,i,j,k+1)*wp1 - 
     >                          u(4,i,j,k-1)*wm1 +
     >                          (u(5,i,j,k+1) - square(i,j,k+1) - 
     >                           u(5,i,j,k-1) + square(i,j,k-1))
     >                          *c2)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dz5tz1 * 
     >                   (u(5,i,j,k+1) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j,k-1)) +
     >                   zzcon3 * (qs(i,j,k+1) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j,k-1)) +
     >                   zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk + 
     >                             wm1*wm1) +
     >                   zzcon5 * (u(5,i,j,k+1)*rho_i(i,j,k+1) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j,k-1)*rho_i(i,j,k-1)) -
     >                   tz2 * ( (c1*u(5,i,j,k+1) - 
     >                            c2*square(i,j,k+1))*wp1 -
     >                           (c1*u(5,i,j,k-1) - 
     >                            c2*square(i,j,k-1))*wm1)
             end do
          end do
       end do

c---------------------------------------------------------------------
c      add fourth order zeta-direction dissipation                
c---------------------------------------------------------------------

       k = 1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) +
     >                            u(m,i,j,k+2))
             end do
          end do
       end do

       k = 2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j,k-1) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j,k+1) + u(m,i,j,k+2))
             end do
          end do
       end do

       do     k = 3, nz2-2
          do     j = 1, ny2
             do     i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) + 
     >                         u(m,i,j,k+2) )
                end do
             end do
          end do
       end do
 
       k = nz2-1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) )
             end do
          end do
       end do

       k = nz2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.d0*u(m,i,j,k-1) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsz)

       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                do    m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) * dt
                end do
             end do
          end do
       end do
    
       return
       end

invivorhs_compute_rhs__273

Real time: 6.486764e+09 -- Predicted time: 6.425747e+09

Callcount: 402

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
78	### 1	403.282463847128	1.593361e+07	1.619263e+07	1.5

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

       subroutine  txinvr

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

c---------------------------------------------------------------------
c block-diagonal matrix-vector multiplication                  
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k
       double precision t1, t2, t3, ac, ru1, uu, vv, ww, r1, r2, r3, 
     >                  r4, r5, ac2inv


       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2

                ru1 = rho_i(i,j,k)
                uu = us(i,j,k)
                vv = vs(i,j,k)
                ww = ws(i,j,k)
                ac = speed(i,j,k)
c                ac2inv = ainv(i,j,k)*ainv(i,j,k)
                ac2inv = 1.0d0 / ( ac*ac )

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)

                t1 = c2 * ac2inv * ( qs(i,j,k)*r1 - uu*r2  - 
     >                  vv*r3 - ww*r4 + r5 )
                t2 = bt * ru1 * ( uu * r1 - r2 )
                t3 = ( bt * ru1 * ac ) * t1

                rhs(1,i,j,k) = r1 - t1
                rhs(2,i,j,k) = - ru1 * ( ww*r1 - r4 )
                rhs(3,i,j,k) =   ru1 * ( vv*r1 - r3 )
                rhs(4,i,j,k) = - t2 + t3
                rhs(5,i,j,k) =   t2 + t3

             end do
          end do
       end do

       return
       end

invivotxinvr_txinvr__21

Real time: 2.554633e+09 -- Predicted time: 2.359392e+09

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
23	### 1	397.989488499392	5.883671e+06	6.363316e+06	7.5
133	### 2	3.00483108029631	5.908012e+06	6.543464e+06	9.7

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

       subroutine  tzetar

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

c---------------------------------------------------------------------
c   block-diagonal matrix-vector multiplication                       
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k
       double precision  t1, t2, t3, ac, xvel, yvel, zvel, r1, r2, r3, 
     >                   r4, r5, btuz, ac2u, uzik1
c     >                   , acinv


       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2

                xvel = us(i,j,k)
                yvel = vs(i,j,k)
                zvel = ws(i,j,k)
                ac   = speed(i,j,k)
c                acinv = ainv(i,j,k)

                ac2u = ac*ac

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)      

                uzik1 = u(1,i,j,k)
                btuz  = bt * uzik1

c                t1 = btuz*acinv * (r4 + r5)
                t1 = btuz/ac * (r4 + r5)
                t2 = r3 + t1
                t3 = btuz * (r4 - r5)

                rhs(1,i,j,k) = t2
                rhs(2,i,j,k) = -uzik1*r2 + xvel*t2
                rhs(3,i,j,k) =  uzik1*r1 + yvel*t2
                rhs(4,i,j,k) =  zvel*t2  + t3
                rhs(5,i,j,k) =  uzik1*(-xvel*r2 + yvel*r1) + 
     >                    qs(i,j,k)*t2 + c2iv*ac2u*t1 + zvel*t3

             end do
          end do
       end do

       return
       end

invivotzetar_tzetar__22

Real time: 2.982277e+09 -- Predicted time: 3.128143e+09

Callcount: 401

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
58	### 1	400.153952326549	7.817349e+06	7.418972e+06	5.0

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

       subroutine compute_rhs

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

       include 'header.h'

       integer i, j, k, m
       double precision aux, rho_inv, uijk, up1, um1, vijk, vp1, vm1,
     >                  wijk, wp1, wm1


c---------------------------------------------------------------------
c      compute the reciprocal of density, and the kinetic energy, 
c      and the speed of sound. 
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                rho_inv = 1.0d0/u(1,i,j,k)
                rho_i(i,j,k) = rho_inv
                us(i,j,k) = u(2,i,j,k) * rho_inv
                vs(i,j,k) = u(3,i,j,k) * rho_inv
                ws(i,j,k) = u(4,i,j,k) * rho_inv
                square(i,j,k)     = 0.5d0* (
     >                        u(2,i,j,k)*u(2,i,j,k) + 
     >                        u(3,i,j,k)*u(3,i,j,k) +
     >                        u(4,i,j,k)*u(4,i,j,k) ) * rho_inv
                qs(i,j,k) = square(i,j,k) * rho_inv
c---------------------------------------------------------------------
c               (don't need speed and ainx until the lhs computation)
c---------------------------------------------------------------------
                aux = c1c2*rho_inv* (u(5,i,j,k) - square(i,j,k))
                speed(i,j,k) = dsqrt(aux)
c                ainv(i,j,k)  = 1.0d0/speed(i,j,k)
             end do
          end do
       end do

c---------------------------------------------------------------------
c copy the exact forcing term to the right hand side;  because 
c this forcing term is known, we can store it on the whole grid
c including the boundary                   
c---------------------------------------------------------------------

       do    k = 0, grid_points(3)-1
          do    j = 0, grid_points(2)-1
             do    i = 0, grid_points(1)-1
                do    m = 1, 5
                   rhs(m,i,j,k) = forcing(m,i,j,k)
                end do
             end do
          end do
       end do

c---------------------------------------------------------------------
c      compute xi-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsx)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                uijk = us(i,j,k)
                up1  = us(i+1,j,k)
                um1  = us(i-1,j,k)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dx1tx1 * 
     >                    (u(1,i+1,j,k) - 2.0d0*u(1,i,j,k) + 
     >                     u(1,i-1,j,k)) -
     >                    tx2 * (u(2,i+1,j,k) - u(2,i-1,j,k))

                rhs(2,i,j,k) = rhs(2,i,j,k) + dx2tx1 * 
     >                    (u(2,i+1,j,k) - 2.0d0*u(2,i,j,k) + 
     >                     u(2,i-1,j,k)) +
     >                    xxcon2*con43 * (up1 - 2.0d0*uijk + um1) -
     >                    tx2 * (u(2,i+1,j,k)*up1 - 
     >                           u(2,i-1,j,k)*um1 +
     >                           (u(5,i+1,j,k)- square(i+1,j,k)-
     >                            u(5,i-1,j,k)+ square(i-1,j,k))*
     >                            c2)

                rhs(3,i,j,k) = rhs(3,i,j,k) + dx3tx1 * 
     >                    (u(3,i+1,j,k) - 2.0d0*u(3,i,j,k) +
     >                     u(3,i-1,j,k)) +
     >                    xxcon2 * (vs(i+1,j,k) - 2.0d0*vs(i,j,k) +
     >                              vs(i-1,j,k)) -
     >                    tx2 * (u(3,i+1,j,k)*up1 - 
     >                           u(3,i-1,j,k)*um1)

                rhs(4,i,j,k) = rhs(4,i,j,k) + dx4tx1 * 
     >                    (u(4,i+1,j,k) - 2.0d0*u(4,i,j,k) +
     >                     u(4,i-1,j,k)) +
     >                    xxcon2 * (ws(i+1,j,k) - 2.0d0*ws(i,j,k) +
     >                              ws(i-1,j,k)) -
     >                    tx2 * (u(4,i+1,j,k)*up1 - 
     >                           u(4,i-1,j,k)*um1)

                rhs(5,i,j,k) = rhs(5,i,j,k) + dx5tx1 * 
     >                    (u(5,i+1,j,k) - 2.0d0*u(5,i,j,k) +
     >                     u(5,i-1,j,k)) +
     >                    xxcon3 * (qs(i+1,j,k) - 2.0d0*qs(i,j,k) +
     >                              qs(i-1,j,k)) +
     >                    xxcon4 * (up1*up1 -       2.0d0*uijk*uijk + 
     >                              um1*um1) +
     >                    xxcon5 * (u(5,i+1,j,k)*rho_i(i+1,j,k) - 
     >                              2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                              u(5,i-1,j,k)*rho_i(i-1,j,k)) -
     >                    tx2 * ( (c1*u(5,i+1,j,k) - 
     >                             c2*square(i+1,j,k))*up1 -
     >                            (c1*u(5,i-1,j,k) - 
     >                             c2*square(i-1,j,k))*um1 )
             end do

c---------------------------------------------------------------------
c      add fourth order xi-direction dissipation               
c---------------------------------------------------------------------

             i = 1
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) +
     >                            u(m,i+2,j,k))
             end do

             i = 2
             do    m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i-1,j,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i+1,j,k) + u(m,i+2,j,k))
             end do

             do  i = 3, nx2-2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) + 
     >                         u(m,i+2,j,k) )
                end do
             end do

             i = nx2-1
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.0d0*u(m,i-1,j,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i+1,j,k) )
             end do

             i = nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i-2,j,k) - 4.d0*u(m,i-1,j,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsx)

c---------------------------------------------------------------------
c      compute eta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsy)
       do     k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                vijk = vs(i,j,k)
                vp1  = vs(i,j+1,k)
                vm1  = vs(i,j-1,k)
                rhs(1,i,j,k) = rhs(1,i,j,k) + dy1ty1 * 
     >                   (u(1,i,j+1,k) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j-1,k)) -
     >                   ty2 * (u(3,i,j+1,k) - u(3,i,j-1,k))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dy2ty1 * 
     >                   (u(2,i,j+1,k) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j-1,k)) +
     >                   yycon2 * (us(i,j+1,k) - 2.0d0*us(i,j,k) + 
     >                             us(i,j-1,k)) -
     >                   ty2 * (u(2,i,j+1,k)*vp1 - 
     >                          u(2,i,j-1,k)*vm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dy3ty1 * 
     >                   (u(3,i,j+1,k) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j-1,k)) +
     >                   yycon2*con43 * (vp1 - 2.0d0*vijk + vm1) -
     >                   ty2 * (u(3,i,j+1,k)*vp1 - 
     >                          u(3,i,j-1,k)*vm1 +
     >                          (u(5,i,j+1,k) - square(i,j+1,k) - 
     >                           u(5,i,j-1,k) + square(i,j-1,k))
     >                          *c2)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dy4ty1 * 
     >                   (u(4,i,j+1,k) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j-1,k)) +
     >                   yycon2 * (ws(i,j+1,k) - 2.0d0*ws(i,j,k) + 
     >                             ws(i,j-1,k)) -
     >                   ty2 * (u(4,i,j+1,k)*vp1 - 
     >                          u(4,i,j-1,k)*vm1)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dy5ty1 * 
     >                   (u(5,i,j+1,k) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j-1,k)) +
     >                   yycon3 * (qs(i,j+1,k) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j-1,k)) +
     >                   yycon4 * (vp1*vp1       - 2.0d0*vijk*vijk + 
     >                             vm1*vm1) +
     >                   yycon5 * (u(5,i,j+1,k)*rho_i(i,j+1,k) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j-1,k)*rho_i(i,j-1,k)) -
     >                   ty2 * ((c1*u(5,i,j+1,k) - 
     >                           c2*square(i,j+1,k)) * vp1 -
     >                          (c1*u(5,i,j-1,k) - 
     >                           c2*square(i,j-1,k)) * vm1)
             end do
          end do

c---------------------------------------------------------------------
c      add fourth order eta-direction dissipation         
c---------------------------------------------------------------------

          j = 1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) +
     >                            u(m,i,j+2,k))
             end do
          end do

          j = 2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j-1,k) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j+1,k) + u(m,i,j+2,k))
             end do
          end do

          do    j = 3, ny2-2
             do  i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) + 
     >                         u(m,i,j+2,k) )
                end do
             end do
          end do
 
          j = ny2-1
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.0d0*u(m,i,j-1,k) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j+1,k) )
             end do
          end do

          j = ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j-2,k) - 4.d0*u(m,i,j-1,k) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsy)

c---------------------------------------------------------------------
c      compute zeta-direction fluxes 
c---------------------------------------------------------------------
       if (timeron) call timer_start(t_rhsz)
       do    k = 1, nz2
          do     j = 1, ny2
             do     i = 1, nx2
                wijk = ws(i,j,k)
                wp1  = ws(i,j,k+1)
                wm1  = ws(i,j,k-1)

                rhs(1,i,j,k) = rhs(1,i,j,k) + dz1tz1 * 
     >                   (u(1,i,j,k+1) - 2.0d0*u(1,i,j,k) + 
     >                    u(1,i,j,k-1)) -
     >                   tz2 * (u(4,i,j,k+1) - u(4,i,j,k-1))
                rhs(2,i,j,k) = rhs(2,i,j,k) + dz2tz1 * 
     >                   (u(2,i,j,k+1) - 2.0d0*u(2,i,j,k) + 
     >                    u(2,i,j,k-1)) +
     >                   zzcon2 * (us(i,j,k+1) - 2.0d0*us(i,j,k) + 
     >                             us(i,j,k-1)) -
     >                   tz2 * (u(2,i,j,k+1)*wp1 - 
     >                          u(2,i,j,k-1)*wm1)
                rhs(3,i,j,k) = rhs(3,i,j,k) + dz3tz1 * 
     >                   (u(3,i,j,k+1) - 2.0d0*u(3,i,j,k) + 
     >                    u(3,i,j,k-1)) +
     >                   zzcon2 * (vs(i,j,k+1) - 2.0d0*vs(i,j,k) + 
     >                             vs(i,j,k-1)) -
     >                   tz2 * (u(3,i,j,k+1)*wp1 - 
     >                          u(3,i,j,k-1)*wm1)
                rhs(4,i,j,k) = rhs(4,i,j,k) + dz4tz1 * 
     >                   (u(4,i,j,k+1) - 2.0d0*u(4,i,j,k) + 
     >                    u(4,i,j,k-1)) +
     >                   zzcon2*con43 * (wp1 - 2.0d0*wijk + wm1) -
     >                   tz2 * (u(4,i,j,k+1)*wp1 - 
     >                          u(4,i,j,k-1)*wm1 +
     >                          (u(5,i,j,k+1) - square(i,j,k+1) - 
     >                           u(5,i,j,k-1) + square(i,j,k-1))
     >                          *c2)
                rhs(5,i,j,k) = rhs(5,i,j,k) + dz5tz1 * 
     >                   (u(5,i,j,k+1) - 2.0d0*u(5,i,j,k) + 
     >                    u(5,i,j,k-1)) +
     >                   zzcon3 * (qs(i,j,k+1) - 2.0d0*qs(i,j,k) + 
     >                             qs(i,j,k-1)) +
     >                   zzcon4 * (wp1*wp1 - 2.0d0*wijk*wijk + 
     >                             wm1*wm1) +
     >                   zzcon5 * (u(5,i,j,k+1)*rho_i(i,j,k+1) - 
     >                             2.0d0*u(5,i,j,k)*rho_i(i,j,k) +
     >                             u(5,i,j,k-1)*rho_i(i,j,k-1)) -
     >                   tz2 * ( (c1*u(5,i,j,k+1) - 
     >                            c2*square(i,j,k+1))*wp1 -
     >                           (c1*u(5,i,j,k-1) - 
     >                            c2*square(i,j,k-1))*wm1)
             end do
          end do
       end do

c---------------------------------------------------------------------
c      add fourth order zeta-direction dissipation                
c---------------------------------------------------------------------

       k = 1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k)- dssp * 
     >                    ( 5.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) +
     >                            u(m,i,j,k+2))
             end do
          end do
       end do

       k = 2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (-4.0d0*u(m,i,j,k-1) + 6.0d0*u(m,i,j,k) -
     >                      4.0d0*u(m,i,j,k+1) + u(m,i,j,k+2))
             end do
          end do
       end do

       do     k = 3, nz2-2
          do     j = 1, ny2
             do     i = 1,nx2
                do     m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) - dssp * 
     >                    (  u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                     6.0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) + 
     >                         u(m,i,j,k+2) )
                end do
             end do
          end do
       end do
 
       k = nz2-1
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.0d0*u(m,i,j,k-1) + 
     >                      6.0d0*u(m,i,j,k) - 4.0d0*u(m,i,j,k+1) )
             end do
          end do
       end do

       k = nz2
       do     j = 1, ny2
          do     i = 1, nx2
             do     m = 1, 5
                rhs(m,i,j,k) = rhs(m,i,j,k) - dssp *
     >                    ( u(m,i,j,k-2) - 4.d0*u(m,i,j,k-1) +
     >                      5.d0*u(m,i,j,k) )
             end do
          end do
       end do
       if (timeron) call timer_stop(t_rhsz)

       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2
                do    m = 1, 5
                   rhs(m,i,j,k) = rhs(m,i,j,k) * dt
                end do
             end do
          end do
       end do
    
       return
       end

invivorhs_compute_rhs__21

Real time: 3.234310e+09 -- Predicted time: 2.880587e+09

Callcount: 402

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	8.068562e+06	1.793138e+07	55.0
68	### 2	397.082756010758	7.197631e+06	8.106964e+06	11.2
54	### 3	1	7.224746e+06	8.326252e+06	13.2
233	### 4	1	7.238414e+06	9.239280e+06	21.6

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

       subroutine pinvr

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

c---------------------------------------------------------------------
c   block-diagonal matrix-vector multiplication                       
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k
       double precision r1, r2, r3, r4, r5, t1, t2

       do   k = 1, nz2
          do   j = 1, ny2
             do   i = 1, nx2

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)

                t1 = bt * r1
                t2 = 0.5d0 * ( r4 + r5 )

                rhs(1,i,j,k) =  bt * ( r4 - r5 )
                rhs(2,i,j,k) = -r3
                rhs(3,i,j,k) =  r2
                rhs(4,i,j,k) = -t1 + t2
                rhs(5,i,j,k) =  t1 + t2
             end do
          end do
       end do

       return
       end

SP

Exec Time : 1.012064e+11 cycles

Matching : 94.5%

__invivo__rhs_compute_rhs__50

Real time: 2.248011e+09 -- Predicted time: 2.079801e+09

Callcount: 402

__invivo__rhs_compute_rhs__166

Real time: 8.819109e+09 -- Predicted time: 9.420095e+09

Callcount: 402

__invivo__z_solve_z_solve__27

Real time: 2.145082e+10 -- Predicted time: 2.303774e+10

Callcount: 401

__invivo__add_add__18

Real time: 2.046508e+09 -- Predicted time: 1.931876e+09

Callcount: 401

__invivo__y_solve_y_solve__27

Real time: 2.131032e+10 -- Predicted time: 2.150313e+10

Callcount: 401

__invivo__ninvr_ninvr__19

Real time: 1.303880e+09 -- Predicted time: 1.229155e+09

Callcount: 401

__invivo__rhs_compute_rhs__273

Real time: 6.486764e+09 -- Predicted time: 6.425747e+09

Callcount: 402

__invivo__txinvr_txinvr__21

Real time: 2.554633e+09 -- Predicted time: 2.359392e+09

Callcount: 401

__invivo__tzetar_tzetar__22

Real time: 2.982277e+09 -- Predicted time: 3.128143e+09

Callcount: 401

__invivo__rhs_compute_rhs__21

Real time: 3.234310e+09 -- Predicted time: 2.880587e+09

Callcount: 402

invivorhs_compute_rhs__50

invivorhs_compute_rhs__166

invivoz_solve_z_solve__27

invivoadd_add__18

invivoy_solve_y_solve__27

invivoninvr_ninvr__19

invivorhs_compute_rhs__273

invivotxinvr_txinvr__21

invivotzetar_tzetar__22

invivorhs_compute_rhs__21