File: N:\mfix\model\leqsol_mod.f

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