LU Report

LU
Exec Time (%)	Codelet Name	Error (%)
19.5	buts_buts__69	9.3
13.3	jacu_jacu__39	18.7
21.2	blts_blts__72	4.3
2.3	blts_blts__56	2.1
2.2	buts_buts__55	14.9
13.4	jacld_jacld__39	21.2
1.9	rhs_rhs__34	2.9
5.6	rhs_rhs__55	2.0
9.4	rhs_rhs__177	5.1
9.3	rhs_rhs__302	6.1

Exec Time : 1.153398e+11 cycles

Matching : 71.4%

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

      subroutine rhs

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

c---------------------------------------------------------------------
c   compute the right hand sides
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m
      double precision  q
      double precision  tmp
      double precision  u21, u31, u41
      double precision  u21i, u31i, u41i, u51i
      double precision  u21j, u31j, u41j, u51j
      double precision  u21k, u31k, u41k, u51k
      double precision  u21im1, u31im1, u41im1, u51im1
      double precision  u21jm1, u31jm1, u41jm1, u51jm1
      double precision  u21km1, u31km1, u41km1, u51km1



      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
               tmp = 1.0d+00 / u(1,i,j,k)
               rho_i(i,j,k) = tmp
               qs(i,j,k) = 0.50d+00 * (  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) )
     >                      * tmp
            end do
         end do
      end do

      if (timeron) call timer_start(t_rhsx)
c---------------------------------------------------------------------
c   xi-direction flux differences
c---------------------------------------------------------------------

      do k = 2, nz - 1
         do j = jst, jend
            do i = 1, nx
               flux(1,i) = u(2,i,j,k)
               u21 = u(2,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,i) = u(2,i,j,k) * u21 + c2 * 
     >                        ( u(5,i,j,k) - q )
               flux(3,i) = u(3,i,j,k) * u21
               flux(4,i) = u(4,i,j,k) * u21
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
            end do

            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                 - tx2 * ( flux(m,i+1) - flux(m,i-1) )
               end do
            end do

            do i = ist, nx
               tmp = rho_i(i,j,k)

               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)

               tmp = rho_i(i-1,j,k)

               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)

               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tx3 * ( u21i**2 - u21im1**2 )
     >              + c1 * c5 * tx3 * ( u51i - u51im1 )
            end do

            do i = ist, iend
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dx1 * tx1 * (            u(1,i-1,j,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i+1,j,k) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     >              + dx2 * tx1 * (            u(2,i-1,j,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i+1,j,k) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     >              + dx3 * tx1 * (            u(3,i-1,j,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i+1,j,k) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     >              + dx4 * tx1 * (            u(4,i-1,j,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i+1,j,k) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     >              + dx5 * tx1 * (            u(5,i-1,j,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i+1,j,k) )
            end do

c---------------------------------------------------------------------
c   Fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,2,j,k) = rsd(m,2,j,k)
     >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     >                      - 4.0d+00 * u(m,3,j,k)
     >                      +           u(m,4,j,k) )
               rsd(m,3,j,k) = rsd(m,3,j,k)
     >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     >                      + 6.0d+00 * u(m,3,j,k)
     >                      - 4.0d+00 * u(m,4,j,k)
     >                      +           u(m,5,j,k) )
            end do

            do i = 4, nx - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i-2,j,k)
     >                        - 4.0d+00 * u(m,i-1,j,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i+1,j,k)
     >                        +           u(m,i+2,j,k) )
               end do
            end do


            do m = 1, 5
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     >           - dssp * (             u(m,nx-4,j,k)
     >                      - 4.0d+00 * u(m,nx-3,j,k)
     >                      + 6.0d+00 * u(m,nx-2,j,k)
     >                      - 4.0d+00 * u(m,nx-1,j,k)  )
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     >           - dssp * (             u(m,nx-3,j,k)
     >                      - 4.0d+00 * u(m,nx-2,j,k)
     >                      + 5.0d+00 * u(m,nx-1,j,k) )
            end do

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

      if (timeron) call timer_start(t_rhsy)
c---------------------------------------------------------------------
c   eta-direction flux differences
c---------------------------------------------------------------------
      do k = 2, nz - 1
         do i = ist, iend
            do j = 1, ny
               flux(1,j) = u(3,i,j,k)
               u31 = u(3,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,j) = u(2,i,j,k) * u31 
               flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
               flux(4,j) = u(4,i,j,k) * u31
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
            end do

            do j = jst, jend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                   - ty2 * ( flux(m,j+1) - flux(m,j-1) )
               end do
            end do

            do j = jst, ny
               tmp = rho_i(i,j,k)

               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)

               tmp = rho_i(i,j-1,k)
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)

               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * ty3 * ( u31j**2 - u31jm1**2 )
     >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
            end do

            do j = jst, jend

               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dy1 * ty1 * (            u(1,i,j-1,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j+1,k) )

               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     >              + dy2 * ty1 * (            u(2,i,j-1,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j+1,k) )

               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     >              + dy3 * ty1 * (            u(3,i,j-1,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j+1,k) )

               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     >              + dy4 * ty1 * (            u(4,i,j-1,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j+1,k) )

               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     >              + dy5 * ty1 * (            u(5,i,j-1,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j+1,k) )

            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,2,k) = rsd(m,i,2,k)
     >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     >                      - 4.0d+00 * u(m,i,3,k)
     >                      +           u(m,i,4,k) )
               rsd(m,i,3,k) = rsd(m,i,3,k)
     >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     >                      + 6.0d+00 * u(m,i,3,k)
     >                      - 4.0d+00 * u(m,i,4,k)
     >                      +           u(m,i,5,k) )
            end do

            do j = 4, ny - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j-2,k)
     >                        - 4.0d+00 * u(m,i,j-1,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j+1,k)
     >                        +           u(m,i,j+2,k) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     >           - dssp * (             u(m,i,ny-4,k)
     >                      - 4.0d+00 * u(m,i,ny-3,k)
     >                      + 6.0d+00 * u(m,i,ny-2,k)
     >                      - 4.0d+00 * u(m,i,ny-1,k)  )
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     >           - dssp * (             u(m,i,ny-3,k)
     >                      - 4.0d+00 * u(m,i,ny-2,k)
     >                      + 5.0d+00 * u(m,i,ny-1,k) )
            end do

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

      if (timeron) call timer_start(t_rhsz)
c---------------------------------------------------------------------
c   zeta-direction flux differences
c---------------------------------------------------------------------
      do j = jst, jend
         do i = ist, iend
      	    do k = 1, nz
               flux(1,k) = u(4,i,j,k)
               u41 = u(4,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,k) = u(2,i,j,k) * u41 
               flux(3,k) = u(3,i,j,k) * u41 
               flux(4,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
               flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
            end do

            do k = 2, nz - 1
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                - tz2 * ( flux(m,k+1) - flux(m,k-1) )
               end do
            end do

            do k = 2, nz
               tmp = rho_i(i,j,k)

               u21k = tmp * u(2,i,j,k)
               u31k = tmp * u(3,i,j,k)
               u41k = tmp * u(4,i,j,k)
               u51k = tmp * u(5,i,j,k)

               tmp = rho_i(i,j,k-1)

               u21km1 = tmp * u(2,i,j,k-1)
               u31km1 = tmp * u(3,i,j,k-1)
               u41km1 = tmp * u(4,i,j,k-1)
               u51km1 = tmp * u(5,i,j,k-1)

               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tz3 * ( u41k**2 - u41km1**2 )
     >              + c1 * c5 * tz3 * ( u51k - u51km1 )
            end do

            do k = 2, nz - 1
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dz1 * tz1 * (            u(1,i,j,k-1)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j,k+1) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     >              + dz2 * tz1 * (            u(2,i,j,k-1)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j,k+1) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     >              + dz3 * tz1 * (            u(3,i,j,k-1)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j,k+1) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     >              + dz4 * tz1 * (            u(4,i,j,k-1)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j,k+1) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     >              + dz5 * tz1 * (            u(5,i,j,k-1)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j,k+1) )
            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,j,2) = rsd(m,i,j,2)
     >           - dssp * ( + 5.0d+00 * u(m,i,j,2)
     >                      - 4.0d+00 * u(m,i,j,3)
     >                      +           u(m,i,j,4) )
               rsd(m,i,j,3) = rsd(m,i,j,3)
     >           - dssp * ( - 4.0d+00 * u(m,i,j,2)
     >                      + 6.0d+00 * u(m,i,j,3)
     >                      - 4.0d+00 * u(m,i,j,4)
     >                      +           u(m,i,j,5) )
            end do

            do k = 4, nz - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j,k-2)
     >                        - 4.0d+00 * u(m,i,j,k-1)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j,k+1)
     >                        +           u(m,i,j,k+2) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
     >           - dssp * (             u(m,i,j,nz-4)
     >                      - 4.0d+00 * u(m,i,j,nz-3)
     >                      + 6.0d+00 * u(m,i,j,nz-2)
     >                      - 4.0d+00 * u(m,i,j,nz-1)  )
               rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
     >           - dssp * (             u(m,i,j,nz-3)
     >                      - 4.0d+00 * u(m,i,j,nz-2)
     >                      + 5.0d+00 * u(m,i,j,nz-1) )
            end do
         end do
      end do
      if (timeron) call timer_stop(t_rhsz)

      return
      end

invivorhs_rhs__302

Real time: 1.078415e+10 -- Predicted time: 1.148506e+10

Callcount: 251

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
201	### 1	251.134123944069	4.573276e+07	4.350898e+07	4.8

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

      subroutine blts ( ldmx, ldmy, ldmz,
     >                  nx, ny, nz, k,
     >                  omega,
     >                  v, tv,
     >                  ldz, ldy, ldx, d,
     >                  ist, iend, jst, jend,
     >                  nx0, ny0 )

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

c---------------------------------------------------------------------
c
c   compute the regular-sparse, block lower triangular solution:
c
c                     v <-- ( L-inv ) * v
c
c---------------------------------------------------------------------

      implicit none

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer ldmx, ldmy, ldmz
      integer nx, ny, nz
      integer k
      double precision  omega
c---------------------------------------------------------------------
c   To improve cache performance, second two dimensions padded by 1 
c   for even number sizes only.  Only needed in v.
c---------------------------------------------------------------------
      double precision  v( 5, ldmx/2*2+1, ldmy/2*2+1, *),
     >        tv( 5, ldmx/2*2+1, ldmy),
     >        ldz( 5, 5, ldmx/2*2+1, ldmy),
     >        ldy( 5, 5, ldmx/2*2+1, ldmy),
     >        ldx( 5, 5, ldmx/2*2+1, ldmy),
     >        d( 5, 5, ldmx/2*2+1, ldmy)
      integer ist, iend
      integer jst, jend
      integer nx0, ny0

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, m
      double precision  tmp, tmp1
      double precision  tmat(5,5)




      do j = jst, jend
         do i = ist, iend
            do m = 1, 5

                  tv( m, i, j ) =  v( m, i, j, k )
     >    - omega * (  ldz( m, 1, i, j ) * v( 1, i, j, k-1 )
     >               + ldz( m, 2, i, j ) * v( 2, i, j, k-1 )
     >               + ldz( m, 3, i, j ) * v( 3, i, j, k-1 )
     >               + ldz( m, 4, i, j ) * v( 4, i, j, k-1 )
     >               + ldz( m, 5, i, j ) * v( 5, i, j, k-1 )  )

            end do
         end do
      end do


      do j=jst,jend
        do i = ist, iend

            do m = 1, 5

                  tv( m, i, j ) =  tv( m, i, j )
     > - omega * ( ldy( m, 1, i, j ) * v( 1, i, j-1, k )
     >           + ldx( m, 1, i, j ) * v( 1, i-1, j, k )
     >           + ldy( m, 2, i, j ) * v( 2, i, j-1, k )
     >           + ldx( m, 2, i, j ) * v( 2, i-1, j, k )
     >           + ldy( m, 3, i, j ) * v( 3, i, j-1, k )
     >           + ldx( m, 3, i, j ) * v( 3, i-1, j, k )
     >           + ldy( m, 4, i, j ) * v( 4, i, j-1, k )
     >           + ldx( m, 4, i, j ) * v( 4, i-1, j, k )
     >           + ldy( m, 5, i, j ) * v( 5, i, j-1, k )
     >           + ldx( m, 5, i, j ) * v( 5, i-1, j, k ) )

            end do
       
c---------------------------------------------------------------------
c   diagonal block inversion
c
c   forward elimination
c---------------------------------------------------------------------
            do m = 1, 5
               tmat( m, 1 ) = d( m, 1, i, j )
               tmat( m, 2 ) = d( m, 2, i, j )
               tmat( m, 3 ) = d( m, 3, i, j )
               tmat( m, 4 ) = d( m, 4, i, j )
               tmat( m, 5 ) = d( m, 5, i, j )
            end do

            tmp1 = 1.0d+00 / tmat( 1, 1 )
            tmp = tmp1 * tmat( 2, 1 )
            tmat( 2, 2 ) =  tmat( 2, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 2, 3 ) =  tmat( 2, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 2, 4 ) =  tmat( 2, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 2, 5 ) =  tmat( 2, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 2, i, j ) = tv( 2, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 3, 1 )
            tmat( 3, 2 ) =  tmat( 3, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 1 )
            tmat( 4, 2 ) =  tmat( 4, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 1 )
            tmat( 5, 2 ) =  tmat( 5, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 1, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 2, 2 )
            tmp = tmp1 * tmat( 3, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 2, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 3, 3 )
            tmp = tmp1 * tmat( 4, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 3, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 3, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 4, 4 )
            tmp = tmp1 * tmat( 5, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 4, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 4, i, j ) * tmp

c---------------------------------------------------------------------
c   back substitution
c---------------------------------------------------------------------
            v( 5, i, j, k ) = tv( 5, i, j )
     >                      / tmat( 5, 5 )

            tv( 4, i, j ) = tv( 4, i, j )
     >           - tmat( 4, 5 ) * v( 5, i, j, k )
            v( 4, i, j, k ) = tv( 4, i, j )
     >                      / tmat( 4, 4 )

            tv( 3, i, j ) = tv( 3, i, j )
     >           - tmat( 3, 4 ) * v( 4, i, j, k )
     >           - tmat( 3, 5 ) * v( 5, i, j, k )
            v( 3, i, j, k ) = tv( 3, i, j )
     >                      / tmat( 3, 3 )

            tv( 2, i, j ) = tv( 2, i, j )
     >           - tmat( 2, 3 ) * v( 3, i, j, k )
     >           - tmat( 2, 4 ) * v( 4, i, j, k )
     >           - tmat( 2, 5 ) * v( 5, i, j, k )
            v( 2, i, j, k ) = tv( 2, i, j )
     >                      / tmat( 2, 2 )

            tv( 1, i, j ) = tv( 1, i, j )
     >           - tmat( 1, 2 ) * v( 2, i, j, k )
     >           - tmat( 1, 3 ) * v( 3, i, j, k )
     >           - tmat( 1, 4 ) * v( 4, i, j, k )
     >           - tmat( 1, 5 ) * v( 5, i, j, k )
            v( 1, i, j, k ) = tv( 1, i, j )
     >                      / tmat( 1, 1 )


        enddo
      enddo


      return
      end

invivoblts_blts__72

Real time: 2.450386e+10 -- Predicted time: 2.561577e+10

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
248	### 1	15524.3213228324	1.650041e+06	1.573956e+06	4.6

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

      subroutine jacu(k)

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

c---------------------------------------------------------------------
c   compute the upper triangular part of the jacobian matrix
c---------------------------------------------------------------------


      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer k

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j
      double precision  r43
      double precision  c1345
      double precision  c34
      double precision  tmp1, tmp2, tmp3



      r43 = ( 4.0d+00 / 3.0d+00 )
      c1345 = c1 * c3 * c4 * c5
      c34 = c3 * c4

         do j = jst, jend
            do i = ist, iend

c---------------------------------------------------------------------
c   form the block daigonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               d(1,1,i,j) =  1.0d+00
     >                       + dt * 2.0d+00 * (   tx1 * dx1
     >                                          + ty1 * dy1
     >                                          + tz1 * dz1 )
               d(1,2,i,j) =  0.0d+00
               d(1,3,i,j) =  0.0d+00
               d(1,4,i,j) =  0.0d+00
               d(1,5,i,j) =  0.0d+00

               d(2,1,i,j) =  dt * 2.0d+00
     >           * ( - tx1 * r43 - ty1 - tz1 )
     >           * ( c34 * tmp2 * u(2,i,j,k) )
               d(2,2,i,j) =  1.0d+00
     >          + dt * 2.0d+00 * c34 * tmp1 
     >          * (  tx1 * r43 + ty1 + tz1 )
     >          + dt * 2.0d+00 * (   tx1 * dx2
     >                             + ty1 * dy2
     >                             + tz1 * dz2  )
               d(2,3,i,j) = 0.0d+00
               d(2,4,i,j) = 0.0d+00
               d(2,5,i,j) = 0.0d+00

               d(3,1,i,j) = dt * 2.0d+00
     >           * ( - tx1 - ty1 * r43 - tz1 )
     >           * ( c34 * tmp2 * u(3,i,j,k) )
               d(3,2,i,j) = 0.0d+00
               d(3,3,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 * r43 + tz1 )
     >         + dt * 2.0d+00 * (  tx1 * dx3
     >                           + ty1 * dy3
     >                           + tz1 * dz3 )
               d(3,4,i,j) = 0.0d+00
               d(3,5,i,j) = 0.0d+00

               d(4,1,i,j) = dt * 2.0d+00
     >           * ( - tx1 - ty1 - tz1 * r43 )
     >           * ( c34 * tmp2 * u(4,i,j,k) )
               d(4,2,i,j) = 0.0d+00
               d(4,3,i,j) = 0.0d+00
               d(4,4,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 + tz1 * r43 )
     >         + dt * 2.0d+00 * (  tx1 * dx4
     >                           + ty1 * dy4
     >                           + tz1 * dz4 )
               d(4,5,i,j) = 0.0d+00

               d(5,1,i,j) = -dt * 2.0d+00
     >  * ( ( ( tx1 * ( r43*c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(2,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( r43*c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(3,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( r43*c34 - c1345 ) ) * ( u(4,i,j,k) ** 2 )
     >      ) * tmp3
     >   + ( tx1 + ty1 + tz1 ) * c1345 * tmp2 * u(5,i,j,k) )

               d(5,2,i,j) = dt * 2.0d+00
     > * ( tx1 * ( r43*c34 - c1345 )
     >   + ty1 * (     c34 - c1345 )
     >   + tz1 * (     c34 - c1345 ) ) * tmp2 * u(2,i,j,k)
               d(5,3,i,j) = dt * 2.0d+00
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( r43*c34 -c1345 )
     >   + tz1 * ( c34 - c1345 ) ) * tmp2 * u(3,i,j,k)
               d(5,4,i,j) = dt * 2.0d+00
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( c34 - c1345 )
     >   + tz1 * ( r43*c34 - c1345 ) ) * tmp2 * u(4,i,j,k)
               d(5,5,i,j) = 1.0d+00
     >   + dt * 2.0d+00 * ( tx1 + ty1 + tz1 ) * c1345 * tmp1
     >   + dt * 2.0d+00 * (  tx1 * dx5
     >                    +  ty1 * dy5
     >                    +  tz1 * dz5 )

c---------------------------------------------------------------------
c   form the first block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i+1,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               a(1,1,i,j) = - dt * tx1 * dx1
               a(1,2,i,j) =   dt * tx2
               a(1,3,i,j) =   0.0d+00
               a(1,4,i,j) =   0.0d+00
               a(1,5,i,j) =   0.0d+00

               a(2,1,i,j) =  dt * tx2
     >          * ( - ( u(2,i+1,j,k) * tmp1 ) ** 2
     >     + c2 * qs(i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i+1,j,k) )
               a(2,2,i,j) =  dt * tx2
     >          * ( ( 2.0d+00 - c2 ) * ( u(2,i+1,j,k) * tmp1 ) )
     >          - dt * tx1 * ( r43 * c34 * tmp1 )
     >          - dt * tx1 * dx2
               a(2,3,i,j) =  dt * tx2
     >              * ( - c2 * ( u(3,i+1,j,k) * tmp1 ) )
               a(2,4,i,j) =  dt * tx2
     >              * ( - c2 * ( u(4,i+1,j,k) * tmp1 ) )
               a(2,5,i,j) =  dt * tx2 * c2 

               a(3,1,i,j) =  dt * tx2
     >              * ( - ( u(2,i+1,j,k) * u(3,i+1,j,k) ) * tmp2 )
     >         - dt * tx1 * ( - c34 * tmp2 * u(3,i+1,j,k) )
               a(3,2,i,j) =  dt * tx2 * ( u(3,i+1,j,k) * tmp1 )
               a(3,3,i,j) =  dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx3
               a(3,4,i,j) = 0.0d+00
               a(3,5,i,j) = 0.0d+00

               a(4,1,i,j) = dt * tx2
     >          * ( - ( u(2,i+1,j,k)*u(4,i+1,j,k) ) * tmp2 )
     >          - dt * tx1 * ( - c34 * tmp2 * u(4,i+1,j,k) )
               a(4,2,i,j) = dt * tx2 * ( u(4,i+1,j,k) * tmp1 )
               a(4,3,i,j) = 0.0d+00
               a(4,4,i,j) = dt * tx2 * ( u(2,i+1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx4
               a(4,5,i,j) = 0.0d+00

               a(5,1,i,j) = dt * tx2
     >          * ( ( c2 * 2.0d0 * qs(i+1,j,k)
     >              - c1 * u(5,i+1,j,k) )
     >          * ( u(2,i+1,j,k) * tmp2 ) )
     >          - dt * tx1
     >          * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i+1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(3,i+1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(4,i+1,j,k)**2 )
     >              - c1345 * tmp2 * u(5,i+1,j,k) )
               a(5,2,i,j) = dt * tx2
     >          * ( c1 * ( u(5,i+1,j,k) * tmp1 )
     >             - c2
     >             * (  u(2,i+1,j,k)*u(2,i+1,j,k) * tmp2
     >                  + qs(i+1,j,k) * tmp1 ) )
     >           - dt * tx1
     >           * ( r43*c34 - c1345 ) * tmp2 * u(2,i+1,j,k)
               a(5,3,i,j) = dt * tx2
     >           * ( - c2 * ( u(3,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(3,i+1,j,k)
               a(5,4,i,j) = dt * tx2
     >           * ( - c2 * ( u(4,i+1,j,k)*u(2,i+1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(4,i+1,j,k)
               a(5,5,i,j) = dt * tx2
     >           * ( c1 * ( u(2,i+1,j,k) * tmp1 ) )
     >           - dt * tx1 * c1345 * tmp1
     >           - dt * tx1 * dx5

c---------------------------------------------------------------------
c   form the second block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j+1,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               b(1,1,i,j) = - dt * ty1 * dy1
               b(1,2,i,j) =   0.0d+00
               b(1,3,i,j) =  dt * ty2
               b(1,4,i,j) =   0.0d+00
               b(1,5,i,j) =   0.0d+00

               b(2,1,i,j) =  dt * ty2
     >           * ( - ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
     >           - dt * ty1 * ( - c34 * tmp2 * u(2,i,j+1,k) )
               b(2,2,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
     >          - dt * ty1 * ( c34 * tmp1 )
     >          - dt * ty1 * dy2
               b(2,3,i,j) =  dt * ty2 * ( u(2,i,j+1,k) * tmp1 )
               b(2,4,i,j) = 0.0d+00
               b(2,5,i,j) = 0.0d+00

               b(3,1,i,j) =  dt * ty2
     >           * ( - ( u(3,i,j+1,k) * tmp1 ) ** 2
     >      + c2 * ( qs(i,j+1,k) * tmp1 ) )
     >       - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j+1,k) )
               b(3,2,i,j) =  dt * ty2
     >                   * ( - c2 * ( u(2,i,j+1,k) * tmp1 ) )
               b(3,3,i,j) =  dt * ty2 * ( ( 2.0d+00 - c2 )
     >                   * ( u(3,i,j+1,k) * tmp1 ) )
     >       - dt * ty1 * ( r43 * c34 * tmp1 )
     >       - dt * ty1 * dy3
               b(3,4,i,j) =  dt * ty2
     >                   * ( - c2 * ( u(4,i,j+1,k) * tmp1 ) )
               b(3,5,i,j) =  dt * ty2 * c2

               b(4,1,i,j) =  dt * ty2
     >              * ( - ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
     >       - dt * ty1 * ( - c34 * tmp2 * u(4,i,j+1,k) )
               b(4,2,i,j) = 0.0d+00
               b(4,3,i,j) =  dt * ty2 * ( u(4,i,j+1,k) * tmp1 )
               b(4,4,i,j) =  dt * ty2 * ( u(3,i,j+1,k) * tmp1 )
     >                        - dt * ty1 * ( c34 * tmp1 )
     >                        - dt * ty1 * dy4
               b(4,5,i,j) = 0.0d+00

               b(5,1,i,j) =  dt * ty2
     >          * ( ( c2 * 2.0d0 * qs(i,j+1,k)
     >               - c1 * u(5,i,j+1,k) )
     >          * ( u(3,i,j+1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( - (     c34 - c1345 )*tmp3*(u(2,i,j+1,k)**2)
     >              - ( r43*c34 - c1345 )*tmp3*(u(3,i,j+1,k)**2)
     >              - (     c34 - c1345 )*tmp3*(u(4,i,j+1,k)**2)
     >              - c1345*tmp2*u(5,i,j+1,k) )
               b(5,2,i,j) =  dt * ty2
     >          * ( - c2 * ( u(2,i,j+1,k)*u(3,i,j+1,k) ) * tmp2 )
     >          - dt * ty1
     >          * ( c34 - c1345 ) * tmp2 * u(2,i,j+1,k)
               b(5,3,i,j) =  dt * ty2
     >          * ( c1 * ( u(5,i,j+1,k) * tmp1 )
     >          - c2 
     >          * ( qs(i,j+1,k) * tmp1
     >               + u(3,i,j+1,k)*u(3,i,j+1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j+1,k)
               b(5,4,i,j) =  dt * ty2
     >          * ( - c2 * ( u(3,i,j+1,k)*u(4,i,j+1,k) ) * tmp2 )
     >          - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j+1,k)
               b(5,5,i,j) =  dt * ty2
     >          * ( c1 * ( u(3,i,j+1,k) * tmp1 ) )
     >          - dt * ty1 * c1345 * tmp1
     >          - dt * ty1 * dy5

c---------------------------------------------------------------------
c   form the third block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k+1)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               c(1,1,i,j) = - dt * tz1 * dz1
               c(1,2,i,j) =   0.0d+00
               c(1,3,i,j) =   0.0d+00
               c(1,4,i,j) = dt * tz2
               c(1,5,i,j) =   0.0d+00

               c(2,1,i,j) = dt * tz2
     >           * ( - ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k+1) )
               c(2,2,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
     >           - dt * tz1 * c34 * tmp1
     >           - dt * tz1 * dz2 
               c(2,3,i,j) = 0.0d+00
               c(2,4,i,j) = dt * tz2 * ( u(2,i,j,k+1) * tmp1 )
               c(2,5,i,j) = 0.0d+00

               c(3,1,i,j) = dt * tz2
     >           * ( - ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k+1) )
               c(3,2,i,j) = 0.0d+00
               c(3,3,i,j) = dt * tz2 * ( u(4,i,j,k+1) * tmp1 )
     >           - dt * tz1 * ( c34 * tmp1 )
     >           - dt * tz1 * dz3
               c(3,4,i,j) = dt * tz2 * ( u(3,i,j,k+1) * tmp1 )
               c(3,5,i,j) = 0.0d+00

               c(4,1,i,j) = dt * tz2
     >        * ( - ( u(4,i,j,k+1) * tmp1 ) ** 2
     >            + c2 * ( qs(i,j,k+1) * tmp1 ) )
     >        - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k+1) )
               c(4,2,i,j) = dt * tz2
     >             * ( - c2 * ( u(2,i,j,k+1) * tmp1 ) )
               c(4,3,i,j) = dt * tz2
     >             * ( - c2 * ( u(3,i,j,k+1) * tmp1 ) )
               c(4,4,i,j) = dt * tz2 * ( 2.0d+00 - c2 )
     >             * ( u(4,i,j,k+1) * tmp1 )
     >             - dt * tz1 * ( r43 * c34 * tmp1 )
     >             - dt * tz1 * dz4
               c(4,5,i,j) = dt * tz2 * c2

               c(5,1,i,j) = dt * tz2
     >     * ( ( c2 * 2.0d0 * qs(i,j,k+1)
     >       - c1 * u(5,i,j,k+1) )
     >            * ( u(4,i,j,k+1) * tmp2 ) )
     >       - dt * tz1
     >       * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k+1)**2)
     >           - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k+1)**2)
     >           - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k+1)**2)
     >          - c1345 * tmp2 * u(5,i,j,k+1) )
               c(5,2,i,j) = dt * tz2
     >       * ( - c2 * ( u(2,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k+1)
               c(5,3,i,j) = dt * tz2
     >       * ( - c2 * ( u(3,i,j,k+1)*u(4,i,j,k+1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k+1)
               c(5,4,i,j) = dt * tz2
     >       * ( c1 * ( u(5,i,j,k+1) * tmp1 )
     >       - c2
     >       * ( qs(i,j,k+1) * tmp1
     >            + u(4,i,j,k+1)*u(4,i,j,k+1) * tmp2 ) )
     >       - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k+1)
               c(5,5,i,j) = dt * tz2
     >       * ( c1 * ( u(4,i,j,k+1) * tmp1 ) )
     >       - dt * tz1 * c1345 * tmp1
     >       - dt * tz1 * dz5

            end do
         end do

      return
      end

invivojacu_jacu__39

Real time: 1.539620e+10 -- Predicted time: 1.894338e+10

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
545	### 1	11030.0924812814	1.220802e+06	9.968720e+05	18.3
69	### 2	3951.75289592521	1.214169e+06	9.858680e+05	18.8
8679	### 3	557.594332547675	1.219035e+06	1.045408e+06	14.2

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

      subroutine blts ( ldmx, ldmy, ldmz,
     >                  nx, ny, nz, k,
     >                  omega,
     >                  v, tv,
     >                  ldz, ldy, ldx, d,
     >                  ist, iend, jst, jend,
     >                  nx0, ny0 )

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

c---------------------------------------------------------------------
c
c   compute the regular-sparse, block lower triangular solution:
c
c                     v <-- ( L-inv ) * v
c
c---------------------------------------------------------------------

      implicit none

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer ldmx, ldmy, ldmz
      integer nx, ny, nz
      integer k
      double precision  omega
c---------------------------------------------------------------------
c   To improve cache performance, second two dimensions padded by 1 
c   for even number sizes only.  Only needed in v.
c---------------------------------------------------------------------
      double precision  v( 5, ldmx/2*2+1, ldmy/2*2+1, *),
     >        tv( 5, ldmx/2*2+1, ldmy),
     >        ldz( 5, 5, ldmx/2*2+1, ldmy),
     >        ldy( 5, 5, ldmx/2*2+1, ldmy),
     >        ldx( 5, 5, ldmx/2*2+1, ldmy),
     >        d( 5, 5, ldmx/2*2+1, ldmy)
      integer ist, iend
      integer jst, jend
      integer nx0, ny0

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, m
      double precision  tmp, tmp1
      double precision  tmat(5,5)




      do j = jst, jend
         do i = ist, iend
            do m = 1, 5

                  tv( m, i, j ) =  v( m, i, j, k )
     >    - omega * (  ldz( m, 1, i, j ) * v( 1, i, j, k-1 )
     >               + ldz( m, 2, i, j ) * v( 2, i, j, k-1 )
     >               + ldz( m, 3, i, j ) * v( 3, i, j, k-1 )
     >               + ldz( m, 4, i, j ) * v( 4, i, j, k-1 )
     >               + ldz( m, 5, i, j ) * v( 5, i, j, k-1 )  )

            end do
         end do
      end do


      do j=jst,jend
        do i = ist, iend

            do m = 1, 5

                  tv( m, i, j ) =  tv( m, i, j )
     > - omega * ( ldy( m, 1, i, j ) * v( 1, i, j-1, k )
     >           + ldx( m, 1, i, j ) * v( 1, i-1, j, k )
     >           + ldy( m, 2, i, j ) * v( 2, i, j-1, k )
     >           + ldx( m, 2, i, j ) * v( 2, i-1, j, k )
     >           + ldy( m, 3, i, j ) * v( 3, i, j-1, k )
     >           + ldx( m, 3, i, j ) * v( 3, i-1, j, k )
     >           + ldy( m, 4, i, j ) * v( 4, i, j-1, k )
     >           + ldx( m, 4, i, j ) * v( 4, i-1, j, k )
     >           + ldy( m, 5, i, j ) * v( 5, i, j-1, k )
     >           + ldx( m, 5, i, j ) * v( 5, i-1, j, k ) )

            end do
       
c---------------------------------------------------------------------
c   diagonal block inversion
c
c   forward elimination
c---------------------------------------------------------------------
            do m = 1, 5
               tmat( m, 1 ) = d( m, 1, i, j )
               tmat( m, 2 ) = d( m, 2, i, j )
               tmat( m, 3 ) = d( m, 3, i, j )
               tmat( m, 4 ) = d( m, 4, i, j )
               tmat( m, 5 ) = d( m, 5, i, j )
            end do

            tmp1 = 1.0d+00 / tmat( 1, 1 )
            tmp = tmp1 * tmat( 2, 1 )
            tmat( 2, 2 ) =  tmat( 2, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 2, 3 ) =  tmat( 2, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 2, 4 ) =  tmat( 2, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 2, 5 ) =  tmat( 2, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 2, i, j ) = tv( 2, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 3, 1 )
            tmat( 3, 2 ) =  tmat( 3, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 1 )
            tmat( 4, 2 ) =  tmat( 4, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 1 )
            tmat( 5, 2 ) =  tmat( 5, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 1, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 2, 2 )
            tmp = tmp1 * tmat( 3, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 2, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 3, 3 )
            tmp = tmp1 * tmat( 4, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 3, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 3, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 4, 4 )
            tmp = tmp1 * tmat( 5, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 4, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 4, i, j ) * tmp

c---------------------------------------------------------------------
c   back substitution
c---------------------------------------------------------------------
            v( 5, i, j, k ) = tv( 5, i, j )
     >                      / tmat( 5, 5 )

            tv( 4, i, j ) = tv( 4, i, j )
     >           - tmat( 4, 5 ) * v( 5, i, j, k )
            v( 4, i, j, k ) = tv( 4, i, j )
     >                      / tmat( 4, 4 )

            tv( 3, i, j ) = tv( 3, i, j )
     >           - tmat( 3, 4 ) * v( 4, i, j, k )
     >           - tmat( 3, 5 ) * v( 5, i, j, k )
            v( 3, i, j, k ) = tv( 3, i, j )
     >                      / tmat( 3, 3 )

            tv( 2, i, j ) = tv( 2, i, j )
     >           - tmat( 2, 3 ) * v( 3, i, j, k )
     >           - tmat( 2, 4 ) * v( 4, i, j, k )
     >           - tmat( 2, 5 ) * v( 5, i, j, k )
            v( 2, i, j, k ) = tv( 2, i, j )
     >                      / tmat( 2, 2 )

            tv( 1, i, j ) = tv( 1, i, j )
     >           - tmat( 1, 2 ) * v( 2, i, j, k )
     >           - tmat( 1, 3 ) * v( 3, i, j, k )
     >           - tmat( 1, 4 ) * v( 4, i, j, k )
     >           - tmat( 1, 5 ) * v( 5, i, j, k )
            v( 1, i, j, k ) = tv( 1, i, j )
     >                      / tmat( 1, 1 )


        enddo
      enddo


      return
      end

invivoblts_blts__56

Real time: 2.757386e+09 -- Predicted time: 2.697139e+09

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
2791	### 1	264.980632129278	2.355170e+05	2.019840e+05	14.2
472	### 2	15260.8010393098	1.726470e+05	1.608760e+05	6.8

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

      subroutine rhs

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

c---------------------------------------------------------------------
c   compute the right hand sides
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m
      double precision  q
      double precision  tmp
      double precision  u21, u31, u41
      double precision  u21i, u31i, u41i, u51i
      double precision  u21j, u31j, u41j, u51j
      double precision  u21k, u31k, u41k, u51k
      double precision  u21im1, u31im1, u41im1, u51im1
      double precision  u21jm1, u31jm1, u41jm1, u51jm1
      double precision  u21km1, u31km1, u41km1, u51km1



      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
               tmp = 1.0d+00 / u(1,i,j,k)
               rho_i(i,j,k) = tmp
               qs(i,j,k) = 0.50d+00 * (  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) )
     >                      * tmp
            end do
         end do
      end do

      if (timeron) call timer_start(t_rhsx)
c---------------------------------------------------------------------
c   xi-direction flux differences
c---------------------------------------------------------------------

      do k = 2, nz - 1
         do j = jst, jend
            do i = 1, nx
               flux(1,i) = u(2,i,j,k)
               u21 = u(2,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,i) = u(2,i,j,k) * u21 + c2 * 
     >                        ( u(5,i,j,k) - q )
               flux(3,i) = u(3,i,j,k) * u21
               flux(4,i) = u(4,i,j,k) * u21
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
            end do

            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                 - tx2 * ( flux(m,i+1) - flux(m,i-1) )
               end do
            end do

            do i = ist, nx
               tmp = rho_i(i,j,k)

               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)

               tmp = rho_i(i-1,j,k)

               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)

               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tx3 * ( u21i**2 - u21im1**2 )
     >              + c1 * c5 * tx3 * ( u51i - u51im1 )
            end do

            do i = ist, iend
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dx1 * tx1 * (            u(1,i-1,j,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i+1,j,k) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     >              + dx2 * tx1 * (            u(2,i-1,j,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i+1,j,k) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     >              + dx3 * tx1 * (            u(3,i-1,j,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i+1,j,k) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     >              + dx4 * tx1 * (            u(4,i-1,j,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i+1,j,k) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     >              + dx5 * tx1 * (            u(5,i-1,j,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i+1,j,k) )
            end do

c---------------------------------------------------------------------
c   Fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,2,j,k) = rsd(m,2,j,k)
     >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     >                      - 4.0d+00 * u(m,3,j,k)
     >                      +           u(m,4,j,k) )
               rsd(m,3,j,k) = rsd(m,3,j,k)
     >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     >                      + 6.0d+00 * u(m,3,j,k)
     >                      - 4.0d+00 * u(m,4,j,k)
     >                      +           u(m,5,j,k) )
            end do

            do i = 4, nx - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i-2,j,k)
     >                        - 4.0d+00 * u(m,i-1,j,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i+1,j,k)
     >                        +           u(m,i+2,j,k) )
               end do
            end do


            do m = 1, 5
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     >           - dssp * (             u(m,nx-4,j,k)
     >                      - 4.0d+00 * u(m,nx-3,j,k)
     >                      + 6.0d+00 * u(m,nx-2,j,k)
     >                      - 4.0d+00 * u(m,nx-1,j,k)  )
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     >           - dssp * (             u(m,nx-3,j,k)
     >                      - 4.0d+00 * u(m,nx-2,j,k)
     >                      + 5.0d+00 * u(m,nx-1,j,k) )
            end do

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

      if (timeron) call timer_start(t_rhsy)
c---------------------------------------------------------------------
c   eta-direction flux differences
c---------------------------------------------------------------------
      do k = 2, nz - 1
         do i = ist, iend
            do j = 1, ny
               flux(1,j) = u(3,i,j,k)
               u31 = u(3,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,j) = u(2,i,j,k) * u31 
               flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
               flux(4,j) = u(4,i,j,k) * u31
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
            end do

            do j = jst, jend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                   - ty2 * ( flux(m,j+1) - flux(m,j-1) )
               end do
            end do

            do j = jst, ny
               tmp = rho_i(i,j,k)

               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)

               tmp = rho_i(i,j-1,k)
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)

               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * ty3 * ( u31j**2 - u31jm1**2 )
     >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
            end do

            do j = jst, jend

               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dy1 * ty1 * (            u(1,i,j-1,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j+1,k) )

               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     >              + dy2 * ty1 * (            u(2,i,j-1,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j+1,k) )

               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     >              + dy3 * ty1 * (            u(3,i,j-1,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j+1,k) )

               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     >              + dy4 * ty1 * (            u(4,i,j-1,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j+1,k) )

               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     >              + dy5 * ty1 * (            u(5,i,j-1,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j+1,k) )

            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,2,k) = rsd(m,i,2,k)
     >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     >                      - 4.0d+00 * u(m,i,3,k)
     >                      +           u(m,i,4,k) )
               rsd(m,i,3,k) = rsd(m,i,3,k)
     >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     >                      + 6.0d+00 * u(m,i,3,k)
     >                      - 4.0d+00 * u(m,i,4,k)
     >                      +           u(m,i,5,k) )
            end do

            do j = 4, ny - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j-2,k)
     >                        - 4.0d+00 * u(m,i,j-1,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j+1,k)
     >                        +           u(m,i,j+2,k) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     >           - dssp * (             u(m,i,ny-4,k)
     >                      - 4.0d+00 * u(m,i,ny-3,k)
     >                      + 6.0d+00 * u(m,i,ny-2,k)
     >                      - 4.0d+00 * u(m,i,ny-1,k)  )
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     >           - dssp * (             u(m,i,ny-3,k)
     >                      - 4.0d+00 * u(m,i,ny-2,k)
     >                      + 5.0d+00 * u(m,i,ny-1,k) )
            end do

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

      if (timeron) call timer_start(t_rhsz)
c---------------------------------------------------------------------
c   zeta-direction flux differences
c---------------------------------------------------------------------
      do j = jst, jend
         do i = ist, iend
      	    do k = 1, nz
               flux(1,k) = u(4,i,j,k)
               u41 = u(4,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,k) = u(2,i,j,k) * u41 
               flux(3,k) = u(3,i,j,k) * u41 
               flux(4,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
               flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
            end do

            do k = 2, nz - 1
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                - tz2 * ( flux(m,k+1) - flux(m,k-1) )
               end do
            end do

            do k = 2, nz
               tmp = rho_i(i,j,k)

               u21k = tmp * u(2,i,j,k)
               u31k = tmp * u(3,i,j,k)
               u41k = tmp * u(4,i,j,k)
               u51k = tmp * u(5,i,j,k)

               tmp = rho_i(i,j,k-1)

               u21km1 = tmp * u(2,i,j,k-1)
               u31km1 = tmp * u(3,i,j,k-1)
               u41km1 = tmp * u(4,i,j,k-1)
               u51km1 = tmp * u(5,i,j,k-1)

               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tz3 * ( u41k**2 - u41km1**2 )
     >              + c1 * c5 * tz3 * ( u51k - u51km1 )
            end do

            do k = 2, nz - 1
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dz1 * tz1 * (            u(1,i,j,k-1)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j,k+1) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     >              + dz2 * tz1 * (            u(2,i,j,k-1)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j,k+1) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     >              + dz3 * tz1 * (            u(3,i,j,k-1)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j,k+1) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     >              + dz4 * tz1 * (            u(4,i,j,k-1)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j,k+1) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     >              + dz5 * tz1 * (            u(5,i,j,k-1)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j,k+1) )
            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,j,2) = rsd(m,i,j,2)
     >           - dssp * ( + 5.0d+00 * u(m,i,j,2)
     >                      - 4.0d+00 * u(m,i,j,3)
     >                      +           u(m,i,j,4) )
               rsd(m,i,j,3) = rsd(m,i,j,3)
     >           - dssp * ( - 4.0d+00 * u(m,i,j,2)
     >                      + 6.0d+00 * u(m,i,j,3)
     >                      - 4.0d+00 * u(m,i,j,4)
     >                      +           u(m,i,j,5) )
            end do

            do k = 4, nz - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j,k-2)
     >                        - 4.0d+00 * u(m,i,j,k-1)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j,k+1)
     >                        +           u(m,i,j,k+2) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
     >           - dssp * (             u(m,i,j,nz-4)
     >                      - 4.0d+00 * u(m,i,j,nz-3)
     >                      + 6.0d+00 * u(m,i,j,nz-2)
     >                      - 4.0d+00 * u(m,i,j,nz-1)  )
               rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
     >           - dssp * (             u(m,i,j,nz-3)
     >                      - 4.0d+00 * u(m,i,j,nz-2)
     >                      + 5.0d+00 * u(m,i,j,nz-1) )
            end do
         end do
      end do
      if (timeron) call timer_stop(t_rhsz)

      return
      end

invivorhs_rhs__34

Real time: 2.288832e+09 -- Predicted time: 2.220453e+09

Callcount: 251

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	9.728824e+06	1.113942e+07	12.6
144	### 2	249.98773158285	8.843331e+06	9.225640e+06	4.1

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

      subroutine rhs

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

c---------------------------------------------------------------------
c   compute the right hand sides
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m
      double precision  q
      double precision  tmp
      double precision  u21, u31, u41
      double precision  u21i, u31i, u41i, u51i
      double precision  u21j, u31j, u41j, u51j
      double precision  u21k, u31k, u41k, u51k
      double precision  u21im1, u31im1, u41im1, u51im1
      double precision  u21jm1, u31jm1, u41jm1, u51jm1
      double precision  u21km1, u31km1, u41km1, u51km1



      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
               tmp = 1.0d+00 / u(1,i,j,k)
               rho_i(i,j,k) = tmp
               qs(i,j,k) = 0.50d+00 * (  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) )
     >                      * tmp
            end do
         end do
      end do

      if (timeron) call timer_start(t_rhsx)
c---------------------------------------------------------------------
c   xi-direction flux differences
c---------------------------------------------------------------------

      do k = 2, nz - 1
         do j = jst, jend
            do i = 1, nx
               flux(1,i) = u(2,i,j,k)
               u21 = u(2,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,i) = u(2,i,j,k) * u21 + c2 * 
     >                        ( u(5,i,j,k) - q )
               flux(3,i) = u(3,i,j,k) * u21
               flux(4,i) = u(4,i,j,k) * u21
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
            end do

            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                 - tx2 * ( flux(m,i+1) - flux(m,i-1) )
               end do
            end do

            do i = ist, nx
               tmp = rho_i(i,j,k)

               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)

               tmp = rho_i(i-1,j,k)

               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)

               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tx3 * ( u21i**2 - u21im1**2 )
     >              + c1 * c5 * tx3 * ( u51i - u51im1 )
            end do

            do i = ist, iend
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dx1 * tx1 * (            u(1,i-1,j,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i+1,j,k) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     >              + dx2 * tx1 * (            u(2,i-1,j,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i+1,j,k) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     >              + dx3 * tx1 * (            u(3,i-1,j,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i+1,j,k) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     >              + dx4 * tx1 * (            u(4,i-1,j,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i+1,j,k) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     >              + dx5 * tx1 * (            u(5,i-1,j,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i+1,j,k) )
            end do

c---------------------------------------------------------------------
c   Fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,2,j,k) = rsd(m,2,j,k)
     >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     >                      - 4.0d+00 * u(m,3,j,k)
     >                      +           u(m,4,j,k) )
               rsd(m,3,j,k) = rsd(m,3,j,k)
     >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     >                      + 6.0d+00 * u(m,3,j,k)
     >                      - 4.0d+00 * u(m,4,j,k)
     >                      +           u(m,5,j,k) )
            end do

            do i = 4, nx - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i-2,j,k)
     >                        - 4.0d+00 * u(m,i-1,j,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i+1,j,k)
     >                        +           u(m,i+2,j,k) )
               end do
            end do


            do m = 1, 5
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     >           - dssp * (             u(m,nx-4,j,k)
     >                      - 4.0d+00 * u(m,nx-3,j,k)
     >                      + 6.0d+00 * u(m,nx-2,j,k)
     >                      - 4.0d+00 * u(m,nx-1,j,k)  )
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     >           - dssp * (             u(m,nx-3,j,k)
     >                      - 4.0d+00 * u(m,nx-2,j,k)
     >                      + 5.0d+00 * u(m,nx-1,j,k) )
            end do

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

      if (timeron) call timer_start(t_rhsy)
c---------------------------------------------------------------------
c   eta-direction flux differences
c---------------------------------------------------------------------
      do k = 2, nz - 1
         do i = ist, iend
            do j = 1, ny
               flux(1,j) = u(3,i,j,k)
               u31 = u(3,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,j) = u(2,i,j,k) * u31 
               flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
               flux(4,j) = u(4,i,j,k) * u31
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
            end do

            do j = jst, jend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                   - ty2 * ( flux(m,j+1) - flux(m,j-1) )
               end do
            end do

            do j = jst, ny
               tmp = rho_i(i,j,k)

               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)

               tmp = rho_i(i,j-1,k)
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)

               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * ty3 * ( u31j**2 - u31jm1**2 )
     >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
            end do

            do j = jst, jend

               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dy1 * ty1 * (            u(1,i,j-1,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j+1,k) )

               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     >              + dy2 * ty1 * (            u(2,i,j-1,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j+1,k) )

               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     >              + dy3 * ty1 * (            u(3,i,j-1,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j+1,k) )

               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     >              + dy4 * ty1 * (            u(4,i,j-1,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j+1,k) )

               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     >              + dy5 * ty1 * (            u(5,i,j-1,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j+1,k) )

            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,2,k) = rsd(m,i,2,k)
     >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     >                      - 4.0d+00 * u(m,i,3,k)
     >                      +           u(m,i,4,k) )
               rsd(m,i,3,k) = rsd(m,i,3,k)
     >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     >                      + 6.0d+00 * u(m,i,3,k)
     >                      - 4.0d+00 * u(m,i,4,k)
     >                      +           u(m,i,5,k) )
            end do

            do j = 4, ny - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j-2,k)
     >                        - 4.0d+00 * u(m,i,j-1,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j+1,k)
     >                        +           u(m,i,j+2,k) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     >           - dssp * (             u(m,i,ny-4,k)
     >                      - 4.0d+00 * u(m,i,ny-3,k)
     >                      + 6.0d+00 * u(m,i,ny-2,k)
     >                      - 4.0d+00 * u(m,i,ny-1,k)  )
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     >           - dssp * (             u(m,i,ny-3,k)
     >                      - 4.0d+00 * u(m,i,ny-2,k)
     >                      + 5.0d+00 * u(m,i,ny-1,k) )
            end do

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

      if (timeron) call timer_start(t_rhsz)
c---------------------------------------------------------------------
c   zeta-direction flux differences
c---------------------------------------------------------------------
      do j = jst, jend
         do i = ist, iend
      	    do k = 1, nz
               flux(1,k) = u(4,i,j,k)
               u41 = u(4,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,k) = u(2,i,j,k) * u41 
               flux(3,k) = u(3,i,j,k) * u41 
               flux(4,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
               flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
            end do

            do k = 2, nz - 1
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                - tz2 * ( flux(m,k+1) - flux(m,k-1) )
               end do
            end do

            do k = 2, nz
               tmp = rho_i(i,j,k)

               u21k = tmp * u(2,i,j,k)
               u31k = tmp * u(3,i,j,k)
               u41k = tmp * u(4,i,j,k)
               u51k = tmp * u(5,i,j,k)

               tmp = rho_i(i,j,k-1)

               u21km1 = tmp * u(2,i,j,k-1)
               u31km1 = tmp * u(3,i,j,k-1)
               u41km1 = tmp * u(4,i,j,k-1)
               u51km1 = tmp * u(5,i,j,k-1)

               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tz3 * ( u41k**2 - u41km1**2 )
     >              + c1 * c5 * tz3 * ( u51k - u51km1 )
            end do

            do k = 2, nz - 1
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dz1 * tz1 * (            u(1,i,j,k-1)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j,k+1) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     >              + dz2 * tz1 * (            u(2,i,j,k-1)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j,k+1) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     >              + dz3 * tz1 * (            u(3,i,j,k-1)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j,k+1) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     >              + dz4 * tz1 * (            u(4,i,j,k-1)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j,k+1) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     >              + dz5 * tz1 * (            u(5,i,j,k-1)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j,k+1) )
            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,j,2) = rsd(m,i,j,2)
     >           - dssp * ( + 5.0d+00 * u(m,i,j,2)
     >                      - 4.0d+00 * u(m,i,j,3)
     >                      +           u(m,i,j,4) )
               rsd(m,i,j,3) = rsd(m,i,j,3)
     >           - dssp * ( - 4.0d+00 * u(m,i,j,2)
     >                      + 6.0d+00 * u(m,i,j,3)
     >                      - 4.0d+00 * u(m,i,j,4)
     >                      +           u(m,i,j,5) )
            end do

            do k = 4, nz - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j,k-2)
     >                        - 4.0d+00 * u(m,i,j,k-1)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j,k+1)
     >                        +           u(m,i,j,k+2) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
     >           - dssp * (             u(m,i,j,nz-4)
     >                      - 4.0d+00 * u(m,i,j,nz-3)
     >                      + 6.0d+00 * u(m,i,j,nz-2)
     >                      - 4.0d+00 * u(m,i,j,nz-1)  )
               rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
     >           - dssp * (             u(m,i,j,nz-3)
     >                      - 4.0d+00 * u(m,i,j,nz-2)
     >                      + 5.0d+00 * u(m,i,j,nz-1) )
            end do
         end do
      end do
      if (timeron) call timer_stop(t_rhsz)

      return
      end

invivorhs_rhs__55

Real time: 6.484518e+09 -- Predicted time: 6.622537e+09

Callcount: 251

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
77	### 1	251.162278410786	2.636756e+07	2.514834e+07	4.6

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

      subroutine rhs

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

c---------------------------------------------------------------------
c   compute the right hand sides
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, k, m
      double precision  q
      double precision  tmp
      double precision  u21, u31, u41
      double precision  u21i, u31i, u41i, u51i
      double precision  u21j, u31j, u41j, u51j
      double precision  u21k, u31k, u41k, u51k
      double precision  u21im1, u31im1, u41im1, u51im1
      double precision  u21jm1, u31jm1, u41jm1, u51jm1
      double precision  u21km1, u31km1, u41km1, u51km1



      do k = 1, nz
         do j = 1, ny
            do i = 1, nx
               do m = 1, 5
                  rsd(m,i,j,k) = - frct(m,i,j,k)
               end do
               tmp = 1.0d+00 / u(1,i,j,k)
               rho_i(i,j,k) = tmp
               qs(i,j,k) = 0.50d+00 * (  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) )
     >                      * tmp
            end do
         end do
      end do

      if (timeron) call timer_start(t_rhsx)
c---------------------------------------------------------------------
c   xi-direction flux differences
c---------------------------------------------------------------------

      do k = 2, nz - 1
         do j = jst, jend
            do i = 1, nx
               flux(1,i) = u(2,i,j,k)
               u21 = u(2,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,i) = u(2,i,j,k) * u21 + c2 * 
     >                        ( u(5,i,j,k) - q )
               flux(3,i) = u(3,i,j,k) * u21
               flux(4,i) = u(4,i,j,k) * u21
               flux(5,i) = ( c1 * u(5,i,j,k) - c2 * q ) * u21
            end do

            do i = ist, iend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                 - tx2 * ( flux(m,i+1) - flux(m,i-1) )
               end do
            end do

            do i = ist, nx
               tmp = rho_i(i,j,k)

               u21i = tmp * u(2,i,j,k)
               u31i = tmp * u(3,i,j,k)
               u41i = tmp * u(4,i,j,k)
               u51i = tmp * u(5,i,j,k)

               tmp = rho_i(i-1,j,k)

               u21im1 = tmp * u(2,i-1,j,k)
               u31im1 = tmp * u(3,i-1,j,k)
               u41im1 = tmp * u(4,i-1,j,k)
               u51im1 = tmp * u(5,i-1,j,k)

               flux(2,i) = (4.0d+00/3.0d+00) * tx3 * (u21i-u21im1)
               flux(3,i) = tx3 * ( u31i - u31im1 )
               flux(4,i) = tx3 * ( u41i - u41im1 )
               flux(5,i) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tx3 * ( ( u21i  **2 + u31i  **2 + u41i  **2 )
     >                      - ( u21im1**2 + u31im1**2 + u41im1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tx3 * ( u21i**2 - u21im1**2 )
     >              + c1 * c5 * tx3 * ( u51i - u51im1 )
            end do

            do i = ist, iend
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dx1 * tx1 * (            u(1,i-1,j,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i+1,j,k) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(2,i+1) - flux(2,i) )
     >              + dx2 * tx1 * (            u(2,i-1,j,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i+1,j,k) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(3,i+1) - flux(3,i) )
     >              + dx3 * tx1 * (            u(3,i-1,j,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i+1,j,k) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(4,i+1) - flux(4,i) )
     >              + dx4 * tx1 * (            u(4,i-1,j,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i+1,j,k) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tx3 * c3 * c4 * ( flux(5,i+1) - flux(5,i) )
     >              + dx5 * tx1 * (            u(5,i-1,j,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i+1,j,k) )
            end do

c---------------------------------------------------------------------
c   Fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,2,j,k) = rsd(m,2,j,k)
     >           - dssp * ( + 5.0d+00 * u(m,2,j,k)
     >                      - 4.0d+00 * u(m,3,j,k)
     >                      +           u(m,4,j,k) )
               rsd(m,3,j,k) = rsd(m,3,j,k)
     >           - dssp * ( - 4.0d+00 * u(m,2,j,k)
     >                      + 6.0d+00 * u(m,3,j,k)
     >                      - 4.0d+00 * u(m,4,j,k)
     >                      +           u(m,5,j,k) )
            end do

            do i = 4, nx - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i-2,j,k)
     >                        - 4.0d+00 * u(m,i-1,j,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i+1,j,k)
     >                        +           u(m,i+2,j,k) )
               end do
            end do


            do m = 1, 5
               rsd(m,nx-2,j,k) = rsd(m,nx-2,j,k)
     >           - dssp * (             u(m,nx-4,j,k)
     >                      - 4.0d+00 * u(m,nx-3,j,k)
     >                      + 6.0d+00 * u(m,nx-2,j,k)
     >                      - 4.0d+00 * u(m,nx-1,j,k)  )
               rsd(m,nx-1,j,k) = rsd(m,nx-1,j,k)
     >           - dssp * (             u(m,nx-3,j,k)
     >                      - 4.0d+00 * u(m,nx-2,j,k)
     >                      + 5.0d+00 * u(m,nx-1,j,k) )
            end do

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

      if (timeron) call timer_start(t_rhsy)
c---------------------------------------------------------------------
c   eta-direction flux differences
c---------------------------------------------------------------------
      do k = 2, nz - 1
         do i = ist, iend
            do j = 1, ny
               flux(1,j) = u(3,i,j,k)
               u31 = u(3,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,j) = u(2,i,j,k) * u31 
               flux(3,j) = u(3,i,j,k) * u31 + c2 * (u(5,i,j,k)-q)
               flux(4,j) = u(4,i,j,k) * u31
               flux(5,j) = ( c1 * u(5,i,j,k) - c2 * q ) * u31
            end do

            do j = jst, jend
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                   - ty2 * ( flux(m,j+1) - flux(m,j-1) )
               end do
            end do

            do j = jst, ny
               tmp = rho_i(i,j,k)

               u21j = tmp * u(2,i,j,k)
               u31j = tmp * u(3,i,j,k)
               u41j = tmp * u(4,i,j,k)
               u51j = tmp * u(5,i,j,k)

               tmp = rho_i(i,j-1,k)
               u21jm1 = tmp * u(2,i,j-1,k)
               u31jm1 = tmp * u(3,i,j-1,k)
               u41jm1 = tmp * u(4,i,j-1,k)
               u51jm1 = tmp * u(5,i,j-1,k)

               flux(2,j) = ty3 * ( u21j - u21jm1 )
               flux(3,j) = (4.0d+00/3.0d+00) * ty3 * (u31j-u31jm1)
               flux(4,j) = ty3 * ( u41j - u41jm1 )
               flux(5,j) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * ty3 * ( ( u21j  **2 + u31j  **2 + u41j  **2 )
     >                      - ( u21jm1**2 + u31jm1**2 + u41jm1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * ty3 * ( u31j**2 - u31jm1**2 )
     >              + c1 * c5 * ty3 * ( u51j - u51jm1 )
            end do

            do j = jst, jend

               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dy1 * ty1 * (            u(1,i,j-1,k)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j+1,k) )

               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(2,j+1) - flux(2,j) )
     >              + dy2 * ty1 * (            u(2,i,j-1,k)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j+1,k) )

               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(3,j+1) - flux(3,j) )
     >              + dy3 * ty1 * (            u(3,i,j-1,k)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j+1,k) )

               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(4,j+1) - flux(4,j) )
     >              + dy4 * ty1 * (            u(4,i,j-1,k)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j+1,k) )

               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + ty3 * c3 * c4 * ( flux(5,j+1) - flux(5,j) )
     >              + dy5 * ty1 * (            u(5,i,j-1,k)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j+1,k) )

            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,2,k) = rsd(m,i,2,k)
     >           - dssp * ( + 5.0d+00 * u(m,i,2,k)
     >                      - 4.0d+00 * u(m,i,3,k)
     >                      +           u(m,i,4,k) )
               rsd(m,i,3,k) = rsd(m,i,3,k)
     >           - dssp * ( - 4.0d+00 * u(m,i,2,k)
     >                      + 6.0d+00 * u(m,i,3,k)
     >                      - 4.0d+00 * u(m,i,4,k)
     >                      +           u(m,i,5,k) )
            end do

            do j = 4, ny - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j-2,k)
     >                        - 4.0d+00 * u(m,i,j-1,k)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j+1,k)
     >                        +           u(m,i,j+2,k) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,ny-2,k) = rsd(m,i,ny-2,k)
     >           - dssp * (             u(m,i,ny-4,k)
     >                      - 4.0d+00 * u(m,i,ny-3,k)
     >                      + 6.0d+00 * u(m,i,ny-2,k)
     >                      - 4.0d+00 * u(m,i,ny-1,k)  )
               rsd(m,i,ny-1,k) = rsd(m,i,ny-1,k)
     >           - dssp * (             u(m,i,ny-3,k)
     >                      - 4.0d+00 * u(m,i,ny-2,k)
     >                      + 5.0d+00 * u(m,i,ny-1,k) )
            end do

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

      if (timeron) call timer_start(t_rhsz)
c---------------------------------------------------------------------
c   zeta-direction flux differences
c---------------------------------------------------------------------
      do j = jst, jend
         do i = ist, iend
      	    do k = 1, nz
               flux(1,k) = u(4,i,j,k)
               u41 = u(4,i,j,k) * rho_i(i,j,k)

               q = qs(i,j,k)

               flux(2,k) = u(2,i,j,k) * u41 
               flux(3,k) = u(3,i,j,k) * u41 
               flux(4,k) = u(4,i,j,k) * u41 + c2 * (u(5,i,j,k)-q)
               flux(5,k) = ( c1 * u(5,i,j,k) - c2 * q ) * u41
            end do

            do k = 2, nz - 1
               do m = 1, 5
                  rsd(m,i,j,k) =  rsd(m,i,j,k)
     >                - tz2 * ( flux(m,k+1) - flux(m,k-1) )
               end do
            end do

            do k = 2, nz
               tmp = rho_i(i,j,k)

               u21k = tmp * u(2,i,j,k)
               u31k = tmp * u(3,i,j,k)
               u41k = tmp * u(4,i,j,k)
               u51k = tmp * u(5,i,j,k)

               tmp = rho_i(i,j,k-1)

               u21km1 = tmp * u(2,i,j,k-1)
               u31km1 = tmp * u(3,i,j,k-1)
               u41km1 = tmp * u(4,i,j,k-1)
               u51km1 = tmp * u(5,i,j,k-1)

               flux(2,k) = tz3 * ( u21k - u21km1 )
               flux(3,k) = tz3 * ( u31k - u31km1 )
               flux(4,k) = (4.0d+00/3.0d+00) * tz3 * (u41k-u41km1)
               flux(5,k) = 0.50d+00 * ( 1.0d+00 - c1*c5 )
     >              * tz3 * ( ( u21k  **2 + u31k  **2 + u41k  **2 )
     >                      - ( u21km1**2 + u31km1**2 + u41km1**2 ) )
     >              + (1.0d+00/6.0d+00)
     >              * tz3 * ( u41k**2 - u41km1**2 )
     >              + c1 * c5 * tz3 * ( u51k - u51km1 )
            end do

            do k = 2, nz - 1
               rsd(1,i,j,k) = rsd(1,i,j,k)
     >              + dz1 * tz1 * (            u(1,i,j,k-1)
     >                             - 2.0d+00 * u(1,i,j,k)
     >                             +           u(1,i,j,k+1) )
               rsd(2,i,j,k) = rsd(2,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(2,k+1) - flux(2,k) )
     >              + dz2 * tz1 * (            u(2,i,j,k-1)
     >                             - 2.0d+00 * u(2,i,j,k)
     >                             +           u(2,i,j,k+1) )
               rsd(3,i,j,k) = rsd(3,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(3,k+1) - flux(3,k) )
     >              + dz3 * tz1 * (            u(3,i,j,k-1)
     >                             - 2.0d+00 * u(3,i,j,k)
     >                             +           u(3,i,j,k+1) )
               rsd(4,i,j,k) = rsd(4,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(4,k+1) - flux(4,k) )
     >              + dz4 * tz1 * (            u(4,i,j,k-1)
     >                             - 2.0d+00 * u(4,i,j,k)
     >                             +           u(4,i,j,k+1) )
               rsd(5,i,j,k) = rsd(5,i,j,k)
     >          + tz3 * c3 * c4 * ( flux(5,k+1) - flux(5,k) )
     >              + dz5 * tz1 * (            u(5,i,j,k-1)
     >                             - 2.0d+00 * u(5,i,j,k)
     >                             +           u(5,i,j,k+1) )
            end do

c---------------------------------------------------------------------
c   fourth-order dissipation
c---------------------------------------------------------------------
            do m = 1, 5
               rsd(m,i,j,2) = rsd(m,i,j,2)
     >           - dssp * ( + 5.0d+00 * u(m,i,j,2)
     >                      - 4.0d+00 * u(m,i,j,3)
     >                      +           u(m,i,j,4) )
               rsd(m,i,j,3) = rsd(m,i,j,3)
     >           - dssp * ( - 4.0d+00 * u(m,i,j,2)
     >                      + 6.0d+00 * u(m,i,j,3)
     >                      - 4.0d+00 * u(m,i,j,4)
     >                      +           u(m,i,j,5) )
            end do

            do k = 4, nz - 3
               do m = 1, 5
                  rsd(m,i,j,k) = rsd(m,i,j,k)
     >              - dssp * (            u(m,i,j,k-2)
     >                        - 4.0d+00 * u(m,i,j,k-1)
     >                        + 6.0d+00 * u(m,i,j,k)
     >                        - 4.0d+00 * u(m,i,j,k+1)
     >                        +           u(m,i,j,k+2) )
               end do
            end do

            do m = 1, 5
               rsd(m,i,j,nz-2) = rsd(m,i,j,nz-2)
     >           - dssp * (             u(m,i,j,nz-4)
     >                      - 4.0d+00 * u(m,i,j,nz-3)
     >                      + 6.0d+00 * u(m,i,j,nz-2)
     >                      - 4.0d+00 * u(m,i,j,nz-1)  )
               rsd(m,i,j,nz-1) = rsd(m,i,j,nz-1)
     >           - dssp * (             u(m,i,j,nz-3)
     >                      - 4.0d+00 * u(m,i,j,nz-2)
     >                      + 5.0d+00 * u(m,i,j,nz-1) )
            end do
         end do
      end do
      if (timeron) call timer_stop(t_rhsz)

      return
      end

invivorhs_rhs__177

Real time: 1.091503e+10 -- Predicted time: 1.151312e+10

Callcount: 251

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
159	### 1	249.926441429168	4.589533e+07	4.363086e+07	4.9
123	### 2	1	4.266215e+07	4.440110e+07	3.9

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

      subroutine buts( ldmx, ldmy, ldmz,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 v, tv,
     >                 d, udx, udy, udz,
     >                 ist, iend, jst, jend,
     >                 nx0, ny0 )

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

c---------------------------------------------------------------------
c
c   compute the regular-sparse, block upper triangular solution:
c
c                     v <-- ( U-inv ) * v
c
c---------------------------------------------------------------------

      implicit none

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer ldmx, ldmy, ldmz
      integer nx, ny, nz
      integer k
      double precision  omega
c---------------------------------------------------------------------
c   To improve cache performance, second two dimensions padded by 1 
c   for even number sizes only.  Only needed in v.
c---------------------------------------------------------------------
      double precision  v( 5,ldmx/2*2+1, ldmy/2*2+1, *), 
     >        tv(5, ldmx/2*2+1, ldmy),
     >        d( 5, 5, ldmx/2*2+1, ldmy),
     >        udx( 5, 5, ldmx/2*2+1, ldmy),
     >        udy( 5, 5, ldmx/2*2+1, ldmy),
     >        udz( 5, 5, ldmx/2*2+1, ldmy )
      integer ist, iend
      integer jst, jend
      integer nx0, ny0

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, m
      double precision  tmp, tmp1
      double precision  tmat(5,5)



      do j = jend, jst, -1
         do i = iend, ist, -1
            do m = 1, 5
                  tv( m, i, j ) = 
     >      omega * (  udz( m, 1, i, j ) * v( 1, i, j, k+1 )
     >               + udz( m, 2, i, j ) * v( 2, i, j, k+1 )
     >               + udz( m, 3, i, j ) * v( 3, i, j, k+1 )
     >               + udz( m, 4, i, j ) * v( 4, i, j, k+1 )
     >               + udz( m, 5, i, j ) * v( 5, i, j, k+1 ) )
            end do
         end do
      end do


      do j = jend,jst,-1
        do i = iend,ist,-1

            do m = 1, 5
                  tv( m, i, j ) = tv( m, i, j )
     > + omega * ( udy( m, 1, i, j ) * v( 1, i, j+1, k )
     >           + udx( m, 1, i, j ) * v( 1, i+1, j, k )
     >           + udy( m, 2, i, j ) * v( 2, i, j+1, k )
     >           + udx( m, 2, i, j ) * v( 2, i+1, j, k )
     >           + udy( m, 3, i, j ) * v( 3, i, j+1, k )
     >           + udx( m, 3, i, j ) * v( 3, i+1, j, k )
     >           + udy( m, 4, i, j ) * v( 4, i, j+1, k )
     >           + udx( m, 4, i, j ) * v( 4, i+1, j, k )
     >           + udy( m, 5, i, j ) * v( 5, i, j+1, k )
     >           + udx( m, 5, i, j ) * v( 5, i+1, j, k ) )
            end do

c---------------------------------------------------------------------
c   diagonal block inversion
c---------------------------------------------------------------------
            do m = 1, 5
               tmat( m, 1 ) = d( m, 1, i, j )
               tmat( m, 2 ) = d( m, 2, i, j )
               tmat( m, 3 ) = d( m, 3, i, j )
               tmat( m, 4 ) = d( m, 4, i, j )
               tmat( m, 5 ) = d( m, 5, i, j )
            end do

            tmp1 = 1.0d+00 / tmat( 1, 1 )
            tmp = tmp1 * tmat( 2, 1 )
            tmat( 2, 2 ) =  tmat( 2, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 2, 3 ) =  tmat( 2, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 2, 4 ) =  tmat( 2, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 2, 5 ) =  tmat( 2, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 2, i, j ) = tv( 2, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 3, 1 )
            tmat( 3, 2 ) =  tmat( 3, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 1 )
            tmat( 4, 2 ) =  tmat( 4, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 1 )
            tmat( 5, 2 ) =  tmat( 5, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 1, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 2, 2 )
            tmp = tmp1 * tmat( 3, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 2, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 3, 3 )
            tmp = tmp1 * tmat( 4, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 3, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 3, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 4, 4 )
            tmp = tmp1 * tmat( 5, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 4, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 4, i, j ) * tmp

c---------------------------------------------------------------------
c   back substitution
c---------------------------------------------------------------------
            tv( 5, i, j ) = tv( 5, i, j )
     >                      / tmat( 5, 5 )

            tv( 4, i, j ) = tv( 4, i, j )
     >           - tmat( 4, 5 ) * tv( 5, i, j )
            tv( 4, i, j ) = tv( 4, i, j )
     >                      / tmat( 4, 4 )

            tv( 3, i, j ) = tv( 3, i, j )
     >           - tmat( 3, 4 ) * tv( 4, i, j )
     >           - tmat( 3, 5 ) * tv( 5, i, j )
            tv( 3, i, j ) = tv( 3, i, j )
     >                      / tmat( 3, 3 )

            tv( 2, i, j ) = tv( 2, i, j )
     >           - tmat( 2, 3 ) * tv( 3, i, j )
     >           - tmat( 2, 4 ) * tv( 4, i, j )
     >           - tmat( 2, 5 ) * tv( 5, i, j )
            tv( 2, i, j ) = tv( 2, i, j )
     >                      / tmat( 2, 2 )

            tv( 1, i, j ) = tv( 1, i, j )
     >           - tmat( 1, 2 ) * tv( 2, i, j )
     >           - tmat( 1, 3 ) * tv( 3, i, j )
     >           - tmat( 1, 4 ) * tv( 4, i, j )
     >           - tmat( 1, 5 ) * tv( 5, i, j )
            tv( 1, i, j ) = tv( 1, i, j )
     >                      / tmat( 1, 1 )

            v( 1, i, j, k ) = v( 1, i, j, k ) - tv( 1, i, j )
            v( 2, i, j, k ) = v( 2, i, j, k ) - tv( 2, i, j )
            v( 3, i, j, k ) = v( 3, i, j, k ) - tv( 3, i, j )
            v( 4, i, j, k ) = v( 4, i, j, k ) - tv( 4, i, j )
            v( 5, i, j, k ) = v( 5, i, j, k ) - tv( 5, i, j )

        enddo
      end do

 
      return
      end

invivobuts_buts__55

Real time: 2.604616e+09 -- Predicted time: 3.062381e+09

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
3174	### 1	433.247890219704	1.962320e+05	1.673160e+05	14.7
31	### 2	15076.8642127648	1.974790e+05	1.523560e+05	22.8

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

      subroutine jacld(k)

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


c---------------------------------------------------------------------
c   compute the lower triangular part of the jacobian matrix
c---------------------------------------------------------------------

      implicit none

      include 'applu.incl'

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer k

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j
      double precision  r43
      double precision  c1345
      double precision  c34
      double precision  tmp1, tmp2, tmp3



      r43 = ( 4.0d+00 / 3.0d+00 )
      c1345 = c1 * c3 * c4 * c5
      c34 = c3 * c4

         do j = jst, jend
            do i = ist, iend

c---------------------------------------------------------------------
c   form the block daigonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               d(1,1,i,j) =  1.0d+00
     >                       + dt * 2.0d+00 * (   tx1 * dx1
     >                                          + ty1 * dy1
     >                                          + tz1 * dz1 )
               d(1,2,i,j) =  0.0d+00
               d(1,3,i,j) =  0.0d+00
               d(1,4,i,j) =  0.0d+00
               d(1,5,i,j) =  0.0d+00

               d(2,1,i,j) = -dt * 2.0d+00
     >          * (  tx1 * r43 + ty1 + tz1  )
     >          * c34 * tmp2 * u(2,i,j,k)
               d(2,2,i,j) =  1.0d+00
     >          + dt * 2.0d+00 * c34 * tmp1 
     >          * (  tx1 * r43 + ty1 + tz1 )
     >          + dt * 2.0d+00 * (   tx1 * dx2
     >                             + ty1 * dy2
     >                             + tz1 * dz2  )
               d(2,3,i,j) = 0.0d+00
               d(2,4,i,j) = 0.0d+00
               d(2,5,i,j) = 0.0d+00

               d(3,1,i,j) = -dt * 2.0d+00
     >           * (  tx1 + ty1 * r43 + tz1  )
     >           * c34 * tmp2 * u(3,i,j,k)
               d(3,2,i,j) = 0.0d+00
               d(3,3,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 * r43 + tz1 )
     >         + dt * 2.0d+00 * (  tx1 * dx3
     >                           + ty1 * dy3
     >                           + tz1 * dz3 )
               d(3,4,i,j) = 0.0d+00
               d(3,5,i,j) = 0.0d+00

               d(4,1,i,j) = -dt * 2.0d+00
     >           * (  tx1 + ty1 + tz1 * r43  )
     >           * c34 * tmp2 * u(4,i,j,k)
               d(4,2,i,j) = 0.0d+00
               d(4,3,i,j) = 0.0d+00
               d(4,4,i,j) = 1.0d+00
     >         + dt * 2.0d+00 * c34 * tmp1
     >              * (  tx1 + ty1 + tz1 * r43 )
     >         + dt * 2.0d+00 * (  tx1 * dx4
     >                           + ty1 * dy4
     >                           + tz1 * dz4 )
               d(4,5,i,j) = 0.0d+00

               d(5,1,i,j) = -dt * 2.0d+00
     >  * ( ( ( tx1 * ( r43*c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(2,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( r43*c34 - c1345 )
     >     + tz1 * ( c34 - c1345 ) ) * ( u(3,i,j,k) ** 2 )
     >   + ( tx1 * ( c34 - c1345 )
     >     + ty1 * ( c34 - c1345 )
     >     + tz1 * ( r43*c34 - c1345 ) ) * ( u(4,i,j,k) ** 2 )
     >      ) * tmp3
     >   + ( tx1 + ty1 + tz1 ) * c1345 * tmp2 * u(5,i,j,k) )

               d(5,2,i,j) = dt * 2.0d+00 * tmp2 * u(2,i,j,k)
     > * ( tx1 * ( r43*c34 - c1345 )
     >   + ty1 * (     c34 - c1345 )
     >   + tz1 * (     c34 - c1345 ) )
               d(5,3,i,j) = dt * 2.0d+00 * tmp2 * u(3,i,j,k)
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( r43*c34 -c1345 )
     >   + tz1 * ( c34 - c1345 ) )
               d(5,4,i,j) = dt * 2.0d+00 * tmp2 * u(4,i,j,k)
     > * ( tx1 * ( c34 - c1345 )
     >   + ty1 * ( c34 - c1345 )
     >   + tz1 * ( r43*c34 - c1345 ) )
               d(5,5,i,j) = 1.0d+00
     >   + dt * 2.0d+00 * ( tx1  + ty1 + tz1 ) * c1345 * tmp1
     >   + dt * 2.0d+00 * (  tx1 * dx5
     >                    +  ty1 * dy5
     >                    +  tz1 * dz5 )

c---------------------------------------------------------------------
c   form the first block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j,k-1)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               a(1,1,i,j) = - dt * tz1 * dz1
               a(1,2,i,j) =   0.0d+00
               a(1,3,i,j) =   0.0d+00
               a(1,4,i,j) = - dt * tz2
               a(1,5,i,j) =   0.0d+00

               a(2,1,i,j) = - dt * tz2
     >           * ( - ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(2,i,j,k-1) )
               a(2,2,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
     >           - dt * tz1 * c34 * tmp1
     >           - dt * tz1 * dz2 
               a(2,3,i,j) = 0.0d+00
               a(2,4,i,j) = - dt * tz2 * ( u(2,i,j,k-1) * tmp1 )
               a(2,5,i,j) = 0.0d+00

               a(3,1,i,j) = - dt * tz2
     >           * ( - ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
     >           - dt * tz1 * ( - c34 * tmp2 * u(3,i,j,k-1) )
               a(3,2,i,j) = 0.0d+00
               a(3,3,i,j) = - dt * tz2 * ( u(4,i,j,k-1) * tmp1 )
     >           - dt * tz1 * ( c34 * tmp1 )
     >           - dt * tz1 * dz3
               a(3,4,i,j) = - dt * tz2 * ( u(3,i,j,k-1) * tmp1 )
               a(3,5,i,j) = 0.0d+00

               a(4,1,i,j) = - dt * tz2
     >        * ( - ( u(4,i,j,k-1) * tmp1 ) ** 2
     >            + c2 * qs(i,j,k-1) * tmp1 )
     >        - dt * tz1 * ( - r43 * c34 * tmp2 * u(4,i,j,k-1) )
               a(4,2,i,j) = - dt * tz2
     >             * ( - c2 * ( u(2,i,j,k-1) * tmp1 ) )
               a(4,3,i,j) = - dt * tz2
     >             * ( - c2 * ( u(3,i,j,k-1) * tmp1 ) )
               a(4,4,i,j) = - dt * tz2 * ( 2.0d+00 - c2 )
     >             * ( u(4,i,j,k-1) * tmp1 )
     >             - dt * tz1 * ( r43 * c34 * tmp1 )
     >             - dt * tz1 * dz4
               a(4,5,i,j) = - dt * tz2 * c2

               a(5,1,i,j) = - dt * tz2
     >       * ( ( c2 * 2.0d0 * qs(i,j,k-1)
     >       - c1 * u(5,i,j,k-1) )
     >            * u(4,i,j,k-1) * tmp2 )
     >       - dt * tz1
     >       * ( - ( c34 - c1345 ) * tmp3 * (u(2,i,j,k-1)**2)
     >           - ( c34 - c1345 ) * tmp3 * (u(3,i,j,k-1)**2)
     >           - ( r43*c34 - c1345 )* tmp3 * (u(4,i,j,k-1)**2)
     >          - c1345 * tmp2 * u(5,i,j,k-1) )
               a(5,2,i,j) = - dt * tz2
     >       * ( - c2 * ( u(2,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(2,i,j,k-1)
               a(5,3,i,j) = - dt * tz2
     >       * ( - c2 * ( u(3,i,j,k-1)*u(4,i,j,k-1) ) * tmp2 )
     >       - dt * tz1 * ( c34 - c1345 ) * tmp2 * u(3,i,j,k-1)
               a(5,4,i,j) = - dt * tz2
     >       * ( c1 * ( u(5,i,j,k-1) * tmp1 )
     >       - c2
     >       * ( qs(i,j,k-1) * tmp1
     >            + u(4,i,j,k-1)*u(4,i,j,k-1) * tmp2 ) )
     >       - dt * tz1 * ( r43*c34 - c1345 ) * tmp2 * u(4,i,j,k-1)
               a(5,5,i,j) = - dt * tz2
     >       * ( c1 * ( u(4,i,j,k-1) * tmp1 ) )
     >       - dt * tz1 * c1345 * tmp1
     >       - dt * tz1 * dz5

c---------------------------------------------------------------------
c   form the second block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i,j-1,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               b(1,1,i,j) = - dt * ty1 * dy1
               b(1,2,i,j) =   0.0d+00
               b(1,3,i,j) = - dt * ty2
               b(1,4,i,j) =   0.0d+00
               b(1,5,i,j) =   0.0d+00

               b(2,1,i,j) = - dt * ty2
     >           * ( - ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
     >           - dt * ty1 * ( - c34 * tmp2 * u(2,i,j-1,k) )
               b(2,2,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
     >          - dt * ty1 * ( c34 * tmp1 )
     >          - dt * ty1 * dy2
               b(2,3,i,j) = - dt * ty2 * ( u(2,i,j-1,k) * tmp1 )
               b(2,4,i,j) = 0.0d+00
               b(2,5,i,j) = 0.0d+00

               b(3,1,i,j) = - dt * ty2
     >           * ( - ( u(3,i,j-1,k) * tmp1 ) ** 2
     >       + c2 * ( qs(i,j-1,k) * tmp1 ) )
     >       - dt * ty1 * ( - r43 * c34 * tmp2 * u(3,i,j-1,k) )
               b(3,2,i,j) = - dt * ty2
     >                   * ( - c2 * ( u(2,i,j-1,k) * tmp1 ) )
               b(3,3,i,j) = - dt * ty2 * ( ( 2.0d+00 - c2 )
     >                   * ( u(3,i,j-1,k) * tmp1 ) )
     >       - dt * ty1 * ( r43 * c34 * tmp1 )
     >       - dt * ty1 * dy3
               b(3,4,i,j) = - dt * ty2
     >                   * ( - c2 * ( u(4,i,j-1,k) * tmp1 ) )
               b(3,5,i,j) = - dt * ty2 * c2

               b(4,1,i,j) = - dt * ty2
     >              * ( - ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
     >       - dt * ty1 * ( - c34 * tmp2 * u(4,i,j-1,k) )
               b(4,2,i,j) = 0.0d+00
               b(4,3,i,j) = - dt * ty2 * ( u(4,i,j-1,k) * tmp1 )
               b(4,4,i,j) = - dt * ty2 * ( u(3,i,j-1,k) * tmp1 )
     >                        - dt * ty1 * ( c34 * tmp1 )
     >                        - dt * ty1 * dy4
               b(4,5,i,j) = 0.0d+00

               b(5,1,i,j) = - dt * ty2
     >          * ( ( c2 * 2.0d0 * qs(i,j-1,k)
     >               - c1 * u(5,i,j-1,k) )
     >          * ( u(3,i,j-1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( - (     c34 - c1345 )*tmp3*(u(2,i,j-1,k)**2)
     >              - ( r43*c34 - c1345 )*tmp3*(u(3,i,j-1,k)**2)
     >              - (     c34 - c1345 )*tmp3*(u(4,i,j-1,k)**2)
     >              - c1345*tmp2*u(5,i,j-1,k) )
               b(5,2,i,j) = - dt * ty2
     >          * ( - c2 * ( u(2,i,j-1,k)*u(3,i,j-1,k) ) * tmp2 )
     >          - dt * ty1
     >          * ( c34 - c1345 ) * tmp2 * u(2,i,j-1,k)
               b(5,3,i,j) = - dt * ty2
     >          * ( c1 * ( u(5,i,j-1,k) * tmp1 )
     >          - c2 
     >          * ( qs(i,j-1,k) * tmp1
     >               + u(3,i,j-1,k)*u(3,i,j-1,k) * tmp2 ) )
     >          - dt * ty1
     >          * ( r43*c34 - c1345 ) * tmp2 * u(3,i,j-1,k)
               b(5,4,i,j) = - dt * ty2
     >          * ( - c2 * ( u(3,i,j-1,k)*u(4,i,j-1,k) ) * tmp2 )
     >          - dt * ty1 * ( c34 - c1345 ) * tmp2 * u(4,i,j-1,k)
               b(5,5,i,j) = - dt * ty2
     >          * ( c1 * ( u(3,i,j-1,k) * tmp1 ) )
     >          - dt * ty1 * c1345 * tmp1
     >          - dt * ty1 * dy5

c---------------------------------------------------------------------
c   form the third block sub-diagonal
c---------------------------------------------------------------------
               tmp1 = rho_i(i-1,j,k)
               tmp2 = tmp1 * tmp1
               tmp3 = tmp1 * tmp2

               c(1,1,i,j) = - dt * tx1 * dx1
               c(1,2,i,j) = - dt * tx2
               c(1,3,i,j) =   0.0d+00
               c(1,4,i,j) =   0.0d+00
               c(1,5,i,j) =   0.0d+00

               c(2,1,i,j) = - dt * tx2
     >          * ( - ( u(2,i-1,j,k) * tmp1 ) ** 2
     >       + c2 * qs(i-1,j,k) * tmp1 )
     >          - dt * tx1 * ( - r43 * c34 * tmp2 * u(2,i-1,j,k) )
               c(2,2,i,j) = - dt * tx2
     >          * ( ( 2.0d+00 - c2 ) * ( u(2,i-1,j,k) * tmp1 ) )
     >          - dt * tx1 * ( r43 * c34 * tmp1 )
     >          - dt * tx1 * dx2
               c(2,3,i,j) = - dt * tx2
     >              * ( - c2 * ( u(3,i-1,j,k) * tmp1 ) )
               c(2,4,i,j) = - dt * tx2
     >              * ( - c2 * ( u(4,i-1,j,k) * tmp1 ) )
               c(2,5,i,j) = - dt * tx2 * c2 

               c(3,1,i,j) = - dt * tx2
     >              * ( - ( u(2,i-1,j,k) * u(3,i-1,j,k) ) * tmp2 )
     >         - dt * tx1 * ( - c34 * tmp2 * u(3,i-1,j,k) )
               c(3,2,i,j) = - dt * tx2 * ( u(3,i-1,j,k) * tmp1 )
               c(3,3,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx3
               c(3,4,i,j) = 0.0d+00
               c(3,5,i,j) = 0.0d+00

               c(4,1,i,j) = - dt * tx2
     >          * ( - ( u(2,i-1,j,k)*u(4,i-1,j,k) ) * tmp2 )
     >          - dt * tx1 * ( - c34 * tmp2 * u(4,i-1,j,k) )
               c(4,2,i,j) = - dt * tx2 * ( u(4,i-1,j,k) * tmp1 )
               c(4,3,i,j) = 0.0d+00
               c(4,4,i,j) = - dt * tx2 * ( u(2,i-1,j,k) * tmp1 )
     >          - dt * tx1 * ( c34 * tmp1 )
     >          - dt * tx1 * dx4
               c(4,5,i,j) = 0.0d+00

               c(5,1,i,j) = - dt * tx2
     >          * ( ( c2 * 2.0d0 * qs(i-1,j,k)
     >              - c1 * u(5,i-1,j,k) )
     >          * u(2,i-1,j,k) * tmp2 )
     >          - dt * tx1
     >          * ( - ( r43*c34 - c1345 ) * tmp3 * ( u(2,i-1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(3,i-1,j,k)**2 )
     >              - (     c34 - c1345 ) * tmp3 * ( u(4,i-1,j,k)**2 )
     >              - c1345 * tmp2 * u(5,i-1,j,k) )
               c(5,2,i,j) = - dt * tx2
     >          * ( c1 * ( u(5,i-1,j,k) * tmp1 )
     >             - c2
     >             * ( u(2,i-1,j,k)*u(2,i-1,j,k) * tmp2
     >                  + qs(i-1,j,k) * tmp1 ) )
     >           - dt * tx1
     >           * ( r43*c34 - c1345 ) * tmp2 * u(2,i-1,j,k)
               c(5,3,i,j) = - dt * tx2
     >           * ( - c2 * ( u(3,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(3,i-1,j,k)
               c(5,4,i,j) = - dt * tx2
     >           * ( - c2 * ( u(4,i-1,j,k)*u(2,i-1,j,k) ) * tmp2 )
     >           - dt * tx1
     >           * (  c34 - c1345 ) * tmp2 * u(4,i-1,j,k)
               c(5,5,i,j) = - dt * tx2
     >           * ( c1 * ( u(2,i-1,j,k) * tmp1 ) )
     >           - dt * tx1 * c1345 * tmp1
     >           - dt * tx1 * dx5

            end do
         end do

      return
      end

invivojacld_jacld__39

Real time: 1.554055e+10 -- Predicted time: 1.974329e+10

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
6325	### 1	250.253528929396	1.271292e+06	1.278008e+06	0.5
2024	### 2	15271.7882375186	1.271963e+06	9.731280e+05	23.4

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

      subroutine buts( ldmx, ldmy, ldmz,
     >                 nx, ny, nz, k,
     >                 omega,
     >                 v, tv,
     >                 d, udx, udy, udz,
     >                 ist, iend, jst, jend,
     >                 nx0, ny0 )

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

c---------------------------------------------------------------------
c
c   compute the regular-sparse, block upper triangular solution:
c
c                     v <-- ( U-inv ) * v
c
c---------------------------------------------------------------------

      implicit none

c---------------------------------------------------------------------
c  input parameters
c---------------------------------------------------------------------
      integer ldmx, ldmy, ldmz
      integer nx, ny, nz
      integer k
      double precision  omega
c---------------------------------------------------------------------
c   To improve cache performance, second two dimensions padded by 1 
c   for even number sizes only.  Only needed in v.
c---------------------------------------------------------------------
      double precision  v( 5,ldmx/2*2+1, ldmy/2*2+1, *), 
     >        tv(5, ldmx/2*2+1, ldmy),
     >        d( 5, 5, ldmx/2*2+1, ldmy),
     >        udx( 5, 5, ldmx/2*2+1, ldmy),
     >        udy( 5, 5, ldmx/2*2+1, ldmy),
     >        udz( 5, 5, ldmx/2*2+1, ldmy )
      integer ist, iend
      integer jst, jend
      integer nx0, ny0

c---------------------------------------------------------------------
c  local variables
c---------------------------------------------------------------------
      integer i, j, m
      double precision  tmp, tmp1
      double precision  tmat(5,5)



      do j = jend, jst, -1
         do i = iend, ist, -1
            do m = 1, 5
                  tv( m, i, j ) = 
     >      omega * (  udz( m, 1, i, j ) * v( 1, i, j, k+1 )
     >               + udz( m, 2, i, j ) * v( 2, i, j, k+1 )
     >               + udz( m, 3, i, j ) * v( 3, i, j, k+1 )
     >               + udz( m, 4, i, j ) * v( 4, i, j, k+1 )
     >               + udz( m, 5, i, j ) * v( 5, i, j, k+1 ) )
            end do
         end do
      end do


      do j = jend,jst,-1
        do i = iend,ist,-1

            do m = 1, 5
                  tv( m, i, j ) = tv( m, i, j )
     > + omega * ( udy( m, 1, i, j ) * v( 1, i, j+1, k )
     >           + udx( m, 1, i, j ) * v( 1, i+1, j, k )
     >           + udy( m, 2, i, j ) * v( 2, i, j+1, k )
     >           + udx( m, 2, i, j ) * v( 2, i+1, j, k )
     >           + udy( m, 3, i, j ) * v( 3, i, j+1, k )
     >           + udx( m, 3, i, j ) * v( 3, i+1, j, k )
     >           + udy( m, 4, i, j ) * v( 4, i, j+1, k )
     >           + udx( m, 4, i, j ) * v( 4, i+1, j, k )
     >           + udy( m, 5, i, j ) * v( 5, i, j+1, k )
     >           + udx( m, 5, i, j ) * v( 5, i+1, j, k ) )
            end do

c---------------------------------------------------------------------
c   diagonal block inversion
c---------------------------------------------------------------------
            do m = 1, 5
               tmat( m, 1 ) = d( m, 1, i, j )
               tmat( m, 2 ) = d( m, 2, i, j )
               tmat( m, 3 ) = d( m, 3, i, j )
               tmat( m, 4 ) = d( m, 4, i, j )
               tmat( m, 5 ) = d( m, 5, i, j )
            end do

            tmp1 = 1.0d+00 / tmat( 1, 1 )
            tmp = tmp1 * tmat( 2, 1 )
            tmat( 2, 2 ) =  tmat( 2, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 2, 3 ) =  tmat( 2, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 2, 4 ) =  tmat( 2, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 2, 5 ) =  tmat( 2, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 2, i, j ) = tv( 2, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 3, 1 )
            tmat( 3, 2 ) =  tmat( 3, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 1 )
            tmat( 4, 2 ) =  tmat( 4, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 1, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 1 )
            tmat( 5, 2 ) =  tmat( 5, 2 )
     >           - tmp * tmat( 1, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 1, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 1, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 1, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 1, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 2, 2 )
            tmp = tmp1 * tmat( 3, 2 )
            tmat( 3, 3 ) =  tmat( 3, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 3, 4 ) =  tmat( 3, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 3, 5 ) =  tmat( 3, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 3, i, j ) = tv( 3, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 4, 2 )
            tmat( 4, 3 ) =  tmat( 4, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 2, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 2 )
            tmat( 5, 3 ) =  tmat( 5, 3 )
     >           - tmp * tmat( 2, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 2, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 2, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 2, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 3, 3 )
            tmp = tmp1 * tmat( 4, 3 )
            tmat( 4, 4 ) =  tmat( 4, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 4, 5 ) =  tmat( 4, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 4, i, j ) = tv( 4, i, j )
     >        - tv( 3, i, j ) * tmp

            tmp = tmp1 * tmat( 5, 3 )
            tmat( 5, 4 ) =  tmat( 5, 4 )
     >           - tmp * tmat( 3, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 3, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 3, i, j ) * tmp



            tmp1 = 1.0d+00 / tmat( 4, 4 )
            tmp = tmp1 * tmat( 5, 4 )
            tmat( 5, 5 ) =  tmat( 5, 5 )
     >           - tmp * tmat( 4, 5 )
            tv( 5, i, j ) = tv( 5, i, j )
     >        - tv( 4, i, j ) * tmp

c---------------------------------------------------------------------
c   back substitution
c---------------------------------------------------------------------
            tv( 5, i, j ) = tv( 5, i, j )
     >                      / tmat( 5, 5 )

            tv( 4, i, j ) = tv( 4, i, j )
     >           - tmat( 4, 5 ) * tv( 5, i, j )
            tv( 4, i, j ) = tv( 4, i, j )
     >                      / tmat( 4, 4 )

            tv( 3, i, j ) = tv( 3, i, j )
     >           - tmat( 3, 4 ) * tv( 4, i, j )
     >           - tmat( 3, 5 ) * tv( 5, i, j )
            tv( 3, i, j ) = tv( 3, i, j )
     >                      / tmat( 3, 3 )

            tv( 2, i, j ) = tv( 2, i, j )
     >           - tmat( 2, 3 ) * tv( 3, i, j )
     >           - tmat( 2, 4 ) * tv( 4, i, j )
     >           - tmat( 2, 5 ) * tv( 5, i, j )
            tv( 2, i, j ) = tv( 2, i, j )
     >                      / tmat( 2, 2 )

            tv( 1, i, j ) = tv( 1, i, j )
     >           - tmat( 1, 2 ) * tv( 2, i, j )
     >           - tmat( 1, 3 ) * tv( 3, i, j )
     >           - tmat( 1, 4 ) * tv( 4, i, j )
     >           - tmat( 1, 5 ) * tv( 5, i, j )
            tv( 1, i, j ) = tv( 1, i, j )
     >                      / tmat( 1, 1 )

            v( 1, i, j, k ) = v( 1, i, j, k ) - tv( 1, i, j )
            v( 2, i, j, k ) = v( 2, i, j, k ) - tv( 2, i, j )
            v( 3, i, j, k ) = v( 3, i, j, k ) - tv( 3, i, j )
            v( 4, i, j, k ) = v( 4, i, j, k ) - tv( 4, i, j )
            v( 5, i, j, k ) = v( 5, i, j, k ) - tv( 5, i, j )

        enddo
      end do

 
      return
      end

invivobuts_buts__69

Real time: 2.254505e+10 -- Predicted time: 2.486108e+10

Callcount: 15500

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
360	### 1	15522.0439121371	1.601663e+06	1.448620e+06	9.5

LU

Exec Time : 1.153398e+11 cycles

Matching : 71.4%

__invivo__rhs_rhs__302

Real time: 1.078415e+10 -- Predicted time: 1.148506e+10

Callcount: 251

__invivo__blts_blts__72

Real time: 2.450386e+10 -- Predicted time: 2.561577e+10

Callcount: 15500

__invivo__jacu_jacu__39

Real time: 1.539620e+10 -- Predicted time: 1.894338e+10

Callcount: 15500

__invivo__blts_blts__56

Real time: 2.757386e+09 -- Predicted time: 2.697139e+09

Callcount: 15500

__invivo__rhs_rhs__34

Real time: 2.288832e+09 -- Predicted time: 2.220453e+09

Callcount: 251

__invivo__rhs_rhs__55

Real time: 6.484518e+09 -- Predicted time: 6.622537e+09

Callcount: 251

__invivo__rhs_rhs__177

Real time: 1.091503e+10 -- Predicted time: 1.151312e+10

Callcount: 251

__invivo__buts_buts__55

Real time: 2.604616e+09 -- Predicted time: 3.062381e+09

Callcount: 15500

__invivo__jacld_jacld__39

Real time: 1.554055e+10 -- Predicted time: 1.974329e+10

Callcount: 15500

__invivo__buts_buts__69

Real time: 2.254505e+10 -- Predicted time: 2.486108e+10

Callcount: 15500

invivorhs_rhs__302

invivoblts_blts__72

invivojacu_jacu__39

invivoblts_blts__56

invivorhs_rhs__34

invivorhs_rhs__55

invivorhs_rhs__177

invivobuts_buts__55

invivojacld_jacld__39

invivobuts_buts__69