Results of parallel program with nested loops differ from serial program

Question

I would like to use OpenMP for this single thread code:

PROGRAM SINGLE
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2

!Just intiial value
  DO I=1, 30000
     SUMGRM(I)=I*3        
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000    
  END DO

!Computing computer's running time (start) : for serial programming
 CALL CPU_TIME(TIME1)

 DO K=1, 50000
    DO I=2, 30000
       HSTEP1X=0.0    
         DO J=SUMGRM(I-1)+1, SUMGRM(I)-1
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X       
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X         
    END DO
 END DO

  PRINT *, 'Results  =', HSTEP1X
  PRINT *, '   '

!Computing computer's running time (finish) : for serial programming
 CALL CPU_TIME(TIME2)
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM SINGLE

As you can see, the main problem is located at the most inner side looping (J) which is also a function of most outer side looping (I). I've tried to parallelize this program like this:

PROGRAM PARALLEL
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2,OMP_GET_WTIME
  INTEGER::Q2,P2

!Just intiial value
  DO I=1, 30000
     SUMGRM(I)=I*3        
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000  
  END DO

!Computing computer's running time (start) : for parallel programming
 TIME1= OMP_GET_WTIME()

 DO K=1, 50000
 !$OMP PARALLEL DO PRIVATE (HSTEP1X,Q2,P2)
    DO I=2, 30000
       HSTEP1X=0.0
       Q2=SUMGRM(I-1)+1
       P2=SUMGRM(I)-1
         DO J=Q2, P2
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X       
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X     
    END DO
 !$OMP END PARALLEL DO
 END DO

 PRINT *, 'Results  =', HSTEP1X
 PRINT *, '   '

!Computing computer's running time (finish) : for parallel programming
 TIME2= OMP_GET_WTIME()
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PARALLEL

I'm using gfortran with -O3 -fopenmp and then export OMP_NUM_THREADS=... The parallel program runs faster but the result is different with the single thread code. By the serial program I got 12.1212 (which it is the correct one) and by parallel I got 0.000 (there must be something wrong).

What did I do wrong?

Use tag ptag:fortran] for all Fortran questions. Add a version tag where necessary to distinguish. Notice the number of people subscribing to [tag:fortran] and to [tag:fortran95] and you want more people to see your question, don't you? — Vladimir F Героям слава, Jul 20 '16 at 16:27
It's almost true to say that parallel programs always produce different results to their serial cousins. Be specific about how the results are different, there is much to be learned from the differences. — High Performance Mark, Jul 20 '16 at 16:28
Please try to use titles specifying your problem, not just the topic. The topic is quite well defined by tags. — Vladimir F Героям слава, Jul 20 '16 at 16:32
@VladimirF: What do you mean? Do you mean I also want to gain popularity by the number of views? No, I never intent to get more views of my questions and I really don't care about it. I just thought that I would get the most proper answer from some experts. I even don't know what the benefit is by gaining the number of views, except I could have the bigger chance to get the answer from some experts. Could you please tell me, whether it is prohibited in this site to tag more like I did? If so, I am willing with pleasure to remove it. Thanks. — bestrong, Jul 20 '16 at 16:40
@HighPerformanceMark: The difference is: by single program I got 12.1212 (which it is the true one) and by parallel I got 0.000 (which there must be something wrong). — bestrong, Jul 20 '16 at 16:43
@VladimirF: Thanks for suggestion about the title. You're absolutely right. — bestrong, Jul 20 '16 at 16:45
bob.bob.bob As you say, you want more people to see your question in order to have a better chance for someone to give you a good answer. @VladimirF was suggesting you should use the more general (and appropriate) tag [tag:fortran] as more people follow this than [tag:fortran95] (at least that was my understanding). — d_1999, Jul 20 '16 at 16:54
`real*8` is not valid Fortran. It was never part of any ISO Fortran standard. — jlokimlin, Jul 20 '16 at 17:59
@d_1999: Oh got it. If so, I'm sorry. I didn't know before. Next time I'll pay more attention about tagging and giving the title. Anyway, all critics and suggestions are really appreciated. Many thanks. — bestrong, Jul 21 '16 at 06:41
@jlokimlin: I'm not quite sure if real*8 is not an ISO part of Fortran, but it works on my case. It is only to define the double precision. Could you please explain me, what ''not valid Fortran'' do you mean? Thanks. — bestrong, Jul 21 '16 at 06:44
@bob.bob.bob `datatype*n` is a common Fortran extension, i.e., not officially part of the language, if used with `datatype` other than `character`. If applied to character type it creates an array of `n` characters (or a string of `n` characters). If applied to another numeric type, it specifies the storage size in bytes. Specifying storage size creates non-portable applications. See: http://stackoverflow.com/questions/838310/fortran-90-kind-parameter — jlokimlin, Jul 21 '16 at 13:58

score 2 · Accepted Answer · edited May 23 '17 at 12:22

Firstly we can note that by default you're likely to find that both j and hstep2x are going to be shared between threads. I don't think this is really what you want as it will lead to some very odd behaviour were multiple threads are using the same iteration index but are trying to loop over different ranges.

Next let's note that your serial code actually just prints the result for the i=30000 iteration as the value of hstep1x is reset to 0 at the start of each iteration. As such to get the "correct" answer in the openmp code we could just focus on reproducing the final iteration -- this completely negates the point of using openmp here I think. I'm guessing this is just a simple case you're trying to use to represent your real problem -- I think you may have missed some of the real problem in producing this.

Nevertheless the below code produces the "correct" answer on my machine. I'm not sure how flexible it is but it works here.

PROGRAM PARALLEL
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2,OMP_GET_WTIME
  INTEGER::Q2,P2

!Just intiial value                                                                                                                                                                                                  
  DO I=1, 30000
     SUMGRM(I)=I*3
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000
  END DO

!Computing computer's running time (start) : for parallel programming                                                                                                                                                
 TIME1= OMP_GET_WTIME()

 DO K=1, 50000
!$OMP PARALLEL DO PRIVATE (Q2,P2,J,HSTEP2X) DEFAULT(SHARED) LASTPRIVATE(HSTEP1X)                                                                                                                                     
    DO I=2, 30000
       HSTEP1X=0.0
       Q2= SUMGRM(I-1)+1
       P2= SUMGRM(I)-1
         DO J=Q2,P2
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X
    END DO
!$OMP END PARALLEL DO                                                                                                                                                                                                
END DO

 PRINT *, 'Results  =', HSTEP1X
 PRINT *, '   '

!Computing computer's running time (finish) : for parallel programming                                                                                                                                               
 TIME2= OMP_GET_WTIME()
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PARALLEL

I've done three things here:

Make sure j and hstep2x are private to each thread.
Explicitly declared the default behaviour to be shared (not needed here but never mind).
Specified that hstep1x is lastprivate. This means that after exiting the parallel region the value of hstep1x is that taken from the thread which executed the last iteration. (see here for details).

Yes you're absolutely true. It is not a real case only an imagination case and for sure [DO I=2, 30000] is absolutely not required. I just added it to understand the parallel code works, and also [DO K=1, 50000] is not required , in order to see the representative of the difference CPU's time. Yeah, I forgot about putting [HSTEP1X] to be on [LASTPRIVATE]. By the way the code also works on my machine. Your explanation was great. — bestrong, Jul 21 '16 at 07:05
Oh yes, I forgot to comment that [J] could be defined optionally as [PRIVATE] since it is inside the looping, but both of [Q2] and [P2] should be kept as [PRIVATE]. — bestrong, Jul 21 '16 at 07:18

score 0 · Answer 2 · answered Jul 22 '16 at 11:26

0

Have you tried using

!$OMP PARALLEL DO DEFAULT(PRIVATE) REDUCTION(+:HSTEP1X)

answered Jul 22 '16 at 11:26

Holmz

51
2

I think it won't work, since the main problem here is that [HSTEP1X] has its own value (zero) before the looping [DO J=Q2,P2]. So, if we're using [DEFAULT(PRIVATE)], the value of [HSTEP1X] will remain unchanged as zero. – bestrong Jul 22 '16 at 16:04
You might be right bob.bob.bob, that gfortran seems a lot harder than ifort. I installed it today, and struggling through the differences. In any case, let's ignor the DEFAULT... the REDUCTION seems like it is needed here? – Holmz Jul 23 '16 at 10:31
Could you please specify what "a lot harder" you mean? Since in my experience I used Gfortran and never had any serious problem, even when I am using the new Intel fortran 16 last week, I got the message like it doesn't support the files with format .F95 and later. And finally I should change the format into .F90 (but with the old one Intel Fortran 13, since the new one still didn't work) and then it worked.... Yes, you could absolutely ignore [DEFAULT(SHARED)] in this case, since it doesn't matter. It just about my habit that prefer to write it completely. – bestrong Jul 23 '16 at 16:57
Well I was looking for the gfortran ewuivilent to intel's -diag-file and -vector-report to get , so a lot harder means "I could not do it in 1/2 hour versus years on Intel." The OpenMP seemed to be relying on the -O# for any vectorisation... And hence the vector messages, and I was getting errors on !$OMP DO SIMD REDUCTION(+:K) – Holmz Jul 24 '16 at 00:24

Results of parallel program with nested loops differ from serial program

2 Answers2