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