diff --git a/base/modules/psblas/psb_d_psblas_mod.F90 b/base/modules/psblas/psb_d_psblas_mod.F90 index b200bc8a..cc791491 100644 --- a/base/modules/psblas/psb_d_psblas_mod.F90 +++ b/base/modules/psblas/psb_d_psblas_mod.F90 @@ -32,6 +32,7 @@ module psb_d_psblas_mod use psb_desc_mod, only : psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_ use psb_d_vect_mod, only : psb_d_vect_type + use psb_d_multivect_mod, only : psb_d_multivect_type use psb_d_mat_mod, only : psb_dspmat_type interface psb_gedot @@ -63,6 +64,26 @@ module psb_d_psblas_mod integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: global end function psb_ddot + function psb_ddot_multivect(x, y, desc_a,info,global) result(res) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_multivect_type, psb_dspmat_type + real(psb_dpk_), dimension(:), allocatable :: res + type(psb_d_multivect_type), intent(inout) :: x, y + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_ddot_multivect + function psb_ddot_mvect_vect(x, y, desc_a,info,global) result(res) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_multivect_type, psb_dspmat_type, & + & psb_d_vect_type + real(psb_dpk_), dimension(:), allocatable :: res + type(psb_d_multivect_type), intent(inout) :: x + type(psb_d_vect_type), intent(inout) :: y + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_ddot_mvect_vect end interface diff --git a/base/modules/serial/psb_d_base_vect_mod.F90 b/base/modules/serial/psb_d_base_vect_mod.F90 index f07b5aed..e32f4e70 100644 --- a/base/modules/serial/psb_d_base_vect_mod.F90 +++ b/base/modules/serial/psb_d_base_vect_mod.F90 @@ -149,7 +149,8 @@ module psb_d_base_vect_mod ! procedure, pass(x) :: dot_v => d_base_dot_v procedure, pass(x) :: dot_a => d_base_dot_a - generic, public :: dot => dot_v, dot_a + procedure, pass(x) :: dot_a2 => d_base_dot_a2 + generic, public :: dot => dot_v, dot_a, dot_a2 procedure, pass(y) :: axpby_v => d_base_axpby_v procedure, pass(y) :: axpby_a => d_base_axpby_a procedure, pass(z) :: axpby_v2 => d_base_axpby_v2 @@ -1017,6 +1018,34 @@ contains end function d_base_dot_a + ! + ! Base workhorse is good old BLAS1 + ! + ! + !> Function base_dot_a + !! \memberof psb_d_base_vect_type + !! \brief Dot product by a normal array + !! \param n Number of entries to be considered + !! \param y(:,:) The matrix to be multiplied by + !! + function d_base_dot_a2(n,x,y) result(res) + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: y(:,:) + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_), allocatable, dimension(:) :: res + real(psb_dpk_), external :: ddot + + ! local + integer(psb_ipk_) :: ncol + + ncol = size(y,2) + allocate(res(ncol)) + + call dgemv('T',n,ncol,done,y,n,x%v,1,dzero,res,1) + + end function d_base_dot_a2 + ! ! AXPBY is invoked via Y, hence the structure below. ! @@ -2312,7 +2341,8 @@ module psb_d_base_multivect_mod ! procedure, pass(x) :: dot_v => d_base_mlv_dot_v procedure, pass(x) :: dot_a => d_base_mlv_dot_a - generic, public :: dot => dot_v, dot_a + procedure, pass(x) :: dot_vect => d_base_mlv_dot_vect + generic, public :: dot => dot_v, dot_a, dot_vect procedure, pass(y) :: axpby_v => d_base_mlv_axpby_v procedure, pass(y) :: axpby_a => d_base_mlv_axpby_a generic, public :: axpby => axpby_v, axpby_a @@ -2883,7 +2913,6 @@ contains integer(psb_ipk_) :: j,nc if (x%is_dev()) call x%sync() - res = dzero ! ! Note: this is the base implementation. ! When we get here, we are sure that X is of @@ -2905,6 +2934,36 @@ contains end function d_base_mlv_dot_v + !> Function d_base_mlv_dot_vect + !! \memberof psb_d_base_multivect_type + !! \brief Dot product by another base_mlv_vector + !! \param n Number of entries to be considered + !! \param y The other (base_vect) to be multiplied by + !! + function d_base_mlv_dot_vect(n,x,y) result(res) + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + class(psb_d_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_), allocatable :: res(:) + real(psb_dpk_), external :: ddot + integer(psb_ipk_) :: j,nc + + if (x%is_dev()) call x%sync() + + select type(yy => y) + type is (psb_d_base_vect_type) + nc = psb_size(x%v,2_psb_ipk_) + allocate(res(nc)) + do j=1,nc + res(j) = ddot(n,x%v(:,j),1,y%v,1) + end do + class default + res = y%dot(n,x%v) + end select + + end function d_base_mlv_dot_vect + ! ! Base workhorse is good old BLAS1 ! diff --git a/base/modules/serial/psb_d_vect_mod.F90 b/base/modules/serial/psb_d_vect_mod.F90 index 2564c714..9406d236 100644 --- a/base/modules/serial/psb_d_vect_mod.F90 +++ b/base/modules/serial/psb_d_vect_mod.F90 @@ -1366,6 +1366,7 @@ end module psb_d_vect_mod module psb_d_multivect_mod use psb_d_base_multivect_mod + use psb_d_vect_mod use psb_const_mod use psb_i_vect_mod @@ -1388,11 +1389,18 @@ module psb_d_multivect_mod procedure, pass(x) :: get_dupl => d_mvect_get_dupl procedure, pass(x) :: set_dupl => d_mvect_set_dupl + procedure, pass(x) :: sync => d_mvect_sync + procedure, pass(x) :: is_host => d_mvect_is_host + procedure, pass(x) :: is_dev => d_mvect_is_dev + procedure, pass(x) :: is_sync => d_mvect_is_sync + procedure, pass(x) :: set_host => d_mvect_set_host + procedure, pass(x) :: set_dev => d_mvect_set_dev + procedure, pass(x) :: set_sync => d_mvect_set_sync + procedure, pass(x) :: all => d_mvect_all procedure, pass(x) :: reall => d_mvect_reall procedure, pass(x) :: zero => d_mvect_zero procedure, pass(x) :: asb => d_mvect_asb - procedure, pass(x) :: sync => d_mvect_sync procedure, pass(x) :: free => d_mvect_free procedure, pass(x) :: ins => d_mvect_ins procedure, pass(x) :: bld_x => d_mvect_bld_x @@ -1411,9 +1419,10 @@ module psb_d_multivect_mod procedure, pass(y) :: sctb => d_mvect_sctb procedure, pass(y) :: sctb_x => d_mvect_sctb_x generic, public :: sct => sctb, sctb_x -!!$ procedure, pass(x) :: dot_v => d_mvect_dot_v -!!$ procedure, pass(x) :: dot_a => d_mvect_dot_a -!!$ generic, public :: dot => dot_v, dot_a + procedure, pass(x) :: dot_v => d_mvect_dot_v + procedure, pass(x) :: dot_a => d_mvect_dot_a + procedure, pass(x) :: dot_vect => d_mvect_dot_vect + generic, public :: dot => dot_v, dot_a, dot_vect !!$ procedure, pass(y) :: axpby_v => d_mvect_axpby_v !!$ procedure, pass(y) :: axpby_a => d_mvect_axpby_a !!$ generic, public :: axpby => axpby_v, axpby_a @@ -1473,6 +1482,75 @@ contains x%dupl = psb_dupl_def_ end if end subroutine d_mvect_set_dupl + + subroutine d_mvect_sync(x) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%sync() + + end subroutine d_mvect_sync + + subroutine d_mvect_set_sync(x) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%set_sync() + + end subroutine d_mvect_set_sync + + subroutine d_mvect_set_host(x) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%set_host() + + end subroutine d_mvect_set_host + + subroutine d_mvect_set_dev(x) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%set_dev() + + end subroutine d_mvect_set_dev + + function d_mvect_is_sync(x) result(res) + implicit none + logical :: res + class(psb_d_multivect_type), intent(inout) :: x + + res = .true. + if (allocated(x%v)) & + & res = x%v%is_sync() + + end function d_mvect_is_sync + + function d_mvect_is_host(x) result(res) + implicit none + logical :: res + class(psb_d_multivect_type), intent(inout) :: x + + res = .true. + if (allocated(x%v)) & + & res = x%v%is_host() + + end function d_mvect_is_host + + function d_mvect_is_dev(x) result(res) + implicit none + logical :: res + class(psb_d_multivect_type), intent(inout) :: x + + res = .false. + if (allocated(x%v)) & + & res = x%v%is_dev() + + end function d_mvect_is_dev function d_mvect_is_remote_build(x) result(res) @@ -1717,15 +1795,6 @@ contains end subroutine d_mvect_asb - subroutine d_mvect_sync(x) - implicit none - class(psb_d_multivect_type), intent(inout) :: x - - if (allocated(x%v)) & - & call x%v%sync() - - end subroutine d_mvect_sync - subroutine d_mvect_gthab(n,idx,alpha,x,beta,y) use psi_serial_mod integer(psb_ipk_) :: n, idx(:) @@ -1840,30 +1909,50 @@ contains end subroutine d_mvect_cnv -!!$ function d_mvect_dot_v(n,x,y) result(res) -!!$ implicit none -!!$ class(psb_d_multivect_type), intent(inout) :: x, y -!!$ integer(psb_ipk_), intent(in) :: n -!!$ real(psb_dpk_) :: res -!!$ -!!$ res = dzero -!!$ if (allocated(x%v).and.allocated(y%v)) & -!!$ & res = x%v%dot(n,y%v) -!!$ -!!$ end function d_mvect_dot_v -!!$ -!!$ function d_mvect_dot_a(n,x,y) result(res) -!!$ implicit none -!!$ class(psb_d_multivect_type), intent(inout) :: x -!!$ real(psb_dpk_), intent(in) :: y(:) -!!$ integer(psb_ipk_), intent(in) :: n -!!$ real(psb_dpk_) :: res -!!$ -!!$ res = dzero -!!$ if (allocated(x%v)) & -!!$ & res = x%v%dot(n,y) -!!$ -!!$ end function d_mvect_dot_a + function d_mvect_dot_v(n,x,y) result(res) + implicit none + class(psb_d_multivect_type), intent(inout) :: x, y + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_), dimension(:), allocatable :: res + + if (allocated(x%v).and.allocated(y%v)) then + res = x%v%dot(n,y%v) + else + allocate(res(1)) + res(1) = psb_err_invalid_vect_state_ + end if + + end function d_mvect_dot_v + + function d_mvect_dot_vect(n,x,y) result(res) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + class(psb_d_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_), dimension(:), allocatable :: res + + if (allocated(x%v).and.allocated(y%v)) then + res = x%v%dot(n,y%v) + else + allocate(res(1)) + res(1) = psb_err_invalid_vect_state_ + end if + + end function d_mvect_dot_vect + + function d_mvect_dot_a(n,x,y) result(res) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: y(:,:) + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_), dimension(:), allocatable :: res + + res = dzero + if (allocated(x%v)) then + res = x%v%dot(n,y) + end if + + end function d_mvect_dot_a !!$ !!$ subroutine d_mvect_axpby_v(m,alpha, x, beta, y, info) !!$ use psi_serial_mod diff --git a/base/psblas/psb_ddot.f90 b/base/psblas/psb_ddot.f90 index 633c7549..ff4e9154 100644 --- a/base/psblas/psb_ddot.f90 +++ b/base/psblas/psb_ddot.f90 @@ -157,6 +157,347 @@ function psb_ddot_vect(x, y, desc_a,info,global) result(res) return end function psb_ddot_vect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! File: psb_ddot.f90 +! +! Function: psb_ddot_vect +! psb_ddot computes the dot product of two distributed vectors, +! +! dot := ( X )**C * ( Y ) +! +! +! Arguments: +! x - type(psb_d_vect_type) The input vector containing the entries of sub( X ). +! y - type(psb_d_vect_type) The input vector containing the entries of sub( Y ). +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global sum, default: .true. +! +! Note: from a functional point of view, X and Y are input, but here +! they are declared INOUT because of the sync() methods. +! +! +function psb_ddot_multivect(x, y, desc_a,info,global) result(res) + use psb_desc_mod + use psb_d_base_mat_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_d_multivect_mod + use psb_d_psblas_mod, psb_protect_name => psb_ddot_multivect + implicit none + real(psb_dpk_), dimension(:), allocatable :: res + type(psb_d_multivect_type), intent(inout) :: x, y + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: np, me, idx, ndm,& + & err_act, iix, jjx, iiy, jjy, i, nr + integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n + logical :: global_ + character(len=20) :: name, ch_err + + name='psb_ddot_multivect' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ctxt=desc_a%get_context() + call psb_info(ctxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + if (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = ione + ijx = ione + + iy = ione + ijy = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,x%get_ncols(),x%get_nrows(),ix,ijx,desc_a,info,iix,jjx) + if (info == psb_success_) & + & call psb_chkvect(m,y%get_ncols(),y%get_nrows(),iy,ijy,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + if (x%get_ncols() /= y%get_ncols()) then + info=psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + else + allocate(res(x%get_ncols()),stat=info) + if (info /= 0) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + end if + + nr = desc_a%get_local_rows() + if(nr > 0) then + res = x%dot(nr,y) + ! FIXME + ! adjust dot_local because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + if (y%is_dev()) call y%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + res(:) = res(:) - (real(ndm-1)/real(ndm))*(x%v%v(idx,:)*y%v%v(idx,:)) + end do + end if + else + res = dzero + end if + + ! compute global sum + if (global_) call psb_sum(ctxt, res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end function psb_ddot_multivect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the following conditions +! are met: +! 1. Redistributions of source code must retain the above copyright +! notice, this list of conditions and the following disclaimer. +! 2. Redistributions in binary form must reproduce the above copyright +! notice, this list of conditions, and the following disclaimer in the +! documentation and/or other materials provided with the distribution. +! 3. The name of the PSBLAS group or the names of its contributors may +! not be used to endorse or promote products derived from this +! software without specific written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +! POSSIBILITY OF SUCH DAMAGE. +! +! +! File: psb_ddot.f90 +! +! Function: psb_ddot_vect +! psb_ddot computes the dot product of two distributed vectors, +! +! dot := ( X )**C * ( Y ) +! +! +! Arguments: +! x - type(psb_d_vect_type) The input vector containing the entries of sub( X ). +! y - type(psb_d_vect_type) The input vector containing the entries of sub( Y ). +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global sum, default: .true. +! +! Note: from a functional point of view, X and Y are input, but here +! they are declared INOUT because of the sync() methods. +! +! +function psb_ddot_mvect_vect(x, y, desc_a,info,global) result(res) + use psb_desc_mod + use psb_d_base_mat_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_d_multivect_mod + use psb_d_vect_mod + use psb_d_psblas_mod, psb_protect_name => psb_ddot_mvect_vect + implicit none + real(psb_dpk_), dimension(:), allocatable :: res + type(psb_d_multivect_type), intent(inout) :: x + type(psb_d_vect_type), intent(inout) :: y + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: np, me, idx, ndm,& + & err_act, iix, jjx, iiy, jjy, i, nr + integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n + logical :: global_ + character(len=20) :: name, ch_err + + name='psb_ddot_mvect_vect' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ctxt=desc_a%get_context() + call psb_info(ctxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + if (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = ione + ijx = ione + + iy = ione + ijy = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,x%get_ncols(),x%get_nrows(),ix,ijx,desc_a,info,iix,jjx) + if (info == psb_success_) & + & call psb_chkvect(m,lone,y%get_nrows(),iy,ijy,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + allocate(res(x%get_ncols()),stat=info) + if (info /= 0) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + + nr = desc_a%get_local_rows() + if(nr > 0) then + res = x%dot(nr,y) + ! FIXME + ! adjust dot_local because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + if (y%is_dev()) call y%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + ! Remove the overlapped elements via dgemv calls + ! res = - (real(ndm-1)/real(ndm))* x(idx,:)^T y(idx) + 1.0 res + call dgemv('T',size(x%v%v,1),size(x%v%v,2),-(real(ndm-1)/real(ndm)), & + & size(x%v%v,1),y%v%v(idx),ione,done,res,ione) + end do + end if + else + res = dzero + end if + + ! compute global sum + if (global_) call psb_sum(ctxt, res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + +end function psb_ddot_mvect_vect ! ! Function: psb_ddot ! psb_ddot computes the dot product of two distributed vectors, diff --git a/test/kernel/Makefile b/test/kernel/Makefile index 7e413092..7dbcb53a 100644 --- a/test/kernel/Makefile +++ b/test/kernel/Makefile @@ -18,7 +18,7 @@ DPGOBJS=pdgenspmv.o DVECOBJS=vecoperation.o EXEDIR=./runs -all: runsd pdgenspmv +all: runsd pdgenspmv vecoperation #d_file_spmv s_file_spmv pdgenspmv vecoperation runsd: diff --git a/test/kernel/vecoperation.f90 b/test/kernel/vecoperation.f90 index 3860a3c3..e13fc37a 100644 --- a/test/kernel/vecoperation.f90 +++ b/test/kernel/vecoperation.f90 @@ -35,15 +35,20 @@ module unittestvector_mod use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& & psb_dspmat_type, psb_d_vect_type, dzero, psb_ctxt_type,& - & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type + & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type,& + & psb_d_base_multivect_type interface psb_gen_const - module procedure psb_d_gen_const + module procedure psb_d_gen_const, psb_d_gen_const_multi end interface psb_gen_const + interface psb_check_ans + module procedure psb_check_ans_v, psb_check_ans_mv + end interface psb_check_ans + contains - function psb_check_ans(v,val,ctxt) result(ans) + function psb_check_ans_v(v,val,ctxt) result(ans) use psb_base_mod implicit none @@ -73,7 +78,39 @@ contains ans = .false. end if - end function psb_check_ans + end function psb_check_ans_v + + function psb_check_ans_mv(v,val,ctxt) result(ans) + use psb_base_mod + + implicit none + + type(psb_d_multivect_type) :: v + real(psb_dpk_) :: val + type(psb_ctxt_type) :: ctxt + logical :: ans + + ! Local variables + integer(psb_ipk_) :: np, iam, info + real(psb_dpk_) :: check + real(psb_dpk_), allocatable :: va(:,:) + + call psb_info(ctxt,iam,np) + + va = v%get_vect() + va = va - val; + + check = maxval(va); + + call psb_sum(ctxt,check) + + if(check == 0.d0) then + ans = .true. + else + ans = .false. + end if + + end function psb_check_ans_mv ! ! subroutine to fill a vector with constant entries ! @@ -160,6 +197,93 @@ contains return end subroutine psb_d_gen_const + ! + ! subroutine to fill a multivectorvector with constant entries + ! + subroutine psb_d_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info) + use psb_base_mod + implicit none + + type(psb_d_multivect_type) :: v + type(psb_desc_type) :: desc_a + integer(psb_lpk_) :: idim + integer(psb_ipk_) :: jdim !! Number of columns of the multivector + type(psb_ctxt_type) :: ctxt + integer(psb_ipk_) :: info + real(psb_dpk_) :: val + + ! Local variables + integer(psb_ipk_), parameter :: nb=20 + real(psb_dpk_) :: zt(nb,jdim) ! Temporary array to fill the vector + character(len=40) :: name, ch_err + integer(psb_ipk_) :: np, iam, nr, nt + integer(psb_ipk_) :: n,nlr,ib,ii + integer(psb_ipk_) :: err_act + integer(psb_lpk_), allocatable :: myidx(:) + + + info = psb_success_ + name = 'create_constant_multivector' + call psb_erractionsave(err_act) + + call psb_info(ctxt, iam, np) + + n = idim*np ! The global dimension is the number of process times + ! the input size + + ! We use a simple minded block distribution + nt = (n+np-1)/np + nr = max(0,min(nt,n-(iam*nt))) + nt = nr + + call psb_sum(ctxt,nt) + if (nt /= n) then + write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n + info = -1 + call psb_barrier(ctxt) + call psb_abort(ctxt) + return + end if + ! Allocate the descriptor with simple minded data distribution + call psb_cdall(ctxt,desc_a,info,nl=nr) + ! Allocate the vector on the recently build descriptor + if (info == psb_success_) call psb_geall(v,desc_a,info,n=jdim) + ! Check that allocation has gone good + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + myidx = desc_a%get_global_indices() + nlr = size(myidx) + + do ii=1,nlr,nb + ib = min(nb,nlr-ii+1) + zt(:,:) = val + call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib,1:jdim),v,desc_a,info) + if(info /= psb_success_) exit + end do + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='insert rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! Assembly of communicator and vector + call psb_cdasb(desc_a,info) + if (info == psb_success_) call psb_geasb(v,desc_a,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ctxt,err_act) + + return + end subroutine psb_d_gen_const_multi end module unittestvector_mod @@ -171,7 +295,8 @@ program vecoperation implicit none ! input parameters - integer(psb_lpk_) :: idim = 100 + integer(psb_lpk_) :: idim = 100 ! Local vector size + integer(psb_ipk_) :: nmv = 10 ! Number of columns of the multivector ! miscellaneous real(psb_dpk_), parameter :: one = 1.d0 @@ -184,13 +309,17 @@ program vecoperation type(psb_desc_type) :: desc_a ! vector type(psb_d_vect_type) :: x,y,z + ! multivector + type(psb_d_multivect_type) :: mv1, mv2 ! blacs parameters type(psb_ctxt_type) :: ctxt integer(psb_ipk_) :: iam, np ! auxiliary parameters + integer(psb_ipk_) :: ii integer(psb_ipk_) :: info character(len=20) :: name,ch_err,readinput real(psb_dpk_) :: ans + real(psb_dpk_), allocatable :: ansmv(:) logical :: hasitnotfailed integer(psb_lpk_), allocatable :: myidx(:) integer(psb_ipk_) :: ib = 1 @@ -363,9 +492,54 @@ program vecoperation if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm")') end if + + ! + ! Test of multivector operation + ! + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'("Multivector Operations")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! X = 1 + call psb_d_gen_const_multi(mv1,one,idim,nmv,ctxt,desc_a,info) + hasitnotfailed = psb_check_ans(mv1,one,ctxt) + if (iam == psb_root_) then + if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector ")') + if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector ")') + end if + + ! + ! Multivector to field operation + ! + if (iam == psb_root_) write(psb_out_unit,'(" ")') + if (iam == psb_root_) write(psb_out_unit,'("Multivector to Field Operations")') + if (iam == psb_root_) write(psb_out_unit,'(" ")') + + ! Dot product: multivector vs multivector + call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info) + call psb_d_gen_const_multi(mv2,onehalf,idim,nmv,ctxt,desc_a,info) + ansmv = psb_gedot(mv1,mv2,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")') + end if + ! Dot product: multivector vs vector + call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info) + call psb_d_gen_const(x,onehalf,idim,ctxt,desc_a,info) + ansmv = psb_gedot(mv1,x,desc_a,info) + if (iam == psb_root_) then + ! write ansmv to check + if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")') + if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")') + end if + + + call psb_gefree(x,desc_a,info) call psb_gefree(y,desc_a,info) call psb_gefree(z,desc_a,info) + call psb_gefree(mv1,desc_a,info) + call psb_gefree(mv2,desc_a,info) call psb_cdfree(desc_a,info) if(info /= psb_success_) then info=psb_err_from_subroutine_