1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
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.
15 c Make sure we treat elements zero to cell_size in the direction
17 c---------------------------------------------------------------------
24 > first, last, recv_id, error, r_status(MPI_STATUS_SIZE),
25 > isize,jsize,ksize,send_id
29 c---------------------------------------------------------------------
30 c in our terminology stage is the number of the cell in the y-direct
31 c i.e. stage = 1 means the start of the line stage=ncells means end
32 c---------------------------------------------------------------------
35 isize = cell_size(1,c) - 1
36 jsize = cell_size(2,c) - 1
37 ksize = cell_size(3,c) - 1
39 c---------------------------------------------------------------------
41 c---------------------------------------------------------------------
42 if (stage .eq. ncells) then
48 if (stage .eq. 1) then
49 c---------------------------------------------------------------------
50 c This is the first cell, so solve without receiving data
51 c---------------------------------------------------------------------
54 call y_solve_cell(first,last,c)
56 c---------------------------------------------------------------------
57 c Not the first cell of this line, so receive info from
58 c processor working on preceeding cell
59 c---------------------------------------------------------------------
61 call y_receive_solve_info(recv_id,c)
62 c---------------------------------------------------------------------
63 c overlap computations and communications
64 c---------------------------------------------------------------------
66 c---------------------------------------------------------------------
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)
78 if (last .eq. 0) call y_send_solve_info(send_id,c)
81 c---------------------------------------------------------------------
82 c now perform backsubstitution in reverse direction
83 c---------------------------------------------------------------------
84 do stage = ncells, 1, -1
88 if (stage .eq. 1) first = 1
89 if (stage .eq. ncells) then
91 c---------------------------------------------------------------------
92 c last cell, so perform back substitute without waiting
93 c---------------------------------------------------------------------
94 call y_backsubstitute(first, last,c)
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)
102 if (first .eq. 0) call y_send_backsub_info(send_id,c)
109 c---------------------------------------------------------------------
110 c---------------------------------------------------------------------
112 subroutine y_unpack_solve_info(c)
114 c---------------------------------------------------------------------
115 c---------------------------------------------------------------------
117 c---------------------------------------------------------------------
118 c unpack C'(-1) and rhs'(-1) for
120 c---------------------------------------------------------------------
124 integer i,k,m,n,ptr,c,jstart
132 lhsc(m,n,i,jstart-1,k,c) = out_buffer(ptr+n)
137 rhs(n,i,jstart-1,k,c) = out_buffer(ptr+n)
146 c---------------------------------------------------------------------
147 c---------------------------------------------------------------------
149 subroutine y_send_solve_info(send_id,c)
151 c---------------------------------------------------------------------
152 c---------------------------------------------------------------------
154 c---------------------------------------------------------------------
155 c pack up and send C'(jend) and rhs'(jend) for
157 c---------------------------------------------------------------------
162 integer i,k,m,n,jsize,ptr,c,ip,kp
163 integer error,send_id,buffer_size
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)
171 c---------------------------------------------------------------------
173 c---------------------------------------------------------------------
179 in_buffer(ptr+n) = lhsc(m,n,i,jsize,k,c)
184 in_buffer(ptr+n) = rhs(n,i,jsize,k,c)
190 c---------------------------------------------------------------------
192 c---------------------------------------------------------------------
193 call mpi_isend(in_buffer, buffer_size,
194 > dp_type, successor(2),
195 > SOUTH+ip+kp*NCELLS, comm_solve,
201 c---------------------------------------------------------------------
202 c---------------------------------------------------------------------
204 subroutine y_send_backsub_info(send_id,c)
206 c---------------------------------------------------------------------
207 c---------------------------------------------------------------------
209 c---------------------------------------------------------------------
210 c pack up and send U(jstart) for all i and k
211 c---------------------------------------------------------------------
216 integer i,k,n,ptr,c,jstart,ip,kp
217 integer error,send_id,buffer_size
219 c---------------------------------------------------------------------
220 c Send element 0 to previous processor
221 c---------------------------------------------------------------------
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
230 in_buffer(ptr+n) = rhs(n,i,jstart,k,c)
235 call mpi_isend(in_buffer, buffer_size,
236 > dp_type, predecessor(2),
237 > NORTH+ip+kp*NCELLS, comm_solve,
243 c---------------------------------------------------------------------
244 c---------------------------------------------------------------------
246 subroutine y_unpack_backsub_info(c)
248 c---------------------------------------------------------------------
249 c---------------------------------------------------------------------
251 c---------------------------------------------------------------------
252 c unpack U(jsize) for all i and k
253 c---------------------------------------------------------------------
263 backsub_info(n,i,k,c) = out_buffer(ptr+n)
272 c---------------------------------------------------------------------
273 c---------------------------------------------------------------------
275 subroutine y_receive_backsub_info(recv_id,c)
277 c---------------------------------------------------------------------
278 c---------------------------------------------------------------------
280 c---------------------------------------------------------------------
282 c---------------------------------------------------------------------
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,
298 c---------------------------------------------------------------------
299 c---------------------------------------------------------------------
301 subroutine y_receive_solve_info(recv_id,c)
303 c---------------------------------------------------------------------
304 c---------------------------------------------------------------------
306 c---------------------------------------------------------------------
308 c---------------------------------------------------------------------
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,
326 c---------------------------------------------------------------------
327 c---------------------------------------------------------------------
329 subroutine y_backsubstitute(first, last, c)
331 c---------------------------------------------------------------------
332 c---------------------------------------------------------------------
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
338 c after call u(jstart) will be sent to next cell
339 c---------------------------------------------------------------------
343 integer first, last, c, i, k
344 integer m,n,j,jsize,isize,ksize,jstart
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---------------------------------------------------------------------
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)
366 do k=start(3,c),ksize
367 do j=jsize-1,jstart,-1
368 do i=start(1,c),isize
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)
382 c---------------------------------------------------------------------
383 c---------------------------------------------------------------------
385 subroutine y_solve_cell(first,last,c)
387 c---------------------------------------------------------------------
388 c---------------------------------------------------------------------
390 c---------------------------------------------------------------------
391 c performs guaussian elimination on this cell.
393 c assumes that unpacking routines for non-first cells
394 c preload C' and rhs' from previous cell.
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---------------------------------------------------------------------
401 include 'work_lhs_vec.h'
404 integer i,j,k,m,n,isize,ksize,jsize,jstart
407 isize = cell_size(1,c)-end(1,c)-1
408 jsize = cell_size(2,c)-1
409 ksize = cell_size(3,c)-end(3,c)-1
411 c---------------------------------------------------------------------
412 c zero the left hand side for starters
413 c set diagonal values to 1. This is overkill, but convenient
414 c---------------------------------------------------------------------
418 lhsa(m,n,i,0) = 0.0d0
419 lhsb(m,n,i,0) = 0.0d0
420 lhsa(m,n,i,jsize) = 0.0d0
421 lhsb(m,n,i,jsize) = 0.0d0
423 lhsb(m,m,i,0) = 1.0d0
424 lhsb(m,m,i,jsize) = 1.0d0
428 do k=start(3,c),ksize
430 c---------------------------------------------------------------------
431 c This function computes the left hand side for the three y-factors
432 c---------------------------------------------------------------------
434 c---------------------------------------------------------------------
435 c Compute the indices for storing the tri-diagonal matrix;
436 c determine a (labeled f) and n jacobians for cell c
437 c---------------------------------------------------------------------
439 do j = start(2,c)-1, cell_size(2,c)-end(2,c)
440 do i=start(1,c),isize
442 tmp1 = 1.0d0 / u(1,i,j,k,c)
446 fjac(1,1,i,j) = 0.0d+00
447 fjac(1,2,i,j) = 0.0d+00
448 fjac(1,3,i,j) = 1.0d+00
449 fjac(1,4,i,j) = 0.0d+00
450 fjac(1,5,i,j) = 0.0d+00
452 fjac(2,1,i,j) = - ( u(2,i,j,k,c)*u(3,i,j,k,c) )
454 fjac(2,2,i,j) = u(3,i,j,k,c) * tmp1
455 fjac(2,3,i,j) = u(2,i,j,k,c) * tmp1
456 fjac(2,4,i,j) = 0.0d+00
457 fjac(2,5,i,j) = 0.0d+00
459 fjac(3,1,i,j) = - ( u(3,i,j,k,c)*u(3,i,j,k,c)*tmp2)
461 fjac(3,2,i,j) = - c2 * u(2,i,j,k,c) * tmp1
462 fjac(3,3,i,j) = ( 2.0d+00 - c2 )
463 > * u(3,i,j,k,c) * tmp1
464 fjac(3,4,i,j) = - c2 * u(4,i,j,k,c) * tmp1
467 fjac(4,1,i,j) = - ( u(3,i,j,k,c)*u(4,i,j,k,c) )
469 fjac(4,2,i,j) = 0.0d+00
470 fjac(4,3,i,j) = u(4,i,j,k,c) * tmp1
471 fjac(4,4,i,j) = u(3,i,j,k,c) * tmp1
472 fjac(4,5,i,j) = 0.0d+00
474 fjac(5,1,i,j) = ( c2 * 2.0d0 * qs(i,j,k,c)
475 > - c1 * u(5,i,j,k,c) * tmp1 )
476 > * u(3,i,j,k,c) * tmp1
477 fjac(5,2,i,j) = - c2 * u(2,i,j,k,c)*u(3,i,j,k,c)
479 fjac(5,3,i,j) = c1 * u(5,i,j,k,c) * tmp1
480 > - c2 * ( qs(i,j,k,c)
481 > + u(3,i,j,k,c)*u(3,i,j,k,c) * tmp2 )
482 fjac(5,4,i,j) = - c2 * ( u(3,i,j,k,c)*u(4,i,j,k,c) )
484 fjac(5,5,i,j) = c1 * u(3,i,j,k,c) * tmp1
486 njac(1,1,i,j) = 0.0d+00
487 njac(1,2,i,j) = 0.0d+00
488 njac(1,3,i,j) = 0.0d+00
489 njac(1,4,i,j) = 0.0d+00
490 njac(1,5,i,j) = 0.0d+00
492 njac(2,1,i,j) = - c3c4 * tmp2 * u(2,i,j,k,c)
493 njac(2,2,i,j) = c3c4 * tmp1
494 njac(2,3,i,j) = 0.0d+00
495 njac(2,4,i,j) = 0.0d+00
496 njac(2,5,i,j) = 0.0d+00
498 njac(3,1,i,j) = - con43 * c3c4 * tmp2 * u(3,i,j,k,c)
499 njac(3,2,i,j) = 0.0d+00
500 njac(3,3,i,j) = con43 * c3c4 * tmp1
501 njac(3,4,i,j) = 0.0d+00
502 njac(3,5,i,j) = 0.0d+00
504 njac(4,1,i,j) = - c3c4 * tmp2 * u(4,i,j,k,c)
505 njac(4,2,i,j) = 0.0d+00
506 njac(4,3,i,j) = 0.0d+00
507 njac(4,4,i,j) = c3c4 * tmp1
508 njac(4,5,i,j) = 0.0d+00
510 njac(5,1,i,j) = - ( c3c4
511 > - c1345 ) * tmp3 * (u(2,i,j,k,c)**2)
513 > - c1345 ) * tmp3 * (u(3,i,j,k,c)**2)
514 > - ( c3c4 - c1345 ) * tmp3 * (u(4,i,j,k,c)**2)
515 > - c1345 * tmp2 * u(5,i,j,k,c)
517 njac(5,2,i,j) = ( c3c4 - c1345 ) * tmp2 * u(2,i,j,k,c)
518 njac(5,3,i,j) = ( con43 * c3c4
519 > - c1345 ) * tmp2 * u(3,i,j,k,c)
520 njac(5,4,i,j) = ( c3c4 - c1345 ) * tmp2 * u(4,i,j,k,c)
521 njac(5,5,i,j) = ( c1345 ) * tmp1
526 c---------------------------------------------------------------------
527 c now joacobians set, so form left hand side in y direction
528 c---------------------------------------------------------------------
529 do j = start(2,c), jsize-end(2,c)
530 do i=start(1,c),isize
535 lhsa(1,1,i,j) = - tmp2 * fjac(1,1,i,j-1)
536 > - tmp1 * njac(1,1,i,j-1)
538 lhsa(1,2,i,j) = - tmp2 * fjac(1,2,i,j-1)
539 > - tmp1 * njac(1,2,i,j-1)
540 lhsa(1,3,i,j) = - tmp2 * fjac(1,3,i,j-1)
541 > - tmp1 * njac(1,3,i,j-1)
542 lhsa(1,4,i,j) = - tmp2 * fjac(1,4,i,j-1)
543 > - tmp1 * njac(1,4,i,j-1)
544 lhsa(1,5,i,j) = - tmp2 * fjac(1,5,i,j-1)
545 > - tmp1 * njac(1,5,i,j-1)
547 lhsa(2,1,i,j) = - tmp2 * fjac(2,1,i,j-1)
548 > - tmp1 * njac(2,1,i,j-1)
549 lhsa(2,2,i,j) = - tmp2 * fjac(2,2,i,j-1)
550 > - tmp1 * njac(2,2,i,j-1)
552 lhsa(2,3,i,j) = - tmp2 * fjac(2,3,i,j-1)
553 > - tmp1 * njac(2,3,i,j-1)
554 lhsa(2,4,i,j) = - tmp2 * fjac(2,4,i,j-1)
555 > - tmp1 * njac(2,4,i,j-1)
556 lhsa(2,5,i,j) = - tmp2 * fjac(2,5,i,j-1)
557 > - tmp1 * njac(2,5,i,j-1)
559 lhsa(3,1,i,j) = - tmp2 * fjac(3,1,i,j-1)
560 > - tmp1 * njac(3,1,i,j-1)
561 lhsa(3,2,i,j) = - tmp2 * fjac(3,2,i,j-1)
562 > - tmp1 * njac(3,2,i,j-1)
563 lhsa(3,3,i,j) = - tmp2 * fjac(3,3,i,j-1)
564 > - tmp1 * njac(3,3,i,j-1)
566 lhsa(3,4,i,j) = - tmp2 * fjac(3,4,i,j-1)
567 > - tmp1 * njac(3,4,i,j-1)
568 lhsa(3,5,i,j) = - tmp2 * fjac(3,5,i,j-1)
569 > - tmp1 * njac(3,5,i,j-1)
571 lhsa(4,1,i,j) = - tmp2 * fjac(4,1,i,j-1)
572 > - tmp1 * njac(4,1,i,j-1)
573 lhsa(4,2,i,j) = - tmp2 * fjac(4,2,i,j-1)
574 > - tmp1 * njac(4,2,i,j-1)
575 lhsa(4,3,i,j) = - tmp2 * fjac(4,3,i,j-1)
576 > - tmp1 * njac(4,3,i,j-1)
577 lhsa(4,4,i,j) = - tmp2 * fjac(4,4,i,j-1)
578 > - tmp1 * njac(4,4,i,j-1)
580 lhsa(4,5,i,j) = - tmp2 * fjac(4,5,i,j-1)
581 > - tmp1 * njac(4,5,i,j-1)
583 lhsa(5,1,i,j) = - tmp2 * fjac(5,1,i,j-1)
584 > - tmp1 * njac(5,1,i,j-1)
585 lhsa(5,2,i,j) = - tmp2 * fjac(5,2,i,j-1)
586 > - tmp1 * njac(5,2,i,j-1)
587 lhsa(5,3,i,j) = - tmp2 * fjac(5,3,i,j-1)
588 > - tmp1 * njac(5,3,i,j-1)
589 lhsa(5,4,i,j) = - tmp2 * fjac(5,4,i,j-1)
590 > - tmp1 * njac(5,4,i,j-1)
591 lhsa(5,5,i,j) = - tmp2 * fjac(5,5,i,j-1)
592 > - tmp1 * njac(5,5,i,j-1)
595 lhsb(1,1,i,j) = 1.0d+00
596 > + tmp1 * 2.0d+00 * njac(1,1,i,j)
597 > + tmp1 * 2.0d+00 * dy1
598 lhsb(1,2,i,j) = tmp1 * 2.0d+00 * njac(1,2,i,j)
599 lhsb(1,3,i,j) = tmp1 * 2.0d+00 * njac(1,3,i,j)
600 lhsb(1,4,i,j) = tmp1 * 2.0d+00 * njac(1,4,i,j)
601 lhsb(1,5,i,j) = tmp1 * 2.0d+00 * njac(1,5,i,j)
603 lhsb(2,1,i,j) = tmp1 * 2.0d+00 * njac(2,1,i,j)
604 lhsb(2,2,i,j) = 1.0d+00
605 > + tmp1 * 2.0d+00 * njac(2,2,i,j)
606 > + tmp1 * 2.0d+00 * dy2
607 lhsb(2,3,i,j) = tmp1 * 2.0d+00 * njac(2,3,i,j)
608 lhsb(2,4,i,j) = tmp1 * 2.0d+00 * njac(2,4,i,j)
609 lhsb(2,5,i,j) = tmp1 * 2.0d+00 * njac(2,5,i,j)
611 lhsb(3,1,i,j) = tmp1 * 2.0d+00 * njac(3,1,i,j)
612 lhsb(3,2,i,j) = tmp1 * 2.0d+00 * njac(3,2,i,j)
613 lhsb(3,3,i,j) = 1.0d+00
614 > + tmp1 * 2.0d+00 * njac(3,3,i,j)
615 > + tmp1 * 2.0d+00 * dy3
616 lhsb(3,4,i,j) = tmp1 * 2.0d+00 * njac(3,4,i,j)
617 lhsb(3,5,i,j) = tmp1 * 2.0d+00 * njac(3,5,i,j)
619 lhsb(4,1,i,j) = tmp1 * 2.0d+00 * njac(4,1,i,j)
620 lhsb(4,2,i,j) = tmp1 * 2.0d+00 * njac(4,2,i,j)
621 lhsb(4,3,i,j) = tmp1 * 2.0d+00 * njac(4,3,i,j)
622 lhsb(4,4,i,j) = 1.0d+00
623 > + tmp1 * 2.0d+00 * njac(4,4,i,j)
624 > + tmp1 * 2.0d+00 * dy4
625 lhsb(4,5,i,j) = tmp1 * 2.0d+00 * njac(4,5,i,j)
627 lhsb(5,1,i,j) = tmp1 * 2.0d+00 * njac(5,1,i,j)
628 lhsb(5,2,i,j) = tmp1 * 2.0d+00 * njac(5,2,i,j)
629 lhsb(5,3,i,j) = tmp1 * 2.0d+00 * njac(5,3,i,j)
630 lhsb(5,4,i,j) = tmp1 * 2.0d+00 * njac(5,4,i,j)
631 lhsb(5,5,i,j) = 1.0d+00
632 > + tmp1 * 2.0d+00 * njac(5,5,i,j)
633 > + tmp1 * 2.0d+00 * dy5
635 lhsc(1,1,i,j,k,c) = tmp2 * fjac(1,1,i,j+1)
636 > - tmp1 * njac(1,1,i,j+1)
638 lhsc(1,2,i,j,k,c) = tmp2 * fjac(1,2,i,j+1)
639 > - tmp1 * njac(1,2,i,j+1)
640 lhsc(1,3,i,j,k,c) = tmp2 * fjac(1,3,i,j+1)
641 > - tmp1 * njac(1,3,i,j+1)
642 lhsc(1,4,i,j,k,c) = tmp2 * fjac(1,4,i,j+1)
643 > - tmp1 * njac(1,4,i,j+1)
644 lhsc(1,5,i,j,k,c) = tmp2 * fjac(1,5,i,j+1)
645 > - tmp1 * njac(1,5,i,j+1)
647 lhsc(2,1,i,j,k,c) = tmp2 * fjac(2,1,i,j+1)
648 > - tmp1 * njac(2,1,i,j+1)
649 lhsc(2,2,i,j,k,c) = tmp2 * fjac(2,2,i,j+1)
650 > - tmp1 * njac(2,2,i,j+1)
652 lhsc(2,3,i,j,k,c) = tmp2 * fjac(2,3,i,j+1)
653 > - tmp1 * njac(2,3,i,j+1)
654 lhsc(2,4,i,j,k,c) = tmp2 * fjac(2,4,i,j+1)
655 > - tmp1 * njac(2,4,i,j+1)
656 lhsc(2,5,i,j,k,c) = tmp2 * fjac(2,5,i,j+1)
657 > - tmp1 * njac(2,5,i,j+1)
659 lhsc(3,1,i,j,k,c) = tmp2 * fjac(3,1,i,j+1)
660 > - tmp1 * njac(3,1,i,j+1)
661 lhsc(3,2,i,j,k,c) = tmp2 * fjac(3,2,i,j+1)
662 > - tmp1 * njac(3,2,i,j+1)
663 lhsc(3,3,i,j,k,c) = tmp2 * fjac(3,3,i,j+1)
664 > - tmp1 * njac(3,3,i,j+1)
666 lhsc(3,4,i,j,k,c) = tmp2 * fjac(3,4,i,j+1)
667 > - tmp1 * njac(3,4,i,j+1)
668 lhsc(3,5,i,j,k,c) = tmp2 * fjac(3,5,i,j+1)
669 > - tmp1 * njac(3,5,i,j+1)
671 lhsc(4,1,i,j,k,c) = tmp2 * fjac(4,1,i,j+1)
672 > - tmp1 * njac(4,1,i,j+1)
673 lhsc(4,2,i,j,k,c) = tmp2 * fjac(4,2,i,j+1)
674 > - tmp1 * njac(4,2,i,j+1)
675 lhsc(4,3,i,j,k,c) = tmp2 * fjac(4,3,i,j+1)
676 > - tmp1 * njac(4,3,i,j+1)
677 lhsc(4,4,i,j,k,c) = tmp2 * fjac(4,4,i,j+1)
678 > - tmp1 * njac(4,4,i,j+1)
680 lhsc(4,5,i,j,k,c) = tmp2 * fjac(4,5,i,j+1)
681 > - tmp1 * njac(4,5,i,j+1)
683 lhsc(5,1,i,j,k,c) = tmp2 * fjac(5,1,i,j+1)
684 > - tmp1 * njac(5,1,i,j+1)
685 lhsc(5,2,i,j,k,c) = tmp2 * fjac(5,2,i,j+1)
686 > - tmp1 * njac(5,2,i,j+1)
687 lhsc(5,3,i,j,k,c) = tmp2 * fjac(5,3,i,j+1)
688 > - tmp1 * njac(5,3,i,j+1)
689 lhsc(5,4,i,j,k,c) = tmp2 * fjac(5,4,i,j+1)
690 > - tmp1 * njac(5,4,i,j+1)
691 lhsc(5,5,i,j,k,c) = tmp2 * fjac(5,5,i,j+1)
692 > - tmp1 * njac(5,5,i,j+1)
699 c---------------------------------------------------------------------
700 c outer most do loops - sweeping in i direction
701 c---------------------------------------------------------------------
702 if (first .eq. 1) then
704 c---------------------------------------------------------------------
705 c multiply c(i,jstart,k) by b_inverse and copy back to c
706 c multiply rhs(jstart) by b_inverse(jstart) and copy to rhs
707 c---------------------------------------------------------------------
709 do i=start(1,c),isize
710 call binvcrhs( lhsb(1,1,i,jstart),
711 > lhsc(1,1,i,jstart,k,c),
712 > rhs(1,i,jstart,k,c) )
717 c---------------------------------------------------------------------
718 c begin inner most do loop
719 c do all the elements of the cell unless last
720 c---------------------------------------------------------------------
721 do j=jstart+first,jsize-last
723 do i=start(1,c),isize
725 c---------------------------------------------------------------------
726 c subtract A*lhs_vector(j-1) from lhs_vector(j)
728 c rhs(j) = rhs(j) - A*rhs(j-1)
729 c---------------------------------------------------------------------
730 call matvec_sub(lhsa(1,1,i,j),
731 > rhs(1,i,j-1,k,c),rhs(1,i,j,k,c))
733 c---------------------------------------------------------------------
734 c B(j) = B(j) - C(j-1)*A(j)
735 c---------------------------------------------------------------------
736 call matmul_sub(lhsa(1,1,i,j),
737 > lhsc(1,1,i,j-1,k,c),
740 c---------------------------------------------------------------------
741 c multiply c(i,j,k) by b_inverse and copy back to c
742 c multiply rhs(i,1,k) by b_inverse(i,1,k) and copy to rhs
743 c---------------------------------------------------------------------
744 call binvcrhs( lhsb(1,1,i,j),
751 c---------------------------------------------------------------------
752 c Now finish up special cases for last cell
753 c---------------------------------------------------------------------
754 if (last .eq. 1) then
757 do i=start(1,c),isize
758 c---------------------------------------------------------------------
759 c rhs(jsize) = rhs(jsize) - A*rhs(jsize-1)
760 c---------------------------------------------------------------------
761 call matvec_sub(lhsa(1,1,i,jsize),
762 > rhs(1,i,jsize-1,k,c),rhs(1,i,jsize,k,c))
764 c---------------------------------------------------------------------
765 c B(jsize) = B(jsize) - C(jsize-1)*A(jsize)
766 c call matmul_sub(aa,i,jsize,k,c,
767 c $ cc,i,jsize-1,k,c,bb,i,jsize,k,c)
768 c---------------------------------------------------------------------
769 call matmul_sub(lhsa(1,1,i,jsize),
770 > lhsc(1,1,i,jsize-1,k,c),
773 c---------------------------------------------------------------------
774 c multiply rhs(jsize) by b_inverse(jsize) and copy to rhs
775 c---------------------------------------------------------------------
776 call binvrhs( lhsb(1,1,i,jsize),
777 > rhs(1,i,jsize,k,c) )