blob: b86e31301e78d25351678e3b91bb8fe81d567848 [file] [log] [blame]
c---------------------------------------------------------------------
c---------------------------------------------------------------------
subroutine buts( ldmx, ldmy, ldmz,
> nx, ny, nz, k,
> omega,
> v, tv,
> d, udx, udy, udz,
> ist, iend, jst, jend,
> nx0, ny0, ipt, jpt )
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
double precision v( 5, -1:ldmx+2, -1:ldmy+2, *),
> tv(5, ldmx, ldmy),
> d( 5, 5, ldmx, ldmy),
> udx( 5, 5, ldmx, ldmy),
> udy( 5, 5, ldmx, ldmy),
> udz( 5, 5, ldmx, ldmy )
integer ist, iend
integer jst, jend
integer nx0, ny0
integer ipt, jpt
include 'timing.h'
c---------------------------------------------------------------------
c local variables
c---------------------------------------------------------------------
integer i, j, m
integer iex
double precision tmp, tmp1
double precision tmat(5,5)
c---------------------------------------------------------------------
c receive data from south and east
c---------------------------------------------------------------------
if (timeron) call timer_start(t_ucomm)
iex = 1
call exchange_1( v,k,iex )
if (timeron) call timer_stop(t_ucomm)
if (timeron) call timer_start(t_buts)
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
if (timeron) call timer_stop(t_buts)
c---------------------------------------------------------------------
c send data to north and west
c---------------------------------------------------------------------
if (timeron) call timer_start(t_ucomm)
iex = 3
call exchange_1( v,k,iex )
if (timeron) call timer_stop(t_ucomm)
return
end