examples/smpi/NAS/LU/blts.f

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