File: N:\mfix\model\leq_bicgst.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: LEQ_BICGSt(Vname, Var, A_m, B_m,                       C
4     !                         cmethod, TOL, ITMAX, IER )
5     !  Purpose: Compute residual of linear system                          C
6     !                                                                      C
7     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !  Variables modified:                                                 C
15     !                                                                      C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     !
20           SUBROUTINE LEQ_BICGSt(VNAME, VNO, VAR, A_M, B_m,  cmethod, TOL, PC, ITMAX,IER)
21     
22     !-----------------------------------------------
23     !   M o d u l e s
24     !-----------------------------------------------
25           USE compar
26           USE exit, only: mfix_exit
27           USE funits
28           USE geometry
29           USE indices
30           USE leqsol
31           USE param
32           USE param1
33           IMPLICIT NONE
34     !-----------------------------------------------
35     !   D u m m y   A r g u m e n t s
36     !-----------------------------------------------
37     !
38     !                      Error indicator
39           INTEGER ::          IER
40     !                      maximum number of iterations
41           INTEGER ::          ITMAX
42     !                      variable number
43           INTEGER ::          VNO
44     !                      convergence tolerance
45           DOUBLE PRECISION ::  TOL
46     !                      Preconditioner
47           CHARACTER(LEN=4)   ::  PC
48     !                      Septadiagonal matrix A_m
49           DOUBLE PRECISION, DIMENSION(-3:3,ijkstart3:ijkend3) :: A_m
50     !                      Vector b_m
51           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: B_m
52     !                      Variable name
53           CHARACTER(LEN=*) ::    Vname
54     !                      Variable
55           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: Var
56     !                    sweep direction
57           CHARACTER(LEN=*) :: CMETHOD
58     !
59     !-------------------------------------------------
60           EXTERNAL LEQ_MATVECt, LEQ_MSOLVEt, LEQ_MSOLVE0t, LEQ_MSOLVE1t
61     !--------------------------------------------------
62     
63           if(PC.eq.'LINE') then
64              call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m,                        &
65              cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVEt, IER )
66           elseif(PC.eq.'DIAG') then
67              call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m,                        &
68              cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVE1t, IER )
69           elseif(PC.eq.'NONE') then
70              call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m,                        &
71              cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVE0t, IER )
72           else
73              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'preconditioner option not found - check mfix.dat and readme'
74              call mfix_exit(myPE)
75           endif
76     
77           return
78           END SUBROUTINE LEQ_BICGSt
79     
80     
81     
82     
83     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
84     !                                                                      C
85     !  Module name: LEQ_BICGS0(Vname, Var, A_m, B_m,                       C
86     !                         cmethod, TOL, ITMAX, MATVEC, MSOLVE, IER )   C
87     !  Purpose: Compute residual of linear system                          C
88     !                                                                      C
89     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
90     !  Reviewer:                                          Date:            C
91     !                                                                      C
92     !                                                                      C
93     !  Literature/Document References:                                     C
94     !                                                                      C
95     !  Variables referenced:                                               C
96     !  Variables modified:                                                 C
97     !                                                                      C
98     !  Local variables:                                                    C
99     !                                                                      C
100     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
101     !
102           SUBROUTINE LEQ_BICGS0t(VNAME, VNO, VAR, A_M, B_m,  cmethod, TOL, ITMAX,  &
103                                 MATVECt, MSOLVEt, IER )
104     !-----------------------------------------------
105     !   M o d u l e s
106     !-----------------------------------------------
107           USE param
108           USE param1
109           USE parallel
110           USE geometry
111           USE compar
112           USE mpi_utility
113           USE sendrecv
114           USE indices
115           USE leqsol
116           USE functions
117           IMPLICIT NONE
118     !-----------------------------------------------
119     !   D u m m y   A r g u m e n t s
120     !-----------------------------------------------
121     !
122     !
123     !                      Error indicator
124           INTEGER ::          IER
125     !                      maximum number of iterations
126           INTEGER ::          ITMAX
127     !                      variable number
128           INTEGER ::          VNO
129     !                      convergence tolerance
130           DOUBLE PRECISION ::  TOL
131     !                      Septadiagonal matrix A_m
132           DOUBLE PRECISION, DIMENSION(-3:3,ijkstart3:ijkend3) :: A_m
133     !                      Vector b_m
134           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: B_m
135     !                      Variable name
136           CHARACTER(LEN=*) ::    Vname
137     !                      Variable
138           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: Var
139     !                    sweep direction
140           CHARACTER(LEN=*) :: CMETHOD
141     !
142     !-----------------------------------------------
143     !   L o c a l   V a r i a b l e s
144     !-----------------------------------------------
145           INTEGER, PARAMETER :: idebugl = 1
146           DOUBLE PRECISION :: ratiotol = 0.2
147     
148           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) ::                       &
149                                     R,Rtilde, P,Phat, Svec, Shat, Tvec,V
150           DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: alpha,beta,omega,rho
151           DOUBLE PRECISION :: TxS, TxT, oam,RtildexV,                   &
152                           RtildexR, aijmax, Rnorm=0, Rnorm0, Snorm, TOLMIN
153           LOGICAL :: isconverged
154           INTEGER :: i, j, k, iter
155     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
156           INTEGER :: ijk2
157           DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
158     
159     !-----------------------------------------------
160     !   E x t e r n a l   F u n c t i o n s
161     !-----------------------------------------------
162           EXTERNAL  MATVECt, MSOLVEt
163     
164           logical, parameter :: do_unit_scaling = .true.
165     
166           is_serial = numPEs.eq.1.and.is_serial
167     
168           alpha(:)  = zero
169           beta(:)   = zero
170           omega(:)  = zero
171           rho(:)    = zero
172     
173     !
174     !     ---------------------------------------------
175     !     zero out R,Rtilde, P,Phat, Svec, Shat, Tvec,V
176     !     ---------------------------------------------
177           if (use_doloop) then
178     
179     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
180     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 201)
181     !!$omp  parallel do private(ijk2)
182              do ijk2=ijkstart3,ijkend3
183                 R(ijk2) = zero
184                 Rtilde(ijk2) = zero
185                 P(ijk2) = zero
186                 Phat(ijk2) = zero
187                 Svec(ijk2) = zero
188                 Shat(ijk2) = zero
189                 Tvec(ijk2) = zero
190                 V(ijk2) = zero
191              enddo
192     
193           else
194     
195              R(:) = zero
196              Rtilde(:) = zero
197              P(:) = zero
198              Phat(:) = zero
199              Svec(:) = zero
200              Shat(:) = zero
201              Tvec(:) = zero
202              V(:) = zero
203     
204           endif
205     
206     
207           TOLMIN = EPSILON( one )
208     
209           if (do_unit_scaling) then
210     !
211     !     Scale matrix to have unit diagonal
212     !
213     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
214     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 233)
215     !!$omp parallel do private(ijk2,i,j,k,oam,aijmax) collapse (3)
216              do k = kstart2,kend2
217                 do i = istart2,iend2
218                    do j = jstart2,jend2
219                       IJK2 = funijk(i,j,k)
220     
221                       aijmax = maxval(abs(A_M(:,ijk2)) )
222     
223                       OAM = one/aijmax
224     
225                       A_M(:,IJK2) = A_M(:,IJK2)*OAM
226     
227                       B_M(IJK2) = B_M(IJK2)*OAM
228     
229                    enddo
230                 enddo
231              enddo
232           endif
233     
234     !
235     !    Compute initial residual, assume initial guess in Var
236     !    r = b - A*x
237     !    rtilde = r
238     
239     
240           call MATVECt( Vname, Var, A_M, R )
241     
242     
243           if (use_doloop) then
244     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
245     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 264)
246     !!!$omp   parallel do private(ijk2)
247              do ijk2=ijkstart3,ijkend3
248                 R(ijk2) = B_m(ijk2) - R(ijk2)
249              enddo
250           else
251              R(:) = B_m(:) - R(:)
252           endif
253     
254           if(is_serial) then
255              Rnorm0 = zero
256              if (use_doloop) then
257     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
258     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 276)
259     !!!$omp          parallel do private(ijk2) reduction(+:Rnorm0)
260                 do ijk2=ijkstart3,ijkend3
261                    Rnorm0 = Rnorm0 + R(ijk2)*R(ijk2)
262                 enddo
263              else
264                 Rnorm0 = dot_product(R,R)
265              endif
266              Rnorm0 = sqrt( Rnorm0 )
267           else
268              Rnorm0 = sqrt( dot_product_par( R, R ) )
269           endif
270     
271           call random_number(Rtilde(:))
272     
273           if (use_doloop) then
274     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
275     !!!$omp   parallel do private(ijk2)
276              do ijk2=ijkstart3,ijkend3
277                 Rtilde(ijk2) = R(ijk2) + (2.0d0*Rtilde(ijk2)-1.0d0)*1.0d-6*Rnorm0
278              enddo
279           else
280              Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
281           endif
282     
283           if (idebugl >= 1) then
284              if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
285           endif
286     !
287     !     Main loop
288     !
289           iter = 1
290           do i=1,itmax
291     
292              if(is_serial) then
293                 if (use_doloop) then
294                    RtildexR = zero
295     !AIKE PFUPGRADE 091409
296     !!!$omp        parallel do private(ijk2) reduction(+:RtildexR)
297                    do ijk2=ijkstart3,ijkend3
298                       RtildexR = RtildexR + Rtilde(ijk2) * R(ijk2)
299                    enddo
300                    rho(i-1) = RtildexR
301                 else
302                    rho(i-1) = dot_product( Rtilde, R )
303                 endif
304              else
305                 rho(i-1) = dot_product_par( Rtilde, R )
306              endif ! is_serial
307     
308     !     print*,'leq_bicgs, initial: ', Vname,' rho(i-1) ', rho(i-1)
309     
310              if (rho(i-1) .eq. zero) then
311                 if(i /= 1)then
312     !     ------------
313     !     Method fails
314     !     ------------
315     !     print*, 'leq_bicgs,',Vname,': rho(i-1) == 0 '
316                    ier = -2
317                 else
318     !     ------------
319     !     converged.  residual is already zero
320     !     ------------
321                    ier = 0
322                 endif
323                 call send_recv(var,2)
324                 return
325              endif ! rho(i-1).eq.0
326     
327              if (i .eq. 1) then
328                 if (use_doloop) then
329     !AIKE PFUPGRADE 091409
330     !!!$omp        parallel do private(ijk2)
331                    do ijk2=ijkstart3,ijkend3
332                       P(ijk2) = R(ijk2)
333                    enddo
334                 else
335                    P(:) = R(:)
336                 endif
337              else
338                 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
339                 if (use_doloop) then
340     !AIKE PFUPGRADE 091409
341     !!!$omp        parallel do private(ijk2)
342                    do ijk2=ijkstart3,ijkend3
343                       P(ijk2) = R(ijk2) + beta(i-1)*( P(ijk2) - omega(i-1)*V(ijk2) )
344                    enddo
345                 else
346                    P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
347                 endif
348              endif ! i.eq.1
349     
350     !
351     !     Solve M Phat(:) = P(:)
352     !     V(:) = A*Phat(:)
353     !
354     
355              call MSOLVEt( Vname, P, A_m, Phat, CMETHOD)
356     
357              call MATVECt( Vname, Phat, A_m, V )
358     
359              if(is_serial) then
360                 if (use_doloop) then
361                    RtildexV = zero
362     !AIKE PFUPGRADE 091409
363     !!!$omp         parallel do private(ijk2) reduction(+:RtildexV)
364                    do ijk2=ijkstart3,ijkend3
365                       RtildexV = RtildexV + Rtilde(ijk2) * V(ijk2)
366                    enddo
367                 else
368                    RtildexV = dot_product( Rtilde, V )
369                 endif
370              else
371                 RtildexV = dot_product_par( Rtilde, V )
372              endif ! is_serial
373     
374     !     print*,'leq_bicgs, initial: ', Vname,' RtildexV ', RtildexV
375     
376              alpha(i) = rho(i-1) / RtildexV
377     
378              if (use_doloop) then
379     !AIKE PFUPGRADE 091409
380     !!!$omp     parallel do private(ijk2)
381                 do ijk2=ijkstart3,ijkend3
382                    Svec(ijk2) = R(ijk2) - alpha(i) * V(ijk2)
383                 enddo
384              else
385                 Svec(:) = R(:) - alpha(i) * V(:)
386              endif ! use_doloop
387     
388              if(.not.minimize_dotproducts) then
389     !
390     !     Check norm of Svec(:); if small enough:
391     !     set X(:) = X(:) + alpha(i)*Phat(:) and stop
392     !
393                 if(is_serial) then
394                    if (use_doloop) then
395                       Snorm = zero
396     !AIKE PFUPGRADE 091409
397     !!!$omp       parallel do private(ijk2) reduction(+:Snorm)
398                       do ijk2=ijkstart3,ijkend3
399                          Snorm = Snorm + Svec(ijk2) * Svec(ijk2)
400                       enddo
401                    else
402                       Snorm = dot_product( Svec, Svec )
403                    endif
404                    Snorm = sqrt( Snorm )
405                 else
406                    Snorm = sqrt( dot_product_par( Svec, Svec ) )
407                 endif               ! is_serial
408     !     print*,'leq_bicgs, initial: ', Vname,' Snorm ', real(Snorm)
409     
410     
411                 if (Snorm <= TOLMIN) then
412                    if (use_doloop) then
413     !AIKE PFUPGRADE 091409
414     !!!$omp          parallel do private(ijk2)
415                       do ijk2=ijkstart3,ijkend3
416                          Var(ijk2) = Var(ijk2) + alpha(i)*Phat(ijk2)
417                       enddo
418                    else
419                       Var(:) = Var(:) + alpha(i)*Phat(:)
420                    endif            ! use_doloop
421     
422                    if (idebugl >= 1) then
423     !
424     !     Recompute residual norm
425     !
426                       call MATVECt( Vname, Var, A_m, R )
427     
428     !     Rnorm = sqrt( dot_product_par( Var, Var ) )
429     !     print*,'leq_bicgs, initial: ', Vname,' Vnorm ', Rnorm
430     
431                       if (use_doloop) then
432     !AIKE PFUPGRADE 091409
433     !!!$omp          parallel do private(ijk2)
434                          do ijk2=ijkstart3,ijkend3
435                             R(ijk2) = B_m(ijk2) - R(ijk2)
436                          enddo
437                       else
438                          R(:) = B_m(:) - R(:)
439                       endif
440     
441                       if(is_serial) then
442                          if (use_doloop) then
443                             Rnorm = zero
444     !AIKE PFUPGRADE 091409
445     !!!$omp            parallel do private(ijk2) reduction(+:Rnorm)
446                             do ijk2=ijkstart3,ijkend3
447                                Rnorm = Rnorm + R(ijk2)*R(ijk2)
448                             enddo
449                          else
450                             Rnorm =  dot_product( R, R )
451                          endif
452                          Rnorm = sqrt( Rnorm )
453                       else
454                          Rnorm = sqrt( dot_product_par( R, R ) )
455                       endif
456     !     print*,'leq_bicgs, initial: ', Vname,' Rnorm ', Rnorm
457                    endif            ! idebugl >= 1
458     
459                    EXIT
460                 endif               ! Snorm <= TOLMIN
461     
462              endif                  ! .not.minimize_dotproducts
463     
464     !
465     !     Solve M Shat(:) = Svec(:)
466     !     Tvec(:) = A * Shat(:)
467     !
468              call MSOLVEt( Vname, Svec, A_m, Shat, CMETHOD)
469     
470              call MATVECt( Vname, Shat, A_m, Tvec )
471     
472              if(is_serial) then
473                 if (use_doloop) then
474                    TxS = zero
475                    TxT = zero
476     !AIKE PFUPGRADE 091409
477     !!!$omp  parallel do private(ijk2) reduction(+:TxS,TxT)
478                    do ijk2=ijkstart3,ijkend3
479                       TxS = TxS + Tvec(ijk2)  * Svec(ijk2)
480                       TxT = TxT + Tvec(ijk2)  * Tvec(ijk2)
481                    enddo
482                 else
483                    TxS = dot_product( Tvec, Svec )
484                    TxT = dot_product( Tvec, Tvec )
485                 endif
486              else
487                 if(.not.minimize_dotproducts) then
488                    TxS = dot_product_par( Tvec, Svec )
489                    TxT = dot_product_par( Tvec, Tvec )
490                 else
491                    TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
492                    TxS = TxS_TxT(1)
493                    TxT = TxS_TxT(2)
494                 endif
495              endif
496              IF(TxT.eq.Zero) TxT = SMALL_NUMBER
497              omega(i) = TxS / TxT
498     
499     
500              if (use_doloop) then
501     !AIKE PFUPGRADE 091409
502     !!!$omp    parallel do private(ijk2)
503                 do ijk2=ijkstart3,ijkend3
504                    Var(ijk2) = Var(ijk2) +                           &
505                    alpha(i)*Phat(ijk2) + omega(i)*Shat(ijk2)
506                    R(ijk2) = Svec(ijk2) - omega(i)*Tvec(ijk2)
507                 enddo
508              else
509                 Var(:) = Var(:) +                           &
510                 alpha(i)*Phat(:) + omega(i)*Shat(:)
511                 R(:) = Svec(:) - omega(i)*Tvec(:)
512              endif
513     
514              if(.not.minimize_dotproducts.or.(mod(iter,5).eq.0)) then
515                 if(is_serial) then
516                    if (use_doloop) then
517                       Rnorm = zero
518     !AIKE PFUPGRADE 091409
519     !!!$omp       parallel do private(ijk2) reduction(+:Rnorm)
520                       do ijk2=ijkstart3,ijkend3
521                          Rnorm = Rnorm + R(ijk2) * R(ijk2)
522                       enddo
523                    else
524                       Rnorm =  dot_product(R, R )
525                    endif
526                    Rnorm = sqrt( Rnorm )
527                 else
528                    Rnorm = sqrt( dot_product_par(R, R) )
529                 endif               ! is_serial
530     
531                 if (idebugl.ge.1) then
532                    if (myPE.eq.PE_IO) then
533                       print*,'iter, Rnorm ', iter, Rnorm, Snorm
534                       print*,'alpha(i), omega(i) ', alpha(i), omega(i)
535                       print*,'TxS, TxT ', TxS, TxT
536                       print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
537                    endif
538                 endif
539     
540     !     call mfix_exit(myPE)
541     
542     !     Check convergence; continue if necessary
543     !     for continuation, it is necessary that omega(i) .ne. 0
544     !
545                 isconverged = (Rnorm <= TOL*Rnorm0)
546     
547                 if (isconverged) then
548                    iter_tot(vno) = iter_tot(vno) + iter + 1
549                    EXIT
550                 endif
551     
552              endif                  ! .not.minimize_dotproducts
553     
554     !     Advance the iteration count
555              iter = iter + 1
556     
557           enddo
558     
559           if (idebugl >= 1) then
560              call MATVECt( Vname, Var, A_m, R )
561              if (use_doloop) then
562     !AIKE PFUPGRADE 091409
563     !!!$omp  parallel do private(ijk2)
564                 do ijk2=ijkstart3,ijkend3
565                    R(ijk2) = R(ijk2) - B_m(ijk2)
566                 enddo
567              else
568                 R(:) = R(:) - B_m(:)
569              endif
570     
571              if(is_serial) then
572                 if (use_doloop) then
573                    Rnorm = zero
574     !AIKE PFUPGRADE 091409
575     !!!$omp         parallel do private(ijk2) reduction(+:Rnorm)
576                    do ijk2=ijkstart3,ijkend3
577                       Rnorm = Rnorm + R(ijk2) * R(ijk2)
578                    enddo
579                 else
580                    Rnorm = dot_product( R,R)
581                 endif
582                 Rnorm = sqrt( Rnorm )
583              else
584                 Rnorm = sqrt( dot_product_par( R,R) )
585              endif
586     
587              if(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
588     
589              if(myPE.eq.0)  print*,'leq_bicgs ratio : ', Vname,' ',iter,     &
590              ' L-2', Rnorm/Rnorm0
591           endif
592     
593           isconverged = (real(Rnorm) <= TOL*Rnorm0);
594     !     write(*,*) '***',iter, isconverged, Rnorm, TOL, Rnorm0, myPE
595           IER = 0
596           if (.not.isconverged) then
597              iter_tot(vno) = iter_tot(vno) + iter
598              IER = -1
599              if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
600                 IER = -2
601              endif
602           endif
603     
604           call send_recv(var,2)
605     
606           return
607           end subroutine LEQ_BICGS0t
608     
609     
610     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
611     !                                                                      C
612     !  Module name: LEQ_ISWEEPt(I, Vname, Var, A_m, B_m )                  C
613     !  Purpose: Perform line sweep at coordiante I                         C
614     !                                                                      C
615     !                                                                      C
616     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
617     !  Reviewer:                                          Date:            C
618     !                                                                      C
619     !                                                                      C
620     !  Literature/Document References:                                     C
621     !                                                                      C
622     !  Variables referenced:                                               C
623     !  Variables modified:                                                 C
624     !                                                                      C
625     !  Local variables:                                                    C
626     !                                                                      C
627     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
628           SUBROUTINE LEQ_ISWEEPt(I,Vname, VAR, A_M, B_M)
629     
630     !-----------------------------------------------
631     !   M o d u l e s
632     !-----------------------------------------------
633           USE param
634           USE param1
635           USE geometry
636           USE compar
637           USE indices
638           USE funits
639           USE sendrecv
640           USE mpi_utility
641           USE functions
642           IMPLICIT NONE
643     !-----------------------------------------------
644     !   G l o b a l   P a r a m e t e r s
645     !-----------------------------------------------
646     !-----------------------------------------------
647     !   D u m m y   A r g u m e n t s
648     !-----------------------------------------------
649     !                      Line position
650           INTEGER          I
651     !
652     !
653     !                      Septadiagonal matrix A_m
654           DOUBLE PRECISION A_m(-3:3,ijkstart3:ijkend3)
655     !
656     !                      Vector b_m
657           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
658     !
659     !                      Variable name
660           CHARACTER(LEN=*)    Vname
661     
662     !
663     !                      Variable
664           DOUBLE PRECISION Var(ijkstart3:ijkend3)
665     
666     !-----------------------------------------------
667     !   L o c a l   P a r a m e t e r s
668     !-----------------------------------------------
669     !
670     
671     
672     !-----------------------------------------------
673     !   L o c a l   V a r i a b l e s
674     !-----------------------------------------------
675     
676           INTEGER :: NSTART, NEND
677           DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC,DD,EE,BB
678           INTEGER :: INFO, IJK, J, K, IM1JK, IP1JK
679     
680           NEND = JEND
681           NSTART = JSTART
682           K = 1
683     
684           DO J=NSTART, NEND
685     !     IJK = FUNIJK(IMAP_C(I),JMAP_C(J),KMAP_C(K))
686              IJK = FUNIJK(I,J,K)
687              IM1JK = IM_OF(IJK)
688              IP1JK = IP_OF(IJK)
689     
690              DD(J) = A_M(0, IJK)
691              CC(J) = A_M(-2, IJK)
692              EE(J) = A_M(2, IJK)
693              BB(J) = B_M(IJK) -  A_M(-1, IJK) * Var( IM1JK )         &
694              -  A_M(1, IJK) * Var( IP1JK )
695     
696           END DO
697     
698           CC(NSTART) = ZERO
699           EE(NEND) = ZERO
700           INFO = 0
701     
702     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
703           CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB,  JEND-JSTART+1, INFO )
704     
705           IF (INFO.NE.0) THEN
706              RETURN
707           ENDIF
708     
709           DO J=NSTART, NEND
710              IJK = FUNIJK(I,J,K)
711              Var(IJK) =  BB(J)
712           ENDDO
713     
714           RETURN
715           END SUBROUTINE  LEQ_ISWEEPt
716     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
717     !                                                                      C
718     !  Module name: LEQ_IKSWEEPt(I, K, Vname, Var, A_m, B_m )              C
719     !  Purpose: Perform line sweep at coordiante I,K                       C
720     !                                                                      C
721     !                                                                      C
722     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
723     !  Reviewer:                                          Date:            C
724     !                                                                      C
725     !                                                                      C
726     !  Literature/Document References:                                     C
727     !                                                                      C
728     !  Variables referenced:                                               C
729     !  Variables modified:                                                 C
730     !                                                                      C
731     !  Local variables:                                                    C
732     !                                                                      C
733     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
734           SUBROUTINE LEQ_IKSWEEPt(I,K,Vname, VAR, A_M, B_M )
735     
736     !-----------------------------------------------
737     !   M o d u l e s
738     !-----------------------------------------------
739           USE param
740           USE param1
741           USE geometry
742           USE compar
743           USE funits
744           USE indices
745           USE sendrecv
746           USE mpi_utility
747           USE functions
748           IMPLICIT NONE
749     !-----------------------------------------------
750     !   G l o b a l   P a r a m e t e r s
751     !-----------------------------------------------
752     !-----------------------------------------------
753     !   D u m m y   A r g u m e n t s
754     !-----------------------------------------------
755     !                      Line position
756           INTEGER          I,K
757     !
758     !                      Septadiagonal matrix A_m
759           DOUBLE PRECISION A_m(-3:3,ijkstart3:ijkend3)
760     !
761     !                      Vector b_m
762           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
763     !
764     !                      Variable name
765           CHARACTER(LEN=*)    Vname
766     
767     !
768     !                      Variable
769           DOUBLE PRECISION Var(ijkstart3:ijkend3)
770     !
771     
772     !-----------------------------------------------
773     !   L o c a l   P a r a m e t e r s
774     !-----------------------------------------------
775     !
776     
777     
778     !-----------------------------------------------
779     !     L o c a l   V a r i a b l e s
780     !-----------------------------------------------
781     
782           DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC,DD,EE, BB
783           INTEGER :: NSTART, NEND, INFO, IJK, J,  IM1JK, IP1JK, IJKM1, IJKP1
784     
785           NEND = JEND
786           NSTART = JSTART
787     
788     !!!$omp parallel do private(j,ijk,im1jk,ip1jk,ijkm1,ijkp1)
789           DO J=NSTART, NEND
790     
791     !     IJK = FUNIJK(IMAP_C(I),JMAP_C(J),KMAP_C(K))
792              IJK = FUNIJK(I,J,K)
793              IM1JK = IM_OF(IJK)
794              IP1JK = IP_OF(IJK)
795              IJKM1 = KM_OF(IJK)
796              IJKP1 = KP_OF(IJK)
797     
798              DD(J) = A_M(0, IJK)
799              CC(J) = A_M(-2, IJK)
800              EE(J) = A_M(2, IJK)
801              BB(J) = B_M(IJK) -  A_M(-1, IJK) * Var( IM1JK )         &
802              -  A_M(1, IJK) * Var( IP1JK )         &
803              -  A_M(-3, IJK) * Var( IJKM1 )         &
804              -  A_M(3, IJK) * Var( IJKP1 )
805     
806           ENDDO
807     
808           CC(NSTART) = ZERO
809           EE(NEND) = ZERO
810           INFO = 0
811     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
812           CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB,  JEND-JSTART+1, INFO )
813     
814           IF (INFO.NE.0) THEN
815              write(*,*) 'leq_iksweep',INFO, myPE
816              RETURN
817           ENDIF
818     
819           DO J=NSTART, NEND
820     
821              IJK = FUNIJK(I,J,K)
822              Var(IJK) = BB(J)
823     
824           ENDDO
825     
826           RETURN
827           END SUBROUTINE  LEQ_IKSWEEPt
828     
829     
830     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
831     !                                                                      C
832     !  Module name: LEQ_MATVECt(Vname, Var, A_m, B_m )                     C
833     !  Purpose: Compute residual of linear system                          C
834     !                                                                      C
835     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
836     !  Reviewer:                                          Date:            C
837     !                                                                      C
838     !                                                                      C
839     !  Literature/Document References:                                     C
840     !                                                                      C
841     !  Variables referenced:                                               C
842     !  Variables modified:                                                 C
843     !                                                                      C
844     !  Local variables:                                                    C
845     !                                                                      C
846     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
847     !
848           SUBROUTINE LEQ_MATVECt(VNAME, VAR, A_M, Avar )
849     !
850     !-----------------------------------------------
851     !   M o d u l e s
852     !-----------------------------------------------
853           USE param
854           USE param1
855           USE geometry
856           USE compar
857           USE indices
858           USE sendrecv
859           USE mpi_utility
860           USE functions
861           IMPLICIT NONE
862     !-----------------------------------------------
863     !   G l o b a l   P a r a m e t e r s
864     !-----------------------------------------------
865     !-----------------------------------------------
866     !   D u m m y   A r g u m e n t s
867     !-----------------------------------------------
868     !
869     !
870     !                      Septadiagonal matrix A_m
871           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
872     !
873     !                      Vector AVar
874           DOUBLE PRECISION AVar(ijkstart3:ijkend3)
875     !
876     !                      Variable name
877           CHARACTER(LEN=*)    Vname
878     !
879     !                      Variable
880           DOUBLE PRECISION Var(ijkstart3:ijkend3)
881     !                      Variable
882     !
883     !-----------------------------------------------
884     !   L o c a l   P a r a m e t e r s
885     !-----------------------------------------------
886     !
887     !-----------------------------------------------
888     !   L o c a l   V a r i a b l e s
889     !-----------------------------------------------
890     !
891           INTEGER          I,  J, K, IJK2
892     
893           integer :: im1jk,ip1jk, ijm1k,ijp1k, ijkm1,ijkp1
894     
895           if (do_k) then
896     
897     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
898     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 938)
899     !!$omp    parallel  do &
900     !!$omp&   private(     &
901     !!$omp&           ijk2,i,j,k, &
902     !!$omp&           im1jk,ip1jk,ijm1k,ijp1k,ijkm1,ijkp1) collapse (3)
903              do k = kstart,kend
904                 do i = istart,iend
905                    do j = jstart,jend
906     
907                       IJK2 = funijk(i,j,k)
908     
909                       im1jk = im_of(ijk2)
910                       ip1jk = ip_of(ijk2)
911                       ijm1k = jm_of(ijk2)
912                       ijp1k = jp_of(ijk2)
913     !
914                       ijkm1 = km_of(ijk2)
915                       ijkp1 = kp_of(ijk2)
916     
917     
918                       AVar(ijk2) =      A_m(-3, ijk2) * Var(ijkm1)   &
919                       + A_m(-2, ijk2) * Var(ijm1k)   &
920                       + A_m(-1, ijk2) * Var(im1jk)   &
921                       + A_m(0, ijk2) * Var(ijk2)     &
922                       + A_m(1, ijk2) * Var(ip1jk)   &
923                       + A_m(2, ijk2) * Var(ijp1k)   &
924                       + A_m(3, ijk2) * Var(ijkp1)
925     
926                    enddo
927                 enddo
928              enddo
929     
930           else
931              k = 1
932     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
933     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 971)
934     !!$omp parallel do private(i,j,ijk2,im1jk,ip1jk,ijm1k,ijp1k) collapse (2)
935              do i = istart,iend
936                 do j = jstart,jend
937     
938     
939                    IJK2 = funijk(i,j,k)
940     
941     
942                    im1jk = im_of(ijk2)
943                    ip1jk = ip_of(ijk2)
944                    ijm1k = jm_of(ijk2)
945                    ijp1k = jp_of(ijk2)
946                    AVar(ijk2) =      A_m(-2, ijk2) * Var(ijm1k)   &
947                    + A_m(-1, ijk2) * Var(im1jk)   &
948                    + A_m(0, ijk2) * Var(ijk2)     &
949                    + A_m(1, ijk2) * Var(ip1jk)   &
950                    + A_m(2, ijk2) * Var(ijp1k)
951     
952                 enddo
953              enddo
954     
955           endif
956     
957           call send_recv(Avar,nlayers_bicgs)
958           return
959           END SUBROUTINE LEQ_MATVECt
960     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
961     !                                                                      C
962     !  Module name: LEQ_MSOLVE(Vname, B_m, A_m, Var, CMETHOD)
963     !  Purpose: Successive line over-relaxation method -- Cyclic bc        C
964     !                                                                      C
965     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
966     !  Reviewer:                                          Date:            C
967     !                                                                      C
968     !                                                                      C
969     !  Literature/Document References:                                     C
970     !                                                                      C
971     !  Variables referenced:                                               C
972     !  Variables modified:                                                 C
973     !                                                                      C
974     !  Local variables:                                                    C
975     !                                                                      C
976     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
977     !
978           SUBROUTINE LEQ_MSOLVEt(VNAME, B_m, A_M, Var, CMETHOD)
979     !
980     !-----------------------------------------------
981     !   M o d u l e s
982     !-----------------------------------------------
983           USE param
984           USE param1
985           USE geometry
986           USE compar
987           USE indices
988           USE sendrecv
989           USE functions
990           IMPLICIT NONE
991     !-----------------------------------------------
992     !   G l o b a l   P a r a m e t e r s
993     !-----------------------------------------------
994     !-----------------------------------------------
995     !   D u m m y   A r g u m e n t s
996     !-----------------------------------------------
997     !
998     !                      Septadiagonal matrix A_m
999           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1000     !
1001     !                      Vector b_m
1002           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1003     !
1004     !                      Variable name
1005           CHARACTER(LEN=*)    Vname
1006     !
1007     !                      Variable
1008           DOUBLE PRECISION Var(ijkstart3:ijkend3)
1009     
1010     
1011     !-----------------------------------------------
1012     !   L o c a l   P a r a m e t e r s
1013     !-----------------------------------------------
1014     !
1015     !-----------------------------------------------
1016     !   L o c a l   V a r i a b l e s
1017     !-----------------------------------------------
1018     
1019     !
1020           INTEGER ::   I, J, K, ITER, NITER, IJK2
1021           INTEGER ::   I1 , K1 , I2, K2, IK, ISIZE, KSIZE
1022           INTEGER ::   ICASE
1023     
1024     !     CHARACTER(LEN=4), PARAMETER :: CMETHOD = 'II'
1025     !     sweep direction
1026           CHARACTER(LEN=4) :: CMETHOD
1027           CHARACTER :: CH
1028           LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
1029           LOGICAL :: DO_SENDRECV, DO_REDBLACK
1030           LOGICAL, PARAMETER :: USE_IKLOOP = .FALSE.
1031     
1032           LOGICAL, PARAMETER :: SETGUESS = .TRUE.
1033     
1034     !!$      double precision omp_start, omp_end
1035     !!$      double precision omp_get_wtime
1036     !       by Tingwen
1037     !!$      omp_start=omp_get_wtime()
1038     
1039           IF (SETGUESS) THEN
1040     
1041     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
1042     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 1077)
1043     !!$omp   parallel do private(i,j,k,ijk2) collapse (3)
1044              do k = kstart3,kend3
1045                 do i = istart3,iend3
1046                    do j = jstart3,jend3
1047     
1048                       IJK2 = funijk(i,j,k)
1049     
1050                       VAR(IJK2) = B_M(IJK2)
1051     
1052                    enddo
1053                 enddo
1054              enddo
1055     
1056           call send_recv(var,nlayers_bicgs)
1057     
1058           ENDIF
1059     
1060           NITER = LEN( CMETHOD )
1061     
1062           DO ITER=1,NITER
1063     !
1064     !     Perform sweeps
1065     !
1066              CH = CMETHOD( ITER:ITER )
1067              DO_ISWEEP = (CH .EQ. 'I') .OR. (CH .EQ. 'i')
1068              DO_JSWEEP = (CH .EQ. 'J') .OR. (CH .EQ. 'j')
1069              DO_KSWEEP = (CH .EQ. 'K') .OR. (CH .EQ. 'k')
1070              DO_SENDRECV = (CH .EQ. 'S') .OR. (CH .EQ. 's')
1071              DO_REDBLACK = (CH .EQ. 'R') .OR. (CH .EQ. 'r')
1072     
1073              IF (NO_K) THEN
1074     ! 2D run no need to enable openmp parallel
1075                 IF ( DO_ISWEEP ) THEN
1076     !!!$omp   parallel do private(I)
1077                    DO I=istart,iend,1
1078                       CALL LEQ_ISWEEPt( I, Vname, Var, A_m, B_m )
1079                    ENDDO
1080                 ENDIF
1081     
1082              ELSE
1083     
1084                 IF(DO_REDBLACK) THEN
1085     
1086                    i1 = istart
1087                    k1 = kstart
1088                    i2 = iend
1089                    k2 = kend
1090                    isize = i2-i1+1
1091                    ksize = k2-k1+1
1092     
1093                    DO icase = 1, 2
1094     !!$omp   parallel do private(K,I,IK)
1095                       DO IK=icase, ksize*isize, 2
1096                          if (mod(ik,isize).ne.0) then
1097                             k = int( ik/isize ) + k1
1098                          else
1099                             k = int( ik/isize ) + k1 -1
1100                          endif
1101                          i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
1102                          if(i.gt.i2) i=i-i2 + i1 -1
1103                          CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1104                       ENDDO
1105                    ENDDO
1106     
1107                 ENDIF
1108     !  Not sure the purpose of us_ikloop
1109     !  The SMP directives below need review                        !Tingwen Jan 2012
1110                 IF(USE_IKLOOP) THEN
1111     
1112                    i1 = istart
1113                    k1 = kstart
1114                    i2 = iend
1115                    k2 = kend
1116                    isize = i2-i1+1
1117                    ksize = k2-k1+1
1118     
1119                    IF (DO_ISWEEP) THEN
1120     !!!$omp   parallel do private(K,I,IK)
1121                       DO IK=1, ksize*isize
1122                          if (mod(ik,isize).ne.0) then
1123                             k = int( ik/isize ) + k1
1124                          else
1125                             k = int( ik/isize ) + k1 -1
1126                          endif
1127                          i = (ik-1-(k-k1)*isize) + i1
1128                          CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1129                       ENDDO
1130                    ENDIF
1131     
1132                    IF (DO_KSWEEP) THEN
1133     !!!$omp   parallel do private(K,I,IK)
1134                       DO IK=1, ksize*isize
1135                          if (mod(ik,ksize).ne.0) then
1136                             i = int( ik/ksize ) + i1
1137                          else
1138                             i = int( ik/ksize ) + i1 -1
1139                          endif
1140                          k = (ik-1-(i-i1)*ksize) + k1
1141     
1142                          CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1143                       ENDDO
1144                    ENDIF
1145     
1146                 ELSE
1147     !  Not sure the purpose of us_ikloop
1148     !  The SMP directives below need review                        !Tingwen Jan 2012
1149                    IF (DO_ISWEEP) THEN
1150     !!!$omp   parallel do private(K,I)
1151                       DO K=kstart,kend
1152                          DO I=istart,iend
1153                             CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1154                          ENDDO
1155                       ENDDO
1156                    ENDIF
1157     
1158                    IF (DO_KSWEEP) THEN
1159     !!!$omp   parallel do private(K,I)
1160                       DO I=istart,iend
1161                          DO K=kstart,kend
1162                             CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1163                          ENDDO
1164                       ENDDO
1165                    ENDIF
1166     
1167                 ENDIF
1168     
1169                 IF (DO_SENDRECV) call send_recv(var,nlayers_bicgs)
1170     
1171              ENDIF
1172     
1173           ENDDO
1174     !!$      omp_end=omp_get_wtime()
1175     !!$      write(*,*)'leq_msolvet:',omp_end - omp_start
1176           RETURN
1177           END SUBROUTINE LEQ_MSOLVEt
1178     
1179     
1180     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1181     !                                                                      C
1182     !  Module name: LEQ_JKSWEEPt(J, K, Vname, Var, A_m, B_m )              C
1183     !  Purpose: Perform line sweep at coordiante I,K                       C
1184     !                                                                      C
1185     !                                                                      C
1186     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1187     !  Reviewer:                                          Date:            C
1188     !                                                                      C
1189     !                                                                      C
1190     !  Literature/Document References:                                     C
1191     !                                                                      C
1192     !  Variables referenced:                                               C
1193     !  Variables modified:                                                 C
1194     !                                                                      C
1195     !  Local variables:                                                    C
1196     !                                                                      C
1197     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1198           SUBROUTINE LEQ_JKSWEEPt(J,K,Vname, VAR, A_M, B_M )
1199     
1200     !-----------------------------------------------
1201     !   M o d u l e s
1202     !-----------------------------------------------
1203           USE compar
1204           USE exit, only: mfix_exit
1205           USE functions
1206           USE funits
1207           USE geometry
1208           USE indices
1209           USE param
1210           USE param1
1211           IMPLICIT NONE
1212     !-----------------------------------------------
1213     !   G l o b a l   P a r a m e t e r s
1214     !-----------------------------------------------
1215     !-----------------------------------------------
1216     !   D u m m y   A r g u m e n t s
1217     !-----------------------------------------------
1218     !                      Line position
1219           INTEGER          J,K
1220     !
1221     !                      Septadiagonal matrix A_m
1222           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1223     !
1224     !                      Vector b_m
1225           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1226     !
1227     !                      Variable name
1228           CHARACTER(LEN=*)    Vname
1229     !
1230     !                      Variable
1231           DOUBLE PRECISION Var(ijkstart3:ijkend3)
1232     
1233     !-----------------------------------------------
1234     !   L o c a l   P a r a m e t e r s
1235     !-----------------------------------------------
1236     !
1237     
1238     !-----------------------------------------------
1239     !   L o c a l   V a r i a b l e s
1240     !-----------------------------------------------
1241     
1242           DOUBLE PRECISION, DIMENSION (IMAX2) :: CC,DD,EE,BB
1243           INTEGER :: NN, INFO, IJK, I
1244     
1245           NN = IMAX2
1246     
1247           DO I=1,NN
1248              IJK = FUNIJK(I,J,K)
1249     
1250              DD(I) = A_M(0, IJK)
1251              CC(I) = A_M(-1, IJK)
1252              EE(I) = A_M(1, IJK)
1253              BB(I) = B_M(IJK)    -  A_M(-2,IJK) * Var( JM_OF(IJK) )         &
1254              -  A_M(2, IJK) * Var( JP_OF(IJK) )         &
1255              -  A_M(-3, IJK) * Var( KM_OF(IJK) )         &
1256              -  A_M(3, IJK) * Var( KP_OF(IJK) )
1257     
1258           ENDDO
1259     
1260           CC(1) = ZERO
1261           EE(NN) = ZERO
1262     !     DL(1:NEND-1) = CC(2:NEND)
1263           INFO = 0
1264           CALL DGTSL( NN, CC, DD, EE, BB, INFO )
1265     !     CALL DGTSV( JEND-JSTART+1, 1, DL, DD, EE, BB,  JEND-JSTART+1, INFO )
1266     
1267           IF (INFO.NE.0) THEN
1268              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1269              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1270              call mfix_exit(myPE)
1271           ENDIF
1272     
1273           DO I=1,NN
1274              IJK = FUNIJK(I,J,K)
1275              Var(IJK) = BB(I)
1276           ENDDO
1277     
1278           RETURN
1279           END SUBROUTINE  LEQ_JKSWEEPt
1280     
1281     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1282     !                                                                      C
1283     !  Module name: LEQ_IJSWEEPt(I,J, Vname, Var, A_m, B_m )               C
1284     !  Purpose: Perform line sweep at coordiante I,K                       C
1285     !                                                                      C
1286     !                                                                      C
1287     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1288     !  Reviewer:                                          Date:            C
1289     !                                                                      C
1290     !                                                                      C
1291     !  Literature/Document References:                                     C
1292     !                                                                      C
1293     !  Variables referenced:                                               C
1294     !  Variables modified:                                                 C
1295     !                                                                      C
1296     !  Local variables:                                                    C
1297     !                                                                      C
1298     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1299           SUBROUTINE LEQ_IJSWEEPt(I,J,Vname, VAR, A_M, B_M )
1300     
1301     !-----------------------------------------------
1302     !   M o d u l e s
1303     !-----------------------------------------------
1304           USE compar
1305           USE exit, only: mfix_exit
1306           USE functions
1307           USE funits
1308           USE geometry
1309           USE indices
1310           USE param
1311           USE param1
1312           IMPLICIT NONE
1313     !-----------------------------------------------
1314     !   G l o b a l   P a r a m e t e r s
1315     !-----------------------------------------------
1316     !-----------------------------------------------
1317     !   D u m m y   A r g u m e n t s
1318     !-----------------------------------------------
1319     !                      Line position
1320           INTEGER          I,J
1321     !
1322     !                      Septadiagonal matrix A_m
1323           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1324     !
1325     !                      Vector b_m
1326           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1327     !
1328     !                      Variable name
1329           CHARACTER(LEN=*)    Vname
1330     
1331     !
1332     !                      Variable
1333           DOUBLE PRECISION Var(ijkstart3:ijkend3)
1334     
1335     !-----------------------------------------------
1336     !   L o c a l   P a r a m e t e r s
1337     !-----------------------------------------------
1338     !
1339     
1340     !-----------------------------------------------
1341     !   L o c a l   V a r i a b l e s
1342     !-----------------------------------------------
1343     
1344           DOUBLE PRECISION, DIMENSION (KMAX2) :: CC,DD,EE,BB
1345           INTEGER :: NN, INFO, IJK, K
1346     
1347           NN = KMAX2
1348     
1349           DO K=1,NN
1350              IJK = FUNIJK(I,J,K)
1351     
1352              DD(K) = A_M(0,IJK)
1353              CC(K) = A_M(-3,IJK)
1354              EE(K) = A_M(3,IJK)
1355              BB(K) = B_M(IJK)    -  A_M(-2,IJK) * Var( JM_OF(IJK) )         &
1356              -  A_M(2,IJK) * Var( JP_OF(IJK) )         &
1357              -  A_M(-1,IJK) * Var( IM_OF(IJK) )         &
1358              -  A_M(1,IJK) * Var( IP_OF(IJK) )
1359     
1360           ENDDO
1361     
1362           CC(1) = ZERO
1363           EE(NN) = ZERO
1364     !     DL(1:NEND-1) = CC(2:NEND)
1365           INFO = 0
1366           CALL DGTSL( NN, CC, DD, EE, BB, INFO )
1367     !     CALL DGTSV( JEND-JSTART+1, 1, DL, DD, EE, BB,  JEND-JSTART+1, INFO )
1368     
1369           IF (INFO.NE.0) THEN
1370              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1371              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1372              call mfix_exit(myPE)
1373           ENDIF
1374     
1375           DO K=1,NN
1376              IJK = FUNIJK(I,J,K)
1377              Var(IJK) = BB(K)
1378           ENDDO
1379     
1380           RETURN
1381           END SUBROUTINE  LEQ_IJSWEEPt
1382     
1383     !-----------------------------------------------
1384           SUBROUTINE LEQ_MSOLVE0t(VNAME, B_m, A_M, Var, CMETHOD )
1385     !-----------------------------------------------
1386     !
1387     !-----------------------------------------------
1388     !   M o d u l e s
1389     !-----------------------------------------------
1390           USE param
1391           USE param1
1392           USE geometry
1393           USE compar
1394           USE indices
1395           USE sendrecv
1396     
1397           use parallel
1398           IMPLICIT NONE
1399     !-----------------------------------------------
1400     !     G l o b a l   P a r a m e t e r s
1401     !-----------------------------------------------
1402     !-----------------------------------------------
1403     !     D u m m y   A r g u m e n t s
1404     !-----------------------------------------------
1405     !
1406     !     Septadiagonal matrix A_m
1407           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1408     !
1409     !     Vector b_m
1410           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1411     !
1412     !     Variable name
1413           CHARACTER(LEN=*)    Vname
1414     !
1415     !     Variable
1416           DOUBLE PRECISION Var(ijkstart3:ijkend3)
1417           integer :: ijk
1418     
1419           CHARACTER(LEN=4) :: CMETHOD
1420     
1421     !     do nothing or no preconditioning
1422     
1423           if (use_doloop) then
1424     
1425     !!!$omp  parallel do private(ijk)
1426              do ijk=ijkstart3,ijkend3
1427                 var(ijk) = b_m(ijk)
1428              enddo
1429           else
1430              var(:) = b_m(:)
1431           endif
1432           call send_recv(var,nlayers_bicgs)
1433     
1434           return
1435           end subroutine leq_msolve0t
1436     
1437     !-----------------------------------------------
1438           SUBROUTINE LEQ_msolve1t(VNAME, B_m, A_M, Var, CMETHOD )
1439     !-----------------------------------------------
1440     !
1441     !-----------------------------------------------
1442     !   M o d u l e s
1443     !-----------------------------------------------
1444           USE param
1445           USE param1
1446           USE geometry
1447           USE compar
1448           USE indices
1449           USE sendrecv
1450           USE parallel
1451           USE functions
1452     
1453           IMPLICIT NONE
1454     !-----------------------------------------------
1455     !   G l o b a l   P a r a m e t e r s
1456     !-----------------------------------------------
1457     !-----------------------------------------------
1458     !   D u m m y   A r g u m e n t s
1459     !-----------------------------------------------
1460     !
1461     !     Septadiagonal matrix A_m
1462           DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1463     !
1464     !     Vector b_m
1465           DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1466     !
1467     !     Variable name
1468           CHARACTER(LEN=*)    Vname
1469     !
1470     !     Variable
1471           DOUBLE PRECISION Var(ijkstart3:ijkend3)
1472     
1473           CHARACTER(LEN=4) :: CMETHOD
1474     
1475           integer :: i,j,k, ijk2
1476     
1477           if (use_doloop) then
1478     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
1479     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 1516)
1480     !!!$omp    parallel do private(ijk2)
1481              do ijk2=ijkstart3,ijkend3
1482                 var(ijk2) = zero
1483              enddo
1484           else
1485              var(:) = ZERO
1486           endif
1487     
1488     !     diagonal scaling
1489     
1490     !AIKE PFUPGRADE 091409 Modified ijk to ijk2 to avoid compilation error since PF upgrade
1491     ! PGF90-S-0155-ijk may not appear in a PRIVATE clause (leq_bicgst.f: 1526)
1492     !!$omp   parallel do private(i,j,k,ijk2) collapse (3)
1493           do k=kstart2,kend2
1494              do i=istart2,iend2
1495                 do j=jstart2,jend2
1496     
1497                    ijk2 = funijk( i,j,k )
1498                    var(ijk2) = b_m(ijk2)/A_m(0,ijk2)
1499     
1500                 enddo
1501              enddo
1502           enddo
1503     
1504           call send_recv(var,nlayers_bicgs)
1505     
1506           return
1507           end subroutine leq_msolve1t
1508     
1509