File: /nfs/home/0/users/jenkins/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     ! External subroutines
77     !-------------------------------------------------
78     ! These procedures are effectively dummy arguments (procedures as
79     ! arguments within the subroutine leq_bicgs0)
80           EXTERNAL LEQ_MATVEC, LEQ_MSOLVE, LEQ_MSOLVE0, LEQ_MSOLVE1
81     !--------------------------------------------------
82     
83     
84           if(PC.eq.'LINE') then   ! default
85              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,  &
86                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE, IER )
87           elseif(PC.eq.'DIAG') then
88              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
89                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE1, IER )
90           elseif(PC.eq.'NONE') then
91              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
92                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE0, IER )
93           else
94              IF(DMP_LOG)WRITE (UNIT_LOG,*) &
95                'preconditioner option not found - check mfix.dat and readme'
96              call mfix_exit(myPE)
97           endif
98     
99           return
100           END SUBROUTINE LEQ_BICGS
101     
102     
103     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
104     !                                                                      C
105     !  Subroutine: LEQ_BICGS0                                              C
106     !  Purpose: Compute residual of linear system                          C
107     !                                                                      C
108     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
109     !  Reviewer:                                          Date:            C
110     !                                                                      C
111     !  Literature/Document References:                                     C
112     !  Variables referenced:                                               C
113     !  Variables modified:                                                 C
114     !  Local variables:                                                    C
115     !                                                                      C
116     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
117     
118           SUBROUTINE LEQ_BICGS0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
119                                 TOL, ITMAX, MATVEC, MSOLVE, IER )
120     
121     !-----------------------------------------------
122     ! Modules
123     !-----------------------------------------------
124           USE param
125           USE param1
126           USE parallel
127           USE matrix
128           USE geometry
129           USE compar
130           USE mpi_utility
131           USE sendrecv
132           USE indices
133           USE leqsol
134           USE cutcell
135           USE functions
136     
137           IMPLICIT NONE
138     !-----------------------------------------------
139     ! Dummy arguments/procedure
140     !-----------------------------------------------
141     ! variable name
142           CHARACTER(LEN=*), INTENT(IN) :: Vname
143     ! variable number (not really used here-see calling subroutine)
144           INTEGER, INTENT(IN) :: VNO
145     ! variable
146     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
147     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
148     !     e_Turb_G
149     !      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
150           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
151     ! Septadiagonal matrix A_m
152     !      DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
153           DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
154     ! Vector b_m
155     !      DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
156           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
157     ! Sweep direction of leq solver (leq_sweep)
158     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
159           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
160     ! convergence tolerance (generally leq_tol)
161           DOUBLE PRECISION, INTENT(IN) :: TOL
162     ! maximum number of iterations (generally leq_it)
163           INTEGER, INTENT(IN) :: ITMAX
164     ! error indicator
165           INTEGER, INTENT(INOUT) :: IER
166     ! dummy arguments/procedures set as indicated
167     !     matvec->leq_matvec
168     ! for preconditioner (leq_pc)
169     !    'line' msolve->leq_msolve  (default)
170     !    'diag' msolve->leq_msolve1
171     !    'none' msolve->leq_msolve0
172           EXTERNAL  MATVEC, MSOLVE
173     !-----------------------------------------------
174     ! Local parameters
175     !-----------------------------------------------
176           INTEGER, PARAMETER :: idebugl = 0
177           DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
178           LOGICAL, PARAMETER :: do_unit_scaling = .true.
179     !-----------------------------------------------
180     ! Local variables
181     !-----------------------------------------------
182     
183           DOUBLE PRECISION, DIMENSION(:), allocatable :: R,Rtilde, P,Phat, Svec, Shat, Tvec,V
184     
185           DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
186                             alpha, beta, omega, rho
187           DOUBLE PRECISION :: TxS, TxT, RtildexV, RtildexR, &
188                               aijmax, oam
189           DOUBLE PRECISION :: Rnorm, Rnorm0, Snorm, TOLMIN, pnorm
190           LOGICAL :: isconverged
191           INTEGER :: i, j, k, ijk
192           INTEGER :: iter
193           DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
194     !-----------------------------------------------
195     ! External subroutines/functions
196     !-----------------------------------------------
197     
198           INTERFACE
199              DOUBLE PRECISION FUNCTION DOT_PRODUCT_PAR( R1, R2 )
200              use compar
201              USE param
202     !         DOUBLE PRECISION, INTENT(IN), DIMENSION(ijkstart3:ijkend3) :: R1,R2
203              DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_3) :: R1,R2
204              END FUNCTION DOT_PRODUCT_PAR
205           END INTERFACE
206     
207           INTERFACE
208              FUNCTION DOT_PRODUCT_PAR2( R1, R2, R3, R4 )
209              use compar
210              USE param
211              USE functions
212              DOUBLE PRECISION, INTENT(IN), DIMENSION(ijkstart3:ijkend3) :: &
213                                            R1, R2, R3, R4
214              DOUBLE PRECISION, DIMENSION(2) :: DOT_PRODUCT_PAR2
215              END FUNCTION DOT_PRODUCT_PAR2
216           END INTERFACE
217     
218     !-----------------------------------------------
219     
220           allocate(R(DIMENSION_3))
221           allocate(Rtilde(DIMENSION_3))
222           allocate(P(DIMENSION_3))
223           allocate(Phat(DIMENSION_3))
224           allocate(Svec(DIMENSION_3))
225           allocate(Shat(DIMENSION_3))
226           allocate(Tvec(DIMENSION_3))
227           allocate(V(DIMENSION_3))
228     
229           is_serial = numPEs.eq.1.and.is_serial
230     
231     ! these scalars should not be necessary to initialize but done as failsafe
232           rnorm = ZERO
233           rnorm0 = ZERO
234           snorm = ZERO
235           pnorm = ZERO
236     
237     ! initializing
238           alpha(:)  = zero
239           beta(:)   = zero
240           omega(:)  = zero
241           rho(:)    = zero
242     
243               use_doloop = .TRUE.   ! set true by Handan Liu for OpenMP do-loop
244     
245     ! Adding all by Handan Liu
246     ! zero out R, Rtilde, P, Phat, Svec, Shat, Tvec, V
247     ! --------------------------------
248           if (use_doloop) then   ! mfix.dat keyword default=false
249     !$omp  parallel do default(shared) private(ijk)
250              do ijk=ijkstart3,ijkend3
251                 R(ijk) = zero
252                 Rtilde(ijk) = zero
253                 P(ijk) = zero
254                 Phat(ijk) = zero
255                 Svec(ijk) = zero
256                 Shat(ijk) = zero
257                 Tvec(ijk) = zero
258                 V(ijk) = zero
259              enddo
260     
261           else
262              R(:) = zero
263              Rtilde(:) = zero
264              P(:) = zero
265              Phat(:) = zero
266              Svec(:) = zero
267              Shat(:) = zero
268              Tvec(:) = zero
269              V(:) = zero
270           endif
271     
272           TOLMIN = EPSILON( one )
273     
274     ! Scale matrix to have unit diagonal
275     ! ---------------------------------------------------------------->>>
276           if (do_unit_scaling) then
277     
278     
279              IF(RE_INDEXING) THEN  ! Loop only over active cells
280     !$omp parallel do default(shared) private(ijk,oam,aijmax)
281                 DO IJK = IJKSTART3,IJKEND3
282                    aijmax = maxval(abs(A_M(ijk,:)) )
283                    OAM = one/aijmax
284                    A_M(IJK,:) = A_M(IJK,:)*OAM
285                    B_M(IJK) = B_M(IJK)*OAM
286                 ENDDO
287     
288              ELSE
289     
290     !$omp parallel do default(shared) private(ijk,i,j,k,oam,aijmax)
291                 do k = kstart2,kend2
292                    do i = istart2,iend2
293                       do j = jstart2,jend2
294                          IJK = funijk(i,j,k)
295                          aijmax = maxval(abs(A_M(ijk,:)) )
296                          OAM = one/aijmax
297                          A_M(IJK,:) = A_M(IJK,:)*OAM
298                          B_M(IJK) = B_M(IJK)*OAM
299                       enddo
300                    enddo
301                 enddo
302     
303              ENDIF
304     
305           endif
306     ! ----------------------------------------------------------------<<<
307     
308     
309     ! Compute initial residual (R = b-A*x) for Ax=b
310     !    assume initial guess in Var
311     !    rtilde = r
312     ! ---------------------------------------------------------------->>>
313           call MATVEC(Vname, Var, A_M, R)   ! returns R=A*Var
314     
315           if (use_doloop) then   ! mfix.dat keyword default=false
316     !$omp parallel do default(shared) private(ijk)
317              do ijk=ijkstart3,ijkend3
318                 R(ijk) = B_m(ijk) - R(ijk)
319              enddo
320           else
321              R(:) = B_m(:) - R(:)
322           endif
323     
324           call send_recv(R,nlayers_bicgs)
325     
326           if(is_serial) then
327              Rnorm0 = zero
328              if (use_doloop) then   ! mfix.dat keyword default=false
329     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm0)
330                 do ijk=ijkstart3,ijkend3
331                    Rnorm0 = Rnorm0 + R(ijk)*R(ijk)
332                 enddo
333              else
334                 Rnorm0 = dot_product(R,R)
335              endif
336              Rnorm0 = sqrt( Rnorm0 )
337           else
338              Rnorm0 = sqrt( dot_product_par( R, R ) )
339           endif
340     
341     ! determine an initial guess for the residual = residual + small random
342     ! number (so it could be set to anything). note that since random_number
343     ! is used to supply the guess, this line could potentially be the source
344     ! of small differences between runs.  the random number is shifted below
345     ! between -1 and 1 and then scaled by factor 1.0D-6*Rnorm0
346           call random_number(Rtilde(:))
347     
348     ! Shift random number array to be consistent with case when RE_INDEXING is .FALSE.
349            IF(RE_INDEXING) CALL SHIFT_DP_ARRAY(Rtilde)
350     
351     
352           if (use_doloop) then   ! mfix.dat keyword default=false
353     !$omp parallel do default(shared) private(ijk)
354              do ijk=ijkstart3,ijkend3
355                 Rtilde(ijk) = R(ijk) + (2.0d0*Rtilde(ijk)-1.0d0)*1.0d-6*Rnorm0
356              enddo
357           else
358              Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
359           endif
360     
361           if (idebugl >= 1) then
362              if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
363           endif
364     ! ----------------------------------------------------------------<<<
365     
366     
367     ! Main loop
368     ! ---------------------------------------------------------------->>>
369           iter = 1
370           do i=1,itmax
371     
372              if(is_serial) then
373                 if (use_doloop) then   ! mfix.dat keyword default=false
374                    RtildexR = zero
375     !$omp parallel do default(shared) private(ijk) reduction(+:RtildexR)
376                    do ijk=ijkstart3,ijkend3
377                       RtildexR = RtildexR + Rtilde(ijk) * R(ijk)
378                    enddo
379                    rho(i-1) = RtildexR
380                 else
381                    rho(i-1) = dot_product( Rtilde, R )
382                 endif
383              else
384                 rho(i-1) = dot_product_par( Rtilde, R )
385              endif ! is_serial
386     
387     !         print*,'leq_bicgs, initial: ', Vname,' rho(i-1) ', rho(i-1)
388     
389              if (rho(i-1) .eq. zero) then
390                 if(i /= 1)then
391     ! Method fails
392     ! --------------------------------
393     !               print*, 'leq_bicgs,',Vname,': rho(i-1) == 0 '
394                    ier = -2
395                 else
396     ! Method converged.  residual is already zero
397     ! --------------------------------
398                    ier = 0
399                 endif
400                 call send_recv(var,2)
401                 return
402              endif ! rho(i-1).eq.0
403     
404              if (i .eq. 1) then
405                 if (use_doloop) then
406     !!$omp        parallel do private(ijk)
407     !$omp parallel do default(shared) private(ijk)
408                    do ijk=ijkstart3,ijkend3
409                       P(ijk) = R(ijk)
410                    enddo
411                 else
412                    P(:) = R(:)
413                 endif
414              else
415                 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
416                 if (use_doloop) then
417     !$omp parallel do default(shared) private(ijk)
418                    do ijk=ijkstart3,ijkend3
419                       P(ijk) = R(ijk) + beta(i-1)*( P(ijk) - omega(i-1)*V(ijk) )
420                    enddo
421                 else
422                    P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
423                 endif
424              endif ! i.eq.1
425     
426     
427     ! Solve A*Phat(:) = P(:)
428     ! V(:) = A*Phat(:)
429     ! --------------------------------
430              call MSOLVE(Vname, P, A_m, Phat, CMETHOD) ! returns Phat
431              call MATVEC(Vname, Phat, A_m, V)   ! returns V=A*Phat
432     
433              if(is_serial) then
434                 if (use_doloop) then
435                    RtildexV = zero
436     !$omp  parallel do default(shared) private(ijk) reduction(+:RtildexV)
437                    do ijk=ijkstart3,ijkend3
438                       RtildexV = RtildexV + Rtilde(ijk) * V(ijk)
439                    enddo
440                 else
441                    RtildexV = dot_product( Rtilde, V )
442                 endif
443              else
444                 RtildexV = dot_product_par( Rtilde, V )
445              endif ! is_serial
446     
447     !         print*,'leq_bicgs, initial: ', Vname,' RtildexV ', RtildexV
448     
449     ! compute alpha
450     ! --------------------------------
451              alpha(i) = rho(i-1) / RtildexV
452     
453     ! compute Svec
454     ! --------------------------------
455              if (use_doloop) then
456     !!$omp     parallel do private(ijk)
457     !$omp parallel do default(shared) private(ijk)
458                 do ijk=ijkstart3,ijkend3
459                    Svec(ijk) = R(ijk) - alpha(i) * V(ijk)
460                 enddo
461     !!$omp end parallel do
462              else
463                 Svec(:) = R(:) - alpha(i) * V(:)
464              endif ! use_doloop
465     
466     
467     ! Check norm of Svec(:); if small enough:
468     ! set X(:) = X(:) + alpha(i)*Phat(:) and stop
469     ! --------------------------------
470              if(.not.minimize_dotproducts) then
471                 if(is_serial) then
472                    if (use_doloop) then
473                       Snorm = zero
474     !$omp parallel do default(shared) private(ijk) reduction(+:Snorm)
475                       do ijk=ijkstart3,ijkend3
476                          Snorm = Snorm + Svec(ijk) * Svec(ijk)
477                       enddo
478                    else
479                       Snorm = dot_product( Svec, Svec )
480                    endif
481                    Snorm = sqrt( Snorm )
482                 else
483                    Snorm = sqrt( dot_product_par( Svec, Svec ) )
484                 endif               ! is_serial
485     !            print*,'leq_bicgs, initial: ', Vname,' Snorm ', real(Snorm)
486     
487     
488                 if (Snorm <= TOLMIN) then
489                    if (use_doloop) then
490     !$omp parallel do default(shared) private(ijk)
491                       do ijk=ijkstart3,ijkend3
492                          Var(ijk) = Var(ijk) + alpha(i)*Phat(ijk)
493                       enddo
494                    else
495                       Var(:) = Var(:) + alpha(i)*Phat(:)
496                    endif            ! use_doloop
497     
498     ! Recompute residual norm
499     ! --------------------------------
500                    if (idebugl >= 1) then
501                       call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
502     !                  Rnorm = sqrt( dot_product_par( Var, Var ) )
503     !                  print*,'leq_bicgs, initial: ', Vname,' Vnorm ', Rnorm
504     
505                       if (use_doloop) then
506     !$omp parallel do default(shared) private(ijk)
507                          do ijk=ijkstart3,ijkend3
508                             R(ijk) = B_m(ijk) - R(ijk)
509                          enddo
510                       else
511                          R(:) = B_m(:) - R(:)
512                       endif
513     
514                       if(is_serial) then
515                          if (use_doloop) then
516                             Rnorm = zero
517     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
518                             do ijk=ijkstart3,ijkend3
519                                Rnorm = Rnorm + R(ijk)*R(ijk)
520                             enddo
521                          else
522                             Rnorm =  dot_product( R, R )
523                          endif
524                          Rnorm = sqrt( Rnorm )
525                       else
526                          Rnorm = sqrt( dot_product_par( R, R ) )
527                       endif
528     !                  print*,'leq_bicgs, initial: ', Vname,' Rnorm ', Rnorm
529                    endif            ! idebugl >= 1
530                                isConverged = .TRUE.
531                                EXIT
532                 endif               ! end if (Snorm <= TOLMIN)
533              endif                  ! end if (.not.minimize_dotproducts)
534     
535     
536     
537     ! Solve A*Shat(:) = Svec(:)
538     ! Tvec(:) = A*Shat(:)
539     ! --------------------------------
540              call MSOLVE( Vname, Svec, A_m, Shat, CMETHOD)  ! returns Shat
541              call MATVEC( Vname, Shat, A_m, Tvec )   ! returns Tvec=A*Shat
542     
543              if(is_serial) then
544                 if (use_doloop) then
545                    TxS = zero
546                    TxT = zero
547     !$omp  parallel do default(shared) private(ijk) reduction(+:TxS,TxT)
548                    do ijk=ijkstart3,ijkend3
549                       TxS = TxS + Tvec(ijk)  * Svec(ijk)
550                       TxT = TxT + Tvec(ijk)  * Tvec(ijk)
551                    enddo
552                 else
553                    TxS = dot_product( Tvec, Svec )
554                    TxT = dot_product( Tvec, Tvec )
555                 endif
556              else
557                 if(.not.minimize_dotproducts) then
558                    TxS = dot_product_par( Tvec, Svec )
559                    TxT = dot_product_par( Tvec, Tvec )
560                 else
561                    TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
562                    TxS = TxS_TxT(1)
563                    TxT = TxS_TxT(2)
564                 endif
565              endif
566     
567              IF(TxT.eq.Zero) TxT = SMALL_NUMBER
568     
569     ! compute omega
570     ! --------------------------------
571              omega(i) = TxS / TxT
572     
573     ! compute new guess for Var
574     ! --------------------------------
575              if (use_doloop) then
576     !$omp parallel do default(shared) private(ijk)
577                 do ijk=ijkstart3,ijkend3
578                    Var(ijk) = Var(ijk) +  &
579                       alpha(i)*Phat(ijk) + omega(i)*Shat(ijk)
580                       R(ijk) = Svec(ijk) - omega(i)*Tvec(ijk)
581                 enddo
582              else
583                 Var(:) = Var(:) +  &
584                    alpha(i)*Phat(:) + omega(i)*Shat(:)
585                    R(:) = Svec(:) - omega(i)*Tvec(:)
586              endif
587     
588     
589     ! --------------------------------
590              if(.not.minimize_dotproducts.or.(mod(iter,icheck_bicgs).eq.0)) then
591                 if(is_serial) then
592                    if (use_doloop) then
593                       Rnorm = zero
594     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
595                       do ijk=ijkstart3,ijkend3
596                          Rnorm = Rnorm + R(ijk) * R(ijk)
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               ! is_serial
605     
606                 if (idebugl.ge.1) then
607                    if (myPE.eq.PE_IO) then
608                       print*,'iter, Rnorm ', iter, Rnorm, Snorm
609                       print*,'alpha(i), omega(i) ', alpha(i), omega(i)
610                       print*,'TxS, TxT ', TxS, TxT
611                       print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
612                    endif
613                 endif
614     
615     !           call mfix_exit(myPE)
616     
617     ! Check convergence; continue if necessary
618     ! for continuation, it is necessary that omega(i) .ne. 0
619                 isconverged = (Rnorm <= TOL*Rnorm0)
620     
621                 if (isconverged) then
622                    iter_tot(vno) = iter_tot(vno) + iter + 1
623                    EXIT
624                 endif
625              endif                  ! end if(.not.minimize_dotproducts)
626     
627     ! Advance the iteration count
628              iter = iter + 1
629     
630           enddo   ! end do i=1,itmax
631     ! end of linear solver loop
632     ! ----------------------------------------------------------------<<<
633     
634     
635           if (idebugl >= 1) then
636              call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
637              if (use_doloop) then
638     !$omp  parallel do default(shared) private(ijk)
639                 do ijk=ijkstart3,ijkend3
640                    R(ijk) = R(ijk) - B_m(ijk)
641                 enddo
642              else
643                 R(:) = R(:) - B_m(:)
644              endif
645     
646              if(is_serial) then
647                 if (use_doloop) then
648                    Rnorm = zero
649     !$omp parallel do default(shared) private(ijk) reduction(+:Rnorm)
650                    do ijk=ijkstart3,ijkend3
651                       Rnorm = Rnorm + R(ijk) * R(ijk)
652                    enddo
653                 else
654                    Rnorm = dot_product( R,R)
655                 endif
656                 Rnorm = sqrt( Rnorm )
657              else
658                 Rnorm = sqrt( dot_product_par( R,R) )
659              endif
660     
661              if(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
662     
663              if(myPE.eq.0)  print*,'leq_bicgs ratio : ', Vname,' ',iter,  &
664              ' L-2', Rnorm/Rnorm0
665           endif   ! end if(idebugl >=1)
666     
667     !      isconverged = (real(Rnorm) <= TOL*Rnorm0);
668           if(.NOT.isConverged) isconverged = (real(Rnorm) <= TOL*Rnorm0);
669     !     write(*,*) '***',iter, isconverged, Rnorm, TOL, Rnorm0, myPE
670           IER = 0
671           if (.not.isconverged) then
672              IER = -1
673              iter_tot(vno) = iter_tot(vno) + iter
674              if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
675                 IER = -2
676              endif
677           endif
678     
679           call send_recv(var,2)
680     
681           deallocate(R)
682           deallocate(Rtilde)
683           deallocate(P)
684           deallocate(Phat)
685           deallocate(Svec)
686           deallocate(Shat)
687           deallocate(Tvec)
688           deallocate(V)
689     
690           return
691           end subroutine LEQ_BICGS0
692     
693     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
694     !                                                                      C
695     !  Subroutine: LEQ_ISWEEP                                              C
696     !  Purpose: Perform line sweep at coordinate I                         C
697     !                                                                      C
698     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
699     !  Reviewer:                                          Date:            C
700     !                                                                      C
701     !  Literature/Document References:                                     C
702     !  Variables referenced:                                               C
703     !  Variables modified:                                                 C
704     !  Local variables:                                                    C
705     !                                                                      C
706     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
707     
708           SUBROUTINE LEQ_ISWEEP(I, Vname, VAR, A_M, B_M)
709     
710     !-----------------------------------------------
711     ! Modules
712     !-----------------------------------------------
713           USE param
714           USE param1
715           USE matrix
716           USE geometry
717           USE compar
718           USE indices
719           USE funits
720           USE sendrecv
721           USE mpi_utility
722           USE functions
723           IMPLICIT NONE
724     !-----------------------------------------------
725     ! Dummy arguments
726     !-----------------------------------------------
727     !  Line position
728           INTEGER, INTENT(IN) :: I
729     ! Variable name
730           CHARACTER(LEN=*), INTENT(IN) :: Vname
731     ! Variable
732           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
733     ! Septadiagonal matrix A_m
734           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
735     ! Vector b_m
736           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
737     !-----------------------------------------------
738     ! Local variables
739     !-----------------------------------------------
740           DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC, DD, EE, BB
741           INTEGER :: NSTART, NEND, INFO
742           INTEGER :: IJK, J, K, IM1JK, IP1JK
743     !-----------------------------------------------
744     
745           NEND = JEND
746           NSTART = JSTART
747           K = 1
748     
749           DO J=NSTART, NEND
750              IJK = FUNIJK(I,J,K)
751              IM1JK = IM_OF(IJK)
752              IP1JK = IP_OF(IJK)
753              DD(J) = A_M(IJK,  0)
754              CC(J) = A_M(IJK, -2)
755              EE(J) = A_M(IJK,  2)
756              BB(J) = B_M(IJK) -  A_M(IJK,-1) * Var( IM1JK )  &
757                               -  A_M(IJK, 1) * Var( IP1JK )
758           ENDDO
759     
760           CC(NSTART) = ZERO
761           EE(NEND) = ZERO
762           INFO = 0
763     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
764           CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB, JEND-JSTART+1, INFO )
765     
766           IF (INFO.NE.0) THEN
767              RETURN
768           ENDIF
769     
770           DO J=NSTART, NEND
771              IJK = FUNIJK(I,J,K)
772              Var(IJK) =  BB(J)
773           ENDDO
774     
775           RETURN
776           END SUBROUTINE LEQ_ISWEEP
777     
778     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
779     !                                                                      C
780     !  Subroutine: LEQ_IKSWEEP                                             C
781     !  Purpose: Perform line sweep at coordinate I, K                      C
782     !                                                                      C
783     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
784     !  Reviewer:                                          Date:            C
785     !                                                                      C
786     !  Literature/Document References:                                     C
787     !  Variables referenced:                                               C
788     !  Variables modified:                                                 C
789     !  Local variables:                                                    C
790     !                                                                      C
791     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
792     
793           SUBROUTINE LEQ_IKSWEEP(I, K, Vname, VAR, A_M, B_M)
794     
795     !-----------------------------------------------
796     ! Modules
797     !-----------------------------------------------
798           USE param
799           USE param1
800           USE matrix
801           USE geometry
802           USE compar
803           USE funits
804           USE indices
805           USE sendrecv
806           USE mpi_utility
807           USE functions
808           IMPLICIT NONE
809     !-----------------------------------------------
810     ! Dummy arguments
811     !-----------------------------------------------
812     ! Line position
813           INTEGER, INTENT(IN) :: I, K
814     ! Variable name
815           CHARACTER(LEN=*), INTENT(IN) :: Vname
816     ! Variable
817           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
818     ! Septadiagonal matrix A_m
819           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
820     ! Vector b_m
821           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
822     !-----------------------------------------------
823     ! Local variables
824     !-----------------------------------------------
825           DOUBLE PRECISION, DIMENSION(JSTART:JEND) :: CC, DD, EE, BB
826           INTEGER :: NSTART, NEND, INFO
827           INTEGER :: IJK, J, IM1JK, IP1JK, IJKM1, IJKP1
828     !-----------------------------------------------
829     
830           NEND = JEND
831           NSTART = JSTART
832     
833     !!$omp parallel do private(j,ijk,im1jk,ip1jk,ijkm1,ijkp1)
834           DO J=NSTART, NEND
835     !         IJK = FUNIJK(IMAP_C(I),JMAP_C(J),KMAP_C(K))
836              IJK = FUNIJK(I,J,K)
837              IM1JK = IM_OF(IJK)
838              IP1JK = IP_OF(IJK)
839              IJKM1 = KM_OF(IJK)
840              IJKP1 = KP_OF(IJK)
841              DD(J) = A_M(IJK,  0)
842              CC(J) = A_M(IJK, -2)
843              EE(J) = A_M(IJK,  2)
844              BB(J) = B_M(IJK) -  A_M(IJK,-1) * Var( IM1JK ) &
845                               -  A_M(IJK, 1) * Var( IP1JK ) &
846                               -  A_M(IJK,-3) * Var( IJKM1 ) &
847                               -  A_M(IJK, 3) * Var( IJKP1 )
848           ENDDO
849     
850           CC(NSTART) = ZERO
851           EE(NEND) = ZERO
852           INFO = 0
853     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
854           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
855     
856           IF (INFO.NE.0) THEN
857              write(*,*) 'leq_iksweep',INFO, myPE
858              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IKSWEEP'
859              RETURN
860           ENDIF
861     
862           DO J=NSTART, NEND
863              IJK = FUNIJK(I,J,K)
864              Var(IJK) = BB(J)
865           ENDDO
866     
867           RETURN
868           END SUBROUTINE LEQ_IKSWEEP
869     
870     
871     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
872     !                                                                      C
873     !  Subroutine: LEQ_JKSWEEP                                             C
874     !  Purpose: Perform line sweep at coordinate I, K                      C
875     !                                                                      C
876     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
877     !  Reviewer:                                          Date:            C
878     !                                                                      C
879     !  Literature/Document References:                                     C
880     !  Variables referenced:                                               C
881     !  Variables modified:                                                 C
882     !  Local variables:                                                    C
883     !                                                                      C
884     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
885     
886           SUBROUTINE LEQ_JKSWEEP(J, K, Vname, VAR, A_M, B_M)
887     
888     !-----------------------------------------------
889     ! Modules
890     !-----------------------------------------------
891           USE param
892           USE param1
893           USE matrix
894           USE geometry
895           USE funits
896           USE compar
897           USE indices
898           USE functions
899           IMPLICIT NONE
900     !-----------------------------------------------
901     ! Dummy arguments
902     !-----------------------------------------------
903     ! Line position
904           INTEGER, INTENT(IN) :: J, K
905     ! Variable name
906           CHARACTER(LEN=*), INTENT(IN) :: Vname
907     ! Variable
908           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
909     ! Septadiagonal matrix A_m
910           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
911     ! Vector b_m
912           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
913     !-----------------------------------------------
914     ! Local variables
915     !-----------------------------------------------
916           DOUBLE PRECISION, DIMENSION (ISTART:IEND) :: CC, DD, EE, BB
917           INTEGER :: NSTART, NEND, INFO, IJK, I
918     !-----------------------------------------------
919     
920           NEND = IEND
921           NSTART = ISTART
922     
923           DO I=NSTART,NEND
924              IJK = FUNIJK(I,J,K)
925              DD(I) = A_M(IJK,  0)
926              CC(I) = A_M(IJK, -1)
927              EE(I) = A_M(IJK,  1)
928              BB(I) = B_M(IJK)    -  A_M(IJK,-2) * Var( JM_OF(IJK) ) &
929                                  -  A_M(IJK, 2) * Var( JP_OF(IJK) ) &
930                                  -  A_M(IJK,-3) * Var( KM_OF(IJK) ) &
931                                  -  A_M(IJK, 3) * Var( KP_OF(IJK) )
932           ENDDO
933     
934           CC(NSTART) = ZERO
935           EE(NEND) = ZERO
936           INFO = 0
937           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
938     
939           IF (INFO.NE.0) THEN
940              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
941              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' JKSWEEP'
942              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
943              call mfix_exit(myPE)
944           ENDIF
945     
946           DO I=NSTART, NEND
947              IJK = FUNIJK(I,J,K)
948              Var(IJK) = BB(I)
949           ENDDO
950     
951           RETURN
952           END SUBROUTINE LEQ_JKSWEEP
953     
954     
955     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
956     !                                                                      C
957     !  Subroutine: LEQ_IJSWEEP                                             C
958     !  Purpose: Perform line sweep at coordinate I, K                      C
959     !                                                                      C
960     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
961     !  Reviewer:                                          Date:            C
962     !                                                                      C
963     !  Literature/Document References:                                     C
964     !  Variables referenced:                                               C
965     !  Variables modified:                                                 C
966     !  Local variables:                                                    C
967     !                                                                      C
968     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
969     
970           SUBROUTINE LEQ_IJSWEEP(I, J, Vname, VAR, A_M, B_M)
971     
972     !-----------------------------------------------
973     ! Modules
974     !-----------------------------------------------
975           USE param
976           USE param1
977           USE matrix
978           USE geometry
979           USE funits
980           USE compar
981           USE indices
982           USE functions
983           IMPLICIT NONE
984     !-----------------------------------------------
985     ! Dummy arguments
986     !-----------------------------------------------
987     ! Line position
988           INTEGER, INTENT(IN) :: I, J
989     ! Variable name
990           CHARACTER(LEN=*), INTENT(IN) :: Vname
991     ! Variable
992           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
993     ! Septadiagonal matrix A_m
994           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
995     ! Vector b_m
996           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
997     !-----------------------------------------------
998     ! Local variables
999     !-----------------------------------------------
1000           DOUBLE PRECISION, DIMENSION (KSTART:KEND) :: CC, DD, EE, BB
1001           INTEGER :: NEND, NSTART, INFO, IJK, K
1002     !-----------------------------------------------
1003     
1004           NEND = KEND
1005           NSTART = KSTART
1006     
1007           DO K=NSTART, NEND
1008              IJK = FUNIJK(I,J,K)
1009              DD(K) = A_M(IJK,  0)
1010              CC(K) = A_M(IJK, -3)
1011              EE(K) = A_M(IJK,  3)
1012              BB(K) = B_M(IJK)    -  A_M(IJK,-2) * Var( JM_OF(IJK) ) &
1013                                  -  A_M(IJK, 2) * Var( JP_OF(IJK) ) &
1014                                  -  A_M(IJK,-1) * Var( IM_OF(IJK) ) &
1015                                  -  A_M(IJK, 1) * Var( IP_OF(IJK) )
1016           ENDDO
1017     
1018           CC(NSTART) = ZERO
1019           EE(NEND) = ZERO
1020           INFO = 0
1021           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
1022     
1023           IF (INFO.NE.0) THEN
1024              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1025              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IJSWEEP'
1026              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1027              call mfix_exit(myPE)
1028           ENDIF
1029     
1030           DO K=NSTART, NEND
1031              IJK = FUNIJK(I,J,K)
1032              Var(IJK) = BB(K)
1033           ENDDO
1034     
1035           RETURN
1036           END SUBROUTINE LEQ_IJSWEEP
1037     
1038     
1039     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1040     !                                                                      C
1041     !  Subroutine: LEQ_MATVEC                                              C
1042     !  Purpose: Compute matrix vector multiplication                       C
1043     !           (for linear equation Ax=b compute Ax                       C
1044     !                                                                      C
1045     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1046     !  Reviewer:                                          Date:            C
1047     !                                                                      C
1048     !  Literature/Document References:                                     C
1049     !  Variables referenced:                                               C
1050     !  Variables modified:                                                 C
1051     !  Local variables:                                                    C
1052     !                                                                      C
1053     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1054     
1055           SUBROUTINE LEQ_MATVEC(VNAME, VAR, A_M, Avar)
1056     
1057     !-----------------------------------------------
1058     ! Modules
1059     !-----------------------------------------------
1060           USE param
1061           USE param1
1062           USE matrix
1063           USE geometry
1064           USE compar
1065           USE indices
1066           USE sendrecv
1067           USE mpi_utility
1068           USE cutcell
1069           USE functions
1070           IMPLICIT NONE
1071     !-----------------------------------------------
1072     ! Dummy arguments
1073     !-----------------------------------------------
1074     ! Variable name
1075           CHARACTER(LEN=*), INTENT(IN) :: Vname
1076     ! Variable
1077     !      DOUBLE PRECISION, INTENT(IN) :: Var(ijkstart3:ijkend3)
1078           DOUBLE PRECISION, INTENT(IN) :: Var(DIMENSION_3)
1079     ! Septadiagonal matrix A_m
1080     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1081           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1082     ! Vector AVar
1083     !      DOUBLE PRECISION, INTENT(OUT) :: AVar(ijkstart3:ijkend3)
1084           DOUBLE PRECISION, INTENT(OUT) :: AVar(DIMENSION_3)
1085     !-----------------------------------------------
1086     ! Local variables
1087     !-----------------------------------------------
1088     ! Variable
1089           INTEGER :: I, J, K, IJK
1090           integer :: im1jk, ip1jk, ijm1k, ijp1k, ijkm1, ijkp1
1091     !-----------------------------------------------
1092     
1093           IF(RE_INDEXING) THEN
1094     
1095              DO IJK = IJKSTART3, IJKEND3  ! Loop only over active cells
1096     
1097                 im1jk = im_of(ijk)
1098                 ip1jk = ip_of(ijk)
1099                 ijm1k = jm_of(ijk)
1100                 ijp1k = jp_of(ijk)
1101     
1102     
1103                 AVar(ijk) =      A_m(ijk,-2) * Var(ijm1k)   &
1104                                + A_m(ijk,-1) * Var(im1jk)   &
1105                                + A_m(ijk, 0) * Var(ijk)     &
1106                                + A_m(ijk, 1) * Var(ip1jk)   &
1107                                + A_m(ijk, 2) * Var(ijp1k)
1108     
1109                 if (do_k) then
1110                    ijkm1 = km_of(ijk)
1111                    ijkp1 = kp_of(ijk)
1112     
1113     
1114                     AVar(ijk) =   AVar(ijk) + A_m(ijk,-3) * Var(ijkm1)   &
1115                                             + A_m(ijk, 3) * Var(ijkp1)
1116     
1117                 endif
1118     
1119              enddo
1120     
1121     
1122           ELSE
1123     
1124     
1125     
1126              if (do_k) then
1127     !$omp    parallel  do &
1128     !$omp&   private(     &
1129     !$omp&           ijk,i,j,k, &
1130     !$omp&           im1jk,ip1jk,ijm1k,ijp1k,ijkm1,ijkp1) collapse (3)
1131                 do k = kstart,kend
1132                    do i = istart,iend
1133                       do j = jstart,jend
1134                          IJK = funijk(i,j,k)
1135                          im1jk = im_of(ijk)
1136                          ip1jk = ip_of(ijk)
1137                          ijm1k = jm_of(ijk)
1138                          ijp1k = jp_of(ijk)
1139                          ijkm1 = km_of(ijk)
1140                          ijkp1 = kp_of(ijk)
1141                          AVar(ijk) =  A_m(ijk,-3) * Var(ijkm1)   &
1142                                     + A_m(ijk,-2) * Var(ijm1k)   &
1143                                     + A_m(ijk,-1) * Var(im1jk)   &
1144                                     + A_m(ijk, 0) * Var(ijk)     &
1145                                     + A_m(ijk, 1) * Var(ip1jk)   &
1146                                     + A_m(ijk, 2) * Var(ijp1k)   &
1147                                     + A_m(ijk, 3) * Var(ijkp1)
1148                       enddo
1149                    enddo
1150                 enddo
1151     
1152              else
1153                 k = 1
1154     !$omp parallel do private(i,j,ijk,im1jk,ip1jk,ijm1k,ijp1k) collapse (2)
1155                 do i = istart,iend
1156                    do j = jstart,jend
1157                       IJK = funijk(i,j,k)
1158                       im1jk = im_of(ijk)
1159                       ip1jk = ip_of(ijk)
1160                       ijm1k = jm_of(ijk)
1161                       ijp1k = jp_of(ijk)
1162                       AVar(ijk) =  A_m(ijk,-2) * Var(ijm1k)   &
1163                                  + A_m(ijk,-1) * Var(im1jk)   &
1164                                  + A_m(ijk, 0) * Var(ijk)     &
1165                                  + A_m(ijk, 1) * Var(ip1jk)   &
1166                                  + A_m(ijk, 2) * Var(ijp1k)
1167                    enddo
1168                 enddo
1169     
1170              endif
1171     
1172           ENDIF ! RE_INDEXING
1173     
1174           call send_recv(Avar,nlayers_bicgs)
1175           RETURN
1176           END SUBROUTINE LEQ_MATVEC
1177     
1178     
1179     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1180     !                                                                      C
1181     !  Subroutine: LEQ_MSOLVE                                              C
1182     !  Purpose:                                                            C
1183     !  Notes: if leq_method is biggs or cg then this subroutine is         C
1184     !         invoked when leq_pc='line'. if leq_method is gmres then      C
1185     !         this subroutine is invoked (leq_pc setting does not matter)  C
1186     !                                                                      C
1187     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1188     !  Reviewer:                                          Date:            C
1189     !                                                                      C
1190     !  Literature/Document References:                                     C
1191     !  Variables referenced:                                               C
1192     !  Variables modified:                                                 C
1193     !  Local variables:                                                    C
1194     !                                                                      C
1195     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1196     
1197           SUBROUTINE LEQ_MSOLVE(VNAME, B_m, A_M, Var, CMETHOD)
1198     
1199     !-----------------------------------------------
1200     ! Modules
1201     !-----------------------------------------------
1202           USE param
1203           USE param1
1204           USE matrix
1205           USE geometry
1206           USE compar
1207           USE indices
1208           USE sendrecv
1209           USE functions
1210           IMPLICIT NONE
1211     !-----------------------------------------------
1212     ! Dummy arguments
1213     !-----------------------------------------------
1214     ! Variable name
1215           CHARACTER(LEN=*), INTENT(IN) :: Vname
1216     ! Vector b_m
1217     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
1218           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1219     ! Septadiagonal matrix A_m
1220     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1221           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1222     ! Variable
1223     !      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
1224           DOUBLE PRECISION, INTENT(INOUT) :: Var(DIMENSION_3)
1225     ! Sweep direction of leq solver (leq_sweep)
1226     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
1227           CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1228     !-----------------------------------------------
1229     ! Local parameters
1230     !-----------------------------------------------
1231           LOGICAL, PARAMETER :: USE_IKLOOP = .FALSE.
1232           LOGICAL, PARAMETER :: SETGUESS = .TRUE.
1233     !-----------------------------------------------
1234     ! Local variables
1235     !-----------------------------------------------
1236     !
1237           INTEGER :: ITER, NITER
1238           INTEGER :: IJK, I , J, K
1239           INTEGER :: I1, J1, K1, I2, J2, K2, IK, JK, IJ
1240           INTEGER :: ISIZE, JSIZE, KSIZE
1241           INTEGER :: ICASE
1242     
1243     !     CHARACTER(LEN=4), PARAMETER :: CMETHOD = 'II'
1244           CHARACTER :: CH
1245           LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
1246           LOGICAL :: DO_SENDRECV, DO_REDBLACK, DO_ALL
1247     
1248     !-----------------------------------------------
1249     !!$      double precision omp_start, omp_end
1250     !!$      double precision omp_get_wtime
1251     !       by Tingwen
1252     !!$      omp_start=omp_get_wtime()
1253     
1254           IF (SETGUESS) THEN
1255     !$omp   parallel do private(i,j,k,ijk)
1256              do k = kstart3,kend3
1257                 do i = istart3,iend3
1258                    do j = jstart3,jend3
1259                       IJK = funijk(i,j,k)
1260                       VAR(IJK) = B_M(IJK)
1261                    enddo
1262                 enddo
1263              enddo
1264     
1265              call send_recv(var,nlayers_bicgs)
1266           ENDIF
1267     
1268           NITER = LEN( CMETHOD )
1269     
1270           DO ITER=1,NITER
1271     
1272     ! Perform sweeps
1273              CH = CMETHOD( ITER:ITER )
1274              DO_ISWEEP = (CH .EQ. 'I') .OR. (CH .EQ. 'i')
1275              DO_JSWEEP = (CH .EQ. 'J') .OR. (CH .EQ. 'j')
1276              DO_KSWEEP = (CH .EQ. 'K') .OR. (CH .EQ. 'k')
1277              DO_ALL = (CH .EQ. 'A') .OR. (CH .EQ. 'a')
1278              DO_REDBLACK = (CH .EQ. 'R') .OR. (CH .EQ. 'r')
1279              DO_SENDRECV = (CH .EQ. 'S') .OR. (CH .EQ. 's')
1280     
1281              IF (NO_K) THEN   ! two dimensional
1282     ! 2D run no need to enable openmp parallel
1283                 IF ( DO_ISWEEP ) THEN
1284     !!$omp   parallel do private(I)
1285                    DO I=istart,iend,1
1286                       CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1287                    ENDDO
1288                 ENDIF
1289     ! ----------------------------------------------------------------<<<
1290     ! Handan Liu added 2D RSRS sweep and parallelized this loop on Jan 22 2013:
1291                         IF (DO_REDBLACK) THEN
1292     !$omp parallel do private(I)
1293                                   DO I=istart,iend,2
1294                                          CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1295                       ENDDO
1296     !$omp parallel do private(I)
1297                       DO I=istart+1,iend,2
1298                          CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1299                       ENDDO
1300                 ENDIF
1301     ! ---------------------------------------------------------------->>>
1302              ELSE   ! three dimensional
1303     
1304     
1305     ! do_all true only for leq_pc='asas'
1306     ! ---------------------------------------------------------------->>>
1307                 IF(DO_ALL) THEN        ! redblack for all sweeps, not used by default
1308     ! JK Loop
1309     ! --------------------------------
1310                    j1 = jstart
1311                    k1 = kstart
1312                    j2 = jend
1313                    k2 = kend
1314                    jsize = j2-j1+1
1315                    ksize = k2-k1+1
1316                    DO icase = 1, 2
1317     !!$omp   parallel do private(K,J,JK)
1318                       DO JK=icase, ksize*jsize, 2
1319                          if (mod(jk,jsize).ne.0) then
1320                             k = int( jk/jsize ) + k1
1321                          else
1322                             k = int( jk/jsize ) + k1 -1
1323                          endif
1324                          j = (jk-1-(k-k1)*jsize) + j1 + mod(k,2)
1325                          if(j.gt.j2) j=j-j2 + j1 -1
1326                          CALL LEQ_JKSWEEP(J, K, Vname, Var, A_m, B_m)
1327                       ENDDO
1328     
1329                    ENDDO
1330                    call send_recv(var,nlayers_bicgs)
1331     
1332     ! IJ Loop
1333     ! --------------------------------
1334                    i1 = istart
1335                    j1 = jstart
1336                    i2 = iend
1337                    j2 = jend
1338                    isize = i2-i1+1
1339                    jsize = j2-j1+1
1340                    DO icase = 1, 2
1341     !!$omp   parallel do private(J,I,IJ)
1342                       DO IJ=icase, jsize*isize, 2
1343                          if (mod(ij,isize).ne.0) then
1344                             j = int( ij/isize ) + j1
1345                          else
1346                             j = int( ij/isize ) + j1 -1
1347                          endif
1348                          i = (ij-1-(j-j1)*isize) + i1 + mod(j,2)
1349                          if(i.gt.i2) i=i-i2 + i1 -1
1350                          CALL LEQ_IJSWEEP(I, J, Vname, Var, A_m, B_m)
1351                       ENDDO
1352     
1353                    ENDDO
1354                    call send_recv(var,nlayers_bicgs)
1355     
1356     ! IK Loop
1357     ! --------------------------------
1358                    i1 = istart
1359                    k1 = kstart
1360                    i2 = iend
1361                    k2 = kend
1362                    isize = i2-i1+1
1363                    ksize = k2-k1+1
1364     
1365                    DO icase = 1, 2
1366     !!$omp   parallel do private(K,I,IK)
1367                       DO IK=icase, ksize*isize, 2
1368                          if (mod(ik,isize).ne.0) then
1369                             k = int( ik/isize ) + k1
1370                          else
1371                             k = int( ik/isize ) + k1 -1
1372                          endif
1373                          i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
1374                          if(i.gt.i2) i=i-i2 + i1 -1
1375                          CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1376                       ENDDO
1377     
1378                    ENDDO
1379                 ENDIF ! end DO_ALL
1380     ! ----------------------------------------------------------------<<<
1381     
1382     ! do_redblack only true leq_pc='rsrs'
1383     ! ---------------------------------------------------------------->>>
1384                 IF(DO_REDBLACK) THEN
1385                    !i1 = istart
1386                    !k1 = kstart
1387                    !i2 = iend
1388                    !k2 = kend
1389                    !isize = i2-i1+1
1390                    !ksize = k2-k1+1
1391     !               DO icase = 1, 2
1392     !!$omp   parallel do private(K,I,IK)
1393     !                  DO IK=icase, ksize*isize, 2
1394     !                     if (mod(ik,isize).ne.0) then
1395     !                        k = int( ik/isize ) + k1
1396     !                     else
1397     !                        k = int( ik/isize ) + k1 -1
1398     !                     endif
1399     !                     i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
1400     !                     if(i.gt.i2) i=i-i2 + i1 -1
1401     !                     CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1402     !                  ENDDO
1403     !               ENDDO
1404     !             ELSE
1405     ! Handan Liu split above loop for OpenMP at May 22 2013, modified at July 17
1406     !$omp parallel do default(shared) private(I,K) schedule(auto)
1407                    DO k=kstart,kend
1408                                  IF(mod(k,2).ne.0)THEN
1409                                         DO I=istart+1,iend,2
1410                                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1411                                         ENDDO
1412                                  ELSE
1413                                         DO I=istart,iend,2
1414                                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1415                                         ENDDO
1416                                  ENDIF
1417                            ENDDO
1418     !$omp end parallel do
1419     !$omp parallel do default(shared) private(I,K) schedule(auto)
1420                    DO k=kstart,kend
1421                                  IF(mod(k,2).ne.0)THEN
1422                                        DO I=istart,iend,2
1423                                          CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1424                                        ENDDO
1425                                  ELSE
1426                                        DO I=istart+1,iend,2
1427                                          CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1428                                        ENDDO
1429                                  ENDIF
1430                            ENDDO
1431     !$omp end parallel do
1432     
1433                 ENDIF       ! end if(do_redblack)
1434     ! ----------------------------------------------------------------<<<
1435     
1436     !  Not sure the purpose of us_ikloop
1437     !  The SMP directives below need review                        !Tingwen Jan 2012
1438                 IF(USE_IKLOOP) THEN
1439     ! use_ikloop is currently hard-wired to false (so goto else branch)
1440     ! ---------------------------------------------------------------->>>
1441                    i1 = istart
1442                    k1 = kstart
1443                    i2 = iend
1444                    k2 = kend
1445                    isize = i2-i1+1
1446                    ksize = k2-k1+1
1447                    IF (DO_ISWEEP) THEN
1448     !!$omp   parallel do private(K,I,IK)
1449                       DO IK=1, ksize*isize
1450                          if (mod(ik,isize).ne.0) then
1451                             k = int( ik/isize ) + k1
1452                          else
1453                             k = int( ik/isize ) + k1 -1
1454                          endif
1455                          i = (ik-1-(k-k1)*isize) + i1
1456                          CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1457                       ENDDO
1458                    ENDIF
1459                    IF (DO_KSWEEP) THEN
1460     !!$omp   parallel do private(K,I,IK)
1461                       DO IK=1, ksize*isize
1462                          if (mod(ik,ksize).ne.0) then
1463                             i = int( ik/ksize ) + i1
1464                          else
1465                             i = int( ik/ksize ) + i1 -1
1466                          endif
1467                          k = (ik-1-(i-i1)*ksize) + k1
1468                          CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1469                       ENDDO
1470                    ENDIF
1471     ! ----------------------------------------------------------------<<<
1472                 ELSE   ! else branch of if(use_ikloop)
1473     !  Not sure the purpose of us_ikloop
1474     !  The SMP directives below need review                        !Tingwen Jan 2012
1475     ! ---------------------------------------------------------------->>>
1476                    IF (DO_ISWEEP) THEN
1477     !!$omp   parallel do private(K,I)
1478                       DO K=kstart,kend
1479                          DO I=istart,iend
1480                             CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1481                          ENDDO
1482                       ENDDO
1483                    ENDIF
1484                    IF (DO_KSWEEP) THEN
1485     !!$omp   parallel do private(K,I)
1486                       DO I=istart,iend
1487                          DO K=kstart,kend
1488                             CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1489                          ENDDO
1490                       ENDDO
1491                    ENDIF
1492                 ENDIF   ! end if/else (use(ikloop)
1493     ! ----------------------------------------------------------------<<<
1494     
1495              ENDIF   ! end if/else (do_k)
1496     
1497     
1498     ! this is called for all settings of leq_pc
1499              IF (DO_SENDRECV) call send_recv(var,nlayers_bicgs)
1500     
1501     
1502           ENDDO   ! end do iter=1,niter
1503     !!$      omp_end=omp_get_wtime()
1504     !!$      write(*,*)'leq_msolve:',omp_end - omp_start
1505     
1506           RETURN
1507           END SUBROUTINE LEQ_MSOLVE
1508     
1509     
1510     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1511     !                                                                      C
1512     !  Subroutine: LEQ_MSOLVE0                                             C
1513     !  Notes: do nothing or no preconditioning (leq_pc='none')             C
1514     !                                                                      C
1515     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1516     !  Reviewer:                                          Date:            C
1517     !                                                                      C
1518     !  Literature/Document References:                                     C
1519     !  Variables referenced:                                               C
1520     !  Variables modified:                                                 C
1521     !  Local variables:                                                    C
1522     !                                                                      C
1523     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1524     
1525           SUBROUTINE LEQ_MSOLVE0(VNAME, B_m, A_M, Var, CMETHOD)
1526     
1527     !-----------------------------------------------
1528     ! Modules
1529     !-----------------------------------------------
1530           USE param
1531           USE param1
1532           USE matrix
1533           USE geometry
1534           USE compar
1535           USE indices
1536           USE sendrecv
1537           use parallel
1538           IMPLICIT NONE
1539     !-----------------------------------------------
1540     ! Dummy arguments
1541     !-----------------------------------------------
1542     ! Variable name
1543           CHARACTER(LEN=*), INTENT(IN) :: Vname
1544     ! Vector b_m
1545     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
1546           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1547     ! Septadiagonal matrix A_m
1548     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1549           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1550     ! Variable
1551     !      DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
1552           DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
1553     ! sweep direction
1554           CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1555     !-----------------------------------------------
1556     ! Local variables
1557     !-----------------------------------------------
1558           integer :: ijk
1559     !-----------------------------------------------
1560     
1561     ! do nothing or no preconditioning
1562           if (use_doloop) then   ! mfix.dat keyword default=false
1563     !!$omp  parallel do private(ijk)
1564              do ijk=ijkstart3,ijkend3
1565                 var(ijk) = b_m(ijk)
1566              enddo
1567           else
1568              var(:) = b_m(:)
1569           endif
1570           call send_recv(var,nlayers_bicgs)
1571     
1572           return
1573           end subroutine leq_msolve0
1574     
1575     
1576     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1577     !                                                                      C
1578     !  Subroutine: LEQ_MSOLVE1                                             C
1579     !  Notes: diagonal scaling (leq_pc='diag')                             C
1580     !                                                                      C
1581     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1582     !  Reviewer:                                          Date:            C
1583     !                                                                      C
1584     !  Literature/Document References:                                     C
1585     !  Variables referenced:                                               C
1586     !  Variables modified:                                                 C
1587     !  Local variables:                                                    C
1588     !                                                                      C
1589     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1590     
1591           SUBROUTINE LEQ_MSOLVE1(VNAME, B_m, A_M, Var, CMETHOD)
1592     
1593     !-----------------------------------------------
1594     ! Modules
1595     !-----------------------------------------------
1596           USE param
1597           USE param1
1598           USE matrix
1599           USE geometry
1600           USE compar
1601           USE indices
1602           USE sendrecv
1603           use parallel
1604     !      USE cutcell, only: RE_INDEXING,INTERIOR_CELL_AT
1605           USE cutcell
1606           USE functions
1607           IMPLICIT NONE
1608     !-----------------------------------------------
1609     ! Dummy arguments
1610     !-----------------------------------------------
1611     ! Variable name
1612           CHARACTER(LEN=*), INTENT(IN) :: Vname
1613     ! Vector b_m
1614     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
1615           DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1616     ! Septadiagonal matrix A_m
1617     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1618           DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1619     ! Variable
1620     !      DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
1621           DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
1622     ! sweep direction
1623           CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1624     !-----------------------------------------------
1625     ! Local variables
1626     !-----------------------------------------------
1627           integer :: i,j,k, ijk
1628     !-----------------------------------------------
1629     
1630           if (use_doloop) then   ! mfix.dat keyword default=false
1631     !!$omp    parallel do private(ijk)
1632              do ijk=ijkstart3,ijkend3
1633                 var(ijk) = zero
1634              enddo
1635           else
1636              var(:) = ZERO
1637           endif
1638     
1639     ! diagonal scaling
1640           IF(.NOT.RE_INDEXING) THEN
1641     !$omp   parallel do private(i,j,k,ijk)  collapse (3)
1642              do k=kstart2,kend2
1643                 do i=istart2,iend2
1644                    do j=jstart2,jend2
1645                       ijk = funijk( i,j,k )
1646                       var(ijk) = b_m(ijk)/A_m(ijk,0)
1647                    enddo
1648                 enddo
1649              enddo
1650           ELSE
1651     !$omp   parallel do private(ijk)  collapse (1)
1652              DO IJK=IJKSTART3,IJKEND3
1653                 var(ijk) = b_m(ijk)/A_m(ijk,0)
1654              ENDDO
1655           ENDIF
1656     
1657           call send_recv(var,nlayers_bicgs)
1658     
1659           return
1660           end subroutine leq_msolve1
1661     
1662     
1663     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1664     !                                                                      C
1665     !                                                                      C
1666     !                                                                      C
1667     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1668     
1669           double precision function dot_product_par(r1,r2)
1670     
1671     !-----------------------------------------------
1672     ! Modules
1673     !-----------------------------------------------
1674           use mpi_utility
1675           use geometry
1676           use compar
1677           use indices
1678           use cutcell
1679           use functions
1680           implicit none
1681     !-----------------------------------------------
1682     ! Dummy arguments
1683     !-----------------------------------------------
1684     !      double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2
1685           double precision, intent(in), dimension(DIMENSION_3) :: r1,r2
1686     !-----------------------------------------------
1687     ! Local parameters
1688     !-----------------------------------------------
1689           logical, parameter :: do_global_sum = .true.
1690     !-----------------------------------------------
1691     ! Local variables
1692     !-----------------------------------------------
1693           DOUBLE PRECISION, allocatable, Dimension(:) :: r1_g, r2_g
1694           double precision :: prod
1695           integer :: i, j, k, ijk
1696     !-----------------------------------------------
1697     
1698           if(do_global_sum) then
1699              prod = 0.0d0
1700     
1701              IF(RE_INDEXING) THEN
1702     !         IF(.FALSE.) THEN
1703                 ! Somehow, looping in this order leads to smaller time step than k,i,j nested loop below ....
1704                 DO IJK = IJKSTART3,IJKEND3
1705                    IF(INTERIOR_CELL_AT(IJK)) prod = prod + r1(ijk)*r2(ijk)
1706                 ENDDO
1707     
1708                 call global_all_sum(prod, dot_product_par)
1709     
1710     
1711              ELSE
1712     
1713     
1714     
1715     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1716                 do k = kstart1, kend1
1717                    do i = istart1, iend1
1718                       do j = jstart1, jend1
1719                          ijk = funijk_map_c (i,j,k)
1720                          prod = prod + r1(ijk)*r2(ijk)
1721                       enddo
1722                    enddo
1723                 enddo
1724     
1725                 call global_all_sum(prod, dot_product_par)
1726     
1727              ENDIF
1728     
1729           else
1730              if(myPE.eq.root) then
1731                 allocate (r1_g(1:ijkmax3))
1732                 allocate (r2_g(1:ijkmax3))
1733              else
1734                 allocate (r1_g(10))
1735                 allocate (r2_g(10))
1736              endif
1737              call gather(r1,r1_g)
1738              call gather(r2,r2_g)
1739     
1740              if(myPE.eq.root) then
1741                 prod = 0.0d0
1742     
1743     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1744                 do k = kmin1, kmax1
1745                    do i = imin1, imax1
1746                       do j = jmin1, jmax1
1747                          ijk = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1748     !                     ijk = funijk_gl (i,j,k)
1749                          prod = prod + r1_g(ijk)*r2_g(ijk)
1750                       enddo
1751                    enddo
1752                 enddo
1753     
1754              endif
1755              call bcast( prod)
1756     
1757              dot_product_par = prod
1758     
1759              deallocate (r1_g)
1760              deallocate (r2_g)
1761     
1762           endif
1763     
1764           end function dot_product_par
1765     
1766     
1767     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1768     !                                                                      C
1769     !                                                                      C
1770     !                                                                      C
1771     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1772     
1773           function dot_product_par2(r1,r2,r3,r4)
1774     
1775     !-----------------------------------------------
1776     ! Modules
1777     !-----------------------------------------------
1778           use mpi_utility
1779           use geometry
1780           use compar
1781           use indices
1782           use functions
1783           implicit none
1784     !-----------------------------------------------
1785     ! Dummy arguments
1786     !-----------------------------------------------
1787           double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2,r3,r4
1788     !-----------------------------------------------
1789     ! Local parameters
1790     !-----------------------------------------------
1791           logical, parameter :: do_global_sum = .true.
1792     !-----------------------------------------------
1793     ! Local variables
1794     !-----------------------------------------------
1795           DOUBLE PRECISION, allocatable, Dimension(:,:) :: r_temp, rg_temp
1796           double precision, Dimension(2) :: prod, dot_product_par2
1797           integer :: i, j, k, ijk
1798     !-----------------------------------------------
1799     
1800           if(do_global_sum) then
1801     
1802              prod(:) = 0.0d0
1803     
1804     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1805              do k = kstart1, kend1
1806                 do i = istart1, iend1
1807                    do j = jstart1, jend1
1808     
1809                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1810     
1811                       ijk = funijk_map_c (i,j,k)
1812                       prod(1) = prod(1) + r1(ijk)*r2(ijk)
1813                       prod(2) = prod(2) + r3(ijk)*r4(ijk)
1814                    enddo
1815                 enddo
1816              enddo
1817     
1818              call global_all_sum(prod, dot_product_par2)
1819     
1820           else
1821              allocate (r_temp(ijkstart3:ijkend3,4))
1822              r_temp(:,1) = r1
1823              r_temp(:,2) = r2
1824              r_temp(:,3) = r3
1825              r_temp(:,4) = r4
1826     
1827              if(myPE.eq.root) then
1828                 allocate (rg_temp(1:ijkmax3,4))
1829              else
1830                 allocate (rg_temp(10,4))
1831              endif
1832              call gather(r_temp,rg_temp)
1833     
1834              if(myPE.eq.root) then
1835                 prod = 0.0d0
1836     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1837                 do k = kmin1, kmax1
1838                    do i = imin1, imax1
1839                       do j = jmin1, jmax1
1840                          IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1841                          ijk = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1842     !                     ijk = funijk_gl (i,j,k)
1843                          prod(1) = prod(1) + rg_temp(ijk,1)*rg_temp(ijk,2)
1844                          prod(2) = prod(2) + rg_temp(ijk,3)*rg_temp(ijk,4)
1845                       enddo
1846                    enddo
1847                 enddo
1848              endif
1849              call bcast( prod)
1850     
1851              dot_product_par2 = prod
1852     
1853              deallocate (r_temp)
1854              deallocate (rg_temp)
1855     
1856           endif
1857     
1858           end function dot_product_par2
1859     
1860     
1861     
1862