c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_solve c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c c Performs line solves in X direction by first factoring c the block-tridiagonal matrix into an upper triangular matrix, c and then performing back substitution to solve for the unknow c vectors of each line. c c Make sure we treat elements zero to cell_size in the direction c of the sweep. c c--------------------------------------------------------------------- include 'header.h' include 'mpinpb.h' integer c, istart, stage, > first, last, recv_id, error, r_status(MPI_STATUS_SIZE), > isize,jsize,ksize,send_id istart = 0 c--------------------------------------------------------------------- c in our terminology stage is the number of the cell in the x-direction c i.e. stage = 1 means the start of the line stage=ncells means end c--------------------------------------------------------------------- do stage = 1,ncells c = slice(1,stage) isize = cell_size(1,c) - 1 jsize = cell_size(2,c) - 1 ksize = cell_size(3,c) - 1 c--------------------------------------------------------------------- c set last-cell flag c--------------------------------------------------------------------- if (stage .eq. ncells) then last = 1 else last = 0 endif if (stage .eq. 1) then c--------------------------------------------------------------------- c This is the first cell, so solve without receiving data c--------------------------------------------------------------------- first = 1 c call lhsx(c) call x_solve_cell(first,last,c) else c--------------------------------------------------------------------- c Not the first cell of this line, so receive info from c processor working on preceeding cell c--------------------------------------------------------------------- first = 0 call x_receive_solve_info(recv_id,c) c--------------------------------------------------------------------- c overlap computations and communications c--------------------------------------------------------------------- c call lhsx(c) c--------------------------------------------------------------------- c wait for completion c--------------------------------------------------------------------- call mpi_wait(send_id,r_status,error) call mpi_wait(recv_id,r_status,error) c--------------------------------------------------------------------- c install C'(istart) and rhs'(istart) to be used in this cell c--------------------------------------------------------------------- call x_unpack_solve_info(c) call x_solve_cell(first,last,c) endif if (last .eq. 0) call x_send_solve_info(send_id,c) enddo c--------------------------------------------------------------------- c now perform backsubstitution in reverse direction c--------------------------------------------------------------------- do stage = ncells, 1, -1 c = slice(1,stage) first = 0 last = 0 if (stage .eq. 1) first = 1 if (stage .eq. ncells) then last = 1 c--------------------------------------------------------------------- c last cell, so perform back substitute without waiting c--------------------------------------------------------------------- call x_backsubstitute(first, last,c) else call x_receive_backsub_info(recv_id,c) call mpi_wait(send_id,r_status,error) call mpi_wait(recv_id,r_status,error) call x_unpack_backsub_info(c) call x_backsubstitute(first,last,c) endif if (first .eq. 0) call x_send_backsub_info(send_id,c) enddo return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_unpack_solve_info(c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c unpack C'(-1) and rhs'(-1) for c all j and k c--------------------------------------------------------------------- include 'header.h' integer j,k,m,n,ptr,c,istart istart = 0 ptr = 0 do k=0,KMAX-1 do j=0,JMAX-1 do m=1,BLOCK_SIZE do n=1,BLOCK_SIZE lhsc(m,n,istart-1,j,k,c) = out_buffer(ptr+n) enddo ptr = ptr+BLOCK_SIZE enddo do n=1,BLOCK_SIZE rhs(n,istart-1,j,k,c) = out_buffer(ptr+n) enddo ptr = ptr+BLOCK_SIZE enddo enddo return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_send_solve_info(send_id,c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c pack up and send C'(iend) and rhs'(iend) for c all j and k c--------------------------------------------------------------------- include 'header.h' include 'mpinpb.h' integer j,k,m,n,isize,ptr,c,jp,kp integer error,send_id,buffer_size isize = cell_size(1,c)-1 jp = cell_coord(2,c) - 1 kp = cell_coord(3,c) - 1 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM* > (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE) c--------------------------------------------------------------------- c pack up buffer c--------------------------------------------------------------------- ptr = 0 do k=0,KMAX-1 do j=0,JMAX-1 do m=1,BLOCK_SIZE do n=1,BLOCK_SIZE in_buffer(ptr+n) = lhsc(m,n,isize,j,k,c) enddo ptr = ptr+BLOCK_SIZE enddo do n=1,BLOCK_SIZE in_buffer(ptr+n) = rhs(n,isize,j,k,c) enddo ptr = ptr+BLOCK_SIZE enddo enddo c--------------------------------------------------------------------- c send buffer c--------------------------------------------------------------------- call mpi_isend(in_buffer, buffer_size, > dp_type, successor(1), > WEST+jp+kp*NCELLS, comm_solve, > send_id,error) return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_send_backsub_info(send_id,c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c pack up and send U(istart) for all j and k c--------------------------------------------------------------------- include 'header.h' include 'mpinpb.h' integer j,k,n,ptr,c,istart,jp,kp integer error,send_id,buffer_size c--------------------------------------------------------------------- c Send element 0 to previous processor c--------------------------------------------------------------------- istart = 0 jp = cell_coord(2,c)-1 kp = cell_coord(3,c)-1 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE ptr = 0 do k=0,KMAX-1 do j=0,JMAX-1 do n=1,BLOCK_SIZE in_buffer(ptr+n) = rhs(n,istart,j,k,c) enddo ptr = ptr+BLOCK_SIZE enddo enddo call mpi_isend(in_buffer, buffer_size, > dp_type, predecessor(1), > EAST+jp+kp*NCELLS, comm_solve, > send_id,error) return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_unpack_backsub_info(c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c unpack U(isize) for all j and k c--------------------------------------------------------------------- include 'header.h' integer j,k,n,ptr,c ptr = 0 do k=0,KMAX-1 do j=0,JMAX-1 do n=1,BLOCK_SIZE backsub_info(n,j,k,c) = out_buffer(ptr+n) enddo ptr = ptr+BLOCK_SIZE enddo enddo return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_receive_backsub_info(recv_id,c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c post mpi receives c--------------------------------------------------------------------- include 'header.h' include 'mpinpb.h' integer error,recv_id,jp,kp,c,buffer_size jp = cell_coord(2,c) - 1 kp = cell_coord(3,c) - 1 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE call mpi_irecv(out_buffer, buffer_size, > dp_type, successor(1), > EAST+jp+kp*NCELLS, comm_solve, > recv_id, error) return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_receive_solve_info(recv_id,c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c post mpi receives c--------------------------------------------------------------------- include 'header.h' include 'mpinpb.h' integer jp,kp,recv_id,error,c,buffer_size jp = cell_coord(2,c) - 1 kp = cell_coord(3,c) - 1 buffer_size=MAX_CELL_DIM*MAX_CELL_DIM* > (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE) call mpi_irecv(out_buffer, buffer_size, > dp_type, predecessor(1), > WEST+jp+kp*NCELLS, comm_solve, > recv_id, error) return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_backsubstitute(first, last, c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c back solve: if last cell, then generate U(isize)=rhs(isize) c else assume U(isize) is loaded in un pack backsub_info c so just use it c after call u(istart) will be sent to next cell c--------------------------------------------------------------------- include 'header.h' integer first, last, c, i, j, k integer m,n,isize,jsize,ksize,istart istart = 0 isize = cell_size(1,c)-1 jsize = cell_size(2,c)-end(2,c)-1 ksize = cell_size(3,c)-end(3,c)-1 if (last .eq. 0) then do k=start(3,c),ksize do j=start(2,c),jsize c--------------------------------------------------------------------- c U(isize) uses info from previous cell if not last cell c--------------------------------------------------------------------- do m=1,BLOCK_SIZE do n=1,BLOCK_SIZE rhs(m,isize,j,k,c) = rhs(m,isize,j,k,c) > - lhsc(m,n,isize,j,k,c)* > backsub_info(n,j,k,c) c--------------------------------------------------------------------- c rhs(m,isize,j,k,c) = rhs(m,isize,j,k,c) c $ - lhsc(m,n,isize,j,k,c)*rhs(n,isize+1,j,k,c) c--------------------------------------------------------------------- enddo enddo enddo enddo endif do k=start(3,c),ksize do j=start(2,c),jsize do i=isize-1,istart,-1 do m=1,BLOCK_SIZE do n=1,BLOCK_SIZE rhs(m,i,j,k,c) = rhs(m,i,j,k,c) > - lhsc(m,n,i,j,k,c)*rhs(n,i+1,j,k,c) enddo enddo enddo enddo enddo return end c--------------------------------------------------------------------- c--------------------------------------------------------------------- subroutine x_solve_cell(first,last,c) c--------------------------------------------------------------------- c--------------------------------------------------------------------- c--------------------------------------------------------------------- c performs guaussian elimination on this cell. c c assumes that unpacking routines for non-first cells c preload C' and rhs' from previous cell. c c assumed send happens outside this routine, but that c c'(IMAX) and rhs'(IMAX) will be sent to next cell c--------------------------------------------------------------------- include 'header.h' include 'work_lhs.h' integer first,last,c integer i,j,k,isize,ksize,jsize,istart istart = 0 isize = cell_size(1,c)-1 jsize = cell_size(2,c)-end(2,c)-1 ksize = cell_size(3,c)-end(3,c)-1 call lhsabinit(lhsa, lhsb, isize) do k=start(3,c),ksize do j=start(2,c),jsize c--------------------------------------------------------------------- c This function computes the left hand side in the xi-direction c--------------------------------------------------------------------- c--------------------------------------------------------------------- c determine a (labeled f) and n jacobians for cell c c--------------------------------------------------------------------- do i = start(1,c)-1, cell_size(1,c) - end(1,c) tmp1 = rho_i(i,j,k,c) tmp2 = tmp1 * tmp1 tmp3 = tmp1 * tmp2 c--------------------------------------------------------------------- c c--------------------------------------------------------------------- fjac(1,1,i) = 0.0d+00 fjac(1,2,i) = 1.0d+00 fjac(1,3,i) = 0.0d+00 fjac(1,4,i) = 0.0d+00 fjac(1,5,i) = 0.0d+00 fjac(2,1,i) = -(u(2,i,j,k,c) * tmp2 * > u(2,i,j,k,c)) > + c2 * qs(i,j,k,c) fjac(2,2,i) = ( 2.0d+00 - c2 ) > * ( u(2,i,j,k,c) * tmp1 ) fjac(2,3,i) = - c2 * ( u(3,i,j,k,c) * tmp1 ) fjac(2,4,i) = - c2 * ( u(4,i,j,k,c) * tmp1 ) fjac(2,5,i) = c2 fjac(3,1,i) = - ( u(2,i,j,k,c)*u(3,i,j,k,c) ) * tmp2 fjac(3,2,i) = u(3,i,j,k,c) * tmp1 fjac(3,3,i) = u(2,i,j,k,c) * tmp1 fjac(3,4,i) = 0.0d+00 fjac(3,5,i) = 0.0d+00 fjac(4,1,i) = - ( u(2,i,j,k,c)*u(4,i,j,k,c) ) * tmp2 fjac(4,2,i) = u(4,i,j,k,c) * tmp1 fjac(4,3,i) = 0.0d+00 fjac(4,4,i) = u(2,i,j,k,c) * tmp1 fjac(4,5,i) = 0.0d+00 fjac(5,1,i) = ( c2 * 2.0d0 * qs(i,j,k,c) > - c1 * ( u(5,i,j,k,c) * tmp1 ) ) > * ( u(2,i,j,k,c) * tmp1 ) fjac(5,2,i) = c1 * u(5,i,j,k,c) * tmp1 > - c2 > * ( u(2,i,j,k,c)*u(2,i,j,k,c) * tmp2 > + qs(i,j,k,c) ) fjac(5,3,i) = - c2 * ( u(3,i,j,k,c)*u(2,i,j,k,c) ) > * tmp2 fjac(5,4,i) = - c2 * ( u(4,i,j,k,c)*u(2,i,j,k,c) ) > * tmp2 fjac(5,5,i) = c1 * ( u(2,i,j,k,c) * tmp1 ) njac(1,1,i) = 0.0d+00 njac(1,2,i) = 0.0d+00 njac(1,3,i) = 0.0d+00 njac(1,4,i) = 0.0d+00 njac(1,5,i) = 0.0d+00 njac(2,1,i) = - con43 * c3c4 * tmp2 * u(2,i,j,k,c) njac(2,2,i) = con43 * c3c4 * tmp1 njac(2,3,i) = 0.0d+00 njac(2,4,i) = 0.0d+00 njac(2,5,i) = 0.0d+00 njac(3,1,i) = - c3c4 * tmp2 * u(3,i,j,k,c) njac(3,2,i) = 0.0d+00 njac(3,3,i) = c3c4 * tmp1 njac(3,4,i) = 0.0d+00 njac(3,5,i) = 0.0d+00 njac(4,1,i) = - c3c4 * tmp2 * u(4,i,j,k,c) njac(4,2,i) = 0.0d+00 njac(4,3,i) = 0.0d+00 njac(4,4,i) = c3c4 * tmp1 njac(4,5,i) = 0.0d+00 njac(5,1,i) = - ( con43 * c3c4 > - c1345 ) * tmp3 * (u(2,i,j,k,c)**2) > - ( c3c4 - c1345 ) * tmp3 * (u(3,i,j,k,c)**2) > - ( c3c4 - c1345 ) * tmp3 * (u(4,i,j,k,c)**2) > - c1345 * tmp2 * u(5,i,j,k,c) njac(5,2,i) = ( con43 * c3c4 > - c1345 ) * tmp2 * u(2,i,j,k,c) njac(5,3,i) = ( c3c4 - c1345 ) * tmp2 * u(3,i,j,k,c) njac(5,4,i) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k,c) njac(5,5,i) = ( c1345 ) * tmp1 enddo c--------------------------------------------------------------------- c now jacobians set, so form left hand side in x direction c--------------------------------------------------------------------- do i = start(1,c), isize - end(1,c) tmp1 = dt * tx1 tmp2 = dt * tx2 lhsa(1,1,i) = - tmp2 * fjac(1,1,i-1) > - tmp1 * njac(1,1,i-1) > - tmp1 * dx1 lhsa(1,2,i) = - tmp2 * fjac(1,2,i-1) > - tmp1 * njac(1,2,i-1) lhsa(1,3,i) = - tmp2 * fjac(1,3,i-1) > - tmp1 * njac(1,3,i-1) lhsa(1,4,i) = - tmp2 * fjac(1,4,i-1) > - tmp1 * njac(1,4,i-1) lhsa(1,5,i) = - tmp2 * fjac(1,5,i-1) > - tmp1 * njac(1,5,i-1) lhsa(2,1,i) = - tmp2 * fjac(2,1,i-1) > - tmp1 * njac(2,1,i-1) lhsa(2,2,i) = - tmp2 * fjac(2,2,i-1) > - tmp1 * njac(2,2,i-1) > - tmp1 * dx2 lhsa(2,3,i) = - tmp2 * fjac(2,3,i-1) > - tmp1 * njac(2,3,i-1) lhsa(2,4,i) = - tmp2 * fjac(2,4,i-1) > - tmp1 * njac(2,4,i-1) lhsa(2,5,i) = - tmp2 * fjac(2,5,i-1) > - tmp1 * njac(2,5,i-1) lhsa(3,1,i) = - tmp2 * fjac(3,1,i-1) > - tmp1 * njac(3,1,i-1) lhsa(3,2,i) = - tmp2 * fjac(3,2,i-1) > - tmp1 * njac(3,2,i-1) lhsa(3,3,i) = - tmp2 * fjac(3,3,i-1) > - tmp1 * njac(3,3,i-1) > - tmp1 * dx3 lhsa(3,4,i) = - tmp2 * fjac(3,4,i-1) > - tmp1 * njac(3,4,i-1) lhsa(3,5,i) = - tmp2 * fjac(3,5,i-1) > - tmp1 * njac(3,5,i-1) lhsa(4,1,i) = - tmp2 * fjac(4,1,i-1) > - tmp1 * njac(4,1,i-1) lhsa(4,2,i) = - tmp2 * fjac(4,2,i-1) > - tmp1 * njac(4,2,i-1) lhsa(4,3,i) = - tmp2 * fjac(4,3,i-1) > - tmp1 * njac(4,3,i-1) lhsa(4,4,i) = - tmp2 * fjac(4,4,i-1) > - tmp1 * njac(4,4,i-1) > - tmp1 * dx4 lhsa(4,5,i) = - tmp2 * fjac(4,5,i-1) > - tmp1 * njac(4,5,i-1) lhsa(5,1,i) = - tmp2 * fjac(5,1,i-1) > - tmp1 * njac(5,1,i-1) lhsa(5,2,i) = - tmp2 * fjac(5,2,i-1) > - tmp1 * njac(5,2,i-1) lhsa(5,3,i) = - tmp2 * fjac(5,3,i-1) > - tmp1 * njac(5,3,i-1) lhsa(5,4,i) = - tmp2 * fjac(5,4,i-1) > - tmp1 * njac(5,4,i-1) lhsa(5,5,i) = - tmp2 * fjac(5,5,i-1) > - tmp1 * njac(5,5,i-1) > - tmp1 * dx5 lhsb(1,1,i) = 1.0d+00 > + tmp1 * 2.0d+00 * njac(1,1,i) > + tmp1 * 2.0d+00 * dx1 lhsb(1,2,i) = tmp1 * 2.0d+00 * njac(1,2,i) lhsb(1,3,i) = tmp1 * 2.0d+00 * njac(1,3,i) lhsb(1,4,i) = tmp1 * 2.0d+00 * njac(1,4,i) lhsb(1,5,i) = tmp1 * 2.0d+00 * njac(1,5,i) lhsb(2,1,i) = tmp1 * 2.0d+00 * njac(2,1,i) lhsb(2,2,i) = 1.0d+00 > + tmp1 * 2.0d+00 * njac(2,2,i) > + tmp1 * 2.0d+00 * dx2 lhsb(2,3,i) = tmp1 * 2.0d+00 * njac(2,3,i) lhsb(2,4,i) = tmp1 * 2.0d+00 * njac(2,4,i) lhsb(2,5,i) = tmp1 * 2.0d+00 * njac(2,5,i) lhsb(3,1,i) = tmp1 * 2.0d+00 * njac(3,1,i) lhsb(3,2,i) = tmp1 * 2.0d+00 * njac(3,2,i) lhsb(3,3,i) = 1.0d+00 > + tmp1 * 2.0d+00 * njac(3,3,i) > + tmp1 * 2.0d+00 * dx3 lhsb(3,4,i) = tmp1 * 2.0d+00 * njac(3,4,i) lhsb(3,5,i) = tmp1 * 2.0d+00 * njac(3,5,i) lhsb(4,1,i) = tmp1 * 2.0d+00 * njac(4,1,i) lhsb(4,2,i) = tmp1 * 2.0d+00 * njac(4,2,i) lhsb(4,3,i) = tmp1 * 2.0d+00 * njac(4,3,i) lhsb(4,4,i) = 1.0d+00 > + tmp1 * 2.0d+00 * njac(4,4,i) > + tmp1 * 2.0d+00 * dx4 lhsb(4,5,i) = tmp1 * 2.0d+00 * njac(4,5,i) lhsb(5,1,i) = tmp1 * 2.0d+00 * njac(5,1,i) lhsb(5,2,i) = tmp1 * 2.0d+00 * njac(5,2,i) lhsb(5,3,i) = tmp1 * 2.0d+00 * njac(5,3,i) lhsb(5,4,i) = tmp1 * 2.0d+00 * njac(5,4,i) lhsb(5,5,i) = 1.0d+00 > + tmp1 * 2.0d+00 * njac(5,5,i) > + tmp1 * 2.0d+00 * dx5 lhsc(1,1,i,j,k,c) = tmp2 * fjac(1,1,i+1) > - tmp1 * njac(1,1,i+1) > - tmp1 * dx1 lhsc(1,2,i,j,k,c) = tmp2 * fjac(1,2,i+1) > - tmp1 * njac(1,2,i+1) lhsc(1,3,i,j,k,c) = tmp2 * fjac(1,3,i+1) > - tmp1 * njac(1,3,i+1) lhsc(1,4,i,j,k,c) = tmp2 * fjac(1,4,i+1) > - tmp1 * njac(1,4,i+1) lhsc(1,5,i,j,k,c) = tmp2 * fjac(1,5,i+1) > - tmp1 * njac(1,5,i+1) lhsc(2,1,i,j,k,c) = tmp2 * fjac(2,1,i+1) > - tmp1 * njac(2,1,i+1) lhsc(2,2,i,j,k,c) = tmp2 * fjac(2,2,i+1) > - tmp1 * njac(2,2,i+1) > - tmp1 * dx2 lhsc(2,3,i,j,k,c) = tmp2 * fjac(2,3,i+1) > - tmp1 * njac(2,3,i+1) lhsc(2,4,i,j,k,c) = tmp2 * fjac(2,4,i+1) > - tmp1 * njac(2,4,i+1) lhsc(2,5,i,j,k,c) = tmp2 * fjac(2,5,i+1) > - tmp1 * njac(2,5,i+1) lhsc(3,1,i,j,k,c) = tmp2 * fjac(3,1,i+1) > - tmp1 * njac(3,1,i+1) lhsc(3,2,i,j,k,c) = tmp2 * fjac(3,2,i+1) > - tmp1 * njac(3,2,i+1) lhsc(3,3,i,j,k,c) = tmp2 * fjac(3,3,i+1) > - tmp1 * njac(3,3,i+1) > - tmp1 * dx3 lhsc(3,4,i,j,k,c) = tmp2 * fjac(3,4,i+1) > - tmp1 * njac(3,4,i+1) lhsc(3,5,i,j,k,c) = tmp2 * fjac(3,5,i+1) > - tmp1 * njac(3,5,i+1) lhsc(4,1,i,j,k,c) = tmp2 * fjac(4,1,i+1) > - tmp1 * njac(4,1,i+1) lhsc(4,2,i,j,k,c) = tmp2 * fjac(4,2,i+1) > - tmp1 * njac(4,2,i+1) lhsc(4,3,i,j,k,c) = tmp2 * fjac(4,3,i+1) > - tmp1 * njac(4,3,i+1) lhsc(4,4,i,j,k,c) = tmp2 * fjac(4,4,i+1) > - tmp1 * njac(4,4,i+1) > - tmp1 * dx4 lhsc(4,5,i,j,k,c) = tmp2 * fjac(4,5,i+1) > - tmp1 * njac(4,5,i+1) lhsc(5,1,i,j,k,c) = tmp2 * fjac(5,1,i+1) > - tmp1 * njac(5,1,i+1) lhsc(5,2,i,j,k,c) = tmp2 * fjac(5,2,i+1) > - tmp1 * njac(5,2,i+1) lhsc(5,3,i,j,k,c) = tmp2 * fjac(5,3,i+1) > - tmp1 * njac(5,3,i+1) lhsc(5,4,i,j,k,c) = tmp2 * fjac(5,4,i+1) > - tmp1 * njac(5,4,i+1) lhsc(5,5,i,j,k,c) = tmp2 * fjac(5,5,i+1) > - tmp1 * njac(5,5,i+1) > - tmp1 * dx5 enddo c--------------------------------------------------------------------- c outer most do loops - sweeping in i direction c--------------------------------------------------------------------- if (first .eq. 1) then c--------------------------------------------------------------------- c multiply c(istart,j,k) by b_inverse and copy back to c c multiply rhs(istart) by b_inverse(istart) and copy to rhs c--------------------------------------------------------------------- call binvcrhs( lhsb(1,1,istart), > lhsc(1,1,istart,j,k,c), > rhs(1,istart,j,k,c) ) endif c--------------------------------------------------------------------- c begin inner most do loop c do all the elements of the cell unless last c--------------------------------------------------------------------- do i=istart+first,isize-last c--------------------------------------------------------------------- c rhs(i) = rhs(i) - A*rhs(i-1) c--------------------------------------------------------------------- call matvec_sub(lhsa(1,1,i), > rhs(1,i-1,j,k,c),rhs(1,i,j,k,c)) c--------------------------------------------------------------------- c B(i) = B(i) - C(i-1)*A(i) c--------------------------------------------------------------------- call matmul_sub(lhsa(1,1,i), > lhsc(1,1,i-1,j,k,c), > lhsb(1,1,i)) c--------------------------------------------------------------------- c multiply c(i,j,k) by b_inverse and copy back to c c multiply rhs(1,j,k) by b_inverse(1,j,k) and copy to rhs c--------------------------------------------------------------------- call binvcrhs( lhsb(1,1,i), > lhsc(1,1,i,j,k,c), > rhs(1,i,j,k,c) ) enddo c--------------------------------------------------------------------- c Now finish up special cases for last cell c--------------------------------------------------------------------- if (last .eq. 1) then c--------------------------------------------------------------------- c rhs(isize) = rhs(isize) - A*rhs(isize-1) c--------------------------------------------------------------------- call matvec_sub(lhsa(1,1,isize), > rhs(1,isize-1,j,k,c),rhs(1,isize,j,k,c)) c--------------------------------------------------------------------- c B(isize) = B(isize) - C(isize-1)*A(isize) c--------------------------------------------------------------------- call matmul_sub(lhsa(1,1,isize), > lhsc(1,1,isize-1,j,k,c), > lhsb(1,1,isize)) c--------------------------------------------------------------------- c multiply rhs() by b_inverse() and copy to rhs c--------------------------------------------------------------------- call binvrhs( lhsb(1,1,isize), > rhs(1,isize,j,k,c) ) endif enddo enddo return end