File: RELATIVE:/../../../mfix.git/model/leq_bicgs.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine LEQ_BICGS                                                C
4     !  Purpose: Solve system of linear system using BICGS method           C
5     !           Biconjugate gradients stabilized                           C
6     !                                                                      C
7     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !  Variables referenced:                                               C
12     !  Variables modified:                                                 C
13     !  Local variables:                                                    C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     ! Handan Liu wrote below:               !Jan 22 2013
17     ! The modification is as below:
18     !       Adding a loop of 2D RSRS sweep and parallelizing for OpenMP.
19     !       Splitting the existing 3D RSRS loop into two loops for OpenMP
20     !               due to data dependency
21     !       Adding openmp directives in all loops in leq_bicgs0.
22     !       Setting 'use_doloop = ture' to activate using do-loop.
23     !
24           SUBROUTINE LEQ_BICGS(VNAME, VNO, VAR, A_M, B_m, cmethod, &
25                                TOL, PC, ITMAX, IER)
26     
27     !-----------------------------------------------
28     ! Modules
29     !-----------------------------------------------
30           USE param
31           USE param1
32           USE matrix
33           USE geometry
34           USE compar
35           USE indices
36           USE leqsol
37           USE funits
38           IMPLICIT NONE
39     !-----------------------------------------------
40     ! Dummy arguments
41     !-----------------------------------------------
42     ! variable name
43           CHARACTER(LEN=*), INTENT(IN) :: Vname
44     ! variable number (not really used here; see calling subroutine)
45           INTEGER, INTENT(IN) :: VNO
46     ! variable
47     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
48     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
49     !     e_Turb_G
50     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: Var
51           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
52     ! Septadiagonal matrix A_m
53     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3,-3:3), INTENT(INOUT) :: A_m
54          DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
55     
56     ! Vector b_m
57     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: B_m
58           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
59     ! Sweep direction of leq solver (leq_sweep)
60     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
61     ! Note: this setting only seems to matter when leq_pc='line'
62           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
63     ! convergence tolerance (generally leq_tol)
64           DOUBLE PRECISION, INTENT(IN) :: TOL
65     ! preconditioner (leq_pc)
66     !     options = 'line' (default), 'diag', 'none'
67           CHARACTER(LEN=4), INTENT(IN) ::  PC
68     ! maximum number of iterations (generally leq_it)
69           INTEGER, INTENT(IN) :: ITMAX
70     ! error indicator
71           INTEGER, INTENT(INOUT) :: IER
72     !-------------------------------------------------
73     ! Local Variables
74     !-------------------------------------------------
75     
76           if(PC.eq.'LINE') then   ! default
77              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,  &
78                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE, IER )
79           elseif(PC.eq.'DIAG') then
80              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
81                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE1, IER )
82           elseif(PC.eq.'NONE') then
83              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
84                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE0, IER )
85           else
86              IF(DMP_LOG)WRITE (UNIT_LOG,*) &
87                'preconditioner option not found - check mfix.dat and readme'
88              call mfix_exit(myPE)
89           endif
90     
91           return
92           END SUBROUTINE LEQ_BICGS
93     
94     
95     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
96     !                                                                      C
97     !  Subroutine: LEQ_BICGS0                                              C
98     !  Purpose: Compute residual of linear system                          C
99     !                                                                      C
100     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
101     !  Reviewer:                                          Date:            C
102     !                                                                      C
103     !  Literature/Document References:                                     C
104     !  Variables referenced:                                               C
105     !  Variables modified:                                                 C
106     !  Local variables:                                                    C
107     !                                                                      C
108     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
109     
110           SUBROUTINE LEQ_BICGS0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
111                                 TOL, ITMAX, MATVEC, MSOLVE, IER )
112     
113     !-----------------------------------------------
114     ! Modules
115     !-----------------------------------------------
116           USE param
117           USE param1
118           USE parallel
119           USE matrix
120           USE geometry
121           USE compar
122           USE mpi_utility
123           USE sendrecv
124           USE indices
125           USE leqsol
126           USE cutcell
127           USE functions
128     
129           IMPLICIT NONE
130     !-----------------------------------------------
131     ! Dummy arguments/procedure
132     !-----------------------------------------------
133     ! variable name
134           CHARACTER(LEN=*), INTENT(IN) :: Vname
135     ! variable number (not really used here-see calling subroutine)
136           INTEGER, INTENT(IN) :: VNO
137     ! variable
138     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
139     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
140     !     e_Turb_G
141     !      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
142           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
143     ! Septadiagonal matrix A_m
144     !      DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
145           DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
146     ! Vector b_m
147     !      DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
148           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
149     ! Sweep direction of leq solver (leq_sweep)
150     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
151           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
152     ! convergence tolerance (generally leq_tol)
153           DOUBLE PRECISION, INTENT(IN) :: TOL
154     ! maximum number of iterations (generally leq_it)
155           INTEGER, INTENT(IN) :: ITMAX
156     ! error indicator
157           INTEGER, INTENT(INOUT) :: IER
158     ! dummy arguments/procedures set as indicated
159     !     matvec->leq_matvec
160     ! for preconditioner (leq_pc)
161     !    'line' msolve->leq_msolve  (default)
162     !    'diag' msolve->leq_msolve1
163     !    'none' msolve->leq_msolve0
164     !-----------------------------------------------
165     ! Local parameters
166     !-----------------------------------------------
167           INTEGER, PARAMETER :: idebugl = 0
168           DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
169           LOGICAL, PARAMETER :: do_unit_scaling = .true.
170     !-----------------------------------------------
171     ! Local variables
172     !-----------------------------------------------
173     
174           DOUBLE PRECISION, DIMENSION(:), allocatable :: R,Rtilde, P,Phat, Svec, Shat, Tvec,V
175     
176           DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
177                             alpha, beta, omega, rho
178           DOUBLE PRECISION :: TxS, TxT, RtildexV, RtildexR, &
179                               aijmax, oam
180           DOUBLE PRECISION :: Rnorm, Rnorm0, Snorm, TOLMIN, pnorm
181           LOGICAL :: isconverged
182           INTEGER :: i, j, k, ijk
183           INTEGER :: iter
184           DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
185     !-----------------------------------------------
186     
187           allocate(R(DIMENSION_3))
188           allocate(Rtilde(DIMENSION_3))
189           allocate(P(DIMENSION_3))
190           allocate(Phat(DIMENSION_3))
191           allocate(Svec(DIMENSION_3))
192           allocate(Shat(DIMENSION_3))
193           allocate(Tvec(DIMENSION_3))
194           allocate(V(DIMENSION_3))
195     
196           is_serial = numPEs.eq.1.and.is_serial
197     
198     ! these scalars should not be necessary to initialize but done as failsafe
199           rnorm = ZERO
200           rnorm0 = ZERO
201           snorm = ZERO
202           pnorm = ZERO
203     
204     ! initializing
205           alpha(:)  = zero
206           beta(:)   = zero
207           omega(:)  = zero
208           rho(:)    = zero
209     
210               use_doloop = .TRUE.   ! set true by Handan Liu for OpenMP do-loop
211     
212     ! Adding all by Handan Liu
213     ! zero out R, Rtilde, P, Phat, Svec, Shat, Tvec, V
214     ! --------------------------------
215           if (use_doloop) then   ! mfix.dat keyword default=false
216     !$omp  parallel do default(shared) private(ijk)
217              do ijk=ijkstart3,ijkend3
218                 R(ijk) = zero
219                 Rtilde(ijk) = zero
220                 P(ijk) = zero
221                 Phat(ijk) = zero
222                 Svec(ijk) = zero
223                 Shat(ijk) = zero
224                 Tvec(ijk) = zero
225                 V(ijk) = zero
226              enddo
227     
228           else
229              R(:) = zero
230              Rtilde(:) = zero
231              P(:) = zero
232              Phat(:) = zero
233              Svec(:) = zero
234              Shat(:) = zero
235              Tvec(:) = zero
236              V(:) = zero
237           endif
238     
239           TOLMIN = EPSILON( one )
240     
241     ! Scale matrix to have unit diagonal
242     ! ---------------------------------------------------------------->>>
243           if (do_unit_scaling) then
244     
245     
246              IF(RE_INDEXING) THEN  ! Loop only over active cells
247     !$omp parallel do default(shared) private(ijk,oam,aijmax)
248                 DO IJK = IJKSTART3,IJKEND3
249                    aijmax = maxval(abs(A_M(ijk,:)) )
250                    OAM = one/aijmax
251                    A_M(IJK,:) = A_M(IJK,:)*OAM
252                    B_M(IJK) = B_M(IJK)*OAM
253                 ENDDO
254     
255              ELSE
256     
257     !$omp parallel do default(shared) private(ijk,i,j,k,oam,aijmax)
258                 do k = kstart2,kend2
259                    do i = istart2,iend2
260                       do j = jstart2,jend2
261                          IJK = funijk(i,j,k)
262                          aijmax = maxval(abs(A_M(ijk,:)) )
263                          OAM = one/aijmax
264                          A_M(IJK,:) = A_M(IJK,:)*OAM
265                          B_M(IJK) = B_M(IJK)*OAM
266                       enddo
267                    enddo
268                 enddo
269     
270              ENDIF
271     
272           endif
273     ! ----------------------------------------------------------------<<<
274     
275     
276     ! Compute initial residual (R = b-A*x) for Ax=b
277     !    assume initial guess in Var
278     !    rtilde = r
279     ! ---------------------------------------------------------------->>>
280           call MATVEC(Vname, Var, A_M, R)   ! returns R=A*Var
281     
282           if (use_doloop) then   ! mfix.dat keyword default=false
283     !$omp parallel do default(shared) private(ijk)
284              do ijk=ijkstart3,ijkend3
285                 R(ijk) = B_m(ijk) - R(ijk)
286              enddo
287           else
288              R(:) = B_m(:) - R(:)
289           endif
290     
291           call send_recv(R,nlayers_bicgs)
292     
293           if(is_serial) then
294              Rnorm0 = zero
295              if (use_doloop) then   ! mfix.dat keyword default=false
296     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm0)
297                 do ijk=ijkstart3,ijkend3
298                    Rnorm0 = Rnorm0 + R(ijk)*R(ijk)
299                 enddo
300              else
301                 Rnorm0 = dot_product(R,R)
302              endif
303              Rnorm0 = sqrt( Rnorm0 )
304           else
305              Rnorm0 = sqrt( dot_product_par( R, R ) )
306           endif
307     
308     ! determine an initial guess for the residual = residual + small random
309     ! number (so it could be set to anything). note that since random_number
310     ! is used to supply the guess, this line could potentially be the source
311     ! of small differences between runs.  the random number is shifted below
312     ! between -1 and 1 and then scaled by factor 1.0D-6*Rnorm0
313           call random_number(Rtilde(:))
314     
315     ! Shift random number array to be consistent with case when RE_INDEXING is .FALSE.
316            IF(RE_INDEXING) CALL SHIFT_DP_ARRAY(Rtilde)
317     
318     
319           if (use_doloop) then   ! mfix.dat keyword default=false
320     !$omp parallel do default(shared) private(ijk)
321              do ijk=ijkstart3,ijkend3
322                 Rtilde(ijk) = R(ijk) + (2.0d0*Rtilde(ijk)-1.0d0)*1.0d-6*Rnorm0
323              enddo
324           else
325              Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
326           endif
327     
328           if (idebugl >= 1) then
329              if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
330           endif
331     ! ----------------------------------------------------------------<<<
332     
333     
334     ! Main loop
335     ! ---------------------------------------------------------------->>>
336           iter = 1
337           do i=1,itmax
338     
339              if(is_serial) then
340                 if (use_doloop) then   ! mfix.dat keyword default=false
341                    RtildexR = zero
342     !$omp parallel do default(shared) private(ijk) reduction(+:RtildexR)
343                    do ijk=ijkstart3,ijkend3
344                       RtildexR = RtildexR + Rtilde(ijk) * R(ijk)
345                    enddo
346                    rho(i-1) = RtildexR
347                 else
348                    rho(i-1) = dot_product( Rtilde, R )
349                 endif
350              else
351                 rho(i-1) = dot_product_par( Rtilde, R )
352              endif ! is_serial
353     
354     !         print*,'leq_bicgs, initial: ', Vname,' rho(i-1) ', rho(i-1)
355     
356              if (rho(i-1) .eq. zero) then
357                 if(i /= 1)then
358     ! Method fails
359     ! --------------------------------
360     !               print*, 'leq_bicgs,',Vname,': rho(i-1) == 0 '
361                    ier = -2
362                 else
363     ! Method converged.  residual is already zero
364     ! --------------------------------
365                    ier = 0
366                 endif
367                 call send_recv(var,2)
368                 return
369              endif ! rho(i-1).eq.0
370     
371              if (i .eq. 1) then
372                 if (use_doloop) then
373     !!$omp        parallel do private(ijk)
374     !$omp parallel do default(shared) private(ijk)
375                    do ijk=ijkstart3,ijkend3
376                       P(ijk) = R(ijk)
377                    enddo
378                 else
379                    P(:) = R(:)
380                 endif
381              else
382                 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
383                 if (use_doloop) then
384     !$omp parallel do default(shared) private(ijk)
385                    do ijk=ijkstart3,ijkend3
386                       P(ijk) = R(ijk) + beta(i-1)*( P(ijk) - omega(i-1)*V(ijk) )
387                    enddo
388                 else
389                    P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
390                 endif
391              endif ! i.eq.1
392     
393     
394     ! Solve A*Phat(:) = P(:)
395     ! V(:) = A*Phat(:)
396     ! --------------------------------
397              call MSOLVE(Vname, P, A_m, Phat, CMETHOD) ! returns Phat
398              call MATVEC(Vname, Phat, A_m, V)   ! returns V=A*Phat
399     
400              if(is_serial) then
401                 if (use_doloop) then
402                    RtildexV = zero
403     !$omp  parallel do default(shared) private(ijk) reduction(+:RtildexV)
404                    do ijk=ijkstart3,ijkend3
405                       RtildexV = RtildexV + Rtilde(ijk) * V(ijk)
406                    enddo
407                 else
408                    RtildexV = dot_product( Rtilde, V )
409                 endif
410              else
411                 RtildexV = dot_product_par( Rtilde, V )
412              endif ! is_serial
413     
414     !         print*,'leq_bicgs, initial: ', Vname,' RtildexV ', RtildexV
415     
416     ! compute alpha
417     ! --------------------------------
418              alpha(i) = rho(i-1) / RtildexV
419     
420     ! compute Svec
421     ! --------------------------------
422              if (use_doloop) then
423     !!$omp     parallel do private(ijk)
424     !$omp parallel do default(shared) private(ijk)
425                 do ijk=ijkstart3,ijkend3
426                    Svec(ijk) = R(ijk) - alpha(i) * V(ijk)
427                 enddo
428     !!$omp end parallel do
429              else
430                 Svec(:) = R(:) - alpha(i) * V(:)
431              endif ! use_doloop
432     
433     
434     ! Check norm of Svec(:); if small enough:
435     ! set X(:) = X(:) + alpha(i)*Phat(:) and stop
436     ! --------------------------------
437              if(.not.minimize_dotproducts) then
438                 if(is_serial) then
439                    if (use_doloop) then
440                       Snorm = zero
441     !$omp parallel do default(shared) private(ijk) reduction(+:Snorm)
442                       do ijk=ijkstart3,ijkend3
443                          Snorm = Snorm + Svec(ijk) * Svec(ijk)
444                       enddo
445                    else
446                       Snorm = dot_product( Svec, Svec )
447                    endif
448                    Snorm = sqrt( Snorm )
449                 else
450                    Snorm = sqrt( dot_product_par( Svec, Svec ) )
451                 endif               ! is_serial
452     !            print*,'leq_bicgs, initial: ', Vname,' Snorm ', real(Snorm)
453     
454     
455                 if (Snorm <= TOLMIN) then
456                    if (use_doloop) then
457     !$omp parallel do default(shared) private(ijk)
458                       do ijk=ijkstart3,ijkend3
459                          Var(ijk) = Var(ijk) + alpha(i)*Phat(ijk)
460                       enddo
461                    else
462                       Var(:) = Var(:) + alpha(i)*Phat(:)
463                    endif            ! use_doloop
464     
465     ! Recompute residual norm
466     ! --------------------------------
467                    if (idebugl >= 1) then
468                       call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
469     !                  Rnorm = sqrt( dot_product_par( Var, Var ) )
470     !                  print*,'leq_bicgs, initial: ', Vname,' Vnorm ', Rnorm
471     
472                       if (use_doloop) then
473     !$omp parallel do default(shared) private(ijk)
474                          do ijk=ijkstart3,ijkend3
475                             R(ijk) = B_m(ijk) - R(ijk)
476                          enddo
477                       else
478                          R(:) = B_m(:) - R(:)
479                       endif
480     
481                       if(is_serial) then
482                          if (use_doloop) then
483                             Rnorm = zero
484     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
485                             do ijk=ijkstart3,ijkend3
486                                Rnorm = Rnorm + R(ijk)*R(ijk)
487                             enddo
488                          else
489                             Rnorm =  dot_product( R, R )
490                          endif
491                          Rnorm = sqrt( Rnorm )
492                       else
493                          Rnorm = sqrt( dot_product_par( R, R ) )
494                       endif
495     !                  print*,'leq_bicgs, initial: ', Vname,' Rnorm ', Rnorm
496                    endif            ! idebugl >= 1
497                                isConverged = .TRUE.
498                                EXIT
499                 endif               ! end if (Snorm <= TOLMIN)
500              endif                  ! end if (.not.minimize_dotproducts)
501     
502     
503     
504     ! Solve A*Shat(:) = Svec(:)
505     ! Tvec(:) = A*Shat(:)
506     ! --------------------------------
507              call MSOLVE( Vname, Svec, A_m, Shat, CMETHOD)  ! returns Shat
508              call MATVEC( Vname, Shat, A_m, Tvec )   ! returns Tvec=A*Shat
509     
510              if(is_serial) then
511                 if (use_doloop) then
512                    TxS = zero
513                    TxT = zero
514     !$omp  parallel do default(shared) private(ijk) reduction(+:TxS,TxT)
515                    do ijk=ijkstart3,ijkend3
516                       TxS = TxS + Tvec(ijk)  * Svec(ijk)
517                       TxT = TxT + Tvec(ijk)  * Tvec(ijk)
518                    enddo
519                 else
520                    TxS = dot_product( Tvec, Svec )
521                    TxT = dot_product( Tvec, Tvec )
522                 endif
523              else
524                 if(.not.minimize_dotproducts) then
525                    TxS = dot_product_par( Tvec, Svec )
526                    TxT = dot_product_par( Tvec, Tvec )
527                 else
528                    TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
529                    TxS = TxS_TxT(1)
530                    TxT = TxS_TxT(2)
531                 endif
532              endif
533     
534              IF(TxT.eq.Zero) TxT = SMALL_NUMBER
535     
536     ! compute omega
537     ! --------------------------------
538              omega(i) = TxS / TxT
539     
540     ! compute new guess for Var
541     ! --------------------------------
542              if (use_doloop) then
543     !$omp parallel do default(shared) private(ijk)
544                 do ijk=ijkstart3,ijkend3
545                    Var(ijk) = Var(ijk) +  &
546                       alpha(i)*Phat(ijk) + omega(i)*Shat(ijk)
547                       R(ijk) = Svec(ijk) - omega(i)*Tvec(ijk)
548                 enddo
549              else
550                 Var(:) = Var(:) +  &
551                    alpha(i)*Phat(:) + omega(i)*Shat(:)
552                    R(:) = Svec(:) - omega(i)*Tvec(:)
553              endif
554     
555     
556     ! --------------------------------
557              if(.not.minimize_dotproducts.or.(mod(iter,icheck_bicgs).eq.0)) then
558                 if(is_serial) then
559                    if (use_doloop) then
560                       Rnorm = zero
561     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
562                       do ijk=ijkstart3,ijkend3
563                          Rnorm = Rnorm + R(ijk) * R(ijk)
564                       enddo
565                    else
566                       Rnorm =  dot_product(R, R )
567                    endif
568                    Rnorm = sqrt( Rnorm )
569                 else
570                    Rnorm = sqrt( dot_product_par(R, R) )
571                 endif               ! is_serial
572     
573                 if (idebugl.ge.1) then
574                    if (myPE.eq.PE_IO) then
575                       print*,'iter, Rnorm ', iter, Rnorm, Snorm
576                       print*,'alpha(i), omega(i) ', alpha(i), omega(i)
577                       print*,'TxS, TxT ', TxS, TxT
578                       print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
579                    endif
580                 endif
581     
582     !           call mfix_exit(myPE)
583     
584     ! Check convergence; continue if necessary
585     ! for continuation, it is necessary that omega(i) .ne. 0
586                 isconverged = (Rnorm <= TOL*Rnorm0)
587     
588                 if (isconverged) then
589                    iter_tot(vno) = iter_tot(vno) + iter + 1
590                    EXIT
591                 endif
592              endif                  ! end if(.not.minimize_dotproducts)
593     
594     ! Advance the iteration count
595              iter = iter + 1
596     
597           enddo   ! end do i=1,itmax
598     ! end of linear solver loop
599     ! ----------------------------------------------------------------<<<
600     
601     
602           if (idebugl >= 1) then
603              call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
604              if (use_doloop) then
605     !$omp  parallel do default(shared) private(ijk)
606                 do ijk=ijkstart3,ijkend3
607                    R(ijk) = R(ijk) - B_m(ijk)
608                 enddo
609              else
610                 R(:) = R(:) - B_m(:)
611              endif
612     
613              if(is_serial) then
614                 if (use_doloop) then
615                    Rnorm = zero
616     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
617                    do ijk=ijkstart3,ijkend3
618                       Rnorm = Rnorm + R(ijk) * R(ijk)
619                    enddo
620                 else
621                    Rnorm = dot_product( R,R)
622                 endif
623                 Rnorm = sqrt( Rnorm )
624              else
625                 Rnorm = sqrt( dot_product_par( R,R) )
626              endif
627     
628              if(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
629     
630              if(myPE.eq.0)  print*,'leq_bicgs ratio : ', Vname,' ',iter,  &
631              ' L-2', Rnorm/Rnorm0
632           endif   ! end if(idebugl >=1)
633     
634     !      isconverged = (real(Rnorm) <= TOL*Rnorm0);
635           if(.NOT.isConverged) isconverged = (real(Rnorm) <= TOL*Rnorm0);
636     !     write(*,*) '***',iter, isconverged, Rnorm, TOL, Rnorm0, myPE
637           IER = 0
638           if (.not.isconverged) then
639              IER = -1
640              iter_tot(vno) = iter_tot(vno) + iter
641              if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
642                 IER = -2
643              endif
644           endif
645     
646           call send_recv(var,2)
647     
648           deallocate(R)
649           deallocate(Rtilde)
650           deallocate(P)
651           deallocate(Phat)
652           deallocate(Svec)
653           deallocate(Shat)
654           deallocate(Tvec)
655           deallocate(V)
656     
657           return
658           end subroutine LEQ_BICGS0
659