
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine x_solve

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c this function performs the solution of the approximate factorization
c step in the x-direction for all five matrix components
c simultaneously. The Thomas algorithm is employed to solve the
c systems for the x-lines. Boundary conditions are non-periodic
c---------------------------------------------------------------------

       include 'header.h'
       include 'mpinpb.h'


       integer i, j, k, jp, kp, n, iend, jsize, ksize, i1, i2,
     >         buffer_size, c, m, p, istart, stage, error,
     >         requests(2), statuses(MPI_STATUS_SIZE, 2)
       double precision  r1, r2, d, e, s(5), sm1, sm2,
     >                   fac1, fac2


c---------------------------------------------------------------------
c      OK, now we know that there are multiple processors
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c now do a sweep on a layer-by-layer basis, i.e. sweeping through cells
c on this node in the direction of increasing i for the forward sweep,
c and after that reversing the direction for the backsubstitution.
c---------------------------------------------------------------------

       if (timeron) call timer_start(t_xsolve)
c---------------------------------------------------------------------
c                          FORWARD ELIMINATION  
c---------------------------------------------------------------------
       do    stage = 1, ncells
          c         = slice(1,stage)

          istart = 0
          iend   = cell_size(1,c)-1

          jsize     = cell_size(2,c)
          ksize     = cell_size(3,c)
          jp        = cell_coord(2,c)-1
          kp        = cell_coord(3,c)-1

          buffer_size = (jsize-start(2,c)-end(2,c)) * 
     >                  (ksize-start(3,c)-end(3,c))

          if ( stage .ne. 1) then

c---------------------------------------------------------------------
c            if this is not the first processor in this row of cells, 
c            receive data from predecessor containing the right hand
c            sides and the upper diagonal elements of the previous two rows
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_xcomm)
             call mpi_irecv(in_buffer, 22*buffer_size, 
     >                      dp_type, predecessor(1), 
     >                      DEFAULT_TAG,  comm_solve, 
     >                      requests(1), error)
             if (timeron) call timer_stop(t_xcomm)


c---------------------------------------------------------------------
c            communication has already been started. 
c            compute the left hand side while waiting for the msg
c---------------------------------------------------------------------
             call lhsx(c)

c---------------------------------------------------------------------
c            wait for pending communication to complete
c            This waits on the current receive and on the send
c            from the previous stage. They always come in pairs. 
c---------------------------------------------------------------------

             if (timeron) call timer_start(t_xcomm)
             call mpi_waitall(2, requests, statuses, error)
             if (timeron) call timer_stop(t_xcomm)

c---------------------------------------------------------------------
c            unpack the buffer                                 
c---------------------------------------------------------------------
             i  = istart
             i1 = istart + 1
             n = 0

c---------------------------------------------------------------------
c            create a running pointer
c---------------------------------------------------------------------
             p = 0
             do    k = start(3,c), ksize-end(3,c)-1
                do    j = start(2,c), jsize-end(2,c)-1
                   lhs(i,j,k,n+2,c) = lhs(i,j,k,n+2,c) -
     >                       in_buffer(p+1) * lhs(i,j,k,n+1,c)
                   lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) -
     >                       in_buffer(p+2) * lhs(i,j,k,n+1,c)
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                       in_buffer(p+2+m) * lhs(i,j,k,n+1,c)
                   end do
                   d            = in_buffer(p+6)
                   e            = in_buffer(p+7)
                   do    m = 1, 3
                      s(m) = in_buffer(p+7+m)
                   end do
                   r1 = lhs(i,j,k,n+2,c)
                   lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) - d * r1
                   lhs(i,j,k,n+4,c) = lhs(i,j,k,n+4,c) - e * r1
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - s(m) * r1
                   end do
                   r2 = lhs(i1,j,k,n+1,c)
                   lhs(i1,j,k,n+2,c) = lhs(i1,j,k,n+2,c) - d * r2
                   lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) - e * r2
                   do    m = 1, 3
                      rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) - s(m) * r2
                   end do
                   p = p + 10
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    j = start(2,c), jsize-end(2,c)-1
                      lhs(i,j,k,n+2,c) = lhs(i,j,k,n+2,c) -
     >                          in_buffer(p+1) * lhs(i,j,k,n+1,c)
                      lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) -
     >                          in_buffer(p+2) * lhs(i,j,k,n+1,c)
                      rhs(i,j,k,m,c)   = rhs(i,j,k,m,c) -
     >                          in_buffer(p+3) * lhs(i,j,k,n+1,c)
                      d                = in_buffer(p+4)
                      e                = in_buffer(p+5)
                      s(m)             = in_buffer(p+6)
                      r1 = lhs(i,j,k,n+2,c)
                      lhs(i,j,k,n+3,c) = lhs(i,j,k,n+3,c) - d * r1
                      lhs(i,j,k,n+4,c) = lhs(i,j,k,n+4,c) - e * r1
                      rhs(i,j,k,m,c)   = rhs(i,j,k,m,c) - s(m) * r1
                      r2 = lhs(i1,j,k,n+1,c)
                      lhs(i1,j,k,n+2,c) = lhs(i1,j,k,n+2,c) - d * r2
                      lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) - e * r2
                      rhs(i1,j,k,m,c)   = rhs(i1,j,k,m,c) - s(m) * r2
                      p = p + 6
                   end do
                end do
             end do

          else            

c---------------------------------------------------------------------
c            if this IS the first cell, we still compute the lhs
c---------------------------------------------------------------------
             call lhsx(c)
          endif

c---------------------------------------------------------------------
c         perform the Thomas algorithm; first, FORWARD ELIMINATION     
c---------------------------------------------------------------------
          n = 0

          do    k = start(3,c), ksize-end(3,c)-1
             do    j = start(2,c), jsize-end(2,c)-1
                do    i = istart, iend-2
                   i1 = i  + 1
                   i2 = i  + 2
                   fac1               = 1.d0/lhs(i,j,k,n+3,c)
                   lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                   lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                   end do
                   lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) -
     >                         lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+4,c)
                   lhs(i1,j,k,n+4,c) = lhs(i1,j,k,n+4,c) -
     >                         lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) -
     >                         lhs(i1,j,k,n+2,c)*rhs(i,j,k,m,c)
                   end do
                   lhs(i2,j,k,n+2,c) = lhs(i2,j,k,n+2,c) -
     >                         lhs(i2,j,k,n+1,c)*lhs(i,j,k,n+4,c)
                   lhs(i2,j,k,n+3,c) = lhs(i2,j,k,n+3,c) -
     >                         lhs(i2,j,k,n+1,c)*lhs(i,j,k,n+5,c)
                   do    m = 1, 3
                      rhs(i2,j,k,m,c) = rhs(i2,j,k,m,c) -
     >                         lhs(i2,j,k,n+1,c)*rhs(i,j,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         The last two rows in this grid block are a bit different, 
c         since they do not have two more rows available for the
c         elimination of off-diagonal entries
c---------------------------------------------------------------------

          i  = iend - 1
          i1 = iend
          do    k = start(3,c), ksize-end(3,c)-1
             do    j = start(2,c), jsize-end(2,c)-1
                fac1               = 1.d0/lhs(i,j,k,n+3,c)
                lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                do    m = 1, 3
                   rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                end do
                lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) -
     >                      lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+4,c)
                lhs(i1,j,k,n+4,c) = lhs(i1,j,k,n+4,c) -
     >                      lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+5,c)
                do    m = 1, 3
                   rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) -
     >                      lhs(i1,j,k,n+2,c)*rhs(i,j,k,m,c)
                end do
c---------------------------------------------------------------------
c               scale the last row immediately (some of this is
c               overkill in case this is the last cell)
c---------------------------------------------------------------------
                fac2               = 1.d0/lhs(i1,j,k,n+3,c)
                lhs(i1,j,k,n+4,c) = fac2*lhs(i1,j,k,n+4,c)
                lhs(i1,j,k,n+5,c) = fac2*lhs(i1,j,k,n+5,c)  
                do    m = 1, 3
                   rhs(i1,j,k,m,c) = fac2*rhs(i1,j,k,m,c)
                end do
             end do
          end do

c---------------------------------------------------------------------
c         do the u+c and the u-c factors                 
c---------------------------------------------------------------------

          do    m = 4, 5
             n = (m-3)*5
             do    k = start(3,c), ksize-end(3,c)-1
                do    j = start(2,c), jsize-end(2,c)-1
                   do    i = istart, iend-2
                   i1 = i  + 1
                   i2 = i  + 2
                   fac1               = 1.d0/lhs(i,j,k,n+3,c)
                   lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                   lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                   rhs(i,j,k,m,c) = fac1*rhs(i,j,k,m,c)
                   lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) -
     >                         lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+4,c)
                   lhs(i1,j,k,n+4,c) = lhs(i1,j,k,n+4,c) -
     >                         lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+5,c)
                   rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) -
     >                         lhs(i1,j,k,n+2,c)*rhs(i,j,k,m,c)
                   lhs(i2,j,k,n+2,c) = lhs(i2,j,k,n+2,c) -
     >                         lhs(i2,j,k,n+1,c)*lhs(i,j,k,n+4,c)
                   lhs(i2,j,k,n+3,c) = lhs(i2,j,k,n+3,c) -
     >                         lhs(i2,j,k,n+1,c)*lhs(i,j,k,n+5,c)
                   rhs(i2,j,k,m,c) = rhs(i2,j,k,m,c) -
     >                         lhs(i2,j,k,n+1,c)*rhs(i,j,k,m,c)
                end do
             end do
          end do

c---------------------------------------------------------------------
c            And again the last two rows separately
c---------------------------------------------------------------------
             i  = iend - 1
             i1 = iend
             do    k = start(3,c), ksize-end(3,c)-1
                do    j = start(2,c), jsize-end(2,c)-1
                fac1               = 1.d0/lhs(i,j,k,n+3,c)
                lhs(i,j,k,n+4,c)   = fac1*lhs(i,j,k,n+4,c)
                lhs(i,j,k,n+5,c)   = fac1*lhs(i,j,k,n+5,c)
                rhs(i,j,k,m,c)     = fac1*rhs(i,j,k,m,c)
                lhs(i1,j,k,n+3,c) = lhs(i1,j,k,n+3,c) -
     >                      lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+4,c)
                lhs(i1,j,k,n+4,c) = lhs(i1,j,k,n+4,c) -
     >                      lhs(i1,j,k,n+2,c)*lhs(i,j,k,n+5,c)
                rhs(i1,j,k,m,c)   = rhs(i1,j,k,m,c) -
     >                      lhs(i1,j,k,n+2,c)*rhs(i,j,k,m,c)
c---------------------------------------------------------------------
c               Scale the last row immediately (some of this is overkill
c               if this is the last cell)
c---------------------------------------------------------------------
                fac2               = 1.d0/lhs(i1,j,k,n+3,c)
                lhs(i1,j,k,n+4,c) = fac2*lhs(i1,j,k,n+4,c)
                lhs(i1,j,k,n+5,c) = fac2*lhs(i1,j,k,n+5,c)
                rhs(i1,j,k,m,c)   = fac2*rhs(i1,j,k,m,c)

             end do
          end do
       end do

c---------------------------------------------------------------------
c         send information to the next processor, except when this
c         is the last grid block
c---------------------------------------------------------------------
          if (stage .ne. ncells) then

c---------------------------------------------------------------------
c            create a running pointer for the send buffer  
c---------------------------------------------------------------------
             p = 0
             n = 0
             do    k = start(3,c), ksize-end(3,c)-1
                do    j = start(2,c), jsize-end(2,c)-1
                   do    i = iend-1, iend
                      out_buffer(p+1) = lhs(i,j,k,n+4,c)
                      out_buffer(p+2) = lhs(i,j,k,n+5,c)
                      do    m = 1, 3
                         out_buffer(p+2+m) = rhs(i,j,k,m,c)
                      end do
                      p = p+5
                   end do
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    j = start(2,c), jsize-end(2,c)-1
                      do    i = iend-1, iend
                         out_buffer(p+1) = lhs(i,j,k,n+4,c)
                         out_buffer(p+2) = lhs(i,j,k,n+5,c)
                         out_buffer(p+3) = rhs(i,j,k,m,c)
                         p = p + 3
                      end do
                   end do
                end do
             end do

c---------------------------------------------------------------------
c send data to next phase
c can't receive data yet because buffer size will be wrong 
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_xcomm)
             call mpi_isend(out_buffer, 22*buffer_size, 
     >                     dp_type, successor(1), 
     >                     DEFAULT_TAG, comm_solve, 
     >                     requests(2), error)
             if (timeron) call timer_stop(t_xcomm)

          endif
       end do

c---------------------------------------------------------------------
c      now go in the reverse direction                      
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c                         BACKSUBSTITUTION 
c---------------------------------------------------------------------
       do    stage = ncells, 1, -1
          c = slice(1,stage)

          istart = 0
          iend   = cell_size(1,c)-1

          jsize = cell_size(2,c)
          ksize = cell_size(3,c)
          jp    = cell_coord(2,c)-1
          kp    = cell_coord(3,c)-1

          buffer_size = (jsize-start(2,c)-end(2,c)) * 
     >                  (ksize-start(3,c)-end(3,c))

          if (stage .ne. ncells) then

c---------------------------------------------------------------------
c            if this is not the starting cell in this row of cells, 
c            wait for a message to be received, containing the 
c            solution of the previous two stations     
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_xcomm)
             call mpi_irecv(in_buffer, 10*buffer_size, 
     >                      dp_type, successor(1), 
     >                      DEFAULT_TAG, comm_solve, 
     >                      requests(1), error)
             if (timeron) call timer_stop(t_xcomm)


c---------------------------------------------------------------------
c            communication has already been started
c            while waiting, do the block-diagonal inversion for the 
c            cell that was just finished                
c---------------------------------------------------------------------

             call ninvr(slice(1,stage+1))

c---------------------------------------------------------------------
c            wait for pending communication to complete
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_xcomm)
             call mpi_waitall(2, requests, statuses, error)
             if (timeron) call timer_stop(t_xcomm)

c---------------------------------------------------------------------
c            unpack the buffer for the first three factors         
c---------------------------------------------------------------------
             n = 0
             p = 0
             i  = iend
             i1 = i - 1
             do    m = 1, 3
                do   k = start(3,c), ksize-end(3,c)-1
                   do   j = start(2,c), jsize-end(2,c)-1
                      sm1 = in_buffer(p+1)
                      sm2 = in_buffer(p+2)
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                        lhs(i,j,k,n+4,c)*sm1 -
     >                        lhs(i,j,k,n+5,c)*sm2
                      rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) -
     >                        lhs(i1,j,k,n+4,c) * rhs(i,j,k,m,c) - 
     >                        lhs(i1,j,k,n+5,c) * sm1
                      p = p + 2
                   end do
                end do
             end do

c---------------------------------------------------------------------
c            now unpack the buffer for the remaining two factors
c---------------------------------------------------------------------
             do    m = 4, 5
                n = (m-3)*5
                do   k = start(3,c), ksize-end(3,c)-1
                   do   j = start(2,c), jsize-end(2,c)-1
                      sm1 = in_buffer(p+1)
                      sm2 = in_buffer(p+2)
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                        lhs(i,j,k,n+4,c)*sm1 -
     >                        lhs(i,j,k,n+5,c)*sm2
                      rhs(i1,j,k,m,c) = rhs(i1,j,k,m,c) -
     >                        lhs(i1,j,k,n+4,c) * rhs(i,j,k,m,c) - 
     >                        lhs(i1,j,k,n+5,c) * sm1
                      p = p + 2
                   end do
                end do
             end do

          else

c---------------------------------------------------------------------
c            now we know this is the first grid block on the back sweep,
c            so we don't need a message to start the substitution. 
c---------------------------------------------------------------------
             i  = iend-1
             i1 = iend
             n = 0
             do   m = 1, 3
                do   k = start(3,c), ksize-end(3,c)-1
                   do   j = start(2,c), jsize-end(2,c)-1
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                             lhs(i,j,k,n+4,c)*rhs(i1,j,k,m,c)
                   end do
                end do
             end do

             do    m = 4, 5
                n = (m-3)*5
                do   k = start(3,c), ksize-end(3,c)-1
                   do   j = start(2,c), jsize-end(2,c)-1
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) -
     >                             lhs(i,j,k,n+4,c)*rhs(i1,j,k,m,c)
                   end do
                end do
             end do
          endif

c---------------------------------------------------------------------
c         Whether or not this is the last processor, we always have
c         to complete the back-substitution 
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c         The first three factors
c---------------------------------------------------------------------
          n = 0
          do   m = 1, 3
             do   k = start(3,c), ksize-end(3,c)-1
                do   j = start(2,c), jsize-end(2,c)-1
                   do    i = iend-2, istart, -1
                      i1 = i  + 1
                      i2 = i  + 2
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                          lhs(i,j,k,n+4,c)*rhs(i1,j,k,m,c) -
     >                          lhs(i,j,k,n+5,c)*rhs(i2,j,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         And the remaining two
c---------------------------------------------------------------------
          do    m = 4, 5
             n = (m-3)*5
             do   k = start(3,c), ksize-end(3,c)-1
                do   j = start(2,c), jsize-end(2,c)-1
                   do    i = iend-2, istart, -1
                      i1 = i  + 1
                      i2 = i  + 2
                      rhs(i,j,k,m,c) = rhs(i,j,k,m,c) - 
     >                          lhs(i,j,k,n+4,c)*rhs(i1,j,k,m,c) -
     >                          lhs(i,j,k,n+5,c)*rhs(i2,j,k,m,c)
                   end do
                end do
             end do
          end do

c---------------------------------------------------------------------
c         send on information to the previous processor, if needed
c---------------------------------------------------------------------
          if (stage .ne.  1) then
             i  = istart
             i1 = istart+1
             p = 0
             do    m = 1, 5
                do    k = start(3,c), ksize-end(3,c)-1
                   do    j = start(2,c), jsize-end(2,c)-1
                      out_buffer(p+1) = rhs(i,j,k,m,c)
                       out_buffer(p+2) = rhs(i1,j,k,m,c)
                      p = p + 2
                   end do
                end do
             end do

c---------------------------------------------------------------------
c            pack and send the buffer
c---------------------------------------------------------------------
             if (timeron) call timer_start(t_xcomm)
             call mpi_isend(out_buffer, 10*buffer_size, 
     >                     dp_type, predecessor(1), 
     >                     DEFAULT_TAG, comm_solve, 
     >                     requests(2), error)
             if (timeron) call timer_stop(t_xcomm)

          endif

c---------------------------------------------------------------------
c         If this was the last stage, do the block-diagonal inversion          
c---------------------------------------------------------------------
          if (stage .eq. 1) call ninvr(c)

       end do

       if (timeron) call timer_stop(t_xsolve)

       return
       end