fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fix bug in COO impl. OpenMP version.

Browse Source
omp-threadsafe
Salvatore Filippone 3 years ago
parent d0a5ff9893
commit 7ee9001fde
8 changed files with 5072 additions and 1748 deletions
Show Stats Download Patch File Download Diff File

1501
base/serial/impl/psb_c_coo_impl.F90
View File

File diff suppressed because it is too large Load Diff

204
base/serial/impl/psb_c_csr_impl.f90
Unescape Escape View File

@ -2851,6 +2851,13 @@ subroutine psb_c_cp_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='c_cp_csr_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s,j
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem,maxthreads
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -2894,17 +2901,92 @@ subroutine psb_c_cp_csr_from_coo(a,b,info)
endif
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$ maxthreads = omp_get_max_threads()
!$ allocate(sum(maxthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP PARALLEL default(none) &
!$OMP shared(nza,itemp,a,nthreads,sum,nr) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$OMP BARRIER
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
@ -3020,6 +3102,13 @@ subroutine psb_c_mv_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='mv_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -3045,17 +3134,88 @@ subroutine psb_c_mv_csr_from_coo(a,b,info)
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$OMP PARALLEL default(none) &
!$OMP shared(sum,nthreads,nr,a,itemp,nza) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$ allocate(sum(nthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
end subroutine psb_c_mv_csr_from_coo

1501
base/serial/impl/psb_d_coo_impl.F90
View File

File diff suppressed because it is too large Load Diff

204
base/serial/impl/psb_d_csr_impl.f90
Unescape Escape View File

@ -2851,6 +2851,13 @@ subroutine psb_d_cp_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='d_cp_csr_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s,j
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem,maxthreads
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -2894,17 +2901,92 @@ subroutine psb_d_cp_csr_from_coo(a,b,info)
endif
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$ maxthreads = omp_get_max_threads()
!$ allocate(sum(maxthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP PARALLEL default(none) &
!$OMP shared(nza,itemp,a,nthreads,sum,nr) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$OMP BARRIER
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
@ -3020,6 +3102,13 @@ subroutine psb_d_mv_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='mv_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -3045,17 +3134,88 @@ subroutine psb_d_mv_csr_from_coo(a,b,info)
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$OMP PARALLEL default(none) &
!$OMP shared(sum,nthreads,nr,a,itemp,nza) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$ allocate(sum(nthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
end subroutine psb_d_mv_csr_from_coo

1501
base/serial/impl/psb_s_coo_impl.F90
View File

File diff suppressed because it is too large Load Diff

204
base/serial/impl/psb_s_csr_impl.f90
Unescape Escape View File

@ -2851,6 +2851,13 @@ subroutine psb_s_cp_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='s_cp_csr_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s,j
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem,maxthreads
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -2894,17 +2901,92 @@ subroutine psb_s_cp_csr_from_coo(a,b,info)
endif
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$ maxthreads = omp_get_max_threads()
!$ allocate(sum(maxthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP PARALLEL default(none) &
!$OMP shared(nza,itemp,a,nthreads,sum,nr) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$OMP BARRIER
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
@ -3020,6 +3102,13 @@ subroutine psb_s_mv_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='mv_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -3045,17 +3134,88 @@ subroutine psb_s_mv_csr_from_coo(a,b,info)
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$OMP PARALLEL default(none) &
!$OMP shared(sum,nthreads,nr,a,itemp,nza) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$ allocate(sum(nthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
end subroutine psb_s_mv_csr_from_coo

1501
base/serial/impl/psb_z_coo_impl.F90
View File

File diff suppressed because it is too large Load Diff

204
base/serial/impl/psb_z_csr_impl.f90
Unescape Escape View File

@ -2851,6 +2851,13 @@ subroutine psb_z_cp_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='z_cp_csr_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s,j
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem,maxthreads
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -2894,17 +2901,92 @@ subroutine psb_z_cp_csr_from_coo(a,b,info)
endif
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$ maxthreads = omp_get_max_threads()
!$ allocate(sum(maxthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP PARALLEL default(none) &
!$OMP shared(nza,itemp,a,nthreads,sum,nr) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$OMP BARRIER
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
@ -3020,6 +3102,13 @@ subroutine psb_z_mv_csr_from_coo(a,b,info)
integer(psb_ipk_), Parameter :: maxtry=8
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name='mv_from_coo'
logical :: use_openmp = .false.
!$ integer(psb_ipk_), allocatable :: sum(:)
!$ integer(psb_ipk_) :: first_idx,last_idx,work,ithread,nthreads,s
!$ integer(psb_ipk_) :: nxt_val,old_val,saved_elem
!$ use_openmp = .true.
info = psb_success_
debug_unit = psb_get_debug_unit()
@ -3045,17 +3134,88 @@ subroutine psb_z_mv_csr_from_coo(a,b,info)
a%irp(:) = 0
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
if (use_openmp) then
!$OMP PARALLEL default(none) &
!$OMP shared(sum,nthreads,nr,a,itemp,nza) &
!$OMP private(ithread,work,first_idx,last_idx,s,saved_elem,nxt_val,old_val)
!$OMP DO schedule(STATIC) &
!$OMP private(k,i)
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
!$OMP END DO
!$OMP SINGLE
!$ nthreads = omp_get_num_threads()
!$ allocate(sum(nthreads+1))
!$ sum(:) = 0
!$ sum(1) = 1
!$OMP END SINGLE
!$ ithread = omp_get_thread_num()
!$ work = nr/nthreads
!$ if (ithread < MOD(nr,nthreads)) then
!$ work = work + 1
!$ first_idx = ithread*work + 1
!$ else
!$ first_idx = ithread*work + MOD(nr,nthreads) + 1
!$ end if
!$ last_idx = first_idx + work - 1
!$ s = 0
!$ do i=first_idx,last_idx
!$ s = s + a%irp(i)
!$ end do
!$ if (work > 0) then
!$ sum(ithread+2) = s
!$ end if
!$OMP BARRIER
!$OMP SINGLE
!$ do i=2,nthreads+1
!$ sum(i) = sum(i) + sum(i-1)
!$ end do
!$OMP END SINGLE
!$ if (work > 0) then
!$ saved_elem = a%irp(first_idx)
!$ end if
!$ if (ithread == 0) then
!$ a%irp(1) = 1
!$ end if
!$ if (work > 0) then
!$ old_val = a%irp(first_idx+1)
!$ a%irp(first_idx+1) = saved_elem + sum(ithread+1)
!$ end if
!$ do i=first_idx+2,last_idx+1
!$ nxt_val = a%irp(i)
!$ a%irp(i) = a%irp(i-1) + old_val
!$ old_val = nxt_val
!$ end do
!$OMP END PARALLEL
else
do k=1,nza
i = itemp(k)
a%irp(i) = a%irp(i) + 1
end do
ip = 1
do i=1,nr
ncl = a%irp(i)
a%irp(i) = ip
ip = ip + ncl
end do
a%irp(nr+1) = ip
end if
call a%set_host()
end subroutine psb_z_mv_csr_from_coo

Write Preview
Loading…
Cancel
Save