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