File: /nfs/home/0/users/jenkins/mfix.git/model/leq_gmres.f

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