Don't always get expected results

Question

The following is the code of a simple tidal-transport model that I have included OpenMP to parallelize the computation:

!$OMP PARALLEL SHARED(w, u, v, nthreads, chunk) PRIVATE(i, j, tid)
do it = 1, itlast
    !$OMP DO SCHEDULE(DYNAMIC, CHUNK)
    do j = 2, nyw-1
        do i = 2, nxw-1
            w(i,j) = w(i,j) - rx*depth*(u(i,j) - u(i-1,j))                  &
                            - ry*depth*(v(i,j) - v(i,j-1))
        end do
    end do
    !$OMP END DO
    !$OMP SINGLE
    call boudary_condition(it)
    !$OMP END SINGLE
    !$OMP DO SCHEDULE(DYNAMIC, CHUNK)
    do j = 1, nyw
        jv1 = j
        if (jv1 .ge. nyw-1) jv1 = nyw-1
        do i = 1, nxw-1
            u(i,j) = u(i,j) - rxg*(w(i+1,j) - w(i,j))                       &
                        - constant*u(i,j)*sqrt((u(i,j)**2.) + (v(i,jv1)**2.))
        end do
    end do
    !$OMP END DO
    !$OMP DO SCHEDULE(DYNAMIC, CHUNK)
    do j = 1, nyw-1
        do i = 1, nxw
            iu1 = i
            if (iu1 .ge. nxw-1) iu1 = nxw-1
            v(i,j) = v(i,j) - ryg*(w(i,j+1) - w(i,j))                       &
                        - constant*v(i,j)*sqrt((u(iu1,j)**2.) + (v(i,j)**2.))
        end do
    end do
    !$OMP END DO
    call transport_equation(chunk)
    !$OMP MASTER        
    mtprint = it/ntserprint
    if (ntserprint*mtprint .ne. it) goto 20
        call timeseries(it)
20  continue
    !$OMP END MASTER
end do
!$OMP END PARALLEL

The problem is I don't always get the expected results. Using the same input file, I should always get the same results, but sometimes it produces NaN in the output file. I'm not quite understand why this happens. I'm using Intel Visual Fortran Composer XE 2013 on Windows 10 to compile and run the executable file.

Gilles · Accepted Answer · 2015-09-24T12:39:33.937

2

You need at least to have it, jv1, and ui1 private. Try first fixing these, and let us know.

edited Sep 24 '15 at 12:39

answered Sep 24 '15 at 11:10

Gilles

9,269
4
34
53

Hi @Gilles, thank you for your help. I understand why `iu1` and `jv1` have to be private, but I am not quite understand why `it` also has to be private in this case. `it` here is the iteration, i.e. a time step, and my attached code is time-independent. – Wai Kiat Sep 24 '15 at 13:14
In Fortran, I understand the `SINGLE` directive will order the other threads to wait until the selected thread completes executing the code within the `!$OMP SINGLE` block. Doesn't `MASTER` directive do the same? All threads will synchronize at `!$OMP END MASTER`? – Wai Kiat Sep 24 '15 at 13:18
On a general note, whatever variable that is updated inside a `parallel` region should be examine carefully.Here, `it` being a loop counter for a loop that isn't explicitly `omp do` decorated, it **has to** be declared `private`. Otherwise, different threads will update it on a quite random order, and give undefined behaviour. – Gilles Sep 24 '15 at 13:18
`single` implies `barrier` but do not define which thread will execute the section. `master` restricts the code section to the master thread, but doesn't imply any `barrier` upon exit. – Gilles Sep 24 '15 at 13:22
Hi @Gilles, thank you for the quick respond. I understand it, after having `it`, `iu1` and `jv1` to be within `PRIVATE` clause, and running the executable for 50 times, the model so far didn't produce NaN. Thank you so much. – Wai Kiat Sep 24 '15 at 13:24
Okay, in this case, I need the `barrier` there. I will change that. Thank you again for your help. – Wai Kiat Sep 24 '15 at 13:26
Actually, all loop counters are private automatically, they do not have to be defined as private explicitly. In Fortran (unlike C) it surprisingly even concerns the loops that are nested inside the parallelized one so even i and j are private automatically. – Vladimir F Героям слава Sep 24 '15 at 18:05
1

@VladimirF you're correct, but only for within an`omp do` construct. In this case, `it` is the index counter of a `do` loop, that is inside an `omp parallel` construct, but outside of any `omp do` one. So it is neither the index of a "work-shared" loop, nor the one of a loop nested into such a loop. So my reading of lines 7-9, section 2.7.1 p55 of the [4.0 standard](http://www.openmp.org/mp-documents/OpenMP4.0.0.pdf) let me think that `it` isn't privatised by default... Did I miss or miss-read something? – Gilles Sep 25 '15 at 04:33
No, you are right, I missed that the it loop is just within the parallel region. – Vladimir F Героям слава Sep 25 '15 at 07:43
@VladimirF @Gilles So, in the `PRIVATE` clause, I at least need to have `it`, `iu1`, and `jv1`, and `i`, `j` can not to be included in `PRIVATE` clause? Should I put an `!$OMP BARRIER`, after `!$OMP END MASTER` and before `end do` of the `it` `do`-loop to ensure all threads are synced before going to next iteration? – Wai Kiat Sep 25 '15 at 08:53
Yes, I think so. However, I also strongly encourage you to keep both `i` and `j` in the `private` clause, even if they already are implicitly there, since it makes things clearer. – Gilles Sep 25 '15 at 08:57
@Gilles Thank you for your advice. I will keep that in mind. Thanks! – Wai Kiat Sep 25 '15 at 10:33

Don't always get expected results

1 Answers1