examples/smpi/NAS/BT/x_solve.f

   1
   2 c---------------------------------------------------------------------
   3 c---------------------------------------------------------------------
   4
   5       subroutine x_solve
   6
   7 c---------------------------------------------------------------------
   8 c---------------------------------------------------------------------
   9
  10 c---------------------------------------------------------------------
  11 c
  12 c     Performs line solves in X direction by first factoring
  13 c     the block-tridiagonal matrix into an upper triangular matrix,
  14 c     and then performing back substitution to solve for the unknow
  15 c     vectors of each line.
  16 c
  17 c     Make sure we treat elements zero to cell_size in the direction
  18 c     of the sweep.
  19 c
  20 c---------------------------------------------------------------------
  21
  22       include 'header.h'
  23       include 'mpinpb.h'
  24       integer  c, istart, stage,
  25      >     first, last, recv_id, error, r_status(MPI_STATUS_SIZE),
  26      >     isize,jsize,ksize,send_id
  27
  28       istart = 0
  29
  30 c---------------------------------------------------------------------
  31 c     in our terminology stage is the number of the cell in the x-direction
  32 c     i.e. stage = 1 means the start of the line stage=ncells means end
  33 c---------------------------------------------------------------------
  34       do stage = 1,ncells
  35          c = slice(1,stage)
  36          isize = cell_size(1,c) - 1
  37          jsize = cell_size(2,c) - 1
  38          ksize = cell_size(3,c) - 1
  39
  40 c---------------------------------------------------------------------
  41 c     set last-cell flag
  42 c---------------------------------------------------------------------
  43          if (stage .eq. ncells) then
  44             last = 1
  45          else
  46             last = 0
  47          endif
  48
  49          if (stage .eq. 1) then
  50 c---------------------------------------------------------------------
  51 c     This is the first cell, so solve without receiving data
  52 c---------------------------------------------------------------------
  53             first = 1
  54 c            call lhsx(c)
  55             call x_solve_cell(first,last,c)
  56          else
  57 c---------------------------------------------------------------------
  58 c     Not the first cell of this line, so receive info from
  59 c     processor working on preceeding cell
  60 c---------------------------------------------------------------------
  61             first = 0
  62             call x_receive_solve_info(recv_id,c)
  63 c---------------------------------------------------------------------
  64 c     overlap computations and communications
  65 c---------------------------------------------------------------------
  66 c            call lhsx(c)
  67 c---------------------------------------------------------------------
  68 c     wait for completion
  69 c---------------------------------------------------------------------
  70             call mpi_wait(send_id,r_status,error)
  71             call mpi_wait(recv_id,r_status,error)
  72 c---------------------------------------------------------------------
  73 c     install C'(istart) and rhs'(istart) to be used in this cell
  74 c---------------------------------------------------------------------
  75             call x_unpack_solve_info(c)
  76             call x_solve_cell(first,last,c)
  77          endif
  78
  79          if (last .eq. 0) call x_send_solve_info(send_id,c)
  80       enddo
  81
  82 c---------------------------------------------------------------------
  83 c     now perform backsubstitution in reverse direction
  84 c---------------------------------------------------------------------
  85       do stage = ncells, 1, -1
  86          c = slice(1,stage)
  87          first = 0
  88          last = 0
  89          if (stage .eq. 1) first = 1
  90          if (stage .eq. ncells) then
  91             last = 1
  92 c---------------------------------------------------------------------
  93 c     last cell, so perform back substitute without waiting
  94 c---------------------------------------------------------------------
  95             call x_backsubstitute(first, last,c)
  96          else
  97             call x_receive_backsub_info(recv_id,c)
  98             call mpi_wait(send_id,r_status,error)
  99             call mpi_wait(recv_id,r_status,error)
 100             call x_unpack_backsub_info(c)
 101             call x_backsubstitute(first,last,c)
 102          endif
 103          if (first .eq. 0) call x_send_backsub_info(send_id,c)
 104       enddo
 105
 106
 107       return
 108       end
 109
 110
 111 c---------------------------------------------------------------------
 112 c---------------------------------------------------------------------
 113
 114       subroutine x_unpack_solve_info(c)
 115
 116 c---------------------------------------------------------------------
 117 c---------------------------------------------------------------------
 118
 119 c---------------------------------------------------------------------
 120 c     unpack C'(-1) and rhs'(-1) for
 121 c     all j and k
 122 c---------------------------------------------------------------------
 123
 124       include 'header.h'
 125       integer j,k,m,n,ptr,c,istart
 126
 127       istart = 0
 128       ptr = 0
 129       do k=0,KMAX-1
 130          do j=0,JMAX-1
 131             do m=1,BLOCK_SIZE
 132                do n=1,BLOCK_SIZE
 133                   lhsc(m,n,istart-1,j,k,c) = out_buffer(ptr+n)
 134                enddo
 135                ptr = ptr+BLOCK_SIZE
 136             enddo
 137             do n=1,BLOCK_SIZE
 138                rhs(n,istart-1,j,k,c) = out_buffer(ptr+n)
 139             enddo
 140             ptr = ptr+BLOCK_SIZE
 141          enddo
 142       enddo
 143
 144       return
 145       end
 146
 147 c---------------------------------------------------------------------
 148 c---------------------------------------------------------------------
 149
 150       subroutine x_send_solve_info(send_id,c)
 151
 152 c---------------------------------------------------------------------
 153 c---------------------------------------------------------------------
 154
 155 c---------------------------------------------------------------------
 156 c     pack up and send C'(iend) and rhs'(iend) for
 157 c     all j and k
 158 c---------------------------------------------------------------------
 159
 160       include 'header.h'
 161       include 'mpinpb.h'
 162
 163       integer j,k,m,n,isize,ptr,c,jp,kp
 164       integer error,send_id,buffer_size
 165
 166       isize = cell_size(1,c)-1
 167       jp = cell_coord(2,c) - 1
 168       kp = cell_coord(3,c) - 1
 169       buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*
 170      >     (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE)
 171
 172 c---------------------------------------------------------------------
 173 c     pack up buffer
 174 c---------------------------------------------------------------------
 175       ptr = 0
 176       do k=0,KMAX-1
 177          do j=0,JMAX-1
 178             do m=1,BLOCK_SIZE
 179                do n=1,BLOCK_SIZE
 180                   in_buffer(ptr+n) = lhsc(m,n,isize,j,k,c)
 181                enddo
 182                ptr = ptr+BLOCK_SIZE
 183             enddo
 184             do n=1,BLOCK_SIZE
 185                in_buffer(ptr+n) = rhs(n,isize,j,k,c)
 186             enddo
 187             ptr = ptr+BLOCK_SIZE
 188          enddo
 189       enddo
 190
 191 c---------------------------------------------------------------------
 192 c     send buffer
 193 c---------------------------------------------------------------------
 194       call mpi_isend(in_buffer, buffer_size,
 195      >     dp_type, successor(1),
 196      >     WEST+jp+kp*NCELLS, comm_solve,
 197      >     send_id,error)
 198
 199       return
 200       end
 201
 202 c---------------------------------------------------------------------
 203 c---------------------------------------------------------------------
 204
 205       subroutine x_send_backsub_info(send_id,c)
 206
 207 c---------------------------------------------------------------------
 208 c---------------------------------------------------------------------
 209
 210 c---------------------------------------------------------------------
 211 c     pack up and send U(istart) for all j and k
 212 c---------------------------------------------------------------------
 213
 214       include 'header.h'
 215       include 'mpinpb.h'
 216
 217       integer j,k,n,ptr,c,istart,jp,kp
 218       integer error,send_id,buffer_size
 219
 220 c---------------------------------------------------------------------
 221 c     Send element 0 to previous processor
 222 c---------------------------------------------------------------------
 223       istart = 0
 224       jp = cell_coord(2,c)-1
 225       kp = cell_coord(3,c)-1
 226       buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE
 227       ptr = 0
 228       do k=0,KMAX-1
 229          do j=0,JMAX-1
 230             do n=1,BLOCK_SIZE
 231                in_buffer(ptr+n) = rhs(n,istart,j,k,c)
 232             enddo
 233             ptr = ptr+BLOCK_SIZE
 234          enddo
 235       enddo
 236       call mpi_isend(in_buffer, buffer_size,
 237      >     dp_type, predecessor(1),
 238      >     EAST+jp+kp*NCELLS, comm_solve,
 239      >     send_id,error)
 240
 241       return
 242       end
 243
 244 c---------------------------------------------------------------------
 245 c---------------------------------------------------------------------
 246
 247       subroutine x_unpack_backsub_info(c)
 248
 249 c---------------------------------------------------------------------
 250 c---------------------------------------------------------------------
 251
 252 c---------------------------------------------------------------------
 253 c     unpack U(isize) for all j and k
 254 c---------------------------------------------------------------------
 255
 256       include 'header.h'
 257       integer j,k,n,ptr,c
 258
 259       ptr = 0
 260       do k=0,KMAX-1
 261          do j=0,JMAX-1
 262             do n=1,BLOCK_SIZE
 263                backsub_info(n,j,k,c) = out_buffer(ptr+n)
 264             enddo
 265             ptr = ptr+BLOCK_SIZE
 266          enddo
 267       enddo
 268
 269       return
 270       end
 271
 272 c---------------------------------------------------------------------
 273 c---------------------------------------------------------------------
 274
 275       subroutine x_receive_backsub_info(recv_id,c)
 276
 277 c---------------------------------------------------------------------
 278 c---------------------------------------------------------------------
 279
 280 c---------------------------------------------------------------------
 281 c     post mpi receives
 282 c---------------------------------------------------------------------
 283
 284       include 'header.h'
 285       include 'mpinpb.h'
 286
 287       integer error,recv_id,jp,kp,c,buffer_size
 288       jp = cell_coord(2,c) - 1
 289       kp = cell_coord(3,c) - 1
 290       buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*BLOCK_SIZE
 291       call mpi_irecv(out_buffer, buffer_size,
 292      >     dp_type, successor(1),
 293      >     EAST+jp+kp*NCELLS, comm_solve,
 294      >     recv_id, error)
 295
 296       return
 297       end
 298
 299 c---------------------------------------------------------------------
 300 c---------------------------------------------------------------------
 301
 302       subroutine x_receive_solve_info(recv_id,c)
 303
 304 c---------------------------------------------------------------------
 305 c---------------------------------------------------------------------
 306
 307 c---------------------------------------------------------------------
 308 c     post mpi receives
 309 c---------------------------------------------------------------------
 310
 311       include 'header.h'
 312       include 'mpinpb.h'
 313
 314       integer jp,kp,recv_id,error,c,buffer_size
 315       jp = cell_coord(2,c) - 1
 316       kp = cell_coord(3,c) - 1
 317       buffer_size=MAX_CELL_DIM*MAX_CELL_DIM*
 318      >     (BLOCK_SIZE*BLOCK_SIZE + BLOCK_SIZE)
 319       call mpi_irecv(out_buffer, buffer_size,
 320      >     dp_type, predecessor(1),
 321      >     WEST+jp+kp*NCELLS,  comm_solve,
 322      >     recv_id, error)
 323
 324       return
 325       end
 326
 327 c---------------------------------------------------------------------
 328 c---------------------------------------------------------------------
 329
 330       subroutine x_backsubstitute(first, last, c)
 331
 332 c---------------------------------------------------------------------
 333 c---------------------------------------------------------------------
 334
 335 c---------------------------------------------------------------------
 336 c     back solve: if last cell, then generate U(isize)=rhs(isize)
 337 c     else assume U(isize) is loaded in un pack backsub_info
 338 c     so just use it
 339 c     after call u(istart) will be sent to next cell
 340 c---------------------------------------------------------------------
 341
 342       include 'header.h'
 343
 344       integer first, last, c, i, j, k
 345       integer m,n,isize,jsize,ksize,istart
 346
 347       istart = 0
 348       isize = cell_size(1,c)-1
 349       jsize = cell_size(2,c)-end(2,c)-1
 350       ksize = cell_size(3,c)-end(3,c)-1
 351       if (last .eq. 0) then
 352          do k=start(3,c),ksize
 353             do j=start(2,c),jsize
 354 c---------------------------------------------------------------------
 355 c     U(isize) uses info from previous cell if not last cell
 356 c---------------------------------------------------------------------
 357                do m=1,BLOCK_SIZE
 358                   do n=1,BLOCK_SIZE
 359                      rhs(m,isize,j,k,c) = rhs(m,isize,j,k,c)
 360      >                    - lhsc(m,n,isize,j,k,c)*
 361      >                    backsub_info(n,j,k,c)
 362 c---------------------------------------------------------------------
 363 c     rhs(m,isize,j,k,c) = rhs(m,isize,j,k,c)
 364 c     $                    - lhsc(m,n,isize,j,k,c)*rhs(n,isize+1,j,k,c)
 365 c---------------------------------------------------------------------
 366                   enddo
 367                enddo
 368             enddo
 369          enddo
 370       endif
 371       do k=start(3,c),ksize
 372          do j=start(2,c),jsize
 373             do i=isize-1,istart,-1
 374                do m=1,BLOCK_SIZE
 375                   do n=1,BLOCK_SIZE
 376                      rhs(m,i,j,k,c) = rhs(m,i,j,k,c)
 377      >                    - lhsc(m,n,i,j,k,c)*rhs(n,i+1,j,k,c)
 378                   enddo
 379                enddo
 380             enddo
 381          enddo
 382       enddo
 383
 384       return
 385       end
 386
 387
 388 c---------------------------------------------------------------------
 389 c---------------------------------------------------------------------
 390
 391       subroutine x_solve_cell(first,last,c)
 392
 393 c---------------------------------------------------------------------
 394 c---------------------------------------------------------------------
 395
 396 c---------------------------------------------------------------------
 397 c     performs guaussian elimination on this cell.
 398 c
 399 c     assumes that unpacking routines for non-first cells
 400 c     preload C' and rhs' from previous cell.
 401 c
 402 c     assumed send happens outside this routine, but that
 403 c     c'(IMAX) and rhs'(IMAX) will be sent to next cell
 404 c---------------------------------------------------------------------
 405
 406       include 'header.h'
 407       include 'work_lhs.h'
 408
 409       integer first,last,c
 410       integer i,j,k,isize,ksize,jsize,istart
 411
 412       istart = 0
 413       isize = cell_size(1,c)-1
 414       jsize = cell_size(2,c)-end(2,c)-1
 415       ksize = cell_size(3,c)-end(3,c)-1
 416
 417       call lhsabinit(lhsa, lhsb, isize)
 418
 419       do k=start(3,c),ksize
 420          do j=start(2,c),jsize
 421
 422 c---------------------------------------------------------------------
 423 c     This function computes the left hand side in the xi-direction
 424 c---------------------------------------------------------------------
 425
 426 c---------------------------------------------------------------------
 427 c     determine a (labeled f) and n jacobians for cell c
 428 c---------------------------------------------------------------------
 429             do i = start(1,c)-1, cell_size(1,c) - end(1,c)
 430
 431                tmp1 = rho_i(i,j,k,c)
 432                tmp2 = tmp1 * tmp1
 433                tmp3 = tmp1 * tmp2
 434 c---------------------------------------------------------------------
 435 c
 436 c---------------------------------------------------------------------
 437                fjac(1,1,i) = 0.0d+00
 438                fjac(1,2,i) = 1.0d+00
 439                fjac(1,3,i) = 0.0d+00
 440                fjac(1,4,i) = 0.0d+00
 441                fjac(1,5,i) = 0.0d+00
 442
 443                fjac(2,1,i) = -(u(2,i,j,k,c) * tmp2 *
 444      >              u(2,i,j,k,c))
 445      >              + c2 * qs(i,j,k,c)
 446                fjac(2,2,i) = ( 2.0d+00 - c2 )
 447      >              * ( u(2,i,j,k,c) * tmp1 )
 448                fjac(2,3,i) = - c2 * ( u(3,i,j,k,c) * tmp1 )
 449                fjac(2,4,i) = - c2 * ( u(4,i,j,k,c) * tmp1 )
 450                fjac(2,5,i) = c2
 451
 452                fjac(3,1,i) = - ( u(2,i,j,k,c)*u(3,i,j,k,c) ) * tmp2
 453                fjac(3,2,i) = u(3,i,j,k,c) * tmp1
 454                fjac(3,3,i) = u(2,i,j,k,c) * tmp1
 455                fjac(3,4,i) = 0.0d+00
 456                fjac(3,5,i) = 0.0d+00
 457
 458                fjac(4,1,i) = - ( u(2,i,j,k,c)*u(4,i,j,k,c) ) * tmp2
 459                fjac(4,2,i) = u(4,i,j,k,c) * tmp1
 460                fjac(4,3,i) = 0.0d+00
 461                fjac(4,4,i) = u(2,i,j,k,c) * tmp1
 462                fjac(4,5,i) = 0.0d+00
 463
 464                fjac(5,1,i) = ( c2 * 2.0d0 * qs(i,j,k,c)
 465      >              - c1 * ( u(5,i,j,k,c) * tmp1 ) )
 466      >              * ( u(2,i,j,k,c) * tmp1 )
 467                fjac(5,2,i) = c1 *  u(5,i,j,k,c) * tmp1
 468      >              - c2
 469      >              * ( u(2,i,j,k,c)*u(2,i,j,k,c) * tmp2
 470      >              + qs(i,j,k,c) )
 471                fjac(5,3,i) = - c2 * ( u(3,i,j,k,c)*u(2,i,j,k,c) )
 472      >              * tmp2
 473                fjac(5,4,i) = - c2 * ( u(4,i,j,k,c)*u(2,i,j,k,c) )
 474      >              * tmp2
 475                fjac(5,5,i) = c1 * ( u(2,i,j,k,c) * tmp1 )
 476
 477                njac(1,1,i) = 0.0d+00
 478                njac(1,2,i) = 0.0d+00
 479                njac(1,3,i) = 0.0d+00
 480                njac(1,4,i) = 0.0d+00
 481                njac(1,5,i) = 0.0d+00
 482
 483                njac(2,1,i) = - con43 * c3c4 * tmp2 * u(2,i,j,k,c)
 484                njac(2,2,i) =   con43 * c3c4 * tmp1
 485                njac(2,3,i) =   0.0d+00
 486                njac(2,4,i) =   0.0d+00
 487                njac(2,5,i) =   0.0d+00
 488
 489                njac(3,1,i) = - c3c4 * tmp2 * u(3,i,j,k,c)
 490                njac(3,2,i) =   0.0d+00
 491                njac(3,3,i) =   c3c4 * tmp1
 492                njac(3,4,i) =   0.0d+00
 493                njac(3,5,i) =   0.0d+00
 494
 495                njac(4,1,i) = - c3c4 * tmp2 * u(4,i,j,k,c)
 496                njac(4,2,i) =   0.0d+00
 497                njac(4,3,i) =   0.0d+00
 498                njac(4,4,i) =   c3c4 * tmp1
 499                njac(4,5,i) =   0.0d+00
 500
 501                njac(5,1,i) = - ( con43 * c3c4
 502      >              - c1345 ) * tmp3 * (u(2,i,j,k,c)**2)
 503      >              - ( c3c4 - c1345 ) * tmp3 * (u(3,i,j,k,c)**2)
 504      >              - ( c3c4 - c1345 ) * tmp3 * (u(4,i,j,k,c)**2)
 505      >              - c1345 * tmp2 * u(5,i,j,k,c)
 506
 507                njac(5,2,i) = ( con43 * c3c4
 508      >              - c1345 ) * tmp2 * u(2,i,j,k,c)
 509                njac(5,3,i) = ( c3c4 - c1345 ) * tmp2 * u(3,i,j,k,c)
 510                njac(5,4,i) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k,c)
 511                njac(5,5,i) = ( c1345 ) * tmp1
 512
 513             enddo
 514 c---------------------------------------------------------------------
 515 c     now jacobians set, so form left hand side in x direction
 516 c---------------------------------------------------------------------
 517             do i = start(1,c), isize - end(1,c)
 518
 519                tmp1 = dt * tx1
 520                tmp2 = dt * tx2
 521
 522                lhsa(1,1,i) = - tmp2 * fjac(1,1,i-1)
 523      >              - tmp1 * njac(1,1,i-1)
 524      >              - tmp1 * dx1
 525                lhsa(1,2,i) = - tmp2 * fjac(1,2,i-1)
 526      >              - tmp1 * njac(1,2,i-1)
 527                lhsa(1,3,i) = - tmp2 * fjac(1,3,i-1)
 528      >              - tmp1 * njac(1,3,i-1)
 529                lhsa(1,4,i) = - tmp2 * fjac(1,4,i-1)
 530      >              - tmp1 * njac(1,4,i-1)
 531                lhsa(1,5,i) = - tmp2 * fjac(1,5,i-1)
 532      >              - tmp1 * njac(1,5,i-1)
 533
 534                lhsa(2,1,i) = - tmp2 * fjac(2,1,i-1)
 535      >              - tmp1 * njac(2,1,i-1)
 536                lhsa(2,2,i) = - tmp2 * fjac(2,2,i-1)
 537      >              - tmp1 * njac(2,2,i-1)
 538      >              - tmp1 * dx2
 539                lhsa(2,3,i) = - tmp2 * fjac(2,3,i-1)
 540      >              - tmp1 * njac(2,3,i-1)
 541                lhsa(2,4,i) = - tmp2 * fjac(2,4,i-1)
 542      >              - tmp1 * njac(2,4,i-1)
 543                lhsa(2,5,i) = - tmp2 * fjac(2,5,i-1)
 544      >              - tmp1 * njac(2,5,i-1)
 545
 546                lhsa(3,1,i) = - tmp2 * fjac(3,1,i-1)
 547      >              - tmp1 * njac(3,1,i-1)
 548                lhsa(3,2,i) = - tmp2 * fjac(3,2,i-1)
 549      >              - tmp1 * njac(3,2,i-1)
 550                lhsa(3,3,i) = - tmp2 * fjac(3,3,i-1)
 551      >              - tmp1 * njac(3,3,i-1)
 552      >              - tmp1 * dx3
 553                lhsa(3,4,i) = - tmp2 * fjac(3,4,i-1)
 554      >              - tmp1 * njac(3,4,i-1)
 555                lhsa(3,5,i) = - tmp2 * fjac(3,5,i-1)
 556      >              - tmp1 * njac(3,5,i-1)
 557
 558                lhsa(4,1,i) = - tmp2 * fjac(4,1,i-1)
 559      >              - tmp1 * njac(4,1,i-1)
 560                lhsa(4,2,i) = - tmp2 * fjac(4,2,i-1)
 561      >              - tmp1 * njac(4,2,i-1)
 562                lhsa(4,3,i) = - tmp2 * fjac(4,3,i-1)
 563      >              - tmp1 * njac(4,3,i-1)
 564                lhsa(4,4,i) = - tmp2 * fjac(4,4,i-1)
 565      >              - tmp1 * njac(4,4,i-1)
 566      >              - tmp1 * dx4
 567                lhsa(4,5,i) = - tmp2 * fjac(4,5,i-1)
 568      >              - tmp1 * njac(4,5,i-1)
 569
 570                lhsa(5,1,i) = - tmp2 * fjac(5,1,i-1)
 571      >              - tmp1 * njac(5,1,i-1)
 572                lhsa(5,2,i) = - tmp2 * fjac(5,2,i-1)
 573      >              - tmp1 * njac(5,2,i-1)
 574                lhsa(5,3,i) = - tmp2 * fjac(5,3,i-1)
 575      >              - tmp1 * njac(5,3,i-1)
 576                lhsa(5,4,i) = - tmp2 * fjac(5,4,i-1)
 577      >              - tmp1 * njac(5,4,i-1)
 578                lhsa(5,5,i) = - tmp2 * fjac(5,5,i-1)
 579      >              - tmp1 * njac(5,5,i-1)
 580      >              - tmp1 * dx5
 581
 582                lhsb(1,1,i) = 1.0d+00
 583      >              + tmp1 * 2.0d+00 * njac(1,1,i)
 584      >              + tmp1 * 2.0d+00 * dx1
 585                lhsb(1,2,i) = tmp1 * 2.0d+00 * njac(1,2,i)
 586                lhsb(1,3,i) = tmp1 * 2.0d+00 * njac(1,3,i)
 587                lhsb(1,4,i) = tmp1 * 2.0d+00 * njac(1,4,i)
 588                lhsb(1,5,i) = tmp1 * 2.0d+00 * njac(1,5,i)
 589
 590                lhsb(2,1,i) = tmp1 * 2.0d+00 * njac(2,1,i)
 591                lhsb(2,2,i) = 1.0d+00
 592      >              + tmp1 * 2.0d+00 * njac(2,2,i)
 593      >              + tmp1 * 2.0d+00 * dx2
 594                lhsb(2,3,i) = tmp1 * 2.0d+00 * njac(2,3,i)
 595                lhsb(2,4,i) = tmp1 * 2.0d+00 * njac(2,4,i)
 596                lhsb(2,5,i) = tmp1 * 2.0d+00 * njac(2,5,i)
 597
 598                lhsb(3,1,i) = tmp1 * 2.0d+00 * njac(3,1,i)
 599                lhsb(3,2,i) = tmp1 * 2.0d+00 * njac(3,2,i)
 600                lhsb(3,3,i) = 1.0d+00
 601      >              + tmp1 * 2.0d+00 * njac(3,3,i)
 602      >              + tmp1 * 2.0d+00 * dx3
 603                lhsb(3,4,i) = tmp1 * 2.0d+00 * njac(3,4,i)
 604                lhsb(3,5,i) = tmp1 * 2.0d+00 * njac(3,5,i)
 605
 606                lhsb(4,1,i) = tmp1 * 2.0d+00 * njac(4,1,i)
 607                lhsb(4,2,i) = tmp1 * 2.0d+00 * njac(4,2,i)
 608                lhsb(4,3,i) = tmp1 * 2.0d+00 * njac(4,3,i)
 609                lhsb(4,4,i) = 1.0d+00
 610      >              + tmp1 * 2.0d+00 * njac(4,4,i)
 611      >              + tmp1 * 2.0d+00 * dx4
 612                lhsb(4,5,i) = tmp1 * 2.0d+00 * njac(4,5,i)
 613
 614                lhsb(5,1,i) = tmp1 * 2.0d+00 * njac(5,1,i)
 615                lhsb(5,2,i) = tmp1 * 2.0d+00 * njac(5,2,i)
 616                lhsb(5,3,i) = tmp1 * 2.0d+00 * njac(5,3,i)
 617                lhsb(5,4,i) = tmp1 * 2.0d+00 * njac(5,4,i)
 618                lhsb(5,5,i) = 1.0d+00
 619      >              + tmp1 * 2.0d+00 * njac(5,5,i)
 620      >              + tmp1 * 2.0d+00 * dx5
 621
 622                lhsc(1,1,i,j,k,c) =  tmp2 * fjac(1,1,i+1)
 623      >              - tmp1 * njac(1,1,i+1)
 624      >              - tmp1 * dx1
 625                lhsc(1,2,i,j,k,c) =  tmp2 * fjac(1,2,i+1)
 626      >              - tmp1 * njac(1,2,i+1)
 627                lhsc(1,3,i,j,k,c) =  tmp2 * fjac(1,3,i+1)
 628      >              - tmp1 * njac(1,3,i+1)
 629                lhsc(1,4,i,j,k,c) =  tmp2 * fjac(1,4,i+1)
 630      >              - tmp1 * njac(1,4,i+1)
 631                lhsc(1,5,i,j,k,c) =  tmp2 * fjac(1,5,i+1)
 632      >              - tmp1 * njac(1,5,i+1)
 633
 634                lhsc(2,1,i,j,k,c) =  tmp2 * fjac(2,1,i+1)
 635      >              - tmp1 * njac(2,1,i+1)
 636                lhsc(2,2,i,j,k,c) =  tmp2 * fjac(2,2,i+1)
 637      >              - tmp1 * njac(2,2,i+1)
 638      >              - tmp1 * dx2
 639                lhsc(2,3,i,j,k,c) =  tmp2 * fjac(2,3,i+1)
 640      >              - tmp1 * njac(2,3,i+1)
 641                lhsc(2,4,i,j,k,c) =  tmp2 * fjac(2,4,i+1)
 642      >              - tmp1 * njac(2,4,i+1)
 643                lhsc(2,5,i,j,k,c) =  tmp2 * fjac(2,5,i+1)
 644      >              - tmp1 * njac(2,5,i+1)
 645
 646                lhsc(3,1,i,j,k,c) =  tmp2 * fjac(3,1,i+1)
 647      >              - tmp1 * njac(3,1,i+1)
 648                lhsc(3,2,i,j,k,c) =  tmp2 * fjac(3,2,i+1)
 649      >              - tmp1 * njac(3,2,i+1)
 650                lhsc(3,3,i,j,k,c) =  tmp2 * fjac(3,3,i+1)
 651      >              - tmp1 * njac(3,3,i+1)
 652      >              - tmp1 * dx3
 653                lhsc(3,4,i,j,k,c) =  tmp2 * fjac(3,4,i+1)
 654      >              - tmp1 * njac(3,4,i+1)
 655                lhsc(3,5,i,j,k,c) =  tmp2 * fjac(3,5,i+1)
 656      >              - tmp1 * njac(3,5,i+1)
 657
 658                lhsc(4,1,i,j,k,c) =  tmp2 * fjac(4,1,i+1)
 659      >              - tmp1 * njac(4,1,i+1)
 660                lhsc(4,2,i,j,k,c) =  tmp2 * fjac(4,2,i+1)
 661      >              - tmp1 * njac(4,2,i+1)
 662                lhsc(4,3,i,j,k,c) =  tmp2 * fjac(4,3,i+1)
 663      >              - tmp1 * njac(4,3,i+1)
 664                lhsc(4,4,i,j,k,c) =  tmp2 * fjac(4,4,i+1)
 665      >              - tmp1 * njac(4,4,i+1)
 666      >              - tmp1 * dx4
 667                lhsc(4,5,i,j,k,c) =  tmp2 * fjac(4,5,i+1)
 668      >              - tmp1 * njac(4,5,i+1)
 669
 670                lhsc(5,1,i,j,k,c) =  tmp2 * fjac(5,1,i+1)
 671      >              - tmp1 * njac(5,1,i+1)
 672                lhsc(5,2,i,j,k,c) =  tmp2 * fjac(5,2,i+1)
 673      >              - tmp1 * njac(5,2,i+1)
 674                lhsc(5,3,i,j,k,c) =  tmp2 * fjac(5,3,i+1)
 675      >              - tmp1 * njac(5,3,i+1)
 676                lhsc(5,4,i,j,k,c) =  tmp2 * fjac(5,4,i+1)
 677      >              - tmp1 * njac(5,4,i+1)
 678                lhsc(5,5,i,j,k,c) =  tmp2 * fjac(5,5,i+1)
 679      >              - tmp1 * njac(5,5,i+1)
 680      >              - tmp1 * dx5
 681
 682             enddo
 683
 684
 685 c---------------------------------------------------------------------
 686 c     outer most do loops - sweeping in i direction
 687 c---------------------------------------------------------------------
 688             if (first .eq. 1) then
 689
 690 c---------------------------------------------------------------------
 691 c     multiply c(istart,j,k) by b_inverse and copy back to c
 692 c     multiply rhs(istart) by b_inverse(istart) and copy to rhs
 693 c---------------------------------------------------------------------
 694                call binvcrhs( lhsb(1,1,istart),
 695      >                        lhsc(1,1,istart,j,k,c),
 696      >                        rhs(1,istart,j,k,c) )
 697
 698             endif
 699
 700 c---------------------------------------------------------------------
 701 c     begin inner most do loop
 702 c     do all the elements of the cell unless last
 703 c---------------------------------------------------------------------
 704             do i=istart+first,isize-last
 705
 706 c---------------------------------------------------------------------
 707 c     rhs(i) = rhs(i) - A*rhs(i-1)
 708 c---------------------------------------------------------------------
 709                call matvec_sub(lhsa(1,1,i),
 710      >                         rhs(1,i-1,j,k,c),rhs(1,i,j,k,c))
 711
 712 c---------------------------------------------------------------------
 713 c     B(i) = B(i) - C(i-1)*A(i)
 714 c---------------------------------------------------------------------
 715                call matmul_sub(lhsa(1,1,i),
 716      >                         lhsc(1,1,i-1,j,k,c),
 717      >                         lhsb(1,1,i))
 718
 719
 720 c---------------------------------------------------------------------
 721 c     multiply c(i,j,k) by b_inverse and copy back to c
 722 c     multiply rhs(1,j,k) by b_inverse(1,j,k) and copy to rhs
 723 c---------------------------------------------------------------------
 724                call binvcrhs( lhsb(1,1,i),
 725      >                        lhsc(1,1,i,j,k,c),
 726      >                        rhs(1,i,j,k,c) )
 727
 728             enddo
 729
 730 c---------------------------------------------------------------------
 731 c     Now finish up special cases for last cell
 732 c---------------------------------------------------------------------
 733             if (last .eq. 1) then
 734
 735 c---------------------------------------------------------------------
 736 c     rhs(isize) = rhs(isize) - A*rhs(isize-1)
 737 c---------------------------------------------------------------------
 738                call matvec_sub(lhsa(1,1,isize),
 739      >                         rhs(1,isize-1,j,k,c),rhs(1,isize,j,k,c))
 740
 741 c---------------------------------------------------------------------
 742 c     B(isize) = B(isize) - C(isize-1)*A(isize)
 743 c---------------------------------------------------------------------
 744                call matmul_sub(lhsa(1,1,isize),
 745      >                         lhsc(1,1,isize-1,j,k,c),
 746      >                         lhsb(1,1,isize))
 747
 748 c---------------------------------------------------------------------
 749 c     multiply rhs() by b_inverse() and copy to rhs
 750 c---------------------------------------------------------------------
 751                call binvrhs( lhsb(1,1,isize),
 752      >                       rhs(1,isize,j,k,c) )
 753
 754             endif
 755          enddo
 756       enddo
 757
 758
 759       return
 760       end
 761