Define new options for BSRCH, clean interface

omp-threadsafe
sfilippone 2 years ago
parent 40cc78854a
commit 7e5dc20e03

@ -44,7 +44,6 @@ module psb_c_qsort_mod
use psb_const_mod
interface psb_qsort
subroutine psb_cqsort(x,ix,dir,flag)
import

@ -43,14 +43,13 @@
module psb_d_qsort_mod
use psb_const_mod
interface psb_bsrch
function psb_dbsrch(key,n,v) result(ipos)
function psb_dbsrch(key,n,v,dir,find) result(ipos)
import
integer(psb_ipk_) :: ipos, n
real(psb_dpk_) :: key
real(psb_dpk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
end function psb_dbsrch
end interface psb_bsrch

@ -43,14 +43,13 @@
module psb_e_qsort_mod
use psb_const_mod
interface psb_bsrch
function psb_ebsrch(key,n,v) result(ipos)
function psb_ebsrch(key,n,v,dir,find) result(ipos)
import
integer(psb_ipk_) :: ipos, n
integer(psb_epk_) :: key
integer(psb_epk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
end function psb_ebsrch
end interface psb_bsrch

@ -43,14 +43,13 @@
module psb_i2_qsort_mod
use psb_const_mod
interface psb_bsrch
function psb_i2bsrch(key,n,v) result(ipos)
function psb_i2bsrch(key,n,v,dir,find) result(ipos)
import
integer(psb_ipk_) :: ipos, n
integer(psb_i2pk_) :: key
integer(psb_i2pk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
end function psb_i2bsrch
end interface psb_bsrch

@ -43,14 +43,13 @@
module psb_m_qsort_mod
use psb_const_mod
interface psb_bsrch
function psb_mbsrch(key,n,v) result(ipos)
function psb_mbsrch(key,n,v,dir,find) result(ipos)
import
integer(psb_ipk_) :: ipos, n
integer(psb_mpk_) :: key
integer(psb_mpk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
end function psb_mbsrch
end interface psb_bsrch

@ -43,14 +43,13 @@
module psb_s_qsort_mod
use psb_const_mod
interface psb_bsrch
function psb_sbsrch(key,n,v) result(ipos)
function psb_sbsrch(key,n,v,dir,find) result(ipos)
import
integer(psb_ipk_) :: ipos, n
real(psb_spk_) :: key
real(psb_spk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
end function psb_sbsrch
end interface psb_bsrch

@ -44,7 +44,6 @@ module psb_z_qsort_mod
use psb_const_mod
interface psb_qsort
subroutine psb_zqsort(x,ix,dir,flag)
import

@ -185,13 +185,17 @@ module psb_const_mod
! The up/down constant are defined in pairs having
! opposite values. We make use of this fact in the heapsort routine.
!
integer(psb_ipk_), parameter :: psb_sort_up_ = 1, psb_sort_down_ = -1
integer(psb_ipk_), parameter :: psb_lsort_up_ = 2, psb_lsort_down_ = -2
integer(psb_ipk_), parameter :: psb_asort_up_ = 3, psb_asort_down_ = -3
integer(psb_ipk_), parameter :: psb_alsort_up_ = 4, psb_alsort_down_ = -4
integer(psb_ipk_), parameter :: psb_sort_ovw_idx_ = 0, psb_sort_keep_idx_ = 1
integer(psb_ipk_), parameter :: psb_heap_resize = 200
integer(psb_ipk_), parameter :: psb_sort_up_ = 1, psb_sort_down_ = -1
integer(psb_ipk_), parameter :: psb_lsort_up_ = 2, psb_lsort_down_ = -2
integer(psb_ipk_), parameter :: psb_asort_up_ = 3, psb_asort_down_ = -3
integer(psb_ipk_), parameter :: psb_alsort_up_ = 4, psb_alsort_down_ = -4
integer(psb_ipk_), parameter :: psb_sort_ovw_idx_ = 0, psb_sort_keep_idx_ = 1
integer(psb_ipk_), parameter :: psb_heap_resize = 200
integer(psb_ipk_), parameter :: psb_find_any_ = 0
integer(psb_ipk_), parameter :: psb_find_first_ge_ = 1
integer(psb_ipk_), parameter :: psb_find_last_le_ = 2
!
! Sparse matrix constants

@ -3674,6 +3674,7 @@ subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name = 'psb_fixcoo'
logical :: srt_inp, use_buffers
real(psb_dpk_) :: t0, t1
#if defined(OPENMP)
integer(psb_ipk_) :: work,idxstart,idxend,first_elem,last_elem,s,nthreads,ithread
integer(psb_ipk_) :: saved_elem,old_val,nxt_val,err,act_row,act_col,maxthreads
@ -3757,7 +3758,7 @@ subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
! Here, starting from 'iaux', we apply a fixing in order to obtain the starting
! index for each row. We do the same on 'kaux'
!$OMP PARALLEL default(none) &
!$OMP shared(idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP shared(t0,t1,idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP private(s,i,j,k,ithread,idxstart,idxend,work,nxt_val,old_val,saved_elem, &
!$OMP first_elem,last_elem,nzl,iret,act_row) reduction(max: info)
@ -3782,7 +3783,10 @@ subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
!write(0,*) 'fix_coo_inner: trying with exscan'
call psi_exscan(nr+1,iaux,info,shift=ione)
!$OMP BARRIER
!$OMP SINGLE
t0 = omp_get_wtime()
!$OMP END SINGLE
! ------------------ Sorting and buffers -------------------
! Let's use an auxiliary buffer, 'idxaux', to get indices leaving
@ -3803,7 +3807,10 @@ subroutine psb_c_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
end do
!$OMP BARRIER
!$OMP SINGLE
t1 = omp_get_wtime()
write(0,*) 'Srt&Cpy :',t1-t0
!$OMP END SINGLE
! Let's sort column indices and values. After that we will store
! the number of unique values in 'kaux'
do j=idxstart,idxend

@ -3674,6 +3674,7 @@ subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name = 'psb_fixcoo'
logical :: srt_inp, use_buffers
real(psb_dpk_) :: t0, t1
#if defined(OPENMP)
integer(psb_ipk_) :: work,idxstart,idxend,first_elem,last_elem,s,nthreads,ithread
integer(psb_ipk_) :: saved_elem,old_val,nxt_val,err,act_row,act_col,maxthreads
@ -3757,7 +3758,7 @@ subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
! Here, starting from 'iaux', we apply a fixing in order to obtain the starting
! index for each row. We do the same on 'kaux'
!$OMP PARALLEL default(none) &
!$OMP shared(idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP shared(t0,t1,idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP private(s,i,j,k,ithread,idxstart,idxend,work,nxt_val,old_val,saved_elem, &
!$OMP first_elem,last_elem,nzl,iret,act_row) reduction(max: info)
@ -3782,7 +3783,10 @@ subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
!write(0,*) 'fix_coo_inner: trying with exscan'
call psi_exscan(nr+1,iaux,info,shift=ione)
!$OMP BARRIER
!$OMP SINGLE
t0 = omp_get_wtime()
!$OMP END SINGLE
! ------------------ Sorting and buffers -------------------
! Let's use an auxiliary buffer, 'idxaux', to get indices leaving
@ -3803,7 +3807,10 @@ subroutine psb_d_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
end do
!$OMP BARRIER
!$OMP SINGLE
t1 = omp_get_wtime()
write(0,*) 'Srt&Cpy :',t1-t0
!$OMP END SINGLE
! Let's sort column indices and values. After that we will store
! the number of unique values in 'kaux'
do j=idxstart,idxend

@ -3674,6 +3674,7 @@ subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name = 'psb_fixcoo'
logical :: srt_inp, use_buffers
real(psb_dpk_) :: t0, t1
#if defined(OPENMP)
integer(psb_ipk_) :: work,idxstart,idxend,first_elem,last_elem,s,nthreads,ithread
integer(psb_ipk_) :: saved_elem,old_val,nxt_val,err,act_row,act_col,maxthreads
@ -3757,7 +3758,7 @@ subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
! Here, starting from 'iaux', we apply a fixing in order to obtain the starting
! index for each row. We do the same on 'kaux'
!$OMP PARALLEL default(none) &
!$OMP shared(idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP shared(t0,t1,idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP private(s,i,j,k,ithread,idxstart,idxend,work,nxt_val,old_val,saved_elem, &
!$OMP first_elem,last_elem,nzl,iret,act_row) reduction(max: info)
@ -3782,7 +3783,10 @@ subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
!write(0,*) 'fix_coo_inner: trying with exscan'
call psi_exscan(nr+1,iaux,info,shift=ione)
!$OMP BARRIER
!$OMP SINGLE
t0 = omp_get_wtime()
!$OMP END SINGLE
! ------------------ Sorting and buffers -------------------
! Let's use an auxiliary buffer, 'idxaux', to get indices leaving
@ -3803,7 +3807,10 @@ subroutine psb_s_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
end do
!$OMP BARRIER
!$OMP SINGLE
t1 = omp_get_wtime()
write(0,*) 'Srt&Cpy :',t1-t0
!$OMP END SINGLE
! Let's sort column indices and values. After that we will store
! the number of unique values in 'kaux'
do j=idxstart,idxend

@ -3674,6 +3674,7 @@ subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
integer(psb_ipk_) :: debug_level, debug_unit
character(len=20) :: name = 'psb_fixcoo'
logical :: srt_inp, use_buffers
real(psb_dpk_) :: t0, t1
#if defined(OPENMP)
integer(psb_ipk_) :: work,idxstart,idxend,first_elem,last_elem,s,nthreads,ithread
integer(psb_ipk_) :: saved_elem,old_val,nxt_val,err,act_row,act_col,maxthreads
@ -3757,7 +3758,7 @@ subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
! Here, starting from 'iaux', we apply a fixing in order to obtain the starting
! index for each row. We do the same on 'kaux'
!$OMP PARALLEL default(none) &
!$OMP shared(idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP shared(t0,t1,idxaux,ia,ja,val,ias,jas,vs,nthreads,sum,nr,nc,nzin,iaux,kaux,dupl,err) &
!$OMP private(s,i,j,k,ithread,idxstart,idxend,work,nxt_val,old_val,saved_elem, &
!$OMP first_elem,last_elem,nzl,iret,act_row) reduction(max: info)
@ -3782,7 +3783,10 @@ subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
!write(0,*) 'fix_coo_inner: trying with exscan'
call psi_exscan(nr+1,iaux,info,shift=ione)
!$OMP BARRIER
!$OMP SINGLE
t0 = omp_get_wtime()
!$OMP END SINGLE
! ------------------ Sorting and buffers -------------------
! Let's use an auxiliary buffer, 'idxaux', to get indices leaving
@ -3803,7 +3807,10 @@ subroutine psb_z_fix_coo_inner_rowmajor(nr,nc,nzin,dupl,ia,ja,val,iaux,nzout,inf
end do
!$OMP BARRIER
!$OMP SINGLE
t1 = omp_get_wtime()
write(0,*) 'Srt&Cpy :',t1-t0
!$OMP END SINGLE
! Let's sort column indices and values. After that we will store
! the number of unique values in 'kaux'
do j=idxstart,idxend

@ -40,6 +40,7 @@
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_cqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_cqsort
use psb_error_mod

@ -76,64 +76,6 @@ subroutine psb_dmsort_u(x,nout,dir)
return
end subroutine psb_dmsort_u
function psb_dbsrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_dbsrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_dpk_) :: key
real(psb_dpk_) :: v(:)
integer(psb_ipk_) :: lb, ub, m, i
ipos = -1
if (n<5) then
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end if
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
lb = ub + 1
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
return
end function psb_dbsrch
function psb_dssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_dssrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_dpk_) :: key
real(psb_dpk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_dssrch
subroutine psb_dmsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_dmsort
use psb_error_mod

@ -40,6 +40,128 @@
! Data Structures and Algorithms
! Addison-Wesley
!
function psb_dbsrch(key,n,v,dir,find) result(ipos)
use psb_sort_mod, psb_protect_name => psb_dbsrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_dpk_) :: key
real(psb_dpk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
integer(psb_ipk_) :: lb, ub, m, i, k, dir_, find_
if (present(dir)) then
dir_ = dir
else
dir_ = psb_sort_up_
end if
if (present(find)) then
find_ = find
else
find_ = psb_find_any_
end if
ipos = -1
if (dir_ == psb_sort_up_) then
if (n<=5) then
do m=1,n
if (key == v(m)) then
ipos = m
exit
end if
enddo
else
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
exit
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
end if
select case(find_)
case (psb_find_any_ )
! do nothing
case (psb_find_last_le_ )
if ((m>n) .or. (m<1)) then
m = n
do while (m>=1)
if (v(m)<=key) then
ipos = m
exit
end if
m = m - 1
end do
else
do while (m<n)
if (v(m)<=key) then
m=m+1
else
exit
end if
end do
end if
case (psb_find_first_ge_ )
if ((m>n) .or. (m<1)) then
m = 1
do while (m<=n)
if (v(m)>=key) then
ipos = m
exit
end if
m = m + 1
end do
else
do while (m>n)
if (v(m)>=key) then
m=m-1
else
exit
end if
end do
end if
case default
write(0,*) 'Wrong FIND'
end select
else if (dir_ == psb_sort_down_) then
write(0,*) ' bsrch on sort down not implemented'
else
write(0,*) ' bsrch wrong DIR ',dir_,psb_sort_up_,psb_sort_down_
end if
return
end function psb_dbsrch
function psb_dssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_dssrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_dpk_) :: key
real(psb_dpk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_dssrch
subroutine psb_dqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_dqsort
use psb_error_mod

@ -131,64 +131,6 @@ subroutine psb_emsort_u(x,nout,dir)
return
end subroutine psb_emsort_u
function psb_ebsrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_ebsrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_epk_) :: key
integer(psb_epk_) :: v(:)
integer(psb_ipk_) :: lb, ub, m, i
ipos = -1
if (n<5) then
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end if
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
lb = ub + 1
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
return
end function psb_ebsrch
function psb_essrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_essrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_epk_) :: key
integer(psb_epk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_essrch
subroutine psb_emsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_emsort
use psb_error_mod

@ -40,6 +40,128 @@
! Data Structures and Algorithms
! Addison-Wesley
!
function psb_ebsrch(key,n,v,dir,find) result(ipos)
use psb_sort_mod, psb_protect_name => psb_ebsrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_epk_) :: key
integer(psb_epk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
integer(psb_ipk_) :: lb, ub, m, i, k, dir_, find_
if (present(dir)) then
dir_ = dir
else
dir_ = psb_sort_up_
end if
if (present(find)) then
find_ = find
else
find_ = psb_find_any_
end if
ipos = -1
if (dir_ == psb_sort_up_) then
if (n<=5) then
do m=1,n
if (key == v(m)) then
ipos = m
exit
end if
enddo
else
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
exit
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
end if
select case(find_)
case (psb_find_any_ )
! do nothing
case (psb_find_last_le_ )
if ((m>n) .or. (m<1)) then
m = n
do while (m>=1)
if (v(m)<=key) then
ipos = m
exit
end if
m = m - 1
end do
else
do while (m<n)
if (v(m)<=key) then
m=m+1
else
exit
end if
end do
end if
case (psb_find_first_ge_ )
if ((m>n) .or. (m<1)) then
m = 1
do while (m<=n)
if (v(m)>=key) then
ipos = m
exit
end if
m = m + 1
end do
else
do while (m>n)
if (v(m)>=key) then
m=m-1
else
exit
end if
end do
end if
case default
write(0,*) 'Wrong FIND'
end select
else if (dir_ == psb_sort_down_) then
write(0,*) ' bsrch on sort down not implemented'
else
write(0,*) ' bsrch wrong DIR ',dir_,psb_sort_up_,psb_sort_down_
end if
return
end function psb_ebsrch
function psb_essrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_essrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_epk_) :: key
integer(psb_epk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_essrch
subroutine psb_eqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_eqsort
use psb_error_mod

@ -131,64 +131,6 @@ subroutine psb_mmsort_u(x,nout,dir)
return
end subroutine psb_mmsort_u
function psb_mbsrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_mbsrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_mpk_) :: key
integer(psb_mpk_) :: v(:)
integer(psb_ipk_) :: lb, ub, m, i
ipos = -1
if (n<5) then
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end if
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
lb = ub + 1
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
return
end function psb_mbsrch
function psb_mssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_mssrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_mpk_) :: key
integer(psb_mpk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_mssrch
subroutine psb_mmsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_mmsort
use psb_error_mod

@ -40,6 +40,128 @@
! Data Structures and Algorithms
! Addison-Wesley
!
function psb_mbsrch(key,n,v,dir,find) result(ipos)
use psb_sort_mod, psb_protect_name => psb_mbsrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_mpk_) :: key
integer(psb_mpk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
integer(psb_ipk_) :: lb, ub, m, i, k, dir_, find_
if (present(dir)) then
dir_ = dir
else
dir_ = psb_sort_up_
end if
if (present(find)) then
find_ = find
else
find_ = psb_find_any_
end if
ipos = -1
if (dir_ == psb_sort_up_) then
if (n<=5) then
do m=1,n
if (key == v(m)) then
ipos = m
exit
end if
enddo
else
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
exit
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
end if
select case(find_)
case (psb_find_any_ )
! do nothing
case (psb_find_last_le_ )
if ((m>n) .or. (m<1)) then
m = n
do while (m>=1)
if (v(m)<=key) then
ipos = m
exit
end if
m = m - 1
end do
else
do while (m<n)
if (v(m)<=key) then
m=m+1
else
exit
end if
end do
end if
case (psb_find_first_ge_ )
if ((m>n) .or. (m<1)) then
m = 1
do while (m<=n)
if (v(m)>=key) then
ipos = m
exit
end if
m = m + 1
end do
else
do while (m>n)
if (v(m)>=key) then
m=m-1
else
exit
end if
end do
end if
case default
write(0,*) 'Wrong FIND'
end select
else if (dir_ == psb_sort_down_) then
write(0,*) ' bsrch on sort down not implemented'
else
write(0,*) ' bsrch wrong DIR ',dir_,psb_sort_up_,psb_sort_down_
end if
return
end function psb_mbsrch
function psb_mssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_mssrch
implicit none
integer(psb_ipk_) :: ipos, n
integer(psb_mpk_) :: key
integer(psb_mpk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_mssrch
subroutine psb_mqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_mqsort
use psb_error_mod

@ -76,64 +76,6 @@ subroutine psb_smsort_u(x,nout,dir)
return
end subroutine psb_smsort_u
function psb_sbsrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_sbsrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_spk_) :: key
real(psb_spk_) :: v(:)
integer(psb_ipk_) :: lb, ub, m, i
ipos = -1
if (n<5) then
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end if
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
lb = ub + 1
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
return
end function psb_sbsrch
function psb_sssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_sssrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_spk_) :: key
real(psb_spk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_sssrch
subroutine psb_smsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_smsort
use psb_error_mod

@ -40,6 +40,128 @@
! Data Structures and Algorithms
! Addison-Wesley
!
function psb_sbsrch(key,n,v,dir,find) result(ipos)
use psb_sort_mod, psb_protect_name => psb_sbsrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_spk_) :: key
real(psb_spk_) :: v(:)
integer(psb_ipk_), optional :: dir, find
integer(psb_ipk_) :: lb, ub, m, i, k, dir_, find_
if (present(dir)) then
dir_ = dir
else
dir_ = psb_sort_up_
end if
if (present(find)) then
find_ = find
else
find_ = psb_find_any_
end if
ipos = -1
if (dir_ == psb_sort_up_) then
if (n<=5) then
do m=1,n
if (key == v(m)) then
ipos = m
exit
end if
enddo
else
lb = 1
ub = n
do while (lb.le.ub)
m = (lb+ub)/2
if (key.eq.v(m)) then
ipos = m
exit
else if (key < v(m)) then
ub = m-1
else
lb = m + 1
end if
enddo
end if
select case(find_)
case (psb_find_any_ )
! do nothing
case (psb_find_last_le_ )
if ((m>n) .or. (m<1)) then
m = n
do while (m>=1)
if (v(m)<=key) then
ipos = m
exit
end if
m = m - 1
end do
else
do while (m<n)
if (v(m)<=key) then
m=m+1
else
exit
end if
end do
end if
case (psb_find_first_ge_ )
if ((m>n) .or. (m<1)) then
m = 1
do while (m<=n)
if (v(m)>=key) then
ipos = m
exit
end if
m = m + 1
end do
else
do while (m>n)
if (v(m)>=key) then
m=m-1
else
exit
end if
end do
end if
case default
write(0,*) 'Wrong FIND'
end select
else if (dir_ == psb_sort_down_) then
write(0,*) ' bsrch on sort down not implemented'
else
write(0,*) ' bsrch wrong DIR ',dir_,psb_sort_up_,psb_sort_down_
end if
return
end function psb_sbsrch
function psb_sssrch(key,n,v) result(ipos)
use psb_sort_mod, psb_protect_name => psb_sssrch
implicit none
integer(psb_ipk_) :: ipos, n
real(psb_spk_) :: key
real(psb_spk_) :: v(:)
integer(psb_ipk_) :: i
ipos = -1
do i=1,n
if (key.eq.v(i)) then
ipos = i
return
end if
enddo
return
end function psb_sssrch
subroutine psb_sqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_sqsort
use psb_error_mod

@ -40,6 +40,7 @@
! Data Structures and Algorithms
! Addison-Wesley
!
subroutine psb_zqsort(x,ix,dir,flag)
use psb_sort_mod, psb_protect_name => psb_zqsort
use psb_error_mod

@ -135,7 +135,7 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l_ins(ja(1:nz),jla(1:nz),info,&
@ -198,9 +198,18 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
& a_err='allocate',i_err=(/info/))
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
if (local_) then
#if defined(OPENMP)
!$omp workshare
#endif
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
#if defined(OPENMP)
!$omp end workshare
#endif
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
@ -210,7 +219,7 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='a%csput')
goto 9999
!goto 9999
end if
if (a%is_remote_build()) then
nnl = count(ila(1:nz)<0)
@ -229,8 +238,12 @@ subroutine psb_cspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
end if
end if
#if defined(OPENMP)
!$omp end parallel
#endif
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)

@ -135,7 +135,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l_ins(ja(1:nz),jla(1:nz),info,&
@ -198,9 +198,18 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
& a_err='allocate',i_err=(/info/))
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
if (local_) then
#if defined(OPENMP)
!$omp workshare
#endif
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
#if defined(OPENMP)
!$omp end workshare
#endif
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
@ -210,7 +219,7 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='a%csput')
goto 9999
!goto 9999
end if
if (a%is_remote_build()) then
nnl = count(ila(1:nz)<0)
@ -229,8 +238,12 @@ subroutine psb_dspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
end if
end if
#if defined(OPENMP)
!$omp end parallel
#endif
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)

@ -135,7 +135,7 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l_ins(ja(1:nz),jla(1:nz),info,&
@ -198,9 +198,18 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
& a_err='allocate',i_err=(/info/))
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
if (local_) then
#if defined(OPENMP)
!$omp workshare
#endif
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
#if defined(OPENMP)
!$omp end workshare
#endif
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
@ -210,7 +219,7 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='a%csput')
goto 9999
!goto 9999
end if
if (a%is_remote_build()) then
nnl = count(ila(1:nz)<0)
@ -229,8 +238,12 @@ subroutine psb_sspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
end if
end if
#if defined(OPENMP)
!$omp end parallel
#endif
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)

@ -135,7 +135,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l_ins(ja(1:nz),jla(1:nz),info,&
@ -198,9 +198,18 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
& a_err='allocate',i_err=(/info/))
goto 9999
end if
#if defined(OPENMP)
!$omp parallel private(ila,jla,nrow,ncol,nnl,k)
#endif
if (local_) then
#if defined(OPENMP)
!$omp workshare
#endif
ila(1:nz) = ia(1:nz)
jla(1:nz) = ja(1:nz)
#if defined(OPENMP)
!$omp end workshare
#endif
else
call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.)
if (info == 0) call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info,&
@ -210,7 +219,7 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (info /= psb_success_) then
info=psb_err_from_subroutine_
call psb_errpush(info,name,a_err='a%csput')
goto 9999
!goto 9999
end if
if (a%is_remote_build()) then
nnl = count(ila(1:nz)<0)
@ -229,8 +238,12 @@ subroutine psb_zspins(nz,ia,ja,val,a,desc_a,info,rebuild,local)
if (k /= nnl) write(0,*) name,' Wrong conversion?',k,nnl
call a%rmta%csput(nnl,lila,ljla,lval,1_psb_lpk_,desc_a%get_global_rows(),&
& 1_psb_lpk_,desc_a%get_global_rows(),info)
end if
end if
end if
#if defined(OPENMP)
!$omp end parallel
#endif
else
info = psb_err_invalid_cd_state_
call psb_errpush(info,name)

@ -736,6 +736,30 @@ program psb_d_pde3d
write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_
write(*,*) 'This is the ',trim(name),' sample program'
end if
#if 0
block
integer(psb_ipk_), parameter :: ntv=10
integer(psb_ipk_) :: itv(ntv+1),i
itv(:) = 0
do i=1,ntv
itv(i) = 2 + mod(i,2)
end do
write(0,*) 'ITV before : ',itv(:)
call psi_exscan(ntv,itv,info)
write(0,*) 'ITV after : ',itv(:)
itv(:) = 0
do i=1,ntv
itv(i) = 2 + mod(i,2)
end do
write(0,*) 'ITV before 1: ',itv(:)
call psi_exscan(ntv,itv,info,shift=ione)
write(0,*) 'ITV after 1: ',itv(:)
! call a%print('a.mtx',head='Test')
end block
!!$
!!$ call psb_exit(ctxt)
!!$ stop
#endif
!
! get parameters
!
@ -756,6 +780,7 @@ program psb_d_pde3d
end if
if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2
if (iam == psb_root_) write(psb_out_unit,'(" ")')
call a%print('a.mtx',head='Test')
!
! prepare the preconditioner.
!
@ -858,7 +883,6 @@ program psb_d_pde3d
write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt()
end if
!
! cleanup storage and exit
!

@ -2,10 +2,10 @@
BICGSTAB Iterative method BICGSTAB CGS BICG BICGSTABL RGMRES FCG CGR
BJAC Preconditioner NONE DIAG BJAC
CSR Storage format for matrix A: CSR COO
040 Domain size (acutal system is this**3 (pde3d) or **2 (pde2d) )
140 Domain size (acutal system is this**3 (pde3d) or **2 (pde2d) )
3 Partition: 1 BLOCK 3 3D
2 Stopping criterion 1 2
0100 MAXIT
0200 MAXIT
05 ITRACE
002 IRST restart for RGMRES and BiCGSTABL
INVK Block Solver ILU,ILUT,INVK,AINVT,AORTH

Loading…
Cancel
Save