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