examples/smpi/NAS/LU/buts.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4       subroutine buts( ldmx, ldmy, ldmz,
   5      >                 nx, ny, nz, k,
   6      >                 omega,
   7      >                 v, tv,
   8      >                 d, udx, udy, udz,
   9      >                 ist, iend, jst, jend,
  10      >                 nx0, ny0, ipt, jpt )
  11
  12 c---------------------------------------------------------------------
  13 c---------------------------------------------------------------------
  14
  15 c---------------------------------------------------------------------
  16 c
  17 c   compute the regular-sparse, block upper triangular solution:
  18 c
  19 c                     v <-- ( U-inv ) * v
  20 c
  21 c---------------------------------------------------------------------
  22
  23       implicit none
  24
  25 c---------------------------------------------------------------------
  26 c  input parameters
  27 c---------------------------------------------------------------------
  28       integer ldmx, ldmy, ldmz
  29       integer nx, ny, nz
  30       integer k
  31       double precision  omega
  32       double precision  v( 5, -1:ldmx+2, -1:ldmy+2, *),
  33      >        tv(5, ldmx, ldmy),
  34      >        d( 5, 5, ldmx, ldmy),
  35      >        udx( 5, 5, ldmx, ldmy),
  36      >        udy( 5, 5, ldmx, ldmy),
  37      >        udz( 5, 5, ldmx, ldmy )
  38       integer ist, iend
  39       integer jst, jend
  40       integer nx0, ny0
  41       integer ipt, jpt
  42
  43 c---------------------------------------------------------------------
  44 c  local variables
  45 c---------------------------------------------------------------------
  46       integer i, j, m
  47       integer iex
  48       double precision  tmp, tmp1
  49       double precision  tmat(5,5)
  50
  51
  52 c---------------------------------------------------------------------
  53 c   receive data from south and east
  54 c---------------------------------------------------------------------
  55       iex = 1
  56       call exchange_1( v,k,iex )
  57
  58       do j = jend, jst, -1
  59          do i = iend, ist, -1
  60             do m = 1, 5
  61                   tv( m, i, j ) =
  62      >      omega * (  udz( m, 1, i, j ) * v( 1, i, j, k+1 )
  63      >               + udz( m, 2, i, j ) * v( 2, i, j, k+1 )
  64      >               + udz( m, 3, i, j ) * v( 3, i, j, k+1 )
  65      >               + udz( m, 4, i, j ) * v( 4, i, j, k+1 )
  66      >               + udz( m, 5, i, j ) * v( 5, i, j, k+1 ) )
  67             end do
  68          end do
  69       end do
  70
  71
  72       do j = jend,jst,-1
  73         do i = iend,ist,-1
  74
  75             do m = 1, 5
  76                   tv( m, i, j ) = tv( m, i, j )
  77      > + omega * ( udy( m, 1, i, j ) * v( 1, i, j+1, k )
  78      >           + udx( m, 1, i, j ) * v( 1, i+1, j, k )
  79      >           + udy( m, 2, i, j ) * v( 2, i, j+1, k )
  80      >           + udx( m, 2, i, j ) * v( 2, i+1, j, k )
  81      >           + udy( m, 3, i, j ) * v( 3, i, j+1, k )
  82      >           + udx( m, 3, i, j ) * v( 3, i+1, j, k )
  83      >           + udy( m, 4, i, j ) * v( 4, i, j+1, k )
  84      >           + udx( m, 4, i, j ) * v( 4, i+1, j, k )
  85      >           + udy( m, 5, i, j ) * v( 5, i, j+1, k )
  86      >           + udx( m, 5, i, j ) * v( 5, i+1, j, k ) )
  87             end do
  88
  89 c---------------------------------------------------------------------
  90 c   diagonal block inversion
  91 c---------------------------------------------------------------------
  92             do m = 1, 5
  93                tmat( m, 1 ) = d( m, 1, i, j )
  94                tmat( m, 2 ) = d( m, 2, i, j )
  95                tmat( m, 3 ) = d( m, 3, i, j )
  96                tmat( m, 4 ) = d( m, 4, i, j )
  97                tmat( m, 5 ) = d( m, 5, i, j )
  98             end do
  99
 100             tmp1 = 1.0d+00 / tmat( 1, 1 )
 101             tmp = tmp1 * tmat( 2, 1 )
 102             tmat( 2, 2 ) =  tmat( 2, 2 )
 103      >           - tmp * tmat( 1, 2 )
 104             tmat( 2, 3 ) =  tmat( 2, 3 )
 105      >           - tmp * tmat( 1, 3 )
 106             tmat( 2, 4 ) =  tmat( 2, 4 )
 107      >           - tmp * tmat( 1, 4 )
 108             tmat( 2, 5 ) =  tmat( 2, 5 )
 109      >           - tmp * tmat( 1, 5 )
 110             tv( 2, i, j ) = tv( 2, i, j )
 111      >        - tv( 1, i, j ) * tmp
 112
 113             tmp = tmp1 * tmat( 3, 1 )
 114             tmat( 3, 2 ) =  tmat( 3, 2 )
 115      >           - tmp * tmat( 1, 2 )
 116             tmat( 3, 3 ) =  tmat( 3, 3 )
 117      >           - tmp * tmat( 1, 3 )
 118             tmat( 3, 4 ) =  tmat( 3, 4 )
 119      >           - tmp * tmat( 1, 4 )
 120             tmat( 3, 5 ) =  tmat( 3, 5 )
 121      >           - tmp * tmat( 1, 5 )
 122             tv( 3, i, j ) = tv( 3, i, j )
 123      >        - tv( 1, i, j ) * tmp
 124
 125             tmp = tmp1 * tmat( 4, 1 )
 126             tmat( 4, 2 ) =  tmat( 4, 2 )
 127      >           - tmp * tmat( 1, 2 )
 128             tmat( 4, 3 ) =  tmat( 4, 3 )
 129      >           - tmp * tmat( 1, 3 )
 130             tmat( 4, 4 ) =  tmat( 4, 4 )
 131      >           - tmp * tmat( 1, 4 )
 132             tmat( 4, 5 ) =  tmat( 4, 5 )
 133      >           - tmp * tmat( 1, 5 )
 134             tv( 4, i, j ) = tv( 4, i, j )
 135      >        - tv( 1, i, j ) * tmp
 136
 137             tmp = tmp1 * tmat( 5, 1 )
 138             tmat( 5, 2 ) =  tmat( 5, 2 )
 139      >           - tmp * tmat( 1, 2 )
 140             tmat( 5, 3 ) =  tmat( 5, 3 )
 141      >           - tmp * tmat( 1, 3 )
 142             tmat( 5, 4 ) =  tmat( 5, 4 )
 143      >           - tmp * tmat( 1, 4 )
 144             tmat( 5, 5 ) =  tmat( 5, 5 )
 145      >           - tmp * tmat( 1, 5 )
 146             tv( 5, i, j ) = tv( 5, i, j )
 147      >        - tv( 1, i, j ) * tmp
 148
 149
 150
 151             tmp1 = 1.0d+00 / tmat( 2, 2 )
 152             tmp = tmp1 * tmat( 3, 2 )
 153             tmat( 3, 3 ) =  tmat( 3, 3 )
 154      >           - tmp * tmat( 2, 3 )
 155             tmat( 3, 4 ) =  tmat( 3, 4 )
 156      >           - tmp * tmat( 2, 4 )
 157             tmat( 3, 5 ) =  tmat( 3, 5 )
 158      >           - tmp * tmat( 2, 5 )
 159             tv( 3, i, j ) = tv( 3, i, j )
 160      >        - tv( 2, i, j ) * tmp
 161
 162             tmp = tmp1 * tmat( 4, 2 )
 163             tmat( 4, 3 ) =  tmat( 4, 3 )
 164      >           - tmp * tmat( 2, 3 )
 165             tmat( 4, 4 ) =  tmat( 4, 4 )
 166      >           - tmp * tmat( 2, 4 )
 167             tmat( 4, 5 ) =  tmat( 4, 5 )
 168      >           - tmp * tmat( 2, 5 )
 169             tv( 4, i, j ) = tv( 4, i, j )
 170      >        - tv( 2, i, j ) * tmp
 171
 172             tmp = tmp1 * tmat( 5, 2 )
 173             tmat( 5, 3 ) =  tmat( 5, 3 )
 174      >           - tmp * tmat( 2, 3 )
 175             tmat( 5, 4 ) =  tmat( 5, 4 )
 176      >           - tmp * tmat( 2, 4 )
 177             tmat( 5, 5 ) =  tmat( 5, 5 )
 178      >           - tmp * tmat( 2, 5 )
 179             tv( 5, i, j ) = tv( 5, i, j )
 180      >        - tv( 2, i, j ) * tmp
 181
 182
 183
 184             tmp1 = 1.0d+00 / tmat( 3, 3 )
 185             tmp = tmp1 * tmat( 4, 3 )
 186             tmat( 4, 4 ) =  tmat( 4, 4 )
 187      >           - tmp * tmat( 3, 4 )
 188             tmat( 4, 5 ) =  tmat( 4, 5 )
 189      >           - tmp * tmat( 3, 5 )
 190             tv( 4, i, j ) = tv( 4, i, j )
 191      >        - tv( 3, i, j ) * tmp
 192
 193             tmp = tmp1 * tmat( 5, 3 )
 194             tmat( 5, 4 ) =  tmat( 5, 4 )
 195      >           - tmp * tmat( 3, 4 )
 196             tmat( 5, 5 ) =  tmat( 5, 5 )
 197      >           - tmp * tmat( 3, 5 )
 198             tv( 5, i, j ) = tv( 5, i, j )
 199      >        - tv( 3, i, j ) * tmp
 200
 201
 202
 203             tmp1 = 1.0d+00 / tmat( 4, 4 )
 204             tmp = tmp1 * tmat( 5, 4 )
 205             tmat( 5, 5 ) =  tmat( 5, 5 )
 206      >           - tmp * tmat( 4, 5 )
 207             tv( 5, i, j ) = tv( 5, i, j )
 208      >        - tv( 4, i, j ) * tmp
 209
 210 c---------------------------------------------------------------------
 211 c   back substitution
 212 c---------------------------------------------------------------------
 213             tv( 5, i, j ) = tv( 5, i, j )
 214      >                      / tmat( 5, 5 )
 215
 216             tv( 4, i, j ) = tv( 4, i, j )
 217      >           - tmat( 4, 5 ) * tv( 5, i, j )
 218             tv( 4, i, j ) = tv( 4, i, j )
 219      >                      / tmat( 4, 4 )
 220
 221             tv( 3, i, j ) = tv( 3, i, j )
 222      >           - tmat( 3, 4 ) * tv( 4, i, j )
 223      >           - tmat( 3, 5 ) * tv( 5, i, j )
 224             tv( 3, i, j ) = tv( 3, i, j )
 225      >                      / tmat( 3, 3 )
 226
 227             tv( 2, i, j ) = tv( 2, i, j )
 228      >           - tmat( 2, 3 ) * tv( 3, i, j )
 229      >           - tmat( 2, 4 ) * tv( 4, i, j )
 230      >           - tmat( 2, 5 ) * tv( 5, i, j )
 231             tv( 2, i, j ) = tv( 2, i, j )
 232      >                      / tmat( 2, 2 )
 233
 234             tv( 1, i, j ) = tv( 1, i, j )
 235      >           - tmat( 1, 2 ) * tv( 2, i, j )
 236      >           - tmat( 1, 3 ) * tv( 3, i, j )
 237      >           - tmat( 1, 4 ) * tv( 4, i, j )
 238      >           - tmat( 1, 5 ) * tv( 5, i, j )
 239             tv( 1, i, j ) = tv( 1, i, j )
 240      >                      / tmat( 1, 1 )
 241
 242             v( 1, i, j, k ) = v( 1, i, j, k ) - tv( 1, i, j )
 243             v( 2, i, j, k ) = v( 2, i, j, k ) - tv( 2, i, j )
 244             v( 3, i, j, k ) = v( 3, i, j, k ) - tv( 3, i, j )
 245             v( 4, i, j, k ) = v( 4, i, j, k ) - tv( 4, i, j )
 246             v( 5, i, j, k ) = v( 5, i, j, k ) - tv( 5, i, j )
 247
 248
 249         enddo
 250       end do
 251
 252 c---------------------------------------------------------------------
 253 c   send data to north and west
 254 c---------------------------------------------------------------------
 255       iex = 3
 256       call exchange_1( v,k,iex )
 257
 258       return
 259       end