I have a code in which there is a generalization of banded matrices, and
a need to compute Z = Z + A^T B using this representation. A snippet of
the code is below. The non-obvious fields of the derived type are: R1(i)
is the row subscript of the first nonzero element in column i, and R2(i)
is the index in VALUES(:,1) of the last nonzero element in column i (R2
is zero-origin, with R2(0) == 0). Z is just a Fortran rank-2 array.
The code works fine when I compile it without OpenMP. It gets the wrong
answer when I compile it with OpenMP.
I don't have extensive experience with OpenMP, which brings me to the
question: Do you see an obvious blunder, or should I report a bug to
the compiler vendor?
Thanks in advance,
Van Snyder
===========================================================================
do j = 1, zb%ncols ! Columns of Z = columns of YB
yi_1 = yb%r2(j-1)+1 ! index of 1st /=0 el. in this col of yb
yi_n = yb%r2(j) ! index of last /=0 el. in this col of yb
yr_1 = yb%r1(j) ! first row index of yb
yr_n = yr_1 + yi_n - yi_1 ! last row index of yb
mz = zb%nrows
if ( my_upper ) mz = j
!$OMP PARALLEL DO private ( xi_1, xi_n, xr_1, xr_n, cr_1, cr_n, c_n, xd, yd,
xy)
do i = 1, mz ! Rows of Z = columns of XB
! Inner product of column I of XB with column J of YB
xi_1 = xb%r2(i-1)+1 ! index of 1st /=0 el. in this col of xb
xi_n = xb%r2(i) ! index of last /=0 el. in this col of xb
xr_1 = xb%r1(i) ! first row index of xb
xr_n = xr_1 + xi_n - xi_1 ! last row index of xb
! Now work out what they have in common
cr_1 = max ( xr_1, yr_1 )
cr_n = min ( xr_n, yr_n )
c_n = cr_n - cr_1 + 1
if ( c_n <= 0 ) cycle
! Now make xd and yd point to the starts of the common parts
xd = xi_1 + cr_1 - xr_1
yd = yi_1 + cr_1 - yr_1
xy = dot( c_n, xb%values(xd,1), 1, &
& yb%values(yd,1), 1 )
! xy = dot_product( xb%values(xd:xd+c_n-1,1), &
! & yb%values(yd:yd+c_n-1,1) )
z(i,j) = z(i,j) + s * xy
end do ! i
!$OMP END PARALLEL DO
end do ! j
|