File: N:\mfix\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     
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
760