Итеративный решатель Gauss Siedel с использованием OpenMP не сходится

Я хочу написать итеративный решатель Gauss Siedel с использованием OpenMP. Мой решатель не сходится к правильному результату, и я не мог понять почему.

Основное уравнение представляет собой стационарное уравнение тепловыделения: del^2(T)=S, где S — функция источника тепла. S=-35*pi/2*cos(pi*x/2)*cos(3*pi*y/2)*cos(5*pi*z/2)

Я реализую OpenMP внутри цикла dowhile, поскольку он не позволяет мне запускать параллельную форму цикла dowhile. Есть ли способ изменить цикл do while на цикл do?

Результат сходится без параллельных вычислений, но после добавления опенмпа потом уже не сходится.

Вот мой код:

    PROGRAM GAUSS_MP
        USE omp_lib
        IMPLICIT NONE

        REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: S, T
        REAL*8 :: X, Y, Z
        REAL*8, PARAMETER :: PI=2*ASIN(1.0)
        REAL*8 :: DX                     ! STEP SIZE DX=DY=DZ
        REAL*8, PARAMETER :: TOL=1.0E-5  ! TOLERANCE 
        REAL*8 :: DUMAX
        REAL*8 :: T_OLD
        REAL*8 :: T1,T2

        INTEGER, PARAMETER :: N=100        ! GRID NUMBER, START FROM 1
        INTEGER, PARAMETER :: ITERMAX=1E5   ! MAXIMUM ITERATION
        INTEGER :: I, J, ITER, K
        INTEGER :: POINT_NUM
        INTEGER, PARAMETER :: NUM_THREADS=32
    !    INTEGER :: A

        ! INITIALIZE OPENMP THREADS
        CALL OMP_SET_NUM_THREADS(NUM_THREADS)


        ! COMPUTE STEP SIZE
        DX=2.0/REAL(N-1, KIND=8)  
    !    PRINT *, "DX = ", DX
        
        ! INITIALIZE THE T ARRAYS AND S(I)
        ALLOCATE(S(N,N,N), T(N,N,N))
``````````````````````````````````````````````````````````````````````````````            
        ! INITIAL GUESS
        T=1.0
``````````````````````````````````````````````````````````````````````````````            
        ! BOUNDARY CONDITION
        T(1,:,:)=0.0; T(N,:,:)=0.0; T(:,:,1)=0.0; T(:,:,N)=0.0;
        T(:,1,:)=0.0; T(:,N,:)=0.0;
        
        X=0.0D0 ! COORDINATES
        Y=0.0D0

        S=0.0D0 ! SOURCE
``````````````````````````````````````````````````````````````````````````````            
        ! SOURCE MATRIX
        !$OMP PARALLEL DO PRIVATE(I,J,K)
         DO K=2,N-1
            Z=-1.0+(K-1)*DX
            DO I=2,N-1
                Y=-1.0+(I-1)*DX
                DO J=2,N-1
                    X=-1.0+(J-1)*DX
                    S(I,J,K)=-35.0*PI/2.0*COS(PI*X/2.0)*COS(3.0*PI*Y/2.0)&
                            *COS(5.0*PI*Z/2.0)
                END DO
            END DO
        END DO
        !$OMP END PARALLEL DO
``````````````````````````````````````````````````````````````````````````````            
        ! GAUSS-SEIDEL ITERATION
        PRINT *, 'PARALLEL START'
        T1=OMP_GET_WTIME()
      
        ITER=0
        DUMAX=1.0D0 ! UPDATE DIFFERENCE
        DO WHILE(ITER <= ITERMAX .AND. DUMAX >= TOL)
            ITER=ITER+1
            DUMAX=0.0D0
            !$OMP PARALLEL PRIVATE(T_OLD, K, I, J, DUMAX)
            !$OMP DO REDUCTION(MAX:DUMAX)
            DO K=2,N-1
                DO I=2, N-1
                    DO J=2, N-1
                        T_OLD=T(I,J,K)
                        T(I,J,K)=1.0/6.0*(T(I+1,J,K)+T(I-1,J,K)+T(I,J-1,K)+T(I,J+1,K) &
                                    +T(I,J,K+1)+T(I,J,K-1) &
                                    -DX**2*S(I,J,K))
                        DUMAX=MAX(DUMAX, ABS(T_OLD-T(I,J,K)))
                    END DO
                END DO
            END DO
            !$OMP END DO
           !$OMP END PARALLEL
        END DO
        
        T2=OMP_GET_WTIME()

    END PROGRAM GAUSS_MP