2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 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.
17 c Make sure we treat elements zero to cell_size in the direction
20 c---------------------------------------------------------------------
24 integer c, istart, stage,
25 > first, last, recv_id, error, r_status(MPI_STATUS_SIZE),
26 > isize,jsize,ksize,send_id
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---------------------------------------------------------------------
36 isize = cell_size(1,c) - 1
37 jsize = cell_size(2,c) - 1
38 ksize = cell_size(3,c) - 1
40 c---------------------------------------------------------------------
42 c---------------------------------------------------------------------
43 if (stage .eq. ncells) then
49 if (stage .eq. 1) then
50 c---------------------------------------------------------------------
51 c This is the first cell, so solve without receiving data
52 c---------------------------------------------------------------------
55 call x_solve_cell(first,last,c)
57 c---------------------------------------------------------------------
58 c Not the first cell of this line, so receive info from
59 c processor working on preceeding cell
60 c---------------------------------------------------------------------
62 call x_receive_solve_info(recv_id,c)
63 c---------------------------------------------------------------------
64 c overlap computations and communications
65 c---------------------------------------------------------------------
67 c---------------------------------------------------------------------
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)
79 if (last .eq. 0) call x_send_solve_info(send_id,c)
82 c---------------------------------------------------------------------
83 c now perform backsubstitution in reverse direction
84 c---------------------------------------------------------------------
85 do stage = ncells, 1, -1
89 if (stage .eq. 1) first = 1
90 if (stage .eq. ncells) then
92 c---------------------------------------------------------------------
93 c last cell, so perform back substitute without waiting
94 c---------------------------------------------------------------------
95 call x_backsubstitute(first, last,c)
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)
103 if (first .eq. 0) call x_send_backsub_info(send_id,c)
111 c---------------------------------------------------------------------
112 c---------------------------------------------------------------------
114 subroutine x_unpack_solve_info(c)
116 c---------------------------------------------------------------------
117 c---------------------------------------------------------------------
119 c---------------------------------------------------------------------
120 c unpack C'(-1) and rhs'(-1) for
122 c---------------------------------------------------------------------
125 integer j,k,m,n,ptr,c,istart
133 lhsc(m,n,istart-1,j,k,c) = out_buffer(ptr+n)
138 rhs(n,istart-1,j,k,c) = out_buffer(ptr+n)
147 c---------------------------------------------------------------------
148 c---------------------------------------------------------------------
150 subroutine x_send_solve_info(send_id,c)
152 c---------------------------------------------------------------------
153 c---------------------------------------------------------------------
155 c---------------------------------------------------------------------
156 c pack up and send C'(iend) and rhs'(iend) for
158 c---------------------------------------------------------------------
163 integer j,k,m,n,isize,ptr,c,jp,kp
164 integer error,send_id,buffer_size
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)
172 c---------------------------------------------------------------------
174 c---------------------------------------------------------------------
180 in_buffer(ptr+n) = lhsc(m,n,isize,j,k,c)
185 in_buffer(ptr+n) = rhs(n,isize,j,k,c)
191 c---------------------------------------------------------------------
193 c---------------------------------------------------------------------
194 call mpi_isend(in_buffer, buffer_size,
195 > dp_type, successor(1),
196 > WEST+jp+kp*NCELLS, comm_solve,
202 c---------------------------------------------------------------------
203 c---------------------------------------------------------------------
205 subroutine x_send_backsub_info(send_id,c)
207 c---------------------------------------------------------------------
208 c---------------------------------------------------------------------
210 c---------------------------------------------------------------------
211 c pack up and send U(istart) for all j and k
212 c---------------------------------------------------------------------
217 integer j,k,n,ptr,c,istart,jp,kp
218 integer error,send_id,buffer_size
220 c---------------------------------------------------------------------
221 c Send element 0 to previous processor
222 c---------------------------------------------------------------------
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
231 in_buffer(ptr+n) = rhs(n,istart,j,k,c)
236 call mpi_isend(in_buffer, buffer_size,
237 > dp_type, predecessor(1),
238 > EAST+jp+kp*NCELLS, comm_solve,
244 c---------------------------------------------------------------------
245 c---------------------------------------------------------------------
247 subroutine x_unpack_backsub_info(c)
249 c---------------------------------------------------------------------
250 c---------------------------------------------------------------------
252 c---------------------------------------------------------------------
253 c unpack U(isize) for all j and k
254 c---------------------------------------------------------------------
263 backsub_info(n,j,k,c) = out_buffer(ptr+n)
272 c---------------------------------------------------------------------
273 c---------------------------------------------------------------------
275 subroutine x_receive_backsub_info(recv_id,c)
277 c---------------------------------------------------------------------
278 c---------------------------------------------------------------------
280 c---------------------------------------------------------------------
282 c---------------------------------------------------------------------
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,
299 c---------------------------------------------------------------------
300 c---------------------------------------------------------------------
302 subroutine x_receive_solve_info(recv_id,c)
304 c---------------------------------------------------------------------
305 c---------------------------------------------------------------------
307 c---------------------------------------------------------------------
309 c---------------------------------------------------------------------
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,
327 c---------------------------------------------------------------------
328 c---------------------------------------------------------------------
330 subroutine x_backsubstitute(first, last, c)
332 c---------------------------------------------------------------------
333 c---------------------------------------------------------------------
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
339 c after call u(istart) will be sent to next cell
340 c---------------------------------------------------------------------
344 integer first, last, c, i, j, k
345 integer m,n,isize,jsize,ksize,istart
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---------------------------------------------------------------------
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---------------------------------------------------------------------
371 do k=start(3,c),ksize
372 do j=start(2,c),jsize
373 do i=isize-1,istart,-1
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)
388 c---------------------------------------------------------------------
389 c---------------------------------------------------------------------
391 subroutine x_solve_cell(first,last,c)
393 c---------------------------------------------------------------------
394 c---------------------------------------------------------------------
396 c---------------------------------------------------------------------
397 c performs guaussian elimination on this cell.
399 c assumes that unpacking routines for non-first cells
400 c preload C' and rhs' from previous cell.
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---------------------------------------------------------------------
410 integer i,j,k,isize,ksize,jsize,istart
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
417 call lhsabinit(lhsa, lhsb, isize)
419 do k=start(3,c),ksize
420 do j=start(2,c),jsize
422 c---------------------------------------------------------------------
423 c This function computes the left hand side in the xi-direction
424 c---------------------------------------------------------------------
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)
431 tmp1 = rho_i(i,j,k,c)
434 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
443 fjac(2,1,i) = -(u(2,i,j,k,c) * tmp2 *
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 )
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
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
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
469 > * ( u(2,i,j,k,c)*u(2,i,j,k,c) * tmp2
471 fjac(5,3,i) = - c2 * ( u(3,i,j,k,c)*u(2,i,j,k,c) )
473 fjac(5,4,i) = - c2 * ( u(4,i,j,k,c)*u(2,i,j,k,c) )
475 fjac(5,5,i) = c1 * ( u(2,i,j,k,c) * tmp1 )
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
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
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
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
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)
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
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)
522 lhsa(1,1,i) = - tmp2 * fjac(1,1,i-1)
523 > - tmp1 * njac(1,1,i-1)
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)
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)
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)
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)
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)
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)
567 lhsa(4,5,i) = - tmp2 * fjac(4,5,i-1)
568 > - tmp1 * njac(4,5,i-1)
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)
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)
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)
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)
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)
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
622 lhsc(1,1,i,j,k,c) = tmp2 * fjac(1,1,i+1)
623 > - tmp1 * njac(1,1,i+1)
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)
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)
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)
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)
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)
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)
667 lhsc(4,5,i,j,k,c) = tmp2 * fjac(4,5,i+1)
668 > - tmp1 * njac(4,5,i+1)
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)
685 c---------------------------------------------------------------------
686 c outer most do loops - sweeping in i direction
687 c---------------------------------------------------------------------
688 if (first .eq. 1) then
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) )
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
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))
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),
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),
730 c---------------------------------------------------------------------
731 c Now finish up special cases for last cell
732 c---------------------------------------------------------------------
733 if (last .eq. 1) then
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))
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),
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) )