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