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

1     MODULE leqsol
2     
3       use param, only: DIM_EQS
4     
5     ! Maximum number of outer iterations
6       INTEGER :: MAX_NIT
7     
8     ! Automatic adjustment of leq parameters possible (set in iterate after
9     ! the completion of first iteration).
10       LOGICAL :: LEQ_ADJUST
11     
12     ! Maximum number of linear equation solver iterations
13       INTEGER :: LEQ_IT(DIM_EQS)
14     
15     ! Linear equation solver method
16       INTEGER :: LEQ_METHOD(DIM_EQS)
17     
18     ! Total Iterations
19       INTEGER :: ITER_TOT(DIM_EQS) = 0
20     
21     ! Linear equation solver sweep direction
22       CHARACTER(LEN=4) :: LEQ_SWEEP(DIM_EQS)
23     
24     ! Linear equation solver tolerance
25       DOUBLE PRECISION :: LEQ_TOL(DIM_EQS)
26     
27     ! Preconditioner option
28       CHARACTER(LEN=4) :: LEQ_PC(DIM_EQS)
29     
30     ! Option to minimize dot products
31       LOGICAL :: MINIMIZE_DOTPRODUCTS
32     
33     ! Option to transpose A_m
34       LOGICAL :: DO_TRANSPOSE
35     
36     ! Frequency of convergence check in BiCGStab
37       INTEGER :: ICHECK_BICGS
38     
39     ! Optimize for massively parallel machine
40       LOGICAL :: OPT_PARALLEL
41     
42     ! Linear and non-linear solver statistics
43       LOGICAL :: SOLVER_STATISTICS
44     
45     CONTAINS
46     
47     !----------------------------------------------------------------------!
48     !                                                                      !
49     !                                                                      !
50     !                                                                      !
51     !----------------------------------------------------------------------!
52       SUBROUTINE REPORT_SOLVER_STATS(TNIT, STEPS)
53     
54         use error_manager
55     
56         IMPLICIT NONE
57     
58         INTEGER, INTENT(IN) :: TNIT, STEPS
59     
60         INTEGER :: LC
61     
62         WRITE(ERR_MSG,1100) iVal(TNIT), iVal(TNIT/STEPS)
63         CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
64     
65     1100 FORMAT(/2x,'Total number of non-linear iterations: ', A,/2x,&
66              'Average number per time-step: ',A)
67     
68         WRITE(ERR_MSG,1200)
69         CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
70     
71     1200 FORMAT(2x,'|',10('-'),'|',13('-'),'|',14('-'),'|',/&
72              2x,'| Equation |  Number of  |  Avg Solves  |',/&
73              2x,'|  Number  |   Solves    |   for NIT    |',/&
74              2x,'|',10('-'),'|',13('-'),'|',14('-'),'|')
75     
76         DO LC = 1, DIM_EQS
77            WRITE(ERR_MSG,1201) LC, ITER_TOT(LC), ITER_TOT(LC)/TNIT
78            CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
79         ENDDO
80     
81     1201 FORMAT(2x,'|',3x,I3,4x,'|',2x,I9,2x,'|',2x,I10,2x,'|',/ &
82              2x,'|',10('-'),'|',13('-'),'|',14('-'),'|')
83     
84     
85         RETURN
86       END SUBROUTINE REPORT_SOLVER_STATS
87     
88     
89     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
90     !                                                                      C
91     !  Subroutine: LEQ_MATVEC                                              C
92     !  Purpose: Compute matrix vector multiplication                       C
93     !           (for linear equation Ax=b compute Ax                       C
94     !                                                                      C
95     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
96     !  Reviewer:                                          Date:            C
97     !                                                                      C
98     !  Literature/Document References:                                     C
99     !  Variables referenced:                                               C
100     !  Variables modified:                                                 C
101     !  Local variables:                                                    C
102     !                                                                      C
103     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
104     
105       SUBROUTINE LEQ_MATVEC(VNAME, VAR, A_M, Avar)
106     
107     !-----------------------------------------------
108     ! Modules
109     !-----------------------------------------------
110         USE compar, ONLY: istart, iend, jstart, jend, kstart, kend, IJKSTART3, IJKEND3, nlayers_bicgs, c0, c1, c2, mype
111         USE cutcell, ONLY: re_indexing, CARTESIAN_GRID
112         USE geometry, ONLY: do_k, use_corecell_loop, CORE_ISTART, CORE_IEND, CORE_JSTART, CORE_JEND, CORE_KSTART, CORE_KEND
113         USE indices
114         USE param, ONLY: DIMENSION_3
115         USE sendrecv, ONLY: send_recv
116         IMPLICIT NONE
117     !-----------------------------------------------
118     ! Dummy arguments
119     !-----------------------------------------------
120     ! Variable name
121         CHARACTER(LEN=*), INTENT(IN) :: Vname
122     ! Variable
123     !      DOUBLE PRECISION, INTENT(IN) :: Var(ijkstart3:ijkend3)
124         DOUBLE PRECISION, INTENT(IN) :: Var(DIMENSION_3)
125     ! Septadiagonal matrix A_m
126     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
127         DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
128     ! Vector AVar
129     !      DOUBLE PRECISION, INTENT(OUT) :: AVar(ijkstart3:ijkend3)
130         DOUBLE PRECISION, INTENT(OUT) :: AVar(DIMENSION_3)
131     !-----------------------------------------------
132     ! Local variables
133     !-----------------------------------------------
134     ! Variable
135         INTEGER :: I, J, K, IJK
136         integer :: im1jk, ip1jk, ijm1k, ijp1k, ijkm1, ijkp1
137         integer :: class, interval
138         integer :: j_start(2), j_end(2)
139     !-----------------------------------------------
140     
141         IF(RE_INDEXING) THEN
142     
143            DO IJK = IJKSTART3, IJKEND3  ! Loop only over active cells
144     
145               im1jk = im_of(ijk)
146               ip1jk = ip_of(ijk)
147               ijm1k = jm_of(ijk)
148               ijp1k = jp_of(ijk)
149     
150               AVar(ijk) =      A_m(ijk,-2) * Var(ijm1k)   &
151                    + A_m(ijk,-1) * Var(im1jk)   &
152                    + A_m(ijk, 0) * Var(ijk)     &
153                    + A_m(ijk, 1) * Var(ip1jk)   &
154                    + A_m(ijk, 2) * Var(ijp1k)
155     
156               if (do_k) then
157                  ijkm1 = km_of(ijk)
158                  ijkp1 = kp_of(ijk)
159     
160     
161                  AVar(ijk) =   AVar(ijk) + A_m(ijk,-3) * Var(ijkm1)   &
162                       + A_m(ijk, 3) * Var(ijkp1)
163     
164               endif
165     
166            enddo
167     
168         ELSE
169     
170               core_istart = istart+2
171               core_iend = iend-2
172     
173               core_jstart = jstart+2
174               core_jend = jend-2
175     
176               if (do_k) then
177                  core_kstart = kstart+2
178                  core_kend = kend-2
179               else
180                  core_kstart = 1
181                  core_kend = 1
182                  kstart = 1
183                  kend = 1
184               endif
185     
186               if (USE_CORECELL_LOOP) then
187     
188               class = cell_class(funijk(core_istart,core_jstart,core_kstart))
189     
190     !$omp    parallel do default(none) shared(c0,c1,c2,avar,a_m,var,do_k,increment_for_mp,istart,jstart,kstart,iend,jend,kend,cell_class,core_istart,core_jstart,core_kstart,core_iend,core_jend,core_kend,use_corecell_loop,class) &
191     !$omp&   private(ijk,i,j,k) collapse (3)
192                  do k = core_kstart,core_kend
193                     do i = core_istart,core_iend
194                        do j = core_jstart,core_jend
195                           ijk = (j + c0 + i*c1 + k*c2)
196     
197                           AVar(ijk) = &
198                                + A_m(ijk,-2) * Var(ijk+INCREMENT_FOR_MP(3,class))   &
199                                + A_m(ijk,-1) * Var(ijk+INCREMENT_FOR_MP(1,class))   &
200                                + A_m(ijk, 0) * Var(ijk)     &
201                                + A_m(ijk, 1) * Var(ijk+INCREMENT_FOR_MP(2,class))   &
202                                + A_m(ijk, 2) * Var(ijk+INCREMENT_FOR_MP(4,class))
203     
204                           if (do_k) then
205                              AVar(ijk) =  AVar(ijk) + A_m(ijk,-3) * Var(ijk+INCREMENT_FOR_MP(5,class))
206                              AVar(ijk) =  AVar(ijk) + A_m(ijk, 3) * Var(ijk+INCREMENT_FOR_MP(6,class))
207                           endif
208                        enddo
209                     enddo
210                  enddo
211               endif
212     
213               j_start(1) = jstart
214               j_end(1) = jend
215               j_start(2) = 0 ! no iterations
216               j_end(2) = -1  ! no iterations
217     
218     !$omp    parallel do default(none) shared(c0,c1,c2,avar,a_m,var,do_k,increment_for_mp,istart,jstart,kstart,iend,jend,kend,cell_class,core_istart,core_jstart,core_kstart,core_iend,core_jend,core_kend,use_corecell_loop) &
219     !$omp&   private(ijk,i,j,k,class,interval) firstprivate(j_start,j_end) collapse (2)
220               do k = kstart,kend
221                  do i = istart,iend
222     
223                     if  (USE_CORECELL_LOOP) then
224                        if (core_istart<= i .and. i <= core_iend .and. core_kstart <= k .and. k<=core_kend) then
225                           j_start(1) = jstart
226                           j_end(1) = core_jstart-1
227                           j_start(2) = core_jend+1
228                           j_end(2) = jend
229                        else
230                           j_start(1) = jstart
231                           j_end(1) = jend
232                           j_start(2) = 0 ! no iterations
233                           j_end(2) = -1  ! no iterations
234                        endif
235                     endif
236     
237                     do interval=1,2
238                        do j = j_start(interval),j_end(interval)
239                           ijk = (j + c0 + i*c1 + k*c2)
240                           class = cell_class(ijk)
241     
242                           AVar(ijk) = &
243                                + A_m(ijk,-2) * Var(ijk+INCREMENT_FOR_MP(3,class))   &
244                                + A_m(ijk,-1) * Var(ijk+INCREMENT_FOR_MP(1,class))   &
245                                + A_m(ijk, 0) * Var(ijk)     &
246                                + A_m(ijk, 1) * Var(ijk+INCREMENT_FOR_MP(2,class))   &
247                                + A_m(ijk, 2) * Var(ijk+INCREMENT_FOR_MP(4,class))
248     
249                           if (do_k) then
250                              AVar(ijk) =  AVar(ijk) + A_m(ijk,-3) * Var(ijk+INCREMENT_FOR_MP(5,class))
251                              AVar(ijk) =  AVar(ijk) + A_m(ijk, 3) * Var(ijk+INCREMENT_FOR_MP(6,class))
252                           endif
253                        enddo
254                     enddo
255                  enddo
256               enddo
257     
258            ENDIF ! RE_INDEXING
259     
260            call send_recv(Avar,nlayers_bicgs)
261         RETURN
262     
263       CONTAINS
264     
265         INCLUDE 'functions.inc'
266     
267       END SUBROUTINE LEQ_MATVEC
268     
269     
270     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
271     !                                                                      C
272     !  Subroutine: LEQ_MSOLVE                                              C
273     !  Purpose:                                                            C
274     !  Notes: if leq_method is biggs or cg then this subroutine is         C
275     !         invoked when leq_pc='line'. if leq_method is gmres then      C
276     !         this subroutine is invoked (leq_pc setting does not matter)  C
277     !                                                                      C
278     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
279     !  Reviewer:                                          Date:            C
280     !                                                                      C
281     !  Literature/Document References:                                     C
282     !  Variables referenced:                                               C
283     !  Variables modified:                                                 C
284     !  Local variables:                                                    C
285     !                                                                      C
286     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
287     
288       SUBROUTINE LEQ_MSOLVE(VNAME, B_m, A_M, Var, CMETHOD)
289     
290     !-----------------------------------------------
291     ! Modules
292     !-----------------------------------------------
293         USE param
294         USE param1
295         USE matrix
296         USE geometry
297         USE compar
298         USE indices
299         USE sendrecv
300         USE functions
301         IMPLICIT NONE
302     !-----------------------------------------------
303     ! Dummy arguments
304     !-----------------------------------------------
305     ! Variable name
306         CHARACTER(LEN=*), INTENT(IN) :: Vname
307     ! Vector b_m
308     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
309         DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
310     ! Septadiagonal matrix A_m
311     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
312         DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
313     ! Variable
314     !      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
315         DOUBLE PRECISION, INTENT(INOUT) :: Var(DIMENSION_3)
316     ! Sweep direction of leq solver (leq_sweep)
317     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
318         CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
319     !-----------------------------------------------
320     ! Local parameters
321     !-----------------------------------------------
322         LOGICAL, PARAMETER :: USE_IKLOOP = .FALSE.
323         LOGICAL, PARAMETER :: SETGUESS = .TRUE.
324     !-----------------------------------------------
325     ! Local variables
326     !-----------------------------------------------
327     !
328         INTEGER :: ITER, NITER
329         INTEGER :: IJK, I , J, K
330         INTEGER :: I1, J1, K1, I2, J2, K2, IK, JK, IJ
331         INTEGER :: ISIZE, JSIZE, KSIZE
332         INTEGER :: ICASE
333     
334     !     CHARACTER(LEN=4), PARAMETER :: CMETHOD = 'II'
335         CHARACTER :: CH
336         LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
337         LOGICAL :: DO_SENDRECV, DO_REDBLACK, DO_ALL
338     
339     !-----------------------------------------------
340     !!$      double precision omp_start, omp_end
341     !!$      double precision omp_get_wtime
342     !       by Tingwen
343     !!$      omp_start=omp_get_wtime()
344     
345         IF (SETGUESS) THEN
346     !$omp   parallel do private(i,j,k,ijk)
347            do k = kstart3,kend3
348               do i = istart3,iend3
349                  do j = jstart3,jend3
350                     IJK = (J + C0 + I*C1 + K*C2)
351                     VAR(IJK) = B_M(IJK)
352                  enddo
353               enddo
354            enddo
355     
356            call send_recv(var,nlayers_bicgs)
357         ENDIF
358     
359         NITER = LEN( CMETHOD )
360     
361         DO ITER=1,NITER
362     
363     ! Perform sweeps
364            CH = CMETHOD( ITER:ITER )
365            DO_ISWEEP = (CH .EQ. 'I') .OR. (CH .EQ. 'i')
366            DO_JSWEEP = (CH .EQ. 'J') .OR. (CH .EQ. 'j')
367            DO_KSWEEP = (CH .EQ. 'K') .OR. (CH .EQ. 'k')
368            DO_ALL = (CH .EQ. 'A') .OR. (CH .EQ. 'a')
369            DO_REDBLACK = (CH .EQ. 'R') .OR. (CH .EQ. 'r')
370            DO_SENDRECV = (CH .EQ. 'S') .OR. (CH .EQ. 's')
371     
372            IF (NO_K) THEN   ! two dimensional
373     ! 2D run no need to enable openmp parallel
374               IF ( DO_ISWEEP ) THEN
375     !!$omp   parallel do private(I)
376                  DO I=istart,iend,1
377                     CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
378                  ENDDO
379               ENDIF
380     ! ----------------------------------------------------------------<<<
381     ! Handan Liu added 2D RSRS sweep and parallelized this loop on Jan 22 2013:
382               IF (DO_REDBLACK) THEN
383     !$omp parallel do private(I)
384                  DO I=istart,iend,2
385                     CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
386                  ENDDO
387     !$omp parallel do private(I)
388                  DO I=istart+1,iend,2
389                     CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
390                  ENDDO
391               ENDIF
392     ! ---------------------------------------------------------------->>>
393            ELSE   ! three dimensional
394     
395     
396     ! do_all true only for leq_pc='asas'
397     ! ---------------------------------------------------------------->>>
398               IF(DO_ALL) THEN        ! redblack for all sweeps, not used by default
399     ! JK Loop
400     ! --------------------------------
401                  j1 = jstart
402                  k1 = kstart
403                  j2 = jend
404                  k2 = kend
405                  jsize = j2-j1+1
406                  ksize = k2-k1+1
407                  DO icase = 1, 2
408     !!$omp   parallel do private(K,J,JK)
409                     DO JK=icase, ksize*jsize, 2
410                        if (mod(jk,jsize).ne.0) then
411                           k = int( jk/jsize ) + k1
412                        else
413                           k = int( jk/jsize ) + k1 -1
414                        endif
415                        j = (jk-1-(k-k1)*jsize) + j1 + mod(k,2)
416                        if(j.gt.j2) j=j-j2 + j1 -1
417                        CALL LEQ_JKSWEEP(J, K, Vname, Var, A_m, B_m)
418                     ENDDO
419     
420                  ENDDO
421                  call send_recv(var,nlayers_bicgs)
422     
423     ! IJ Loop
424     ! --------------------------------
425                  i1 = istart
426                  j1 = jstart
427                  i2 = iend
428                  j2 = jend
429                  isize = i2-i1+1
430                  jsize = j2-j1+1
431                  DO icase = 1, 2
432     !!$omp   parallel do private(J,I,IJ)
433                     DO IJ=icase, jsize*isize, 2
434                        if (mod(ij,isize).ne.0) then
435                           j = int( ij/isize ) + j1
436                        else
437                           j = int( ij/isize ) + j1 -1
438                        endif
439                        i = (ij-1-(j-j1)*isize) + i1 + mod(j,2)
440                        if(i.gt.i2) i=i-i2 + i1 -1
441                        CALL LEQ_IJSWEEP(I, J, Vname, Var, A_m, B_m)
442                     ENDDO
443     
444                  ENDDO
445                  call send_recv(var,nlayers_bicgs)
446     
447     ! IK Loop
448     ! --------------------------------
449                  i1 = istart
450                  k1 = kstart
451                  i2 = iend
452                  k2 = kend
453                  isize = i2-i1+1
454                  ksize = k2-k1+1
455     
456                  DO icase = 1, 2
457     !!$omp   parallel do private(K,I,IK)
458                     DO IK=icase, ksize*isize, 2
459                        if (mod(ik,isize).ne.0) then
460                           k = int( ik/isize ) + k1
461                        else
462                           k = int( ik/isize ) + k1 -1
463                        endif
464                        i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
465                        if(i.gt.i2) i=i-i2 + i1 -1
466                        CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
467                     ENDDO
468     
469                  ENDDO
470               ENDIF ! end DO_ALL
471     ! ----------------------------------------------------------------<<<
472     
473     ! do_redblack only true leq_pc='rsrs'
474     ! ---------------------------------------------------------------->>>
475               IF(DO_REDBLACK) THEN
476     !i1 = istart
477     !k1 = kstart
478     !i2 = iend
479     !k2 = kend
480     !isize = i2-i1+1
481     !ksize = k2-k1+1
482     !               DO icase = 1, 2
483     !!$omp   parallel do private(K,I,IK)
484     !                  DO IK=icase, ksize*isize, 2
485     !                     if (mod(ik,isize).ne.0) then
486     !                        k = int( ik/isize ) + k1
487     !                     else
488     !                        k = int( ik/isize ) + k1 -1
489     !                     endif
490     !                     i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
491     !                     if(i.gt.i2) i=i-i2 + i1 -1
492     !                     CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
493     !                  ENDDO
494     !               ENDDO
495     !             ELSE
496     ! Handan Liu split above loop for OpenMP at May 22 2013, modified at July 17
497     !$omp parallel do default(shared) private(I,K) schedule(auto)
498                  DO k=kstart,kend
499                     IF(mod(k,2).ne.0)THEN
500                        DO I=istart+1,iend,2
501                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
502                        ENDDO
503                     ELSE
504                        DO I=istart,iend,2
505                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
506                        ENDDO
507                     ENDIF
508                  ENDDO
509     !$omp end parallel do
510     !$omp parallel do default(shared) private(I,K) schedule(auto)
511                  DO k=kstart,kend
512                     IF(mod(k,2).ne.0)THEN
513                        DO I=istart,iend,2
514                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
515                        ENDDO
516                     ELSE
517                        DO I=istart+1,iend,2
518                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
519                        ENDDO
520                     ENDIF
521                  ENDDO
522     !$omp end parallel do
523     
524               ENDIF       ! end if(do_redblack)
525     ! ----------------------------------------------------------------<<<
526     
527     !  Not sure the purpose of us_ikloop
528     !  The SMP directives below need review                        !Tingwen Jan 2012
529               IF(USE_IKLOOP) THEN
530     ! use_ikloop is currently hard-wired to false (so goto else branch)
531     ! ---------------------------------------------------------------->>>
532                  i1 = istart
533                  k1 = kstart
534                  i2 = iend
535                  k2 = kend
536                  isize = i2-i1+1
537                  ksize = k2-k1+1
538                  IF (DO_ISWEEP) THEN
539     !!$omp   parallel do private(K,I,IK)
540                     DO IK=1, ksize*isize
541                        if (mod(ik,isize).ne.0) then
542                           k = int( ik/isize ) + k1
543                        else
544                           k = int( ik/isize ) + k1 -1
545                        endif
546                        i = (ik-1-(k-k1)*isize) + i1
547                        CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
548                     ENDDO
549                  ENDIF
550                  IF (DO_KSWEEP) THEN
551     !!$omp   parallel do private(K,I,IK)
552                     DO IK=1, ksize*isize
553                        if (mod(ik,ksize).ne.0) then
554                           i = int( ik/ksize ) + i1
555                        else
556                           i = int( ik/ksize ) + i1 -1
557                        endif
558                        k = (ik-1-(i-i1)*ksize) + k1
559                        CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
560                     ENDDO
561                  ENDIF
562     ! ----------------------------------------------------------------<<<
563               ELSE   ! else branch of if(use_ikloop)
564     !  Not sure the purpose of us_ikloop
565     !  The SMP directives below need review                        !Tingwen Jan 2012
566     ! ---------------------------------------------------------------->>>
567                  IF (DO_ISWEEP) THEN
568     !!$omp   parallel do private(K,I)
569                     DO K=kstart,kend
570                        DO I=istart,iend
571                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
572                        ENDDO
573                     ENDDO
574                  ENDIF
575                  IF (DO_KSWEEP) THEN
576     !!$omp   parallel do private(K,I)
577                     DO I=istart,iend
578                        DO K=kstart,kend
579                           CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
580                        ENDDO
581                     ENDDO
582                  ENDIF
583               ENDIF   ! end if/else (use(ikloop)
584     ! ----------------------------------------------------------------<<<
585     
586            ENDIF   ! end if/else (do_k)
587     
588     
589     ! this is called for all settings of leq_pc
590            IF (DO_SENDRECV) call send_recv(var,nlayers_bicgs)
591     
592     
593         ENDDO   ! end do iter=1,niter
594     !!$      omp_end=omp_get_wtime()
595     !!$      write(*,*)'leq_msolve:',omp_end - omp_start
596     
597         RETURN
598       END SUBROUTINE LEQ_MSOLVE
599     
600     
601     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
602     !                                                                      C
603     !  Subroutine: LEQ_MSOLVE0                                             C
604     !  Notes: do nothing or no preconditioning (leq_pc='none')             C
605     !                                                                      C
606     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
607     !  Reviewer:                                          Date:            C
608     !                                                                      C
609     !  Literature/Document References:                                     C
610     !  Variables referenced:                                               C
611     !  Variables modified:                                                 C
612     !  Local variables:                                                    C
613     !                                                                      C
614     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
615     
616       SUBROUTINE LEQ_MSOLVE0(VNAME, B_m, A_M, Var, CMETHOD)
617     
618     !-----------------------------------------------
619     ! Modules
620     !-----------------------------------------------
621         USE param
622         USE param1
623         USE matrix
624         USE geometry
625         USE compar
626         USE indices
627         USE sendrecv
628         use parallel
629         IMPLICIT NONE
630     !-----------------------------------------------
631     ! Dummy arguments
632     !-----------------------------------------------
633     ! Variable name
634         CHARACTER(LEN=*), INTENT(IN) :: Vname
635     ! Vector b_m
636     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
637         DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
638     ! Septadiagonal matrix A_m
639     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
640         DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
641     ! Variable
642     !      DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
643         DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
644     ! sweep direction
645         CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
646     !-----------------------------------------------
647     ! Local variables
648     !-----------------------------------------------
649         integer :: ijk
650     !-----------------------------------------------
651     
652     ! do nothing or no preconditioning
653         if (use_doloop) then   ! mfix.dat keyword default=false
654     !!$omp  parallel do private(ijk)
655            do ijk=ijkstart3,ijkend3
656               var(ijk) = b_m(ijk)
657            enddo
658         else
659            var(:) = b_m(:)
660         endif
661         call send_recv(var,nlayers_bicgs)
662     
663         return
664       end subroutine leq_msolve0
665     
666     
667     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
668     !                                                                      C
669     !  Subroutine: LEQ_MSOLVE1                                             C
670     !  Notes: diagonal scaling (leq_pc='diag')                             C
671     !                                                                      C
672     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
673     !  Reviewer:                                          Date:            C
674     !                                                                      C
675     !  Literature/Document References:                                     C
676     !  Variables referenced:                                               C
677     !  Variables modified:                                                 C
678     !  Local variables:                                                    C
679     !                                                                      C
680     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
681     
682       SUBROUTINE LEQ_MSOLVE1(VNAME, B_m, A_M, Var, CMETHOD)
683     
684     !-----------------------------------------------
685     ! Modules
686     !-----------------------------------------------
687         USE param
688         USE param1
689         USE matrix
690         USE geometry
691         USE compar
692         USE indices
693         USE sendrecv
694         use parallel
695     !      USE cutcell, only: RE_INDEXING,INTERIOR_CELL_AT
696         USE cutcell
697         USE functions
698         IMPLICIT NONE
699     !-----------------------------------------------
700     ! Dummy arguments
701     !-----------------------------------------------
702     ! Variable name
703         CHARACTER(LEN=*), INTENT(IN) :: Vname
704     ! Vector b_m
705     !      DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
706         DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
707     ! Septadiagonal matrix A_m
708     !      DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
709         DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
710     ! Variable
711     !      DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
712         DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
713     ! sweep direction
714         CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
715     !-----------------------------------------------
716     ! Local variables
717     !-----------------------------------------------
718         integer :: i,j,k, ijk
719     !-----------------------------------------------
720     
721         if (use_doloop) then   ! mfix.dat keyword default=false
722     !!$omp    parallel do private(ijk)
723            do ijk=ijkstart3,ijkend3
724               var(ijk) = zero
725            enddo
726         else
727            var(:) = ZERO
728         endif
729     
730     ! diagonal scaling
731         IF(.NOT.RE_INDEXING) THEN
732     !$omp   parallel do private(i,j,k,ijk)  collapse (3)
733            do k=kstart2,kend2
734               do i=istart2,iend2
735                  do j=jstart2,jend2
736                     ijk = funijk( i,j,k )
737                     var(ijk) = b_m(ijk)/A_m(ijk,0)
738                  enddo
739               enddo
740            enddo
741         ELSE
742     !$omp   parallel do private(ijk)  collapse (1)
743            DO IJK=IJKSTART3,IJKEND3
744               var(ijk) = b_m(ijk)/A_m(ijk,0)
745            ENDDO
746         ENDIF
747     
748         call send_recv(var,nlayers_bicgs)
749     
750         return
751       end subroutine leq_msolve1
752     
753     
754     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
755     !                                                                      C
756     !  Subroutine: LEQ_ISWEEP                                              C
757     !  Purpose: Perform line sweep at coordinate I                         C
758     !                                                                      C
759     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
760     !  Reviewer:                                          Date:            C
761     !                                                                      C
762     !  Literature/Document References:                                     C
763     !  Variables referenced:                                               C
764     !  Variables modified:                                                 C
765     !  Local variables:                                                    C
766     !                                                                      C
767     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
768     
769           SUBROUTINE LEQ_ISWEEP(I, Vname, VAR, A_M, B_M)
770     
771     !-----------------------------------------------
772     ! Modules
773     !-----------------------------------------------
774           USE param
775           USE param1
776           USE matrix
777           USE geometry
778           USE compar
779           USE indices
780           USE funits
781           USE sendrecv
782           USE mpi_utility
783           USE functions
784           IMPLICIT NONE
785     !-----------------------------------------------
786     ! Dummy arguments
787     !-----------------------------------------------
788     !  Line position
789           INTEGER, INTENT(IN) :: I
790     ! Variable name
791           CHARACTER(LEN=*), INTENT(IN) :: Vname
792     ! Variable
793           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
794     ! Septadiagonal matrix A_m
795           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
796     ! Vector b_m
797           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
798     !-----------------------------------------------
799     ! Local variables
800     !-----------------------------------------------
801           DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC, DD, EE, BB
802           INTEGER :: NSTART, NEND, INFO
803           INTEGER :: IJK, J, K, IM1JK, IP1JK
804     !-----------------------------------------------
805     
806           NEND = JEND
807           NSTART = JSTART
808           K = 1
809     
810           DO J=NSTART, NEND
811              IJK = FUNIJK(I,J,K)
812              IM1JK = IM_OF(IJK)
813              IP1JK = IP_OF(IJK)
814              DD(J) = A_M(IJK,  0)
815              CC(J) = A_M(IJK, -2)
816              EE(J) = A_M(IJK,  2)
817              BB(J) = B_M(IJK) -  A_M(IJK,-1) * Var( IM1JK )  &
818                               -  A_M(IJK, 1) * Var( IP1JK )
819           ENDDO
820     
821           CC(NSTART) = ZERO
822           EE(NEND) = ZERO
823           INFO = 0
824     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
825           CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB, JEND-JSTART+1, INFO )
826     
827           IF (INFO.NE.0) THEN
828              RETURN
829           ENDIF
830     
831           DO J=NSTART, NEND
832              IJK = FUNIJK(I,J,K)
833              Var(IJK) =  BB(J)
834           ENDDO
835     
836           RETURN
837           END SUBROUTINE LEQ_ISWEEP
838     
839     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
840     !                                                                      C
841     !  Subroutine: LEQ_IKSWEEP                                             C
842     !  Purpose: Perform line sweep at coordinate I, K                      C
843     !                                                                      C
844     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
845     !  Reviewer:                                          Date:            C
846     !                                                                      C
847     !  Literature/Document References:                                     C
848     !  Variables referenced:                                               C
849     !  Variables modified:                                                 C
850     !  Local variables:                                                    C
851     !                                                                      C
852     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
853     
854           SUBROUTINE LEQ_IKSWEEP(I, K, Vname, VAR, A_M, B_M)
855     
856     !-----------------------------------------------
857     ! Modules
858     !-----------------------------------------------
859           USE param
860           USE param1
861           USE matrix
862           USE geometry
863           USE compar
864           USE funits
865           USE indices
866           USE sendrecv
867           USE mpi_utility
868           USE functions
869           IMPLICIT NONE
870     !-----------------------------------------------
871     ! Dummy arguments
872     !-----------------------------------------------
873     ! Line position
874           INTEGER, INTENT(IN) :: I, K
875     ! Variable name
876           CHARACTER(LEN=*), INTENT(IN) :: Vname
877     ! Variable
878           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
879     ! Septadiagonal matrix A_m
880           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
881     ! Vector b_m
882           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
883     !-----------------------------------------------
884     ! Local variables
885     !-----------------------------------------------
886           DOUBLE PRECISION, DIMENSION(JSTART:JEND) :: CC, DD, EE, BB
887           INTEGER :: NSTART, NEND, INFO
888           INTEGER :: IJK, J, CLASS
889     !-----------------------------------------------
890     
891           NEND = JEND
892           NSTART = JSTART
893     
894     !!$omp parallel do private(j,ijk,im1jk,ip1jk,ijkm1,ijkp1)
895           DO J=NSTART, NEND
896     !         IJK = FUNIJK(IMAP_C(I),JMAP_C(J),KMAP_C(K))
897              IJK = (J + C0 + I*C1 + K*C2)
898              CLASS = CELL_CLASS(IJK)
899              DD(J) = A_M(IJK,  0)
900              CC(J) = A_M(IJK, -2)
901              EE(J) = A_M(IJK,  2)
902              BB(J) = B_M(IJK) -  A_M(IJK,-1) * Var( IJK+INCREMENT_FOR_MP(1,class) ) &
903                               -  A_M(IJK, 1) * Var( IJK+INCREMENT_FOR_MP(2,class) ) &
904                               -  A_M(IJK,-3) * Var( IJK+INCREMENT_FOR_MP(5,class) ) &
905                               -  A_M(IJK, 3) * Var( IJK+INCREMENT_FOR_MP(6,class) )
906           ENDDO
907     
908           CC(NSTART) = ZERO
909           EE(NEND) = ZERO
910           INFO = 0
911     !     CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
912           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
913     
914           IF (INFO.NE.0) THEN
915              write(*,*) 'leq_iksweep',INFO, myPE
916              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IKSWEEP'
917              RETURN
918           ENDIF
919     
920           DO J=NSTART, NEND
921              Var(J + C0 + I*C1 + K*C2) = BB(J)
922           ENDDO
923     
924           RETURN
925           END SUBROUTINE LEQ_IKSWEEP
926     
927     
928     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
929     !                                                                      C
930     !  Subroutine: LEQ_JKSWEEP                                             C
931     !  Purpose: Perform line sweep at coordinate I, K                      C
932     !                                                                      C
933     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
934     !  Reviewer:                                          Date:            C
935     !                                                                      C
936     !  Literature/Document References:                                     C
937     !  Variables referenced:                                               C
938     !  Variables modified:                                                 C
939     !  Local variables:                                                    C
940     !                                                                      C
941     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
942     
943           SUBROUTINE LEQ_JKSWEEP(J, K, Vname, VAR, A_M, B_M)
944     
945     !-----------------------------------------------
946     ! Modules
947     !-----------------------------------------------
948           USE param
949           USE param1
950           USE matrix
951           USE geometry
952           USE funits
953           USE compar
954           USE indices
955           USE functions
956           IMPLICIT NONE
957     !-----------------------------------------------
958     ! Dummy arguments
959     !-----------------------------------------------
960     ! Line position
961           INTEGER, INTENT(IN) :: J, K
962     ! Variable name
963           CHARACTER(LEN=*), INTENT(IN) :: Vname
964     ! Variable
965           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
966     ! Septadiagonal matrix A_m
967           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
968     ! Vector b_m
969           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
970     !-----------------------------------------------
971     ! Local variables
972     !-----------------------------------------------
973           DOUBLE PRECISION, DIMENSION (ISTART:IEND) :: CC, DD, EE, BB
974           INTEGER :: NSTART, NEND, INFO, IJK, I
975     !-----------------------------------------------
976     
977           NEND = IEND
978           NSTART = ISTART
979     
980           DO I=NSTART,NEND
981              IJK = FUNIJK(I,J,K)
982              DD(I) = A_M(IJK,  0)
983              CC(I) = A_M(IJK, -1)
984              EE(I) = A_M(IJK,  1)
985              BB(I) = B_M(IJK)    -  A_M(IJK,-2) * Var( JM_OF(IJK) ) &
986                                  -  A_M(IJK, 2) * Var( JP_OF(IJK) ) &
987                                  -  A_M(IJK,-3) * Var( KM_OF(IJK) ) &
988                                  -  A_M(IJK, 3) * Var( KP_OF(IJK) )
989           ENDDO
990     
991           CC(NSTART) = ZERO
992           EE(NEND) = ZERO
993           INFO = 0
994           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
995     
996           IF (INFO.NE.0) THEN
997              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
998              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' JKSWEEP'
999              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1000              call mfix_exit(myPE)
1001           ENDIF
1002     
1003           DO I=NSTART, NEND
1004              IJK = FUNIJK(I,J,K)
1005              Var(IJK) = BB(I)
1006           ENDDO
1007     
1008           RETURN
1009           END SUBROUTINE LEQ_JKSWEEP
1010     
1011     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1012     !                                                                      C
1013     !  Subroutine: LEQ_IJSWEEP                                             C
1014     !  Purpose: Perform line sweep at coordinate I, K                      C
1015     !                                                                      C
1016     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
1017     !  Reviewer:                                          Date:            C
1018     !                                                                      C
1019     !  Literature/Document References:                                     C
1020     !  Variables referenced:                                               C
1021     !  Variables modified:                                                 C
1022     !  Local variables:                                                    C
1023     !                                                                      C
1024     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1025     
1026           SUBROUTINE LEQ_IJSWEEP(I, J, Vname, VAR, A_M, B_M)
1027     
1028     !-----------------------------------------------
1029     ! Modules
1030     !-----------------------------------------------
1031           USE param
1032           USE param1
1033           USE matrix
1034           USE geometry
1035           USE funits
1036           USE compar
1037           USE indices
1038           USE functions
1039           IMPLICIT NONE
1040     !-----------------------------------------------
1041     ! Dummy arguments
1042     !-----------------------------------------------
1043     ! Line position
1044           INTEGER, INTENT(IN) :: I, J
1045     ! Variable name
1046           CHARACTER(LEN=*), INTENT(IN) :: Vname
1047     ! Variable
1048           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
1049     ! Septadiagonal matrix A_m
1050           DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1051     ! Vector b_m
1052           DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
1053     !-----------------------------------------------
1054     ! Local variables
1055     !-----------------------------------------------
1056           DOUBLE PRECISION, DIMENSION (KSTART:KEND) :: CC, DD, EE, BB
1057           INTEGER :: NEND, NSTART, INFO, IJK, K
1058     !-----------------------------------------------
1059     
1060           NEND = KEND
1061           NSTART = KSTART
1062     
1063           DO K=NSTART, NEND
1064              IJK = FUNIJK(I,J,K)
1065              DD(K) = A_M(IJK,  0)
1066              CC(K) = A_M(IJK, -3)
1067              EE(K) = A_M(IJK,  3)
1068              BB(K) = B_M(IJK)    -  A_M(IJK,-2) * Var( JM_OF(IJK) ) &
1069                                  -  A_M(IJK, 2) * Var( JP_OF(IJK) ) &
1070                                  -  A_M(IJK,-1) * Var( IM_OF(IJK) ) &
1071                                  -  A_M(IJK, 1) * Var( IP_OF(IJK) )
1072           ENDDO
1073     
1074           CC(NSTART) = ZERO
1075           EE(NEND) = ZERO
1076           INFO = 0
1077           CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
1078     
1079           IF (INFO.NE.0) THEN
1080              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1081              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IJSWEEP'
1082              IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1083              call mfix_exit(myPE)
1084           ENDIF
1085     
1086           DO K=NSTART, NEND
1087              IJK = FUNIJK(I,J,K)
1088              Var(IJK) = BB(K)
1089           ENDDO
1090     
1091           RETURN
1092           END SUBROUTINE LEQ_IJSWEEP
1093     
1094     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1095     !                                                                      C
1096     !                                                                      C
1097     !                                                                      C
1098     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1099     
1100       double precision function dot_product_par(r1,r2)
1101     
1102     !-----------------------------------------------
1103     ! Modules
1104     !-----------------------------------------------
1105         use mpi_utility
1106         use geometry
1107         use compar
1108         use indices
1109         use cutcell
1110         use functions
1111         implicit none
1112     !-----------------------------------------------
1113     ! Dummy arguments
1114     !-----------------------------------------------
1115     !      double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2
1116         double precision, intent(in), dimension(DIMENSION_3) :: r1,r2
1117     !-----------------------------------------------
1118     ! Local parameters
1119     !-----------------------------------------------
1120         logical, parameter :: do_global_sum = .true.
1121     !-----------------------------------------------
1122     ! Local variables
1123     !-----------------------------------------------
1124         DOUBLE PRECISION, allocatable, Dimension(:) :: r1_g, r2_g
1125         double precision :: prod
1126         integer :: i, j, k, ijk
1127     !-----------------------------------------------
1128     
1129         if(do_global_sum) then
1130            prod = 0.0d0
1131     
1132            IF(RE_INDEXING) THEN
1133     !         IF(.FALSE.) THEN
1134     ! Somehow, looping in this order leads to smaller time step than k,i,j nested loop below ....
1135               DO IJK = IJKSTART3,IJKEND3
1136                  IF(INTERIOR_CELL_AT(IJK)) prod = prod + r1(ijk)*r2(ijk)
1137               ENDDO
1138     
1139               call global_all_sum(prod, dot_product_par)
1140     
1141     
1142            ELSE
1143     
1144     
1145     
1146     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1147               do k = kstart1, kend1
1148                  do i = istart1, iend1
1149                     do j = jstart1, jend1
1150                        ijk = funijk_map_c (i,j,k)
1151                        prod = prod + r1(ijk)*r2(ijk)
1152                     enddo
1153                  enddo
1154               enddo
1155     
1156               call global_all_sum(prod, dot_product_par)
1157     
1158            ENDIF
1159     
1160         else
1161            if(myPE.eq.root) then
1162               allocate (r1_g(1:ijkmax3))
1163               allocate (r2_g(1:ijkmax3))
1164            else
1165               allocate (r1_g(10))
1166               allocate (r2_g(10))
1167            endif
1168            call gather(r1,r1_g)
1169            call gather(r2,r2_g)
1170     
1171            if(myPE.eq.root) then
1172               prod = 0.0d0
1173     
1174     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1175               do k = kmin1, kmax1
1176                  do i = imin1, imax1
1177                     do j = jmin1, jmax1
1178                        ijk = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1179     !                     ijk = funijk_gl (i,j,k)
1180                        prod = prod + r1_g(ijk)*r2_g(ijk)
1181                     enddo
1182                  enddo
1183               enddo
1184     
1185            endif
1186            call bcast( prod)
1187     
1188            dot_product_par = prod
1189     
1190            deallocate (r1_g)
1191            deallocate (r2_g)
1192     
1193         endif
1194     
1195       end function dot_product_par
1196     
1197     
1198     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1199     !                                                                      C
1200     !                                                                      C
1201     !                                                                      C
1202     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1203     
1204       function dot_product_par2(r1,r2,r3,r4)
1205     
1206     !-----------------------------------------------
1207     ! Modules
1208     !-----------------------------------------------
1209         use mpi_utility
1210         use geometry
1211         use compar
1212         use indices
1213         use functions
1214         implicit none
1215     !-----------------------------------------------
1216     ! Dummy arguments
1217     !-----------------------------------------------
1218         double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2,r3,r4
1219     !-----------------------------------------------
1220     ! Local parameters
1221     !-----------------------------------------------
1222         logical, parameter :: do_global_sum = .true.
1223     !-----------------------------------------------
1224     ! Local variables
1225     !-----------------------------------------------
1226         DOUBLE PRECISION, allocatable, Dimension(:,:) :: r_temp, rg_temp
1227         double precision, Dimension(2) :: prod, dot_product_par2
1228         integer :: i, j, k, ijk
1229     !-----------------------------------------------
1230     
1231         if(do_global_sum) then
1232     
1233            prod(:) = 0.0d0
1234     
1235     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1236            do k = kstart1, kend1
1237               do i = istart1, iend1
1238                  do j = jstart1, jend1
1239     
1240                     IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1241     
1242                     ijk = funijk_map_c (i,j,k)
1243                     prod(1) = prod(1) + r1(ijk)*r2(ijk)
1244                     prod(2) = prod(2) + r3(ijk)*r4(ijk)
1245                  enddo
1246               enddo
1247            enddo
1248     
1249            call global_all_sum(prod, dot_product_par2)
1250     
1251         else
1252            allocate (r_temp(ijkstart3:ijkend3,4))
1253            r_temp(:,1) = r1
1254            r_temp(:,2) = r2
1255            r_temp(:,3) = r3
1256            r_temp(:,4) = r4
1257     
1258            if(myPE.eq.root) then
1259               allocate (rg_temp(1:ijkmax3,4))
1260            else
1261               allocate (rg_temp(10,4))
1262            endif
1263            call gather(r_temp,rg_temp)
1264     
1265            if(myPE.eq.root) then
1266               prod = 0.0d0
1267     !$omp parallel do private(i,j,k,ijk) reduction(+:prod)  collapse (3)
1268               do k = kmin1, kmax1
1269                  do i = imin1, imax1
1270                     do j = jmin1, jmax1
1271                        IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1272                        ijk = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1273     !                     ijk = funijk_gl (i,j,k)
1274                        prod(1) = prod(1) + rg_temp(ijk,1)*rg_temp(ijk,2)
1275                        prod(2) = prod(2) + rg_temp(ijk,3)*rg_temp(ijk,4)
1276                     enddo
1277                  enddo
1278               enddo
1279            endif
1280            call bcast( prod)
1281     
1282            dot_product_par2 = prod
1283     
1284            deallocate (r_temp)
1285            deallocate (rg_temp)
1286     
1287         endif
1288     
1289       end function dot_product_par2
1290     
1291     END MODULE leqsol
1292