--- /dev/null
+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
+
+c---------------------------------------------------------------------
+c local variables
+c---------------------------------------------------------------------
+ integer i, j, m, l, istp, iendp
+ integer iex
+ double precision tmp, tmp1
+ double precision tmat(5,5)
+
+
+c---------------------------------------------------------------------
+c receive data from south and east
+c---------------------------------------------------------------------
+ iex = 1
+ call exchange_1( v,k,iex )
+
+ 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 l = iend+jend, ist+jst, -1
+ istp = max(l - jend, ist)
+ iendp = min(l - jst, iend)
+
+!dir$ ivdep
+ do i = istp, iendp
+ j = l - i
+
+!!dir$ unroll 5
+! manually unroll the loop
+! do m = 1, 5
+ tv( 1, i, j ) = tv( 1, i, j )
+ > + omega * ( udy( 1, 1, i, j ) * v( 1, i, j+1, k )
+ > + udx( 1, 1, i, j ) * v( 1, i+1, j, k )
+ > + udy( 1, 2, i, j ) * v( 2, i, j+1, k )
+ > + udx( 1, 2, i, j ) * v( 2, i+1, j, k )
+ > + udy( 1, 3, i, j ) * v( 3, i, j+1, k )
+ > + udx( 1, 3, i, j ) * v( 3, i+1, j, k )
+ > + udy( 1, 4, i, j ) * v( 4, i, j+1, k )
+ > + udx( 1, 4, i, j ) * v( 4, i+1, j, k )
+ > + udy( 1, 5, i, j ) * v( 5, i, j+1, k )
+ > + udx( 1, 5, i, j ) * v( 5, i+1, j, k ) )
+ tv( 2, i, j ) = tv( 2, i, j )
+ > + omega * ( udy( 2, 1, i, j ) * v( 1, i, j+1, k )
+ > + udx( 2, 1, i, j ) * v( 1, i+1, j, k )
+ > + udy( 2, 2, i, j ) * v( 2, i, j+1, k )
+ > + udx( 2, 2, i, j ) * v( 2, i+1, j, k )
+ > + udy( 2, 3, i, j ) * v( 3, i, j+1, k )
+ > + udx( 2, 3, i, j ) * v( 3, i+1, j, k )
+ > + udy( 2, 4, i, j ) * v( 4, i, j+1, k )
+ > + udx( 2, 4, i, j ) * v( 4, i+1, j, k )
+ > + udy( 2, 5, i, j ) * v( 5, i, j+1, k )
+ > + udx( 2, 5, i, j ) * v( 5, i+1, j, k ) )
+ tv( 3, i, j ) = tv( 3, i, j )
+ > + omega * ( udy( 3, 1, i, j ) * v( 1, i, j+1, k )
+ > + udx( 3, 1, i, j ) * v( 1, i+1, j, k )
+ > + udy( 3, 2, i, j ) * v( 2, i, j+1, k )
+ > + udx( 3, 2, i, j ) * v( 2, i+1, j, k )
+ > + udy( 3, 3, i, j ) * v( 3, i, j+1, k )
+ > + udx( 3, 3, i, j ) * v( 3, i+1, j, k )
+ > + udy( 3, 4, i, j ) * v( 4, i, j+1, k )
+ > + udx( 3, 4, i, j ) * v( 4, i+1, j, k )
+ > + udy( 3, 5, i, j ) * v( 5, i, j+1, k )
+ > + udx( 3, 5, i, j ) * v( 5, i+1, j, k ) )
+ tv( 4, i, j ) = tv( 4, i, j )
+ > + omega * ( udy( 4, 1, i, j ) * v( 1, i, j+1, k )
+ > + udx( 4, 1, i, j ) * v( 1, i+1, j, k )
+ > + udy( 4, 2, i, j ) * v( 2, i, j+1, k )
+ > + udx( 4, 2, i, j ) * v( 2, i+1, j, k )
+ > + udy( 4, 3, i, j ) * v( 3, i, j+1, k )
+ > + udx( 4, 3, i, j ) * v( 3, i+1, j, k )
+ > + udy( 4, 4, i, j ) * v( 4, i, j+1, k )
+ > + udx( 4, 4, i, j ) * v( 4, i+1, j, k )
+ > + udy( 4, 5, i, j ) * v( 5, i, j+1, k )
+ > + udx( 4, 5, i, j ) * v( 5, i+1, j, k ) )
+ tv( 5, i, j ) = tv( 5, i, j )
+ > + omega * ( udy( 5, 1, i, j ) * v( 1, i, j+1, k )
+ > + udx( 5, 1, i, j ) * v( 1, i+1, j, k )
+ > + udy( 5, 2, i, j ) * v( 2, i, j+1, k )
+ > + udx( 5, 2, i, j ) * v( 2, i+1, j, k )
+ > + udy( 5, 3, i, j ) * v( 3, i, j+1, k )
+ > + udx( 5, 3, i, j ) * v( 3, i+1, j, k )
+ > + udy( 5, 4, i, j ) * v( 4, i, j+1, k )
+ > + udx( 5, 4, i, j ) * v( 4, i+1, j, k )
+ > + udy( 5, 5, i, j ) * v( 5, i, j+1, k )
+ > + udx( 5, 5, i, j ) * v( 5, i+1, j, k ) )
+! end do
+
+c---------------------------------------------------------------------
+c diagonal block inversion
+c---------------------------------------------------------------------
+!!dir$ unroll 5
+! manually unroll the loop
+! do m = 1, 5
+ tmat( 1, 1 ) = d( 1, 1, i, j )
+ tmat( 1, 2 ) = d( 1, 2, i, j )
+ tmat( 1, 3 ) = d( 1, 3, i, j )
+ tmat( 1, 4 ) = d( 1, 4, i, j )
+ tmat( 1, 5 ) = d( 1, 5, i, j )
+ tmat( 2, 1 ) = d( 2, 1, i, j )
+ tmat( 2, 2 ) = d( 2, 2, i, j )
+ tmat( 2, 3 ) = d( 2, 3, i, j )
+ tmat( 2, 4 ) = d( 2, 4, i, j )
+ tmat( 2, 5 ) = d( 2, 5, i, j )
+ tmat( 3, 1 ) = d( 3, 1, i, j )
+ tmat( 3, 2 ) = d( 3, 2, i, j )
+ tmat( 3, 3 ) = d( 3, 3, i, j )
+ tmat( 3, 4 ) = d( 3, 4, i, j )
+ tmat( 3, 5 ) = d( 3, 5, i, j )
+ tmat( 4, 1 ) = d( 4, 1, i, j )
+ tmat( 4, 2 ) = d( 4, 2, i, j )
+ tmat( 4, 3 ) = d( 4, 3, i, j )
+ tmat( 4, 4 ) = d( 4, 4, i, j )
+ tmat( 4, 5 ) = d( 4, 5, i, j )
+ tmat( 5, 1 ) = d( 5, 1, i, j )
+ tmat( 5, 2 ) = d( 5, 2, i, j )
+ tmat( 5, 3 ) = d( 5, 3, i, j )
+ tmat( 5, 4 ) = d( 5, 4, i, j )
+ tmat( 5, 5 ) = d( 5, 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
+
+c---------------------------------------------------------------------
+c send data to north and west
+c---------------------------------------------------------------------
+ iex = 3
+ call exchange_1( v,k,iex )
+
+ return
+ end
