examples/smpi/NAS/BT/y_solve.f

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