diff --git a/base/modules/psblas/psb_c_psblas_mod.F90 b/base/modules/psblas/psb_c_psblas_mod.F90 index 217632fb..a9c37d99 100644 --- a/base/modules/psblas/psb_c_psblas_mod.F90 +++ b/base/modules/psblas/psb_c_psblas_mod.F90 @@ -283,11 +283,22 @@ module psb_c_psblas_mod integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: global end function psb_cnrm2_weight_vect + function psb_cnrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_c_vect_type, psb_cspmat_type + real(psb_spk_) :: res + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: w + type(psb_c_vect_type), intent (inout) :: idv + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_cnrm2_weightmask_vect end interface #if ! defined(HAVE_BUGGY_GENERICS) interface psb_norm2 - procedure psb_cnrm2, psb_cnrm2v, psb_cnrm2_vect, psb_cnrm2_weight_vect + procedure psb_cnrm2, psb_cnrm2v, psb_cnrm2_vect, psb_cnrm2_weight_vect, psb_cnrm2_weightmask_vect end interface #endif @@ -500,5 +511,16 @@ module psb_c_psblas_mod end subroutine psb_cabs_vect end interface + interface psb_gecmp + subroutine psb_ccmp_vect(x,c,z,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_c_vect_type, psb_spk_ + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: z + real(psb_spk_), intent(in) :: c + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_ccmp_vect + end interface end module psb_c_psblas_mod diff --git a/base/modules/psblas/psb_d_psblas_mod.F90 b/base/modules/psblas/psb_d_psblas_mod.F90 index 011ce683..4303f8bc 100644 --- a/base/modules/psblas/psb_d_psblas_mod.F90 +++ b/base/modules/psblas/psb_d_psblas_mod.F90 @@ -283,11 +283,22 @@ module psb_d_psblas_mod integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: global end function psb_dnrm2_weight_vect + function psb_dnrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_vect_type, psb_dspmat_type + real(psb_dpk_) :: res + type(psb_d_vect_type), intent (inout) :: x + type(psb_d_vect_type), intent (inout) :: w + type(psb_d_vect_type), intent (inout) :: idv + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_dnrm2_weightmask_vect end interface #if ! defined(HAVE_BUGGY_GENERICS) interface psb_norm2 - procedure psb_dnrm2, psb_dnrm2v, psb_dnrm2_vect, psb_dnrm2_weight_vect + procedure psb_dnrm2, psb_dnrm2v, psb_dnrm2_vect, psb_dnrm2_weight_vect, psb_dnrm2_weightmask_vect end interface #endif diff --git a/base/modules/psblas/psb_s_psblas_mod.F90 b/base/modules/psblas/psb_s_psblas_mod.F90 index 32296412..04969899 100644 --- a/base/modules/psblas/psb_s_psblas_mod.F90 +++ b/base/modules/psblas/psb_s_psblas_mod.F90 @@ -283,11 +283,22 @@ module psb_s_psblas_mod integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: global end function psb_snrm2_weight_vect + function psb_snrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_s_vect_type, psb_sspmat_type + real(psb_spk_) :: res + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: w + type(psb_s_vect_type), intent (inout) :: idv + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_snrm2_weightmask_vect end interface #if ! defined(HAVE_BUGGY_GENERICS) interface psb_norm2 - procedure psb_snrm2, psb_snrm2v, psb_snrm2_vect, psb_snrm2_weight_vect + procedure psb_snrm2, psb_snrm2v, psb_snrm2_vect, psb_snrm2_weight_vect, psb_snrm2_weightmask_vect end interface #endif diff --git a/base/modules/psblas/psb_z_psblas_mod.F90 b/base/modules/psblas/psb_z_psblas_mod.F90 index 71394e68..dc4d0262 100644 --- a/base/modules/psblas/psb_z_psblas_mod.F90 +++ b/base/modules/psblas/psb_z_psblas_mod.F90 @@ -283,11 +283,22 @@ module psb_z_psblas_mod integer(psb_ipk_), intent(out) :: info logical, intent(in), optional :: global end function psb_znrm2_weight_vect + function psb_znrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_z_vect_type, psb_zspmat_type + real(psb_dpk_) :: res + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: w + type(psb_z_vect_type), intent (inout) :: idv + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + end function psb_znrm2_weightmask_vect end interface #if ! defined(HAVE_BUGGY_GENERICS) interface psb_norm2 - procedure psb_znrm2, psb_znrm2v, psb_znrm2_vect, psb_znrm2_weight_vect + procedure psb_znrm2, psb_znrm2v, psb_znrm2_vect, psb_znrm2_weight_vect, psb_znrm2_weightmask_vect end interface #endif @@ -500,5 +511,16 @@ module psb_z_psblas_mod end subroutine psb_zabs_vect end interface + interface psb_gecmp + subroutine psb_zcmp_vect(x,c,z,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_z_vect_type, psb_dpk_ + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: z + real(psb_dpk_), intent(in) :: c + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_zcmp_vect + end interface end module psb_z_psblas_mod diff --git a/base/modules/serial/psb_c_base_vect_mod.f90 b/base/modules/serial/psb_c_base_vect_mod.f90 index 0eca0832..d629256a 100644 --- a/base/modules/serial/psb_c_base_vect_mod.f90 +++ b/base/modules/serial/psb_c_base_vect_mod.f90 @@ -187,6 +187,12 @@ module psb_c_base_vect_mod procedure, pass(x) :: amax => c_base_amax procedure, pass(x) :: asum => c_base_asum + ! + ! Comparison and mask operation + ! + procedure, pass(z) :: cmp_a2 => c_base_cmp_a2 + procedure, pass(z) :: cmp_v2 => c_base_cmp_v2 + generic, public :: cmp => cmp_a2,cmp_v2 end type psb_c_base_vect_type @@ -1402,6 +1408,60 @@ contains end subroutine c_base_inv_a2_check + + ! + !> Function base_inv_a2_check + !! \memberof psb_c_base_vect_type + !! \brief Compare entry-by-entry the vector x with the scalar c + !! \param x The array to be compared + !! \param z The vector containing in position i 1 if |x(i)| > c, 0 otherwise + !! \param c The comparison term + !! \param info return code + ! + subroutine c_base_cmp_a2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_spk_), intent(in) :: c + complex(psb_spk_), intent(inout) :: x(:) + class(psb_c_base_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + if (z%is_dev()) call z%sync() + + n = size(x) + do i = 1, n, 1 + if ( abs(x(i)).ge.c ) then + z%v(i) = 1_psb_spk_ + else + z%v(i) = 0_psb_spk_ + end if + end do + info = 0 + + end subroutine c_base_cmp_a2 + ! + !> Function base_cmp_v2 + !! \memberof psb_c_base_vect_type + !! \brief Compare entry-by-entry the vector x with the scalar c + !! \param x The vector to be compared + !! \param z The vector containing in position i 1 if |x(i)| > c, 0 otherwise + !! \param c The comparison term + !! \param info return code + ! + subroutine c_base_cmp_v2(x,c,z,info) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + real(psb_spk_), intent(in) :: c + class(psb_c_base_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (x%is_dev()) call x%sync() + call z%cmp(x%v,c,info) + end subroutine c_base_cmp_v2 + ! ! Simple scaling ! diff --git a/base/modules/serial/psb_c_vect_mod.F90 b/base/modules/serial/psb_c_vect_mod.F90 index 8fc78b2c..393cbaa8 100644 --- a/base/modules/serial/psb_c_vect_mod.F90 +++ b/base/modules/serial/psb_c_vect_mod.F90 @@ -110,9 +110,14 @@ module psb_c_vect_mod generic, public :: absval => absval1, absval2 procedure, pass(x) :: nrm2std => c_vect_nrm2 procedure, pass(x) :: nrm2weight => c_vect_nrm2_weight - generic, public :: nrm2 => nrm2std, nrm2weight + procedure, pass(x) :: nrm2weightmask => c_vect_nrm2_weight_mask + generic, public :: nrm2 => nrm2std, nrm2weight, nrm2weightmask procedure, pass(x) :: amax => c_vect_amax procedure, pass(x) :: asum => c_vect_asum + procedure, pass(z) :: cmp_a2 => c_vect_cmp_a2 + procedure, pass(z) :: cmp_v2 => c_vect_cmp_v2 + generic, public :: cmp => cmp_a2, cmp_v2 + end type psb_c_vect_type public :: psb_c_vect @@ -862,6 +867,33 @@ contains end subroutine c_vect_inv_a2_check + subroutine c_vect_cmp_a2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_spk_), intent(in) :: c + complex(psb_spk_), intent(inout) :: x(:) + class(psb_c_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(z%v)) & + & call z%cmp(x,c,info) + + end subroutine c_vect_cmp_a2 + + subroutine c_vect_cmp_v2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_spk_), intent(in) :: c + class(psb_c_vect_type), intent(inout) :: x + class(psb_c_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v).and.allocated(z%v)) & + & call z%v%cmp(x%v,c,info) + + end subroutine c_vect_cmp_v2 subroutine c_vect_scal(alpha, x) use psi_serial_mod @@ -922,6 +954,25 @@ contains end function c_vect_nrm2_weight + function c_vect_nrm2_weight_mask(n,x,w,id) result(res) + implicit none + class(psb_c_vect_type), intent(inout) :: x + class(psb_c_vect_type), intent(inout) :: w + class(psb_c_vect_type), intent(inout) :: id + integer(psb_ipk_), intent(in) :: n + real(psb_spk_) :: res + integer(psb_ipk_) :: info + + if (allocated(x%v).and.allocated(w%v).and.allocated(id%v)) then + call w%v%cmp(id%v,szero,info) + call w%v%mlt(x%v,info) + res = w%v%nrm2(n) + else + res = szero + end if + + end function c_vect_nrm2_weight_mask + function c_vect_amax(n,x) result(res) implicit none class(psb_c_vect_type), intent(inout) :: x diff --git a/base/modules/serial/psb_d_base_vect_mod.f90 b/base/modules/serial/psb_d_base_vect_mod.f90 index dc0c56aa..69241893 100644 --- a/base/modules/serial/psb_d_base_vect_mod.f90 +++ b/base/modules/serial/psb_d_base_vect_mod.f90 @@ -1408,6 +1408,7 @@ contains end subroutine d_base_inv_a2_check + ! !> Function base_inv_a2_check !! \memberof psb_d_base_vect_type @@ -1420,7 +1421,7 @@ contains subroutine d_base_cmp_a2(x,c,z,info) use psi_serial_mod implicit none - real(psb_dpk_), intent(in) :: c + real(psb_dpk_), intent(in) :: c real(psb_dpk_), intent(inout) :: x(:) class(psb_d_base_vect_type), intent(inout) :: z integer(psb_ipk_), intent(out) :: info @@ -1452,7 +1453,7 @@ contains use psi_serial_mod implicit none class(psb_d_base_vect_type), intent(inout) :: x - real(psb_dpk_), intent(in) :: c + real(psb_dpk_), intent(in) :: c class(psb_d_base_vect_type), intent(inout) :: z integer(psb_ipk_), intent(out) :: info @@ -1460,6 +1461,7 @@ contains if (x%is_dev()) call x%sync() call z%cmp(x%v,c,info) end subroutine d_base_cmp_v2 + ! ! Simple scaling ! diff --git a/base/modules/serial/psb_d_vect_mod.F90 b/base/modules/serial/psb_d_vect_mod.F90 index 40c9f9a2..d897e784 100644 --- a/base/modules/serial/psb_d_vect_mod.F90 +++ b/base/modules/serial/psb_d_vect_mod.F90 @@ -110,12 +110,14 @@ module psb_d_vect_mod generic, public :: absval => absval1, absval2 procedure, pass(x) :: nrm2std => d_vect_nrm2 procedure, pass(x) :: nrm2weight => d_vect_nrm2_weight - generic, public :: nrm2 => nrm2std, nrm2weight + procedure, pass(x) :: nrm2weightmask => d_vect_nrm2_weight_mask + generic, public :: nrm2 => nrm2std, nrm2weight, nrm2weightmask procedure, pass(x) :: amax => d_vect_amax procedure, pass(x) :: asum => d_vect_asum procedure, pass(z) :: cmp_a2 => d_vect_cmp_a2 procedure, pass(z) :: cmp_v2 => d_vect_cmp_v2 generic, public :: cmp => cmp_a2, cmp_v2 + end type psb_d_vect_type public :: psb_d_vect @@ -952,6 +954,25 @@ contains end function d_vect_nrm2_weight + function d_vect_nrm2_weight_mask(n,x,w,id) result(res) + implicit none + class(psb_d_vect_type), intent(inout) :: x + class(psb_d_vect_type), intent(inout) :: w + class(psb_d_vect_type), intent(inout) :: id + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_) :: res + integer(psb_ipk_) :: info + + if (allocated(x%v).and.allocated(w%v).and.allocated(id%v)) then + call w%v%cmp(id%v,dzero,info) + call w%v%mlt(x%v,info) + res = w%v%nrm2(n) + else + res = dzero + end if + + end function d_vect_nrm2_weight_mask + function d_vect_amax(n,x) result(res) implicit none class(psb_d_vect_type), intent(inout) :: x diff --git a/base/modules/serial/psb_i_vect_mod.F90 b/base/modules/serial/psb_i_vect_mod.F90 index 7fb233f8..ca857109 100644 --- a/base/modules/serial/psb_i_vect_mod.F90 +++ b/base/modules/serial/psb_i_vect_mod.F90 @@ -79,6 +79,7 @@ module psb_i_vect_mod procedure, pass(x) :: set_dev => i_vect_set_dev procedure, pass(x) :: set_sync => i_vect_set_sync + end type psb_i_vect_type public :: psb_i_vect diff --git a/base/modules/serial/psb_l_vect_mod.F90 b/base/modules/serial/psb_l_vect_mod.F90 index 3a955752..3472e1d8 100644 --- a/base/modules/serial/psb_l_vect_mod.F90 +++ b/base/modules/serial/psb_l_vect_mod.F90 @@ -80,6 +80,7 @@ module psb_l_vect_mod procedure, pass(x) :: set_dev => l_vect_set_dev procedure, pass(x) :: set_sync => l_vect_set_sync + end type psb_l_vect_type public :: psb_l_vect diff --git a/base/modules/serial/psb_s_base_vect_mod.f90 b/base/modules/serial/psb_s_base_vect_mod.f90 index 450fdd29..1c7e7162 100644 --- a/base/modules/serial/psb_s_base_vect_mod.f90 +++ b/base/modules/serial/psb_s_base_vect_mod.f90 @@ -1408,6 +1408,7 @@ contains end subroutine s_base_inv_a2_check + ! !> Function base_inv_a2_check !! \memberof psb_s_base_vect_type @@ -1420,7 +1421,7 @@ contains subroutine s_base_cmp_a2(x,c,z,info) use psi_serial_mod implicit none - real(psb_spk_), intent(in) :: c + real(psb_spk_), intent(in) :: c real(psb_spk_), intent(inout) :: x(:) class(psb_s_base_vect_type), intent(inout) :: z integer(psb_ipk_), intent(out) :: info @@ -1452,7 +1453,7 @@ contains use psi_serial_mod implicit none class(psb_s_base_vect_type), intent(inout) :: x - real(psb_spk_), intent(in) :: c + real(psb_spk_), intent(in) :: c class(psb_s_base_vect_type), intent(inout) :: z integer(psb_ipk_), intent(out) :: info @@ -1460,6 +1461,7 @@ contains if (x%is_dev()) call x%sync() call z%cmp(x%v,c,info) end subroutine s_base_cmp_v2 + ! ! Simple scaling ! diff --git a/base/modules/serial/psb_s_vect_mod.F90 b/base/modules/serial/psb_s_vect_mod.F90 index 26f8a3b1..8ae50c7d 100644 --- a/base/modules/serial/psb_s_vect_mod.F90 +++ b/base/modules/serial/psb_s_vect_mod.F90 @@ -110,12 +110,14 @@ module psb_s_vect_mod generic, public :: absval => absval1, absval2 procedure, pass(x) :: nrm2std => s_vect_nrm2 procedure, pass(x) :: nrm2weight => s_vect_nrm2_weight - generic, public :: nrm2 => nrm2std, nrm2weight + procedure, pass(x) :: nrm2weightmask => s_vect_nrm2_weight_mask + generic, public :: nrm2 => nrm2std, nrm2weight, nrm2weightmask procedure, pass(x) :: amax => s_vect_amax procedure, pass(x) :: asum => s_vect_asum procedure, pass(z) :: cmp_a2 => s_vect_cmp_a2 procedure, pass(z) :: cmp_v2 => s_vect_cmp_v2 generic, public :: cmp => cmp_a2, cmp_v2 + end type psb_s_vect_type public :: psb_s_vect @@ -952,6 +954,25 @@ contains end function s_vect_nrm2_weight + function s_vect_nrm2_weight_mask(n,x,w,id) result(res) + implicit none + class(psb_s_vect_type), intent(inout) :: x + class(psb_s_vect_type), intent(inout) :: w + class(psb_s_vect_type), intent(inout) :: id + integer(psb_ipk_), intent(in) :: n + real(psb_spk_) :: res + integer(psb_ipk_) :: info + + if (allocated(x%v).and.allocated(w%v).and.allocated(id%v)) then + call w%v%cmp(id%v,szero,info) + call w%v%mlt(x%v,info) + res = w%v%nrm2(n) + else + res = szero + end if + + end function s_vect_nrm2_weight_mask + function s_vect_amax(n,x) result(res) implicit none class(psb_s_vect_type), intent(inout) :: x diff --git a/base/modules/serial/psb_z_base_vect_mod.f90 b/base/modules/serial/psb_z_base_vect_mod.f90 index 808d9d4f..f885c149 100644 --- a/base/modules/serial/psb_z_base_vect_mod.f90 +++ b/base/modules/serial/psb_z_base_vect_mod.f90 @@ -187,6 +187,12 @@ module psb_z_base_vect_mod procedure, pass(x) :: amax => z_base_amax procedure, pass(x) :: asum => z_base_asum + ! + ! Comparison and mask operation + ! + procedure, pass(z) :: cmp_a2 => z_base_cmp_a2 + procedure, pass(z) :: cmp_v2 => z_base_cmp_v2 + generic, public :: cmp => cmp_a2,cmp_v2 end type psb_z_base_vect_type @@ -1402,6 +1408,60 @@ contains end subroutine z_base_inv_a2_check + + ! + !> Function base_inv_a2_check + !! \memberof psb_z_base_vect_type + !! \brief Compare entry-by-entry the vector x with the scalar c + !! \param x The array to be compared + !! \param z The vector containing in position i 1 if |x(i)| > c, 0 otherwise + !! \param c The comparison term + !! \param info return code + ! + subroutine z_base_cmp_a2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_dpk_), intent(in) :: c + complex(psb_dpk_), intent(inout) :: x(:) + class(psb_z_base_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + if (z%is_dev()) call z%sync() + + n = size(x) + do i = 1, n, 1 + if ( abs(x(i)).ge.c ) then + z%v(i) = 1_psb_dpk_ + else + z%v(i) = 0_psb_dpk_ + end if + end do + info = 0 + + end subroutine z_base_cmp_a2 + ! + !> Function base_cmp_v2 + !! \memberof psb_z_base_vect_type + !! \brief Compare entry-by-entry the vector x with the scalar c + !! \param x The vector to be compared + !! \param z The vector containing in position i 1 if |x(i)| > c, 0 otherwise + !! \param c The comparison term + !! \param info return code + ! + subroutine z_base_cmp_v2(x,c,z,info) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: c + class(psb_z_base_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (x%is_dev()) call x%sync() + call z%cmp(x%v,c,info) + end subroutine z_base_cmp_v2 + ! ! Simple scaling ! diff --git a/base/modules/serial/psb_z_vect_mod.F90 b/base/modules/serial/psb_z_vect_mod.F90 index 2caf2950..7358fe44 100644 --- a/base/modules/serial/psb_z_vect_mod.F90 +++ b/base/modules/serial/psb_z_vect_mod.F90 @@ -110,9 +110,14 @@ module psb_z_vect_mod generic, public :: absval => absval1, absval2 procedure, pass(x) :: nrm2std => z_vect_nrm2 procedure, pass(x) :: nrm2weight => z_vect_nrm2_weight - generic, public :: nrm2 => nrm2std, nrm2weight + procedure, pass(x) :: nrm2weightmask => z_vect_nrm2_weight_mask + generic, public :: nrm2 => nrm2std, nrm2weight, nrm2weightmask procedure, pass(x) :: amax => z_vect_amax procedure, pass(x) :: asum => z_vect_asum + procedure, pass(z) :: cmp_a2 => z_vect_cmp_a2 + procedure, pass(z) :: cmp_v2 => z_vect_cmp_v2 + generic, public :: cmp => cmp_a2, cmp_v2 + end type psb_z_vect_type public :: psb_z_vect @@ -862,6 +867,33 @@ contains end subroutine z_vect_inv_a2_check + subroutine z_vect_cmp_a2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_dpk_), intent(in) :: c + complex(psb_dpk_), intent(inout) :: x(:) + class(psb_z_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(z%v)) & + & call z%cmp(x,c,info) + + end subroutine z_vect_cmp_a2 + + subroutine z_vect_cmp_v2(x,c,z,info) + use psi_serial_mod + implicit none + real(psb_dpk_), intent(in) :: c + class(psb_z_vect_type), intent(inout) :: x + class(psb_z_vect_type), intent(inout) :: z + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v).and.allocated(z%v)) & + & call z%v%cmp(x%v,c,info) + + end subroutine z_vect_cmp_v2 subroutine z_vect_scal(alpha, x) use psi_serial_mod @@ -922,6 +954,25 @@ contains end function z_vect_nrm2_weight + function z_vect_nrm2_weight_mask(n,x,w,id) result(res) + implicit none + class(psb_z_vect_type), intent(inout) :: x + class(psb_z_vect_type), intent(inout) :: w + class(psb_z_vect_type), intent(inout) :: id + integer(psb_ipk_), intent(in) :: n + real(psb_dpk_) :: res + integer(psb_ipk_) :: info + + if (allocated(x%v).and.allocated(w%v).and.allocated(id%v)) then + call w%v%cmp(id%v,dzero,info) + call w%v%mlt(x%v,info) + res = w%v%nrm2(n) + else + res = dzero + end if + + end function z_vect_nrm2_weight_mask + function z_vect_amax(n,x) result(res) implicit none class(psb_z_vect_type), intent(inout) :: x diff --git a/base/psblas/Makefile b/base/psblas/Makefile index 2e7f89ec..1c54aaa5 100644 --- a/base/psblas/Makefile +++ b/base/psblas/Makefile @@ -13,7 +13,7 @@ OBJS= psb_ddot.o psb_damax.o psb_dasum.o psb_daxpby.o\ psb_cmlt_vect.o psb_dmlt_vect.o psb_zmlt_vect.o psb_smlt_vect.o\ psb_cdiv_vect.o psb_ddiv_vect.o psb_zdiv_vect.o psb_sdiv_vect.o\ psb_cinv_vect.o psb_dinv_vect.o psb_zinv_vect.o psb_sinv_vect.o\ - psb_dcmp_vect.o psb_scmp_vect.o \ + psb_dcmp_vect.o psb_scmp_vect.o psb_ccmp_vect.o psb_zcmp_vect.o\ psb_cabs_vect.o psb_dabs_vect.o psb_sabs_vect.o \ psb_zabs_vect.o diff --git a/base/psblas/psb_ccmp_vect.f90 b/base/psblas/psb_ccmp_vect.f90 new file mode 100644 index 00000000..27d4f989 --- /dev/null +++ b/base/psblas/psb_ccmp_vect.f90 @@ -0,0 +1,105 @@ +! +! 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_ccmp_vect + +subroutine psb_ccmp_vect(x,c,z,desc_a,info) + use psb_base_mod, psb_protect_name => psb_ccmp_vect + implicit none + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: z + real(psb_spk_), intent(in) :: c + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy + integer(psb_lpk_) :: ix, ijx, iy, ijy, m + character(len=20) :: name, ch_err + + name='psb_c_cmp_vect' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, 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(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if(desc_a%get_local_rows() > 0) then + call z%cmp(x,c,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_ccmp_vect diff --git a/base/psblas/psb_cnrm2.f90 b/base/psblas/psb_cnrm2.f90 index c3ad75a4..d357c951 100644 --- a/base/psblas/psb_cnrm2.f90 +++ b/base/psblas/psb_cnrm2.f90 @@ -477,6 +477,115 @@ function psb_cnrm2_weight_vect(x,w, desc_a, info,global) result(res) return end function psb_cnrm2_weight_vect +! Function: psb_cnrm2_weight_vect +! Computes the weighted norm2 of a distributed vector with respect to a mask +! contained in the vector id. +! +! norm2 := sqrt ( (w(id > 0).*X(id > 0))**C * (w(id > 0).*X(id > 0))) +! +! Arguments: +! x - type(psb_c_vect_type) The input vector containing the entries of X. +! w - type(psb_c_vect_type) The input vector containing the entries of W. +! id - type(psb_c_vect_type) The inpute vector containing the mask +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global reduction, default: .true. +! +function psb_cnrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + use psb_desc_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_c_vect_mod + implicit none + + real(psb_spk_) :: res + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: w + type(psb_c_vect_type), intent (inout) :: idv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ndim, i, id, idx, ndm, ldx + integer(psb_lpk_) :: ix, jx, iy, ijy, m + logical :: global_ + real(psb_spk_) :: snrm2, dd + character(len=20) :: name, ch_err + + name='psb_cnrm2v_weightmask' + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -1) 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 (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = 1 + jx = 1 + m = desc_a%get_global_rows() + ldx = x%get_nrows() + call psb_chkvect(m,lone,ldx,ix,jx,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + end if + + if (iix /= 1) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + if (desc_a%get_local_rows() > 0) then + ndim = desc_a%get_local_rows() + res = x%nrm2(ndim,w,idv) + ! adjust because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + dd = dble(ndm-1)/dble(ndm) + res = res - sqrt(cone - dd*(abs(x%v%v(idx))/res)**2) + end do + end if + else + res = szero + end if + + if (global_) call psb_nrm2(ictxt,res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end function psb_cnrm2_weightmask_vect !!$ !!$ Parallel Sparse BLAS version 3.5 diff --git a/base/psblas/psb_dnrm2.f90 b/base/psblas/psb_dnrm2.f90 index 96eded9a..16b18d91 100644 --- a/base/psblas/psb_dnrm2.f90 +++ b/base/psblas/psb_dnrm2.f90 @@ -477,6 +477,115 @@ function psb_dnrm2_weight_vect(x,w, desc_a, info,global) result(res) return end function psb_dnrm2_weight_vect +! Function: psb_dnrm2_weight_vect +! Computes the weighted norm2 of a distributed vector with respect to a mask +! contained in the vector id. +! +! norm2 := sqrt ( (w(id > 0).*X(id > 0))**C * (w(id > 0).*X(id > 0))) +! +! Arguments: +! x - type(psb_d_vect_type) The input vector containing the entries of X. +! w - type(psb_d_vect_type) The input vector containing the entries of W. +! id - type(psb_d_vect_type) The inpute vector containing the mask +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global reduction, default: .true. +! +function psb_dnrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + use psb_desc_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_d_vect_mod + implicit none + + real(psb_dpk_) :: res + type(psb_d_vect_type), intent (inout) :: x + type(psb_d_vect_type), intent (inout) :: w + type(psb_d_vect_type), intent (inout) :: idv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ndim, i, id, idx, ndm, ldx + integer(psb_lpk_) :: ix, jx, iy, ijy, m + logical :: global_ + real(psb_dpk_) :: snrm2, dd + character(len=20) :: name, ch_err + + name='psb_dnrm2v_weightmask' + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -1) 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 (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = 1 + jx = 1 + m = desc_a%get_global_rows() + ldx = x%get_nrows() + call psb_chkvect(m,lone,ldx,ix,jx,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + end if + + if (iix /= 1) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + if (desc_a%get_local_rows() > 0) then + ndim = desc_a%get_local_rows() + res = x%nrm2(ndim,w,idv) + ! adjust because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + dd = dble(ndm-1)/dble(ndm) + res = res - sqrt(done - dd*(abs(x%v%v(idx))/res)**2) + end do + end if + else + res = dzero + end if + + if (global_) call psb_nrm2(ictxt,res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end function psb_dnrm2_weightmask_vect !!$ !!$ Parallel Sparse BLAS version 3.5 diff --git a/base/psblas/psb_snrm2.f90 b/base/psblas/psb_snrm2.f90 index cb2dd507..ab9c56ca 100644 --- a/base/psblas/psb_snrm2.f90 +++ b/base/psblas/psb_snrm2.f90 @@ -477,6 +477,115 @@ function psb_snrm2_weight_vect(x,w, desc_a, info,global) result(res) return end function psb_snrm2_weight_vect +! Function: psb_snrm2_weight_vect +! Computes the weighted norm2 of a distributed vector with respect to a mask +! contained in the vector id. +! +! norm2 := sqrt ( (w(id > 0).*X(id > 0))**C * (w(id > 0).*X(id > 0))) +! +! Arguments: +! x - type(psb_s_vect_type) The input vector containing the entries of X. +! w - type(psb_s_vect_type) The input vector containing the entries of W. +! id - type(psb_s_vect_type) The inpute vector containing the mask +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global reduction, default: .true. +! +function psb_snrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + use psb_desc_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_s_vect_mod + implicit none + + real(psb_spk_) :: res + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: w + type(psb_s_vect_type), intent (inout) :: idv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ndim, i, id, idx, ndm, ldx + integer(psb_lpk_) :: ix, jx, iy, ijy, m + logical :: global_ + real(psb_spk_) :: snrm2, dd + character(len=20) :: name, ch_err + + name='psb_snrm2v_weightmask' + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -1) 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 (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = 1 + jx = 1 + m = desc_a%get_global_rows() + ldx = x%get_nrows() + call psb_chkvect(m,lone,ldx,ix,jx,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + end if + + if (iix /= 1) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + if (desc_a%get_local_rows() > 0) then + ndim = desc_a%get_local_rows() + res = x%nrm2(ndim,w,idv) + ! adjust because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + dd = dble(ndm-1)/dble(ndm) + res = res - sqrt(sone - dd*(abs(x%v%v(idx))/res)**2) + end do + end if + else + res = szero + end if + + if (global_) call psb_nrm2(ictxt,res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end function psb_snrm2_weightmask_vect !!$ !!$ Parallel Sparse BLAS version 3.5 diff --git a/base/psblas/psb_zcmp_vect.f90 b/base/psblas/psb_zcmp_vect.f90 new file mode 100644 index 00000000..fc62d528 --- /dev/null +++ b/base/psblas/psb_zcmp_vect.f90 @@ -0,0 +1,105 @@ +! +! 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_zcmp_vect + +subroutine psb_zcmp_vect(x,c,z,desc_a,info) + use psb_base_mod, psb_protect_name => psb_zcmp_vect + implicit none + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: z + real(psb_dpk_), intent(in) :: c + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy + integer(psb_lpk_) :: ix, ijx, iy, ijy, m + character(len=20) :: name, ch_err + + name='psb_z_cmp_vect' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, 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(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if(desc_a%get_local_rows() > 0) then + call z%cmp(x,c,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_zcmp_vect diff --git a/base/psblas/psb_znrm2.f90 b/base/psblas/psb_znrm2.f90 index 0cde7e21..bfa29d18 100644 --- a/base/psblas/psb_znrm2.f90 +++ b/base/psblas/psb_znrm2.f90 @@ -477,6 +477,115 @@ function psb_znrm2_weight_vect(x,w, desc_a, info,global) result(res) return end function psb_znrm2_weight_vect +! Function: psb_znrm2_weight_vect +! Computes the weighted norm2 of a distributed vector with respect to a mask +! contained in the vector id. +! +! norm2 := sqrt ( (w(id > 0).*X(id > 0))**C * (w(id > 0).*X(id > 0))) +! +! Arguments: +! x - type(psb_z_vect_type) The input vector containing the entries of X. +! w - type(psb_z_vect_type) The input vector containing the entries of W. +! id - type(psb_z_vect_type) The inpute vector containing the mask +! desc_a - type(psb_desc_type). The communication descriptor. +! info - integer. Return code +! global - logical(optional) Whether to perform the global reduction, default: .true. +! +function psb_znrm2_weightmask_vect(x,w,idv, desc_a, info,global) result(res) + use psb_desc_mod + use psb_check_mod + use psb_error_mod + use psb_penv_mod + use psb_z_vect_mod + implicit none + + real(psb_dpk_) :: res + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: w + type(psb_z_vect_type), intent (inout) :: idv + type(psb_desc_type), intent(in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: global + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, ndim, i, id, idx, ndm, ldx + integer(psb_lpk_) :: ix, jx, iy, ijy, m + logical :: global_ + real(psb_dpk_) :: snrm2, dd + character(len=20) :: name, ch_err + + name='psb_znrm2v_weightmask' + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -1) 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 (present(global)) then + global_ = global + else + global_ = .true. + end if + + ix = 1 + jx = 1 + m = desc_a%get_global_rows() + ldx = x%get_nrows() + call psb_chkvect(m,lone,ldx,ix,jx,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect' + call psb_errpush(info,name,a_err=ch_err) + end if + + if (iix /= 1) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + goto 9999 + end if + + if (desc_a%get_local_rows() > 0) then + ndim = desc_a%get_local_rows() + res = x%nrm2(ndim,w,idv) + ! adjust because overlapped elements are computed more than once + if (size(desc_a%ovrlap_elem,1)>0) then + if (x%is_dev()) call x%sync() + do i=1,size(desc_a%ovrlap_elem,1) + idx = desc_a%ovrlap_elem(i,1) + ndm = desc_a%ovrlap_elem(i,2) + dd = dble(ndm-1)/dble(ndm) + res = res - sqrt(zone - dd*(abs(x%v%v(idx))/res)**2) + end do + end if + else + res = dzero + end if + + if (global_) call psb_nrm2(ictxt,res) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end function psb_znrm2_weightmask_vect !!$ !!$ Parallel Sparse BLAS version 3.5 diff --git a/test/kernel/vecoperation.f90 b/test/kernel/vecoperation.f90 index 506a185c..22dd56ae 100644 --- a/test/kernel/vecoperation.f90 +++ b/test/kernel/vecoperation.f90 @@ -53,7 +53,7 @@ program vecoperation integer(psb_ipk_) :: nr, nlr, info, i, ii, ib=1 integer(psb_lpk_) :: nt integer(psb_lpk_), allocatable :: myidx(:) - real(psb_dpk_) :: zt(1), dotresult, norm2, norm1, norminf, norm2w + real(psb_dpk_) :: zt(1), dotresult, norm2, norm1, norminf, norm2w, norm2wc character(len=20) :: name,ch_err,readinput real(psb_dpk_), allocatable :: vx(:), vy(:), vz(:) real(psb_dpk_) :: c @@ -229,11 +229,13 @@ program vecoperation norm1 = psb_norm1(x,desc_a,info) norm2 = psb_norm2(x,desc_a,info) norm2w = psb_norm2(x,absz,desc_a,info) + norm2wc = psb_norm2(x,absz,z,desc_a,info) norminf = psb_normi(x,desc_a,info) - write(psb_out_unit,'("\|x\|_inf : ",es12.5," \|x\|_1 :",es12.5)')norminf,norm1 - write(psb_out_unit,'(" \|x\|_2 : ",es12.5," \|x\|_2,w : ",es12.5)')norm2,norm2w - - + if (iam == psb_root_) then + write(psb_out_unit,'("\|x\|_inf : ",es12.5," \|x\|_1 :",es12.5)')norminf,norm1 + write(psb_out_unit,'(" \|x\|_2 : ",es12.5," \|x\|_2,w : ",es12.5)')norm2,norm2w + write(psb_out_unit,'(" masked \|x\|_2,w : ",es12.5)')norm2wc + end if ! ! cleanup storage and exit