MFIX  2016-1
leq_gmres.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: LEQ_GMRES C
4 ! Purpose: Solve system of linear system using GMRES method C
5 ! generalized minimal residual C
6 ! C
7 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
8 ! Reviewer: Date: C
9 ! C
10 ! Literature/Document References: C
11 ! Variables referenced: C
12 ! Variables modified: C
13 ! Local variables: C
14 ! C
15 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16 
17  SUBROUTINE leq_gmres(VNAME, VNO, VAR, A_M, B_M, &
18  cmethod, tol, itmax, max_it, ier)
19 
20 !-----------------------------------------------
21 ! Modules
22 !-----------------------------------------------
23 
24  USE leqsol, only: leq_msolve, leq_matvec
25  USE mpi_utility, only: global_all_and
26  USE param, only: dimension_3
27  USE param1, only: zero
28 
29  IMPLICIT NONE
30 !-----------------------------------------------
31 ! Dummy arguments
32 !-----------------------------------------------
33 ! variable name
34  CHARACTER(LEN=*), INTENT(IN) :: Vname
35 ! variable number (not really used here; see calling subroutine)
36  INTEGER, INTENT(IN) :: VNO
37 ! variable
38 ! e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
39 ! w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
40 ! e_Turb_G
41  DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
42 ! Septadiagonal matrix A_m
43  DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
44 ! Vector b_m
45  DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
46 ! Sweep direction of leq solver (leq_sweep)
47 ! e.g., options = 'isis', 'rsrs' (default), 'asas'
48  CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
49 ! convergence tolerance (generally leq_tol)
50  DOUBLE PRECISION, INTENT(IN) :: TOL
51 ! maximum number of iterations (generally leq_it)
52  INTEGER, INTENT(IN) :: ITMAX
53 ! maximum number of outer iterations
54 ! (currently set to 1 by calling subroutine-solve_lin_eq)
55  INTEGER, INTENT(IN) :: MAX_IT
56 ! error indicator
57  INTEGER, INTENT(INOUT) :: IER
58 !-------------------------------------------------
59 ! Local Variables
60 !-------------------------------------------------
61  LOGICAL :: IS_BM_ZERO, ALL_IS_BM_ZERO
62 !-------------------------------------------------
63 
64  is_bm_zero = (maxval( abs(b_m(:)) ) .EQ. zero)
65  CALL global_all_and( is_bm_zero, all_is_bm_zero )
66 
67  IF (all_is_bm_zero) THEN
68  var(:) = zero
69  RETURN
70  ENDIF
71 
72  CALL leq_gmres0( vname, vno, var, a_m, b_m, &
73  cmethod, tol, itmax, max_it, leq_matvec, leq_msolve, ier )
74 
75  RETURN
76  END SUBROUTINE leq_gmres
77 
78 
79 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
80 ! C
81 ! Subroutine: LEQ_GMRES0 C
82 ! Purpose: Compute residual of linear system C
83 ! C
84 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
85 ! Reviewer: Date: C
86 ! C
87 ! Literature/Document References: C
88 ! Variables referenced: C
89 ! Variables modified: C
90 ! Local variables: C
91 ! C
92 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
93 
94  SUBROUTINE leq_gmres0(VNAME, VNO, VAR, A_M, B_M, &
95  cmethod, tol, itmax, max_it, &
96  matvec, msolve, ier )
97 
98 !-----------------------------------------------
99 ! Modules
100 !-----------------------------------------------
101 
102  USE debug, only: idebug, write_debug
103  USE functions, only: funijk
104  USE gridmap, only: istart3, iend3, jstart3, jend3, kstart3, kend3
105  USE leqsol, only: dot_product_par
106  USE mpi_utility, only: global_all_and, global_all_or, global_all_max, global_all_min
107  USE param, only: dimension_3
108  USE param1, only: one, zero
109 
110  IMPLICIT NONE
111 !-----------------------------------------------
112 ! Dummy arguments/procedure
113 !-----------------------------------------------
114 ! variable name
115  CHARACTER(LEN=*), INTENT(IN) :: Vname
116 ! variable number (not really used here-see calling subroutine)
117  INTEGER, INTENT(IN) :: VNO
118 ! variable
119 ! e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
120 ! w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
121 ! e_Turb_G
122  DOUBLE PRECISION, INTENT(INOUT) :: Var(dimension_3)
123 ! Septadiagonal matrix A_m
124  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3,-3:3)
125 ! Vector b_m
126  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3)
127 ! Sweep direction of leq solver (leq_sweep)
128 ! e.g., options = 'isis', 'rsrs' (default), 'asas'
129  CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
130 ! convergence tolerance (generally leq_tol)
131  DOUBLE PRECISION, INTENT(IN) :: TOL
132 ! maximum number of iterations (generally leq_it)
133  INTEGER, INTENT(IN) :: ITMAX
134 ! maximum number of outer iterations
135  INTEGER, INTENT(IN) :: MAX_IT
136 ! error indicator
137  INTEGER, INTENT(INOUT) :: IER
138 ! dummy arguments/procedures set as indicated
139 ! matvec->leq_matvec
140 ! for preconditioner (leq_pc)
141 ! msolve->leq_msolve
142 !-----------------------------------------------
143 ! Local parameters
144 !-----------------------------------------------
145  INTEGER JDEBUG
146  parameter(jdebug=0)
147 !-----------------------------------------------
148 ! Local variables
149 !-----------------------------------------------
150 
151  DOUBLE PRECISION, dimension(:,:), allocatable :: V
152  DOUBLE PRECISION H(itmax+1,itmax)
153  DOUBLE PRECISION CS(itmax)
154  DOUBLE PRECISION SN(itmax)
155  DOUBLE PRECISION Y(itmax)
156 
157  DOUBLE PRECISION, dimension(:), allocatable :: R,TEMP,WW
158 
159  DOUBLE PRECISION E1(itmax+2)
160  DOUBLE PRECISION SS(itmax+2)
161 
162  DOUBLE PRECISION BNRM2, ERROR, ERROR0
163  DOUBLE PRECISION NORM_R0, NORM_R, NORM_W
164  DOUBLE PRECISION INV_NORM_R, NORM_S, NORM_Y
165  DOUBLE PRECISION YII
166 
167  INTEGER IJK, II, JJ, KK, I, K
168  INTEGER RESTRT, M, ITER, MDIM
169  DOUBLE PRECISION DTEMP
170  DOUBLE PRECISION MINH, MAXH, CONDH, ALL_MINH, ALL_MAXH
171  DOUBLE PRECISION INV_H_IP1_I
172 
173  LOGICAL IS_BM_ZERO, IS_ERROR, IS_CONVERGED
174  LOGICAL ALL_IS_BM_ZERO, ALL_IS_ERROR, ALL_IS_CONVERGED
175 
176  CHARACTER(LEN=40) :: NAME
177 !-----------------------------------------------
178 
179 
180  ! allocating
181  allocate(v(dimension_3,itmax+1))
182  allocate(r(dimension_3))
183  allocate(temp(dimension_3))
184  allocate(ww(dimension_3))
185 
186 ! initializing
187  name = 'LEQ_GMRES0 ' // trim(vname)
188  restrt = itmax
189  m = restrt
190  iter = 0
191  ier = 0
192 
193 
194 ! clear arrays
195 ! --------------------------------
196  r(:) = zero
197  temp(:) = zero
198 
199  h(:,:) = zero
200  cs(:) = zero
201  e1(:) = zero
202  e1(1) = one
203  ss(:) = zero
204  sn(:) = zero
205  ww(:) = zero
206  v(:,:) = zero
207 
208 
209  bnrm2 = dot_product_par( b_m, b_m )
210  bnrm2 = sqrt( bnrm2 )
211 
212  if (idebug.ge.1) then
213  call write_debug(name, 'bnrm2 = ', bnrm2 )
214  endif
215 
216  is_bm_zero = bnrm2 .EQ. zero
217  CALL global_all_and( is_bm_zero, all_is_bm_zero )
218  IF (all_is_bm_zero) THEN
219  bnrm2 = one
220  ENDIF
221 
222 ! r = M \ (b - A*x)
223 ! error = norm(r) / bnrm2
224 ! --------------------------------
225 
226  CALL matvec(vname, var, a_m, r) ! returns R=A*VAR
227 
228 !!$omp parallel do private(ii,jj,kk,ijk)
229  DO kk=kstart3,kend3
230  DO jj=jstart3,jend3
231  DO ii=istart3,iend3
232  ijk = funijk( ii,jj,kk )
233  temp(ijk) = b_m(ijk) - r(ijk)
234  ENDDO
235  ENDDO
236  ENDDO
237 
238 ! Solve A*R(:) = TEMP(:)
239  CALL msolve(vname, temp, a_m, r, cmethod) ! returns R
240 
241  norm_r = dot_product_par( r, r )
242  norm_r = sqrt( norm_r )
243 
244  error = norm_r/bnrm2
245  error0 = error
246  norm_r0 = norm_r
247 
248  IF (jdebug.GE.1) THEN
249  CALL write_debug( name, ' INITIAL ERROR ', error0 )
250  CALL write_debug( name, ' INITIAL RESIDUAL ', norm_r0)
251  ENDIF
252 
253 
254 ! begin iteration
255  DO iter=1,max_it
256 ! ---------------------------------------------------------------->>>
257 
258 ! r = M \ (b-A*x)
259 ! --------------------------------
260  CALL matvec( vname, var, a_m, r ) ! Returns R=A*VAR
261 ! TEMP(:) = B_M(:) - R(:)
262 
263 !!$omp parallel do private(ii,jj,kk,ijk)
264  DO kk=kstart3,kend3
265  DO jj=jstart3,jend3
266  DO ii=istart3,iend3
267  ijk = funijk( ii,jj,kk )
268  temp(ijk) = b_m(ijk) - r(ijk)
269  ENDDO
270  ENDDO
271  ENDDO
272 
273 ! Solve A*R(:) = TEMP(:)
274  CALL msolve(vname, temp, a_m, r, cmethod) ! returns R
275 
276  norm_r = dot_product_par( r, r )
277  norm_r = sqrt( norm_r )
278  inv_norm_r = one / norm_r
279 
280 ! V(:,1) = R(:) * INV_NORM_R
281 !!$omp parallel do private(ii,jj,kk,ijk)
282  DO kk=kstart3,kend3
283  DO jj=jstart3,jend3
284  DO ii=istart3,iend3
285  ijk = funijk( ii,jj,kk )
286  v(ijk,1) = r(ijk) * inv_norm_r
287  ENDDO
288  ENDDO
289  ENDDO
290 
291  ss(:) = norm_r * e1(:)
292 
293 
294 ! construct orthonormal basis using Gram-Schmidt
295 ! ---------------------------------------------------------------->>>
296  DO i=1,m ! M->restrt->itmax
297 ! w = M \ (A*V(:,i))
298 ! --------------------------------
299  CALL matvec(vname, v(:,i), a_m, temp) ! returns TEMP=A*V
300 ! Solve A*WW(:) = TEMP(:)
301  CALL msolve(vname, temp, a_m, ww, cmethod) ! returns WW
302 
303  DO k=1,i
304  dtemp = dot_product_par( ww(:), v(:,k) )
305  h(k,i) = dtemp
306 ! WW(:) = WW(:) - H(K,I)*V(:,K)
307 
308 !!$omp parallel do private(ii,jj,kk,ijk)
309  DO kk=kstart3,kend3
310  DO jj=jstart3,jend3
311  DO ii=istart3,iend3
312  ijk = funijk( ii,jj,kk )
313  ww(ijk) = ww(ijk) - h(k,i)*v(ijk,k)
314  ENDDO
315  ENDDO
316  ENDDO
317  ENDDO
318 
319  norm_w = dot_product_par( ww(:), ww(:) )
320  norm_w = sqrt( norm_w )
321  h(i+1,i) = norm_w
322 ! V(:,I+1) = WW(:) / H(I+1,I)
323  inv_h_ip1_i = one / h(i+1,i)
324 
325 !!$omp parallel do private(ii,jj,kk,ijk)
326  DO kk=kstart3,kend3
327  DO jj=jstart3,jend3
328  DO ii=istart3,iend3
329  ijk = funijk( ii,jj,kk )
330  v(ijk, i+1) = ww(ijk) * inv_h_ip1_i
331  ENDDO
332  ENDDO
333  ENDDO
334 
335 ! apply Givens rotation
336 ! --------------------------------
337  DO k=1,i-1
338  dtemp = cs(k)*h(k,i) + sn(k)*h(k+1,i)
339  h(k+1,i) = -sn(k)*h(k,i) + cs(k)*h(k+1,i)
340  h(k,i) = dtemp
341  ENDDO
342 
343 ! form i-th rotation matrix approximate residual norm
344 ! --------------------------------
345  CALL rotmat( h(i,i), h(i+1,i), cs(i), sn(i) )
346 
347  dtemp = cs(i)*ss(i)
348  ss(i+1) = -sn(i)*ss(i)
349  ss(i) = dtemp
350  h(i,i) = cs(i)*h(i,i) + sn(i)*h(i+1,i)
351  h(i+1,i) = zero
352  error = abs( ss(i+1) ) / bnrm2
353 
354  is_converged = (error .LE. tol*error0)
355  CALL global_all_and( is_converged, all_is_converged )
356  IF (all_is_converged) THEN
357 ! update approximation and exit
358 ! --------------------------------
359 
360 ! triangular solve with y = H(1:i,1:i) \ s(1:i)
361 ! --------------------------------
362  mdim = i
363  DO ii=1,mdim
364  y(ii) = ss(ii)
365  ENDDO
366 
367  DO ii=mdim,1,-1
368  yii = y(ii)/h(ii,ii)
369  y(ii) = yii
370  DO jj=1,ii-1
371  y(jj) = y(jj) - h(jj,ii)*yii
372  ENDDO
373  ENDDO
374 
375 ! double check
376 ! --------------------------------
377  IF (jdebug.GE.1) THEN
378  maxh = abs(h(1,1))
379  minh = abs(h(1,1))
380  DO ii=1,mdim
381  maxh = max( maxh, abs(h(ii,ii)) )
382  minh = min( minh, abs(h(ii,ii)) )
383  ENDDO
384 
385  CALL global_all_max( maxh, all_maxh )
386  CALL global_all_min( minh, all_minh )
387  maxh = all_maxh
388  minh = all_minh
389  condh = maxh/minh
390 
391  temp(1:mdim) = ss(1:mdim) - &
392  matmul( h(1:mdim,1:mdim ), y(1:mdim) )
393  dtemp = dot_product( temp(1:mdim), temp(1:mdim) )
394  dtemp = sqrt( dtemp )
395 
396  norm_s = dot_product( ss(1:mdim),ss(1:mdim) )
397  norm_s = sqrt( norm_s )
398 
399  norm_y = dot_product( y(1:mdim),y(1:mdim) )
400  norm_y = sqrt( norm_y )
401 
402  is_error = (dtemp .GT. condh*norm_s)
403  CALL global_all_or( is_error, all_is_error )
404  IF (all_is_error) THEN
405  CALL write_debug(name, &
406  'DTEMP, NORM_S ', dtemp, norm_s )
407  CALL write_debug(name, &
408  'CONDH, NORM_Y ', condh, norm_y )
409  CALL write_debug(name, &
410  '** STOP IN LEQ_GMRES ** ')
411  ier = 999
412  RETURN
413  ENDIF
414  ENDIF ! end if(jdebug>=1)
415 
416 
417 ! VAR(:)=VAR(:)+MATMUL(V(:,1:I),Y(1:I))
418 !!$omp parallel do private(ii,jj,kk,ijk)
419  DO kk=kstart3,kend3
420  DO jj=jstart3,jend3
421  DO ii=istart3,iend3
422  ijk = funijk( ii,jj,kk )
423  var(ijk) = var(ijk) + &
424  dot_product( v(ijk,1:i), y(1:i) )
425  ENDDO
426  ENDDO
427  ENDDO
428 
429  EXIT
430  ENDIF !end if(all_is_converged)
431  ENDDO !end do I=1,m (construct orthonormal basis using Gram-Schmidt)
432 ! ----------------------------------------------------------------<<<
433 
434  is_converged = ( error .LE. tol*error0)
435  CALL global_all_and( is_converged, all_is_converged )
436  IF ( all_is_converged ) THEN
437  EXIT
438  ENDIF
439 
440 ! update approximations
441 ! --------------------------------
442 
443 ! y = H(1:m,1:m) \ s(1:m)
444 ! x = x + V(:,1:m)*y
445 ! r = M \ (b-A*x)
446 ! --------------------------------
447  mdim = m
448  DO ii=1,mdim
449  y(ii) = ss(ii)
450  ENDDO
451 
452  DO ii=mdim,1,-1
453  yii = y(ii)/h(ii,ii)
454  y(ii) = yii
455  DO jj=1,ii-1
456  y(jj) = y(jj) - h(jj,ii)*yii
457  ENDDO
458  ENDDO
459 
460 ! double check
461 ! --------------------------------
462  IF (jdebug.GE.1) THEN
463  maxh = abs(h(1,1))
464  minh = abs(h(1,1))
465  DO ii=1,mdim
466  maxh = max( maxh, abs(h(ii,ii)) )
467  minh = min( minh, abs(h(ii,ii)) )
468  ENDDO
469  condh = maxh/minh
470 
471  temp(1:mdim) = ss(1:mdim) - &
472  matmul( h(1:mdim,1:mdim ), y(1:mdim) )
473 
474  dtemp = dot_product( temp(1:mdim), temp(1:mdim) )
475  dtemp = sqrt( dtemp )
476 
477  norm_s = dot_product( ss(1:mdim),ss(1:mdim) )
478  norm_s = sqrt( norm_s )
479 
480  norm_y = dot_product( y(1:mdim),y(1:mdim) )
481  norm_y = sqrt( norm_y )
482 
483  is_error = (dtemp .GT. condh*norm_s)
484  CALL global_all_or( is_error, all_is_error )
485  IF (all_is_error) THEN
486  CALL write_debug(name, &
487  'DTEMP, NORM_S ', dtemp, norm_s )
488  CALL write_debug(name, &
489  'CONDH, NORM_Y ', condh, norm_y )
490  CALL write_debug(name, &
491  '** STOP IN LEQ_GMRES ** ')
492  ier = 999
493  RETURN
494  ENDIF
495  ENDIF ! end if(jdebug>=1)
496 
497 ! VAR(:) = VAR(:) + MATMUL( V(:,1:M), Y(1:M) )
498 !!$omp parallel do private(ii,jj,kk,ijk)
499  DO kk=kstart3,kend3
500  DO jj=jstart3,jend3
501  DO ii=istart3,iend3
502  ijk = funijk( ii,jj,kk )
503  var(ijk) = var(ijk) + &
504  dot_product( v(ijk,1:m), y(1:m) )
505  ENDDO
506  ENDDO
507  ENDDO
508 
509  CALL matvec(vname, var, a_m, r) ! returns R=A*VAR
510 
511 ! TEMP(:) = B_M(:) - R(:)
512 !!$omp parallel do private(ii,jj,kk,ijk)
513  DO kk=kstart3,kend3
514  DO jj=jstart3,jend3
515  DO ii=istart3,iend3
516  ijk = funijk( ii,jj,kk )
517  temp(ijk) = b_m(ijk) - r(ijk)
518  ENDDO
519  ENDDO
520  ENDDO
521 
522 ! Solve A*R(:)=TEMP(:)
523  CALL msolve( vname, temp, a_m, r, cmethod) ! Returns R
524  norm_r = dot_product_par( r, r )
525  norm_r = sqrt( norm_r )
526 
527  ss(i+1) = norm_r
528  error = ss(i+1) / bnrm2
529 
530  IF (jdebug.GE.1) THEN
531  CALL write_debug(name, 'LEQ_GMRES: I, ITER ', i, iter)
532  CALL write_debug(name, 'LEQ_GMRES: ERROR, NORM_R ', &
533  error, norm_r )
534  ENDIF
535 
536  is_converged = (error .LE. tol*error0)
537  CALL global_all_and( is_converged, all_is_converged )
538  IF (all_is_converged) THEN
539  EXIT
540  ENDIF
541 
542  ENDDO ! end do iter=1,max_it
543 ! ----------------------------------------------------------------<<<
544 
545  is_error = (error .GT. tol*error0)
546  CALL global_all_or( is_error, all_is_error )
547  IF (all_is_error) THEN
548  ier = 1
549  ENDIF
550 
551  IF (jdebug.GE.1) THEN
552  CALL matvec(vname, var, a_m, r) ! Returns R=A*VAR
553 
554 ! R(:) = R(:) - B_M(:)
555 !!$omp parallel do private(ii,jj,kk,ijk)
556  DO kk=kstart3,kend3
557  DO jj=jstart3,jend3
558  DO ii=istart3,iend3
559  ijk = funijk( ii,jj,kk )
560  r(ijk) = r(ijk) - b_m(ijk)
561  ENDDO
562  ENDDO
563  ENDDO
564 
565  norm_r=dot_product_par(r,r)
566  norm_r = sqrt( norm_r )
567  CALL write_debug(name, 'ITER, I ', iter, i )
568  CALL write_debug(name, 'VNAME = ' // vname )
569  CALL write_debug(name, 'IER = ', ier )
570  CALL write_debug(name, 'ERROR0 = ', error0 )
571  CALL write_debug(name, 'NORM_R0 = ', norm_r0 )
572  CALL write_debug(name, 'FINAL ERROR = ', error )
573  CALL write_debug(name, 'FINAL RESIDUAL ', norm_r )
574  CALL write_debug(name, 'ERR RATIO ', error/error0 )
575  CALL write_debug(name, 'RESID RATIO ', norm_r/norm_r0 )
576  ENDIF ! end if (jdebug>=1)
577 
578  deallocate(v)
579  deallocate(r)
580  deallocate(temp)
581  deallocate(ww)
582 
583  RETURN
584  END SUBROUTINE leq_gmres0
585 
586 
587 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
588 ! C
589 ! C
590 ! C
591 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
592 
593  SUBROUTINE rotmat( A, B, C, S )
594 
595 !-----------------------------------------------
596 ! Modules
597 !-----------------------------------------------
598  IMPLICIT NONE
599 !-----------------------------------------------
600 ! Dummy arguments
601 !-----------------------------------------------
602  DOUBLE PRECISION A,B,C,S
603 !-----------------------------------------------
604 ! Local parameters
605 !-----------------------------------------------
606  DOUBLE PRECISION ONE,ZERO
607  parameter(one=1.0d0,zero=0.0d0)
608 !-----------------------------------------------
609 ! Local variables
610 !-----------------------------------------------
611  DOUBLE PRECISION TEMP
612 !-----------------------------------------------
613 
614  IF (b.EQ.zero) THEN
615  c = one
616  s = zero
617  ELSEIF (abs(b) .GT. abs(a)) THEN
618  temp = a / b
619  s = one / sqrt( one + temp*temp )
620  c = temp * s
621  ELSE
622  temp = b / a
623  c = one / sqrt( one + temp*temp )
624  s = temp * c
625  ENDIF
626 
627  RETURN
628  END
629 
630 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
631 ! C
632 ! Subroutine: Leq_check C
633 ! Purpose: Verify boundary nodes in A_m are set correctly C
634 ! C
635 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
636 ! Reviewer: Date: C
637 ! C
638 ! Literature/Document References: C
639 ! Variables referenced: C
640 ! Variables modified: C
641 ! Local variables: C
642 ! C
643 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
644 
645  subroutine leq_check( vname, A_m )
646 
647 !-----------------------------------------------
648 ! Modules
649 !-----------------------------------------------
650  USE param
651  USE param1
652  USE geometry
653  USE indices
654  USE debug
655  USE compar
656  USE mpi_utility
657  USE functions
658  IMPLICIT NONE
659 !-----------------------------------------------
660 ! Dummy arguments
661 !-----------------------------------------------
662 ! Septadiagonal matrix A_m
663  DOUBLE PRECISION, INTENT(INOUT) :: A_M(dimension_3, -3:3)
664 ! Variable name
665  CHARACTER(LEN=*), INTENT(IN) :: VNAME
666 !-----------------------------------------------
667 ! Local parameters
668 !-----------------------------------------------
669  logical, parameter :: do_reset = .true.
670 !-----------------------------------------------
671 ! Local variables
672 !-----------------------------------------------
673  integer, dimension(-3:3) :: ijktable
674  integer :: istartl, iendl, jstartl, jendl, kstartl, kendl, elstart, elend
675  integer :: i,j,k,el, ijk,ijk2, ii,jj,kk
676  integer :: nerror, all_nerror
677  logical :: is_in_k, is_in_j, is_in_i, is_in
678  logical :: is_bc_k, is_bc_j, is_bc_i, is_bc, is_ok
679 !-----------------------------------------------
680 
681  kstartl = kstart2
682  kendl = kend2
683  jstartl = jstart2
684  jendl = jend2
685  istartl = istart2
686  iendl = iend2
687 
688  if (no_k) then
689  kstartl = 1
690  kendl = 1
691  endif
692 
693  nerror = 0
694 
695  do k=kstartl,kendl
696  do j=jstartl,jendl
697  do i=istartl,iendl
698  is_in_k = (kstart1 <= k) .and. (k <= kend1)
699  is_in_j = (jstart1 <= j) .and. (j <= jend1)
700  is_in_i = (istart1 <= i) .and. (i <= iend1)
701 
702  is_in = is_in_k .and. is_in_j .and. is_in_i
703  if (is_in) cycle
704 
705  ijk = funijk(i,j,k)
706  ijktable( -2 ) = jm_of(ijk)
707  ijktable( -1 ) = im_of(ijk)
708  ijktable( 0 ) = ijk
709  ijktable( 1 ) = ip_of(ijk)
710  ijktable( 2 ) = jp_of(ijk)
711 
712  if (.not. no_k) then
713  ijktable( -3 ) = km_of(ijk)
714  ijktable( 3 ) = kp_of(ijk)
715  endif
716 
717  elstart = -3
718  elend = 3
719  if (no_k) then
720  elstart = -2
721  elend = 2
722  endif
723 
724  do el=elstart,elend
725  ijk2 = ijktable( el )
726  kk = k_of(ijk2)
727  jj = j_of(ijk2)
728  ii = i_of(ijk2)
729  is_bc_k = (kk < kstart2) .or. (kk > kend2)
730  is_bc_j = (jj < jstart2) .or. (jj > jend2)
731  is_bc_i = (ii < istart2) .or. (ii > iend2)
732  is_bc = is_bc_k .or. is_bc_j .or. is_bc_i
733  if (is_bc) then
734  is_ok = (a_m(ijk,el).eq.zero)
735  if (.not.is_ok) then
736  nerror = nerror + 1
737  if (do_reset) a_m(ijk,el) = zero
738  endif
739  endif
740  enddo ! do el
741  enddo
742  enddo
743  enddo
744 
745  call global_sum( nerror, all_nerror )
746  nerror = all_nerror
747 
748  if ((nerror >= 1) .and. (mype.eq.pe_io)) then
749  if (do_reset) then
750  call write_debug( 'leq_check: ' // trim( vname ), &
751  'A_m is reset. nerror = ', nerror )
752  else
753  call write_debug( 'leq_check: ' // trim( vname ), &
754  'A_m is not reset. nerror = ', nerror )
755  endif
756  endif
757 
758  return
759  end subroutine leq_check
compar::jend2
integer jend2
Definition: compar_mod.f:80
mpi_utility
Definition: mpi_utility_mod.f:6
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
mpi_utility::global_all_and
Definition: mpi_utility_mod.f:126
compar::kend1
integer kend1
Definition: compar_mod.f:80
compar::istart1
integer istart1
Definition: compar_mod.f:80
leqsol::leq_matvec
subroutine leq_matvec(VNAME, VAR, A_M, Avar)
Definition: leqsol_mod.f:104
compar::iend1
integer iend1
Definition: compar_mod.f:80
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
compar::istart2
integer istart2
Definition: compar_mod.f:80
leq_check
subroutine leq_check(vname, A_m)
Definition: leq_gmres.f:646
compar::iend2
integer iend2
Definition: compar_mod.f:80
debug::write_debug
Definition: debug_mod.f:15
leqsol::leq_msolve
subroutine leq_msolve(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:287
compar::kend2
integer kend2
Definition: compar_mod.f:80
compar::kstart2
integer kstart2
Definition: compar_mod.f:80
compar::kstart1
integer kstart1
Definition: compar_mod.f:80
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
mpi_utility::global_all_max
Definition: mpi_utility_mod.f:110
leq_gmres
subroutine leq_gmres(VNAME, VNO, VAR, A_M, B_M, cmethod, TOL, ITMAX, MAX_IT, IER)
Definition: leq_gmres.f:19
leqsol::dot_product_par
double precision function dot_product_par(r1, r2)
Definition: leqsol_mod.f:1095
debug
Definition: debug_mod.f:1
compar::jstart2
integer jstart2
Definition: compar_mod.f:80
mpi_utility::global_all_or
Definition: mpi_utility_mod.f:132
mpi_utility::global_sum
Definition: mpi_utility_mod.f:60
mpi_utility::global_all_min
Definition: mpi_utility_mod.f:88
rotmat
subroutine rotmat(A, B, C, S)
Definition: leq_gmres.f:594
param
Definition: param_mod.f:2
debug::idebug
integer idebug
Definition: debug_mod.f:7
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
compar::mype
integer mype
Definition: compar_mod.f:24
leq_gmres0
subroutine leq_gmres0(VNAME, VNO, VAR, A_M, B_M, CMETHOD, TOL, ITMAX, MAX_IT, MATVEC, MSOLVE, IER)
Definition: leq_gmres.f:97
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
compar::jend1
integer jend1
Definition: compar_mod.f:80
compar::jstart1
integer jstart1
Definition: compar_mod.f:80