examples/smpi/NAS/SP/txinvr.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5        subroutine  txinvr
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c block-diagonal matrix-vector multiplication
  12 c---------------------------------------------------------------------
  13
  14        include 'header.h'
  15
  16        integer c, i, j, k
  17        double precision t1, t2, t3, ac, ru1, uu, vv, ww, r1, r2, r3,
  18      >                  r4, r5, ac2inv
  19
  20 c---------------------------------------------------------------------
  21 c      loop over all cells owned by this node
  22 c---------------------------------------------------------------------
  23        do   c = 1, ncells
  24           do    k = start(3,c), cell_size(3,c)-end(3,c)-1
  25              do    j = start(2,c), cell_size(2,c)-end(2,c)-1
  26                 do    i = start(1,c), cell_size(1,c)-end(1,c)-1
  27
  28                    ru1 = rho_i(i,j,k,c)
  29                    uu = us(i,j,k,c)
  30                    vv = vs(i,j,k,c)
  31                    ww = ws(i,j,k,c)
  32                    ac = speed(i,j,k,c)
  33                    ac2inv = ainv(i,j,k,c)*ainv(i,j,k,c)
  34
  35                    r1 = rhs(i,j,k,1,c)
  36                    r2 = rhs(i,j,k,2,c)
  37                    r3 = rhs(i,j,k,3,c)
  38                    r4 = rhs(i,j,k,4,c)
  39                    r5 = rhs(i,j,k,5,c)
  40
  41                    t1 = c2 * ac2inv * ( qs(i,j,k,c)*r1 - uu*r2  -
  42      >                  vv*r3 - ww*r4 + r5 )
  43                    t2 = bt * ru1 * ( uu * r1 - r2 )
  44                    t3 = ( bt * ru1 * ac ) * t1
  45
  46                    rhs(i,j,k,1,c) = r1 - t1
  47                    rhs(i,j,k,2,c) = - ru1 * ( ww*r1 - r4 )
  48                    rhs(i,j,k,3,c) =   ru1 * ( vv*r1 - r3 )
  49                    rhs(i,j,k,4,c) = - t2 + t3
  50                    rhs(i,j,k,5,c) =   t2 + t3
  51                 end do
  52              end do
  53           end do
  54        end do
  55
  56        return
  57        end
  58
  59