diff --git a/base/modules/auxil/psb_c_qsort_mod.f90 b/base/modules/auxil/psb_c_qsort_mod.f90 index 8b365222..6c4ceb3f 100644 --- a/base/modules/auxil/psb_c_qsort_mod.f90 +++ b/base/modules/auxil/psb_c_qsort_mod.f90 @@ -44,7 +44,6 @@ module psb_c_qsort_mod use psb_const_mod - interface psb_qsort subroutine psb_cqsort(x,ix,dir,flag) import diff --git a/base/modules/auxil/psb_d_qsort_mod.f90 b/base/modules/auxil/psb_d_qsort_mod.f90 index 4e1be1d1..4da0b840 100644 --- a/base/modules/auxil/psb_d_qsort_mod.f90 +++ b/base/modules/auxil/psb_d_qsort_mod.f90 @@ -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 diff --git a/base/modules/auxil/psb_e_qsort_mod.f90 b/base/modules/auxil/psb_e_qsort_mod.f90 index 17943bbf..09f45d45 100644 --- a/base/modules/auxil/psb_e_qsort_mod.f90 +++ b/base/modules/auxil/psb_e_qsort_mod.f90 @@ -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 diff --git a/base/modules/auxil/psb_i2_qsort_mod.f90 b/base/modules/auxil/psb_i2_qsort_mod.f90 index 944a436e..2f192a0a 100644 --- a/base/modules/auxil/psb_i2_qsort_mod.f90 +++ b/base/modules/auxil/psb_i2_qsort_mod.f90 @@ -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 diff --git a/base/modules/auxil/psb_m_qsort_mod.f90 b/base/modules/auxil/psb_m_qsort_mod.f90 index cb4c81c1..bf029065 100644 --- a/base/modules/auxil/psb_m_qsort_mod.f90 +++ b/base/modules/auxil/psb_m_qsort_mod.f90 @@ -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 diff --git a/base/modules/auxil/psb_s_qsort_mod.f90 b/base/modules/auxil/psb_s_qsort_mod.f90 index d4851fd1..a5bdb2d9 100644 --- a/base/modules/auxil/psb_s_qsort_mod.f90 +++ b/base/modules/auxil/psb_s_qsort_mod.f90 @@ -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 diff --git a/base/modules/auxil/psb_z_qsort_mod.f90 b/base/modules/auxil/psb_z_qsort_mod.f90 index 14ee0c57..2fc6baab 100644 --- a/base/modules/auxil/psb_z_qsort_mod.f90 +++ b/base/modules/auxil/psb_z_qsort_mod.f90 @@ -44,7 +44,6 @@ module psb_z_qsort_mod use psb_const_mod - interface psb_qsort subroutine psb_zqsort(x,ix,dir,flag) import diff --git a/base/modules/psb_const_mod.F90 b/base/modules/psb_const_mod.F90 index cc29f049..56134474 100644 --- a/base/modules/psb_const_mod.F90 +++ b/base/modules/psb_const_mod.F90 @@ -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 diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90 index da6e97ca..53bc34f6 100644 --- a/base/serial/impl/psb_c_coo_impl.F90 +++ b/base/serial/impl/psb_c_coo_impl.F90 @@ -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 diff --git a/base/serial/impl/psb_d_coo_impl.F90 b/base/serial/impl/psb_d_coo_impl.F90 index 08276552..030dc867 100644 --- a/base/serial/impl/psb_d_coo_impl.F90 +++ b/base/serial/impl/psb_d_coo_impl.F90 @@ -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 diff --git a/base/serial/impl/psb_s_coo_impl.F90 b/base/serial/impl/psb_s_coo_impl.F90 index 65bc5e10..fb19d523 100644 --- a/base/serial/impl/psb_s_coo_impl.F90 +++ b/base/serial/impl/psb_s_coo_impl.F90 @@ -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 diff --git a/base/serial/impl/psb_z_coo_impl.F90 b/base/serial/impl/psb_z_coo_impl.F90 index 721d2eda..ba14eeba 100644 --- a/base/serial/impl/psb_z_coo_impl.F90 +++ b/base/serial/impl/psb_z_coo_impl.F90 @@ -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 diff --git a/base/serial/sort/psb_c_qsort_impl.f90 b/base/serial/sort/psb_c_qsort_impl.f90 index 712529fc..7f33c099 100644 --- a/base/serial/sort/psb_c_qsort_impl.f90 +++ b/base/serial/sort/psb_c_qsort_impl.f90 @@ -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 diff --git a/base/serial/sort/psb_d_msort_impl.f90 b/base/serial/sort/psb_d_msort_impl.f90 index 11029818..66ad7897 100644 --- a/base/serial/sort/psb_d_msort_impl.f90 +++ b/base/serial/sort/psb_d_msort_impl.f90 @@ -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 diff --git a/base/serial/sort/psb_d_qsort_impl.f90 b/base/serial/sort/psb_d_qsort_impl.f90 index 13328188..4d6918e7 100644 --- a/base/serial/sort/psb_d_qsort_impl.f90 +++ b/base/serial/sort/psb_d_qsort_impl.f90 @@ -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 (mn) .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 diff --git a/base/serial/sort/psb_e_msort_impl.f90 b/base/serial/sort/psb_e_msort_impl.f90 index d8cd6404..b97d448a 100644 --- a/base/serial/sort/psb_e_msort_impl.f90 +++ b/base/serial/sort/psb_e_msort_impl.f90 @@ -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 diff --git a/base/serial/sort/psb_e_qsort_impl.f90 b/base/serial/sort/psb_e_qsort_impl.f90 index 9b95c78e..8be3cd78 100644 --- a/base/serial/sort/psb_e_qsort_impl.f90 +++ b/base/serial/sort/psb_e_qsort_impl.f90 @@ -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 (mn) .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 diff --git a/base/serial/sort/psb_m_msort_impl.f90 b/base/serial/sort/psb_m_msort_impl.f90 index cd99a3c5..437d1069 100644 --- a/base/serial/sort/psb_m_msort_impl.f90 +++ b/base/serial/sort/psb_m_msort_impl.f90 @@ -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 diff --git a/base/serial/sort/psb_m_qsort_impl.f90 b/base/serial/sort/psb_m_qsort_impl.f90 index ac8241f5..460bff43 100644 --- a/base/serial/sort/psb_m_qsort_impl.f90 +++ b/base/serial/sort/psb_m_qsort_impl.f90 @@ -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 (mn) .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 diff --git a/base/serial/sort/psb_s_msort_impl.f90 b/base/serial/sort/psb_s_msort_impl.f90 index dfd7508c..e3382f27 100644 --- a/base/serial/sort/psb_s_msort_impl.f90 +++ b/base/serial/sort/psb_s_msort_impl.f90 @@ -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 diff --git a/base/serial/sort/psb_s_qsort_impl.f90 b/base/serial/sort/psb_s_qsort_impl.f90 index d6e0e66e..44a46e0a 100644 --- a/base/serial/sort/psb_s_qsort_impl.f90 +++ b/base/serial/sort/psb_s_qsort_impl.f90 @@ -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 (mn) .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 diff --git a/base/serial/sort/psb_z_qsort_impl.f90 b/base/serial/sort/psb_z_qsort_impl.f90 index 7b0af1c5..a1cdb193 100644 --- a/base/serial/sort/psb_z_qsort_impl.f90 +++ b/base/serial/sort/psb_z_qsort_impl.f90 @@ -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 diff --git a/base/tools/psb_cspins.F90 b/base/tools/psb_cspins.F90 index 0f5fc9df..f523a529 100644 --- a/base/tools/psb_cspins.F90 +++ b/base/tools/psb_cspins.F90 @@ -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) diff --git a/base/tools/psb_dspins.F90 b/base/tools/psb_dspins.F90 index 06f913b5..3b7e0e80 100644 --- a/base/tools/psb_dspins.F90 +++ b/base/tools/psb_dspins.F90 @@ -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) diff --git a/base/tools/psb_sspins.F90 b/base/tools/psb_sspins.F90 index 56ef9c97..9781eaae 100644 --- a/base/tools/psb_sspins.F90 +++ b/base/tools/psb_sspins.F90 @@ -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) diff --git a/base/tools/psb_zspins.F90 b/base/tools/psb_zspins.F90 index d483f198..36b0b5a5 100644 --- a/base/tools/psb_zspins.F90 +++ b/base/tools/psb_zspins.F90 @@ -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) diff --git a/test/omp/psb_tomp.F90 b/test/omp/psb_tomp.F90 index fda08f4e..79097ca8 100644 --- a/test/omp/psb_tomp.F90 +++ b/test/omp/psb_tomp.F90 @@ -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 ! diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index fb7af68c..40e3358d 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -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