The is a revision of Lecture 10. It discusses race conditions on OpenMP shared variables and gives a few options.
On the p690, with the xlf (or xlf90 compiler) the OpenMP library is specified by the -qsmp=omp flag. To disable automatic parallelization, but recognize OpenMP directives, try the compiler flag -qsmp=omp:noauto. As with the ifc and icc compilers, there appears to be no mechanism for inserting OpenMP directives directly into Fortran or C codes.
In tcsh you could specify use of 4 threads by
>setenv OMP_NUM_THREADS 4
For the p690, the bsub file could be almost the same as on the blades
#!/bin/csh #BSUB -W 10 #BSUB -n 4 #BSUB -o out.%J ./foo
but here we are asking for 4 processors and did not need the #BSUB -R line. The IBM automatic parallelism inserts SMP directives (SMP is the IBM set of directives, which predates OpenMP). OpenMP directives are also translated into SMP. To see a listing of the of the SMP code, compile with
>xlf -qsmp=omp:auto -qreport=smp -c spcol.f
The same spcol.f program for this week's homework compiled to give a spcol.o file and also a file spcol.lst, the listing of which follows.
XL Fortran for AIX Version 08.01.0000.0003 --- spcol.f 09/15/04 12:16:03
>>>>> OPTIONS SECTION <<<<<
*** Options In Effect ***
== On / Off Options ==
32 CCLINES ESCAPE I4
NOLM OBJECT SWAPOMP UNWIND
NOZEROSIZE
== Options Of Integer Type ==
FIXED(72) HOT() MAXMEM(2048)
OPTIMIZE(2) SPILLSIZE(512)
== Options of Integer and Character Type ==
SMP(OMP,AUTO,THRESHOLD(100),SCHEDULE(RUNTIME))
== Options Of Character Type ==
ALIAS(STD,INTPTR) ALIGN(STRUCT(NATURAL))
ARCH(COM) AUTODBL(NONE) DIRECTIVE($OMP)
FLAG(I,I) FLOAT(MAF,FOLD) HALT(S)
HOT(VECTOR) IEEE(NEAR) INTSIZE(4)
LANGLVL(EXTENDED) POSITION(APPENDOLD) REALSIZE(4)
REPORT(SMPLIST) SAVE(ALL) UNROLL(AUTO)
XFLAG() XLF77(NOLEADZERO,GEDIT77,NOBLANKPAD,OLDBOZ,INTARG,INTXOR,PERSISTENT,SOFTEOF)
XLF90(NOSIGNEDZERO,NOAUTODEALLOC)
>>>>> SOURCE SECTION <<<<<
** spcol === End of Compilation 1 ===
>>>>> PARALLELIZATION AND LOOP TRANSFORMATION SECTION <<<<<
PROGRAM spcol ()
! DIR_SCHEDTYPE loopId = -1 schedType = 5
19| #1 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat("input n",7,T_13,7)
CALL _xlfWriteLDChar(%VAL(#1),T_13,7,1)
_xlfEndIO(%VAL(#1))
20| #2 = _xlfBeginIO(5,1,NULL,1024,NULL,0,NULL)
CALL _xlfReadLDInt(%VAL(#2),n,4,4)
_xlfEndIO(%VAL(#2))
21| #3 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat(" input fraction of nonzeros",27,T_14,27)
CALL _xlfWriteLDChar(%VAL(#3),T_14,27,1)
_xlfEndIO(%VAL(#3))
22| #4 = _xlfBeginIO(5,1,NULL,1024,NULL,0,NULL)
CALL _xlfReadLDReal(%VAL(#4),frac,8,8)
_xlfEndIO(%VAL(#4))
23| iseed = n
24| sum = randw(iseed)
32| pretim = ftime()
33| k = 1
34| IF ((n > 0)) THEN
@CIVFINAL1 = int(n)
@CIV1 = 0
Id=1 DO @CIV1 = @CIV1, @CIVFINAL1-1
35| ic((@CIV1 + 1)) = k
36| i = 1
37| lab_2 /* loopid=9 */
IF (.NOT.(i < int((frac * real(n))) + 1)) GOTO lab_3
lab_4 /* loopid=10 */
38| a(k) = ( 1.0000000000000000E+000 / frac) / real(n)
@RET0 = randw(iseed)
39| @CSE1 = real(n)
ir(k) = int((@RET0 * @CSE1 + 1.0000000000000000E+000))
44| @CSE0 = (k - ic((@CIV1 + 1)))
IF ((@CSE0 > 0)) THEN
@CIVFINAL0 = int(@CSE0)
@CIV0 = 0
Id=2 DO @CIV0 = @CIV0, @CIVFINAL0-1
45| IF (ir((ic((@CIV1 + 1)) + @CIV0)) == ir(k)) GOTO &
& lab_2
51| IF ((ir((ic((@CIV1 + 1)) + @CIV0)) > ir(k))) THEN
52| itemp = ir((ic((@CIV1 + 1)) + @CIV0))
53| ir((ic((@CIV1 + 1)) + @CIV0)) = ir(k)
54| ir(k) = itemp
55| ENDIF
56| ENDDO
ENDIF
57| k = (k + 1)
58| i = (i + 1)
59| IF (i < int((frac * @CSE1)) + 1) GOTO lab_4
lab_3
60| ENDDO
ENDIF
61| ic((n + 1)) = k
62| nz = k
63| IF ((n > 0)) THEN
@DCIV0 = 0
@ICM.n0 = n
Id=11 DO @DCIV0 = @DCIV0, int(@ICM.n0)-1
! DIR_INDEPENDENT loopId = 0
64| x((@DCIV0 + 1)) = 1.0000000000000000E+000
63| ENDDO
@DCIV1 = 0
@ICM.n0 = n
Id=12 DO @DCIV1 = @DCIV1, int(@ICM.n0)-1
! DIR_INDEPENDENT loopId = 0
65| b((@DCIV1 + 1)) = 0.0000000000000000E+000
63| ENDDO
ENDIF
66| @RET1 = real(ftime())
67| entim = real((@RET1 - real(pretim)))
68| #5 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat(" initialization time was ",25,T_15,25)
CALL _xlfWriteLDChar(%VAL(#5),T_15,25,1)
CALL _xlfWriteLDReal(%VAL(#5),entim,4,4)
_xlfEndIO(%VAL(#5))
72| its = (100000000 / nz)
73| pretim = ftime()
74| IF ((its > 0)) THEN
@CIV7 = 0
IF ((n > 0)) THEN
@ICM.n0 = n
@ICM.its3 = its
Id=4 DO @CIV7 = @CIV7, int(@ICM.its3)-1
75| @CIV4 = 0
Id=5 DO @CIV4 = @CIV4, int(@ICM.n0)-1
78| @CSE3 = (ic((@CIV4 + 2)) - ic((@CIV4 + 1)))
IF ((@CSE3 > 0)) THEN
@CIV3 = 0
@ICM1 = x((@CIV4 + 1))
@ICM2 = @CSE3
Id=6 DO @CIV3 = @CIV3, int(@ICM2)-1
79| b(ir((ic((@CIV4 + 1)) + @CIV3))) = b(ir((ic((&
& @CIV4 + 1)) + @CIV3))) + a((ic((@CIV4 + 1)) + &
& @CIV3)) * @ICM1
80| ENDDO
ENDIF
81| ENDDO
82| sum = 0.0000000000000000E+000
83| @CIV5 = 0
Id=7 DO @CIV5 = @CIV5, int(@ICM.n0)-1
84| sum = (sum + abs(b((@CIV5 + 1))))
85| ENDDO
86| sum = ( 1.0000000000000000E+000 / sum)
87| @CIV6 = 0
Id=8 DO @CIV6 = @CIV6, int(@ICM.n0)-1
! DIR_INDEPENDENT loopId = 0
88| x((@CIV6 + 1)) = sum * b((@CIV6 + 1))
89| b((@CIV6 + 1)) = 0.0000000000000000E+000
90| ENDDO
91| ENDDO
ELSEIF
@ICM.its3 = its
74|Id=3 DO @CIV7 = @CIV7, int(@ICM.its3)-1
81| sum = 0.0000000000000000E+000
sum = ( 1.0000000000000000E+000 / sum)
91| ENDDO
ENDIF
ENDIF
@RET2 = real(ftime())
92| @CSE2 = (@RET2 - real(pretim))
entim = real(@CSE2)
93| #6 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _xlfWriteLDReal(%VAL(#6),entim,4,4)
CALL _xlfWriteLDInt(%VAL(#6),its,4,4)
_xlfEndIO(%VAL(#6))
94| mflops = ((( 2.0000000000000000E+000 * real(its)) * (real(nz) &
& / real(real(@CSE2)))) / 1.0000000000000000E+006)
95| #7 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat("sum = ",6,T_16,6)
CALL _xlfWriteLDChar(%VAL(#7),T_16,6,1)
T_17 = (sum / real(n))
CALL _xlfWriteLDReal(%VAL(#7),T_17,8,8)
_xlfEndIO(%VAL(#7))
96| #8 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat(" mflops =",9,T_18,9)
CALL _xlfWriteLDChar(%VAL(#8),T_18,9,1)
CALL _xlfWriteLDReal(%VAL(#8),mflops,8,8)
_xlfEndIO(%VAL(#8))
97| #9 = _xlfBeginIO(6,257,NULL,1024,NULL,0,NULL)
CALL _dconcat("nz =",4,T_19,4)
CALL _xlfWriteLDChar(%VAL(#9),T_19,4,1)
CALL _xlfWriteLDInt(%VAL(#9),nz,4,4)
_xlfEndIO(%VAL(#9))
99| CALL _xlfExit(0)
CALL _trap(3)
100| END PROGRAM spcol
Source Source Loop Id Action / Information
File Line
---------- ---------- ---------- ----------------------------------------------------------
0 34 1 Loop cannot be automatically parallelized. Loop
contains a call to "randw" that may have side effects.
0 44 2 Private variables recognized in this loop nest:
"itemp".
0 44 2 Loop cannot be automatically parallelized. A
dependency is carried by variable aliasing or
function call.
0 74 4 Private variables recognized in this loop nest:
"@CIV4", "@CIV3", "@CIV5", "@CIV6", and "@CIV3".
0 74 4 Loop cannot be automatically parallelized. A
dependency is carried by variable "b".
0 75 5 Private variables recognized in this loop nest:
"@CIV3".
0 75 5 Loop cannot be automatically parallelized. A
dependency is carried by variable "b".
0 78 6 Loop cannot be automatically parallelized. A
dependency is carried by variable "b".
>>>>> COMPILATION UNIT EPILOGUE SECTION <<<<<
TOTAL UNRECOVERABLE SEVERE ERROR WARNING INFORMATIONAL
(U) (S) (E) (W) (I)
0 0 0 0 0 0
>>>>> OPTIONS SECTION <<<<<
*** Options In Effect ***
== On / Off Options ==
32 CCLINES ESCAPE I4
NOLM OBJECT SWAPOMP UNWIND
NOZEROSIZE
== Options Of Integer Type ==
FIXED(72) HOT() MAXMEM(2048)
OPTIMIZE(2) SPILLSIZE(512)
== Options of Integer and Character Type ==
SMP(OMP,AUTO,THRESHOLD(100),SCHEDULE(RUNTIME))
== Options Of Character Type ==
ALIAS(STD,INTPTR) ALIGN(STRUCT(NATURAL))
ARCH(COM) AUTODBL(NONE) DIRECTIVE($OMP)
FLAG(I,I) FLOAT(MAF,FOLD) HALT(S)
HOT(VECTOR) IEEE(NEAR) INTSIZE(4)
LANGLVL(EXTENDED) POSITION(APPENDOLD) REALSIZE(4)
REPORT(SMPLIST) SAVE(ALL) UNROLL(AUTO)
XFLAG() XLF77(NOLEADZERO,GEDIT77,NOBLANKPAD,OLDBOZ,INTARG,INTXOR,PERSISTENT,SOFTEOF)
XLF90(NOSIGNEDZERO,NOAUTODEALLOC)
>>>>> SOURCE SECTION <<<<<
** rannor === End of Compilation 2 ===
>>>>> PARALLELIZATION AND LOOP TRANSFORMATION SECTION <<<<<
REAL*8 FUNCTION rannor (ix)
! DIR_SCHEDTYPE loopId = -1 schedType = 5
@RET0 = randw(ix)
110| r = sqrt((-( 2.0000000000000000E+000 * _log(%VAL(@RET0)))))
@RET1 = randw(ix)
113| RETURN
114| END FUNCTION rannor
>>>>> COMPILATION UNIT EPILOGUE SECTION <<<<<
TOTAL UNRECOVERABLE SEVERE ERROR WARNING INFORMATIONAL
(U) (S) (E) (W) (I)
0 0 0 0 0 0
>>>>> OPTIONS SECTION <<<<<
*** Options In Effect ***
== On / Off Options ==
32 CCLINES ESCAPE I4
NOLM OBJECT SWAPOMP UNWIND
NOZEROSIZE
== Options Of Integer Type ==
FIXED(72) HOT() MAXMEM(2048)
OPTIMIZE(2) SPILLSIZE(512)
== Options of Integer and Character Type ==
SMP(OMP,AUTO,THRESHOLD(100),SCHEDULE(RUNTIME))
== Options Of Character Type ==
ALIAS(STD,INTPTR) ALIGN(STRUCT(NATURAL))
ARCH(COM) AUTODBL(NONE) DIRECTIVE($OMP)
FLAG(I,I) FLOAT(MAF,FOLD) HALT(S)
HOT(VECTOR) IEEE(NEAR) INTSIZE(4)
LANGLVL(EXTENDED) POSITION(APPENDOLD) REALSIZE(4)
REPORT(SMPLIST) SAVE(ALL) UNROLL(AUTO)
XFLAG() XLF77(NOLEADZERO,GEDIT77,NOBLANKPAD,OLDBOZ,INTARG,INTXOR,PERSISTENT,SOFTEOF)
XLF90(NOSIGNEDZERO,NOAUTODEALLOC)
>>>>> SOURCE SECTION <<<<<
** randw === End of Compilation 3 ===
>>>>> PARALLELIZATION AND LOOP TRANSFORMATION SECTION <<<<<
REAL*8 FUNCTION randw (ix)
! DIR_SCHEDTYPE loopId = -1 schedType = 5
157| @CSE1 = ix / b16
@CSE4 = (ix - @CSE1 * b16) * a
@CSE0 = @CSE4 / b16
@CSE2 = @CSE1 * a
@CSE3 = (@CSE2 + @CSE0) / b15
ix = ((@CSE4 + @CSE3) + ((@CSE2 + (@CSE0 - @CSE3 * b15)) * &
& b16 - (@CSE0 * b16 + p)))
IF ( ix < 0 ) THEN
ix = ix + p
ELSE
ix = ix
ENDIF
166| RETURN
167| END FUNCTION randw
>>>>> COMPILATION UNIT EPILOGUE SECTION <<<<<
TOTAL UNRECOVERABLE SEVERE ERROR WARNING INFORMATIONAL
(U) (S) (E) (W) (I)
0 0 0 0 0 0
>>>>> FILE TABLE SECTION <<<<<
FILE CREATION FROM
FILE NO FILENAME DATE TIME FILE LINE
0 spcol.f 09/13/04 14:41:32
>>>>> COMPILATION EPILOGUE SECTION <<<<<
FORTRAN Summary of Diagnosed Conditions
TOTAL UNRECOVERABLE SEVERE ERROR WARNING INFORMATIONAL
(U) (S) (E) (W) (I)
0 0 0 0 0 0
Source records read....................................... 168
1501-510 Compilation successful for file spcol.f.
1501-543 Object file created.
Though the loop commands are different than those from OpenMP, we do
get an analaysis of each loop, whether the loop can be parallelized
and if it can be parallelized, what variables should be private or
shared. So we can use this information to aid in detecting dependencies
and in helping to correctly insert OpneMP pragmas.
In the program analyzed above, the loop for sparse matrix multiply of matrices in compressed column storage is
do j=1,n
k1 = ic(j)
k2 = ic(j+1) - 1
do i = k1, k2
b(ir(i)) = b(ir(i)) + a(i)*x(j)
end do
end do
The automatic parallelizer didn't like this loop. Question
would it be correct to parallelize with?
!$omp parallel do
do j=1,n
k1 = ic(j)
k2 = ic(j+1) - 1
do i = k1, k2
b(ir(i)) = b(ir(i)) + a(i)*x(j)
end do
end do
Possibly not. Since elements of the b vector can be updated by
several different threads, there could be a race condition on b,
so that some elements of b could be read by Proc1
and then the location for b could be read by Proc2
before Proc1 has returned its update? Then Proc2 returns an
updated which does not include the updating performed by
Proc1? The "race condition error" may be hard to produce
in the sparse example -- because usually different elements
of b are being updated by different threads.
Let's instead set up a case to see if we can produce a race
condition error.
Consider an alternate code (dense matrix vector multiplication). This was the "optimal" matvec code that got up to 3800 "in-cache" Mflops (with different compiler options than used here. For the fast code we turned on the vectorizor by -xM and set the loop unrolling level to 2. )
do j =1,n,8
do i=1,n
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
end do
end do
Just as in the sparse code, the b vector is updated by each value of j,
so that if we parallelize the loop, processors may encounter a loop
condition. To see what can happen, initialize all values of a
as one and also set the x vector to ones. Then each value of b
should just be the number n of columns.
So the parallelized code is
!$omp parallel do
do j =1,n,8
do i=1,n
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
end do
end do
Setting n to 104 (code only works for n a multiple of 8), compile with
ifc -c openmp calgemv-s.f
and get 1400 Mflops for two processors, but it turns out that some b(i) are 96 instead of 104. .7% of b(i) were miscalculated! So YES THE RACE CONDITION DID CAUSE THE CODE TO FAIL !!
Attempt II.
!$omp parallel do reduction(+:b)
do j =1,n,8
do i=1,n
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
end do
end do
In this case the calculcation proceeds at 777 Mflops (slower) but all
the b(i) are correct. Reservation: it appears that using vectors
as reduction variables is not a part of the standard. So the implementation
may not be portable. Moreover,
the ifc compiler accepts a vector x as a reduction variable only
if x is declared with constant lenght. For example, if x is an
argument to a subroutine
subroutine my(x,n)
double precision x(*)
integer n
!$omp parallel do reduction(+:x)
do j=1,n
do i=1,n
x(j) = x(j) + 1. /(i+j)
end do
end do
return
end
The compiler will give an error.
Attempt III.
!$omp parallel do
do j =1,n,8
do i=1,n
!$omp critical
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
!$omp end critical
end do
end do
The critical pragma ensures that only on processor at a time can
access any variables that are written. Again, the answers
b(i) are all correct.
Unfortunately, the Mflop rate is 47.
Moral: ensuring that only one processor at a time has access to data can be expensive.
Discussion?
To avoid the need to synchronize, we can give each thread its own version of a variable.
!$omp parallel do shared(a) private(b)then if b is updated, other threads don't worry. What happens to b at the end of the do?
If the global b is to the sum of all the local b variables, then we can
!$omp parallel do shared(a) reduction(+:b)
It could also be the case that occasional errors don't bother us. For example, for a parallel sparse matrix vector mutliply in the middle of an iteration, it might be better to just get the iteration done and get on to the next one. Perhaps we'll need a few more iterations but if they run twice as fast? Often it might be better to do more iterations than to wait before performing another.
Here was a use of the critical section code
!$omp parallel do
do j =1,n,8
do i=1,n
!$omp critical
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
!$omp end critical
end do
end do
More efficient in this case might be
do j =1,n,8
!$omp parallel do
do i=1,n
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
end do
end do
No errors. 777 Mflops. But the inner loop in the sparse case is very
short, so might not be as good there.
Another way is with a lock routine
call omp_init_lock(foo)
!$omp parallel do
do j =1,n,8
do i=1,n
call omp_set_lock(foo)
b(i) = b(i) + a(i,j)*x(j) + a(i,j+1)*x(j+1)
+ + a(i,j+2)*x(j+2) + a(i,j+3)*x(j+3) + a(i,j+4)*x(j+4)
+ + a(i,j+5)*x(j+5) + a(i,j+6)*x(j+6) + a(i,j+7)*x(j+7)
call omp_set_unlock(foo)
end do
end do
call omp_destroy_lock(foo)
There are quite a few other OpenMP synchronizations. barrier, master, atomic, flush, ordered do loops. If you seriously want to use OpenMP you will probably end up sorting through these for the best routine. A problem is that synchronization is relatively expensive. From the results below on efficiency the fork and join synchronization of a parallel do loop, run at most 600K times per second. So it takes about 1.6e-6 seconds for a synchronization?
Generally, we expect that forking off a bunch of processes will take some time. Also having a barrier and waiting will take some time. Another common problem is that when several processors want the same bit of data, then "cache-coherence" demands that only one processor have control of the data at a time. Moreover, the data comes in cache lines. So you may not think it is controlled by one processor, but it can happen that data in one cache memory becomes unavailable to the other process.
Remember our small sample program with one parallel do? The following is a bit of an elaboration of that program. The program repeats saxpys long enough to time various parallel loops.
program saxpy
integer n , its ,k ,j, its, Mflops
real*8 a(1000000), b(1000000), alpha, divn
real*4 pretim, entim
external ftime
j = 50
do while(j .lt. 400000)
alpha = 2.d2
divn = 1.0/j
its = 1.e8/j
pretim = ftime()
!$OMP PARALLEL DO
do i=1,j
a(i) = divn
b(i) = i
end do
do k = 1,its
sum = 0.d0
!$OMP PARALLEL DO
do i=1,j
a(i) = alpha* b(i) + a(i)
sum = sum + a(i)
end do
sum = 1./sum
!$OMP PARALLEL DO
do i =1,j
a(i) = a(i)*sum
end do
end do
entim = ftime() - pretim
Mflops = its*4*j/1.e6/entim
print*,' Mflops =', Mflops,j
print*,' ParDos/time =', 2*its/entim
print*,' flops/loop =', j*2,j
j = j*1.5
end do
stop
end
The above code was compiled on the Blade Center by
>ifc -c -tpp7 -xM -openmp sax.f >ifc -o sax -openmp timd.o sax.o
When run, it gave the following output file.