From c1f3b2d9d4a5ac3ba14c06678e2d7d2332163b79 Mon Sep 17 00:00:00 2001 From: cirdans-home Date: Thu, 21 Nov 2019 13:57:11 +0100 Subject: [PATCH] Added entrywise inversion and documentation for vec-vec operation --- base/modules/psblas/psb_c_psblas_mod.F90 | 20 +++ base/modules/psblas/psb_d_psblas_mod.F90 | 20 +++ base/modules/psblas/psb_s_psblas_mod.F90 | 20 +++ base/modules/psblas/psb_z_psblas_mod.F90 | 20 +++ base/modules/serial/psb_c_base_vect_mod.f90 | 111 ++++++++++++ base/modules/serial/psb_c_vect_mod.F90 | 63 +++++++ base/modules/serial/psb_d_base_vect_mod.f90 | 111 ++++++++++++ base/modules/serial/psb_d_vect_mod.F90 | 63 +++++++ base/modules/serial/psb_s_base_vect_mod.f90 | 111 ++++++++++++ base/modules/serial/psb_s_vect_mod.F90 | 63 +++++++ base/modules/serial/psb_z_base_vect_mod.f90 | 111 ++++++++++++ base/modules/serial/psb_z_vect_mod.F90 | 63 +++++++ base/psblas/Makefile | 3 +- base/psblas/psb_cinv_vect.f90 | 179 ++++++++++++++++++++ base/psblas/psb_dinv_vect.f90 | 179 ++++++++++++++++++++ base/psblas/psb_sinv_vect.f90 | 179 ++++++++++++++++++++ base/psblas/psb_zinv_vect.f90 | 179 ++++++++++++++++++++ test/kernel/vecoperation.f90 | 17 +- 18 files changed, 1508 insertions(+), 4 deletions(-) create mode 100644 base/psblas/psb_cinv_vect.f90 create mode 100644 base/psblas/psb_dinv_vect.f90 create mode 100644 base/psblas/psb_sinv_vect.f90 create mode 100644 base/psblas/psb_zinv_vect.f90 diff --git a/base/modules/psblas/psb_c_psblas_mod.F90 b/base/modules/psblas/psb_c_psblas_mod.F90 index 8e9a286a..a4b78538 100644 --- a/base/modules/psblas/psb_c_psblas_mod.F90 +++ b/base/modules/psblas/psb_c_psblas_mod.F90 @@ -459,4 +459,24 @@ module psb_c_psblas_mod end subroutine psb_cdiv_vect_check end interface + interface psb_geinv + subroutine psb_cinv_vect(x,y,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_c_vect_type + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_cinv_vect + subroutine psb_cinv_vect_check(x,y,desc_a,info,flag) + import :: psb_desc_type, psb_ipk_, & + & psb_c_vect_type + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + end subroutine psb_cinv_vect_check + 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 993ac542..c605e9a8 100644 --- a/base/modules/psblas/psb_d_psblas_mod.F90 +++ b/base/modules/psblas/psb_d_psblas_mod.F90 @@ -459,4 +459,24 @@ module psb_d_psblas_mod end subroutine psb_ddiv_vect_check end interface + interface psb_geinv + subroutine psb_dinv_vect(x,y,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_d_vect_type + type(psb_d_vect_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 + end subroutine psb_dinv_vect + subroutine psb_dinv_vect_check(x,y,desc_a,info,flag) + import :: psb_desc_type, psb_ipk_, & + & psb_d_vect_type + type(psb_d_vect_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) :: flag + end subroutine psb_dinv_vect_check + end interface + end module psb_d_psblas_mod diff --git a/base/modules/psblas/psb_s_psblas_mod.F90 b/base/modules/psblas/psb_s_psblas_mod.F90 index 1a2347d1..8abcbaf9 100644 --- a/base/modules/psblas/psb_s_psblas_mod.F90 +++ b/base/modules/psblas/psb_s_psblas_mod.F90 @@ -459,4 +459,24 @@ module psb_s_psblas_mod end subroutine psb_sdiv_vect_check end interface + interface psb_geinv + subroutine psb_sinv_vect(x,y,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_s_vect_type + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_sinv_vect + subroutine psb_sinv_vect_check(x,y,desc_a,info,flag) + import :: psb_desc_type, psb_ipk_, & + & psb_s_vect_type + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + end subroutine psb_sinv_vect_check + end interface + end module psb_s_psblas_mod diff --git a/base/modules/psblas/psb_z_psblas_mod.F90 b/base/modules/psblas/psb_z_psblas_mod.F90 index 7fd96c4c..d284a524 100644 --- a/base/modules/psblas/psb_z_psblas_mod.F90 +++ b/base/modules/psblas/psb_z_psblas_mod.F90 @@ -459,4 +459,24 @@ module psb_z_psblas_mod end subroutine psb_zdiv_vect_check end interface + interface psb_geinv + subroutine psb_zinv_vect(x,y,desc_a,info) + import :: psb_desc_type, psb_ipk_, & + & psb_z_vect_type + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_zinv_vect + subroutine psb_zinv_vect_check(x,y,desc_a,info,flag) + import :: psb_desc_type, psb_ipk_, & + & psb_z_vect_type + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + end subroutine psb_zinv_vect_check + 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 3927a60b..6bf54df0 100644 --- a/base/modules/serial/psb_c_base_vect_mod.f90 +++ b/base/modules/serial/psb_c_base_vect_mod.f90 @@ -171,6 +171,11 @@ module psb_c_base_vect_mod procedure, pass(z) :: div_a2 => c_base_div_a2 procedure, pass(z) :: div_a2_check => c_base_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => c_base_inv_v + procedure, pass(y) :: inv_v_check => c_base_inv_v_check + procedure, pass(y) :: inv_a2 => c_base_inv_a2 + procedure, pass(y) :: inv_a2_check => c_base_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check ! ! Scaling and norms ! @@ -1290,6 +1295,112 @@ contains end subroutine c_base_div_a2_check + ! + !> Function base_inv_v + !! \memberof psb_c_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + subroutine c_base_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + class(psb_c_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info) + + + end subroutine c_base_inv_v + ! + !> Function base_inv_v_check + !! \memberof psb_c_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + subroutine c_base_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + class(psb_c_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info,flag) + + + end subroutine c_base_inv_v_check + ! + !> Function base_inv_a2 + !! \memberof psb_c_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + ! + subroutine c_base_inv_a2(x, y, info) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: y + complex(psb_spk_), intent(in) :: x(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + y%v(i) = 1_psb_spk_/x(i) + end do + + end subroutine c_base_inv_a2 + ! + !> Function base_inv_a2_check + !! \memberof psb_c_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + ! + subroutine c_base_inv_a2_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: y + complex(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + integer(psb_ipk_) :: i, n + + if (flag .eqv. .false.) then + call c_base_inv_a2(x, y, info) + else + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + if (x(i) /= 0) then + y%v(i) = 1_psb_spk_/x(i) + else + info = 1 + exit + end if + end do + end if + + + end subroutine c_base_inv_a2_check ! diff --git a/base/modules/serial/psb_c_vect_mod.F90 b/base/modules/serial/psb_c_vect_mod.F90 index d4eb5a2a..f4fc0873 100644 --- a/base/modules/serial/psb_c_vect_mod.F90 +++ b/base/modules/serial/psb_c_vect_mod.F90 @@ -99,6 +99,11 @@ module psb_c_vect_mod procedure, pass(z) :: div_a2 => c_vect_div_a2 procedure, pass(z) :: div_a2_check => c_vect_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => c_vect_inv_v + procedure, pass(y) :: inv_v_check => c_vect_inv_v_check + procedure, pass(y) :: inv_a2 => c_vect_inv_a2 + procedure, pass(y) :: inv_a2_check => c_vect_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check procedure, pass(x) :: scal => c_vect_scal procedure, pass(x) :: absval1 => c_vect_absval1 procedure, pass(x) :: absval2 => c_vect_absval2 @@ -792,6 +797,64 @@ contains end subroutine c_vect_div_a2_check + subroutine c_vect_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_c_vect_type), intent(inout) :: x + class(psb_c_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info) + + end subroutine c_vect_inv_v + + subroutine c_vect_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_c_vect_type), intent(inout) :: x + class(psb_c_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info,flag) + + end subroutine c_vect_inv_v_check + + subroutine c_vect_inv_a2(x, y, info) + use psi_serial_mod + implicit none + complex(psb_spk_), intent(inout) :: x(:) + class(psb_c_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info) + + end subroutine c_vect_inv_a2 + + subroutine c_vect_inv_a2_check(x, y, info,flag) + use psi_serial_mod + implicit none + complex(psb_spk_), intent(inout) :: x(:) + class(psb_c_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info,flag) + + end subroutine c_vect_inv_a2_check + subroutine c_vect_scal(alpha, x) use psi_serial_mod implicit none diff --git a/base/modules/serial/psb_d_base_vect_mod.f90 b/base/modules/serial/psb_d_base_vect_mod.f90 index 3d00bbef..46ec3413 100644 --- a/base/modules/serial/psb_d_base_vect_mod.f90 +++ b/base/modules/serial/psb_d_base_vect_mod.f90 @@ -171,6 +171,11 @@ module psb_d_base_vect_mod procedure, pass(z) :: div_a2 => d_base_div_a2 procedure, pass(z) :: div_a2_check => d_base_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => d_base_inv_v + procedure, pass(y) :: inv_v_check => d_base_inv_v_check + procedure, pass(y) :: inv_a2 => d_base_inv_a2 + procedure, pass(y) :: inv_a2_check => d_base_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check ! ! Scaling and norms ! @@ -1290,6 +1295,112 @@ contains end subroutine d_base_div_a2_check + ! + !> Function base_inv_v + !! \memberof psb_d_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + subroutine d_base_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + class(psb_d_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info) + + + end subroutine d_base_inv_v + ! + !> Function base_inv_v_check + !! \memberof psb_d_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + subroutine d_base_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + class(psb_d_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info,flag) + + + end subroutine d_base_inv_v_check + ! + !> Function base_inv_a2 + !! \memberof psb_d_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + ! + subroutine d_base_inv_a2(x, y, info) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: y + real(psb_dpk_), intent(in) :: x(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + y%v(i) = 1_psb_dpk_/x(i) + end do + + end subroutine d_base_inv_a2 + ! + !> Function base_inv_a2_check + !! \memberof psb_d_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + ! + subroutine d_base_inv_a2_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: y + real(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + integer(psb_ipk_) :: i, n + + if (flag .eqv. .false.) then + call d_base_inv_a2(x, y, info) + else + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + if (x(i) /= 0) then + y%v(i) = 1_psb_dpk_/x(i) + else + info = 1 + exit + end if + end do + end if + + + end subroutine d_base_inv_a2_check ! diff --git a/base/modules/serial/psb_d_vect_mod.F90 b/base/modules/serial/psb_d_vect_mod.F90 index e7cf6757..4c74b842 100644 --- a/base/modules/serial/psb_d_vect_mod.F90 +++ b/base/modules/serial/psb_d_vect_mod.F90 @@ -99,6 +99,11 @@ module psb_d_vect_mod procedure, pass(z) :: div_a2 => d_vect_div_a2 procedure, pass(z) :: div_a2_check => d_vect_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => d_vect_inv_v + procedure, pass(y) :: inv_v_check => d_vect_inv_v_check + procedure, pass(y) :: inv_a2 => d_vect_inv_a2 + procedure, pass(y) :: inv_a2_check => d_vect_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check procedure, pass(x) :: scal => d_vect_scal procedure, pass(x) :: absval1 => d_vect_absval1 procedure, pass(x) :: absval2 => d_vect_absval2 @@ -792,6 +797,64 @@ contains end subroutine d_vect_div_a2_check + subroutine d_vect_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_d_vect_type), intent(inout) :: x + class(psb_d_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info) + + end subroutine d_vect_inv_v + + subroutine d_vect_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_d_vect_type), intent(inout) :: x + class(psb_d_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info,flag) + + end subroutine d_vect_inv_v_check + + subroutine d_vect_inv_a2(x, y, info) + use psi_serial_mod + implicit none + real(psb_dpk_), intent(inout) :: x(:) + class(psb_d_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info) + + end subroutine d_vect_inv_a2 + + subroutine d_vect_inv_a2_check(x, y, info,flag) + use psi_serial_mod + implicit none + real(psb_dpk_), intent(inout) :: x(:) + class(psb_d_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info,flag) + + end subroutine d_vect_inv_a2_check + subroutine d_vect_scal(alpha, x) use psi_serial_mod implicit none diff --git a/base/modules/serial/psb_s_base_vect_mod.f90 b/base/modules/serial/psb_s_base_vect_mod.f90 index bf229956..3c68d56f 100644 --- a/base/modules/serial/psb_s_base_vect_mod.f90 +++ b/base/modules/serial/psb_s_base_vect_mod.f90 @@ -171,6 +171,11 @@ module psb_s_base_vect_mod procedure, pass(z) :: div_a2 => s_base_div_a2 procedure, pass(z) :: div_a2_check => s_base_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => s_base_inv_v + procedure, pass(y) :: inv_v_check => s_base_inv_v_check + procedure, pass(y) :: inv_a2 => s_base_inv_a2 + procedure, pass(y) :: inv_a2_check => s_base_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check ! ! Scaling and norms ! @@ -1290,6 +1295,112 @@ contains end subroutine s_base_div_a2_check + ! + !> Function base_inv_v + !! \memberof psb_s_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + subroutine s_base_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + class(psb_s_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info) + + + end subroutine s_base_inv_v + ! + !> Function base_inv_v_check + !! \memberof psb_s_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + subroutine s_base_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + class(psb_s_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info,flag) + + + end subroutine s_base_inv_v_check + ! + !> Function base_inv_a2 + !! \memberof psb_s_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + ! + subroutine s_base_inv_a2(x, y, info) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: y + real(psb_spk_), intent(in) :: x(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + y%v(i) = 1_psb_spk_/x(i) + end do + + end subroutine s_base_inv_a2 + ! + !> Function base_inv_a2_check + !! \memberof psb_s_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + ! + subroutine s_base_inv_a2_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: y + real(psb_spk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + integer(psb_ipk_) :: i, n + + if (flag .eqv. .false.) then + call s_base_inv_a2(x, y, info) + else + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + if (x(i) /= 0) then + y%v(i) = 1_psb_spk_/x(i) + else + info = 1 + exit + end if + end do + end if + + + end subroutine s_base_inv_a2_check ! diff --git a/base/modules/serial/psb_s_vect_mod.F90 b/base/modules/serial/psb_s_vect_mod.F90 index 22519456..38478c43 100644 --- a/base/modules/serial/psb_s_vect_mod.F90 +++ b/base/modules/serial/psb_s_vect_mod.F90 @@ -99,6 +99,11 @@ module psb_s_vect_mod procedure, pass(z) :: div_a2 => s_vect_div_a2 procedure, pass(z) :: div_a2_check => s_vect_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => s_vect_inv_v + procedure, pass(y) :: inv_v_check => s_vect_inv_v_check + procedure, pass(y) :: inv_a2 => s_vect_inv_a2 + procedure, pass(y) :: inv_a2_check => s_vect_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check procedure, pass(x) :: scal => s_vect_scal procedure, pass(x) :: absval1 => s_vect_absval1 procedure, pass(x) :: absval2 => s_vect_absval2 @@ -792,6 +797,64 @@ contains end subroutine s_vect_div_a2_check + subroutine s_vect_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_s_vect_type), intent(inout) :: x + class(psb_s_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info) + + end subroutine s_vect_inv_v + + subroutine s_vect_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_s_vect_type), intent(inout) :: x + class(psb_s_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info,flag) + + end subroutine s_vect_inv_v_check + + subroutine s_vect_inv_a2(x, y, info) + use psi_serial_mod + implicit none + real(psb_spk_), intent(inout) :: x(:) + class(psb_s_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info) + + end subroutine s_vect_inv_a2 + + subroutine s_vect_inv_a2_check(x, y, info,flag) + use psi_serial_mod + implicit none + real(psb_spk_), intent(inout) :: x(:) + class(psb_s_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info,flag) + + end subroutine s_vect_inv_a2_check + subroutine s_vect_scal(alpha, x) use psi_serial_mod implicit none diff --git a/base/modules/serial/psb_z_base_vect_mod.f90 b/base/modules/serial/psb_z_base_vect_mod.f90 index 91663600..b87e16a6 100644 --- a/base/modules/serial/psb_z_base_vect_mod.f90 +++ b/base/modules/serial/psb_z_base_vect_mod.f90 @@ -171,6 +171,11 @@ module psb_z_base_vect_mod procedure, pass(z) :: div_a2 => z_base_div_a2 procedure, pass(z) :: div_a2_check => z_base_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => z_base_inv_v + procedure, pass(y) :: inv_v_check => z_base_inv_v_check + procedure, pass(y) :: inv_a2 => z_base_inv_a2 + procedure, pass(y) :: inv_a2_check => z_base_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check ! ! Scaling and norms ! @@ -1290,6 +1295,112 @@ contains end subroutine z_base_div_a2_check + ! + !> Function base_inv_v + !! \memberof psb_z_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + subroutine z_base_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + class(psb_z_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info) + + + end subroutine z_base_inv_v + ! + !> Function base_inv_v_check + !! \memberof psb_z_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x The vector to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + subroutine z_base_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + class(psb_z_base_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (y%is_dev()) call y%sync() + call y%inv(x%v,info,flag) + + + end subroutine z_base_inv_v_check + ! + !> Function base_inv_a2 + !! \memberof psb_z_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + ! + subroutine z_base_inv_a2(x, y, info) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: y + complex(psb_dpk_), intent(in) :: x(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + y%v(i) = 1_psb_dpk_/x(i) + end do + + end subroutine z_base_inv_a2 + ! + !> Function base_inv_a2_check + !! \memberof psb_z_base_vect_type + !! \brief Compute the entry-by-entry inverse of x and put it in y, with 0 check + !! \param x(:) The array to be inverted + !! \param y The vector containing the inverted vector + !! \param info return code + !! \param flag if true does the check, otherwise call base_inv_v + ! + subroutine z_base_inv_a2_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: y + complex(psb_dpk_), intent(inout) :: x(:) + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + integer(psb_ipk_) :: i, n + + if (flag .eqv. .false.) then + call z_base_inv_a2(x, y, info) + else + info = 0 + if (y%is_dev()) call y%sync() + + n = size(x) + do i=1, n + if (x(i) /= 0) then + y%v(i) = 1_psb_dpk_/x(i) + else + info = 1 + exit + end if + end do + end if + + + end subroutine z_base_inv_a2_check ! diff --git a/base/modules/serial/psb_z_vect_mod.F90 b/base/modules/serial/psb_z_vect_mod.F90 index f5e4f31a..be48d0db 100644 --- a/base/modules/serial/psb_z_vect_mod.F90 +++ b/base/modules/serial/psb_z_vect_mod.F90 @@ -99,6 +99,11 @@ module psb_z_vect_mod procedure, pass(z) :: div_a2 => z_vect_div_a2 procedure, pass(z) :: div_a2_check => z_vect_div_a2_check generic, public :: div => div_v, div_v_check, div_a2, div_a2_check + procedure, pass(y) :: inv_v => z_vect_inv_v + procedure, pass(y) :: inv_v_check => z_vect_inv_v_check + procedure, pass(y) :: inv_a2 => z_vect_inv_a2 + procedure, pass(y) :: inv_a2_check => z_vect_inv_a2_check + generic, public :: inv => inv_v, inv_v_check, inv_a2, inv_a2_check procedure, pass(x) :: scal => z_vect_scal procedure, pass(x) :: absval1 => z_vect_absval1 procedure, pass(x) :: absval2 => z_vect_absval2 @@ -792,6 +797,64 @@ contains end subroutine z_vect_div_a2_check + subroutine z_vect_inv_v(x, y, info) + use psi_serial_mod + implicit none + class(psb_z_vect_type), intent(inout) :: x + class(psb_z_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info) + + end subroutine z_vect_inv_v + + subroutine z_vect_inv_v_check(x, y, info, flag) + use psi_serial_mod + implicit none + class(psb_z_vect_type), intent(inout) :: x + class(psb_z_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(x%v).and.allocated(y%v)) & + & call y%v%inv(x%v,info,flag) + + end subroutine z_vect_inv_v_check + + subroutine z_vect_inv_a2(x, y, info) + use psi_serial_mod + implicit none + complex(psb_dpk_), intent(inout) :: x(:) + class(psb_z_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info) + + end subroutine z_vect_inv_a2 + + subroutine z_vect_inv_a2_check(x, y, info,flag) + use psi_serial_mod + implicit none + complex(psb_dpk_), intent(inout) :: x(:) + class(psb_z_vect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i, n + logical, intent(in) :: flag + + info = 0 + if (allocated(y%v)) & + & call y%v%inv(x,info,flag) + + end subroutine z_vect_inv_a2_check + subroutine z_vect_scal(alpha, x) use psi_serial_mod implicit none diff --git a/base/psblas/Makefile b/base/psblas/Makefile index 6e1b1a70..781f454a 100644 --- a/base/psblas/Makefile +++ b/base/psblas/Makefile @@ -11,7 +11,8 @@ OBJS= psb_ddot.o psb_damax.o psb_dasum.o psb_daxpby.o\ psb_camax.o psb_casum.o psb_caxpby.o psb_cdot.o \ psb_cnrm2.o psb_cnrmi.o psb_cspmm.o psb_cspsm.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_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 LIBDIR=.. diff --git a/base/psblas/psb_cinv_vect.f90 b/base/psblas/psb_cinv_vect.f90 new file mode 100644 index 00000000..1df47904 --- /dev/null +++ b/base/psblas/psb_cinv_vect.f90 @@ -0,0 +1,179 @@ +! +! 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_cinv_vect + +subroutine psb_cinv_vect(x,y,desc_a,info) + use psb_base_mod, psb_protect_name => psb_cinv_vect + implicit none + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + 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_inv_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(y%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,y%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 y%div(x,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_cinv_vect + +subroutine psb_cinv_vect_check(x,y,desc_a,info,flag) + use psb_base_mod, psb_protect_name => psb_cinv_vect_check + implicit none + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + + ! 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_inv_vect_check' + 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(y%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,y%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 y%div(x,info,flag) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_cinv_vect_check diff --git a/base/psblas/psb_dinv_vect.f90 b/base/psblas/psb_dinv_vect.f90 new file mode 100644 index 00000000..0ce370b4 --- /dev/null +++ b/base/psblas/psb_dinv_vect.f90 @@ -0,0 +1,179 @@ +! +! 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_dinv_vect + +subroutine psb_dinv_vect(x,y,desc_a,info) + use psb_base_mod, psb_protect_name => psb_dinv_vect + implicit none + type(psb_d_vect_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 + + ! 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_d_inv_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(y%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,y%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 y%div(x,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_dinv_vect + +subroutine psb_dinv_vect_check(x,y,desc_a,info,flag) + use psb_base_mod, psb_protect_name => psb_dinv_vect_check + implicit none + type(psb_d_vect_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) :: flag + + ! 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_d_inv_vect_check' + 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(y%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,y%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 y%div(x,info,flag) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_dinv_vect_check diff --git a/base/psblas/psb_sinv_vect.f90 b/base/psblas/psb_sinv_vect.f90 new file mode 100644 index 00000000..6925d066 --- /dev/null +++ b/base/psblas/psb_sinv_vect.f90 @@ -0,0 +1,179 @@ +! +! 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_sinv_vect + +subroutine psb_sinv_vect(x,y,desc_a,info) + use psb_base_mod, psb_protect_name => psb_sinv_vect + implicit none + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + 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_s_inv_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(y%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,y%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 y%div(x,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_sinv_vect + +subroutine psb_sinv_vect_check(x,y,desc_a,info,flag) + use psb_base_mod, psb_protect_name => psb_sinv_vect_check + implicit none + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + + ! 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_s_inv_vect_check' + 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(y%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,y%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 y%div(x,info,flag) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_sinv_vect_check diff --git a/base/psblas/psb_zinv_vect.f90 b/base/psblas/psb_zinv_vect.f90 new file mode 100644 index 00000000..7b36f6b0 --- /dev/null +++ b/base/psblas/psb_zinv_vect.f90 @@ -0,0 +1,179 @@ +! +! 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_zinv_vect + +subroutine psb_zinv_vect(x,y,desc_a,info) + use psb_base_mod, psb_protect_name => psb_zinv_vect + implicit none + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + 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_inv_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(y%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,y%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 y%div(x,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_zinv_vect + +subroutine psb_zinv_vect_check(x,y,desc_a,info,flag) + use psb_base_mod, psb_protect_name => psb_zinv_vect_check + implicit none + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: flag + + ! 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_inv_vect_check' + 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(y%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,y%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 y%div(x,info,flag) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_zinv_vect_check diff --git a/test/kernel/vecoperation.f90 b/test/kernel/vecoperation.f90 index e9cc88a8..a76049cf 100644 --- a/test/kernel/vecoperation.f90 +++ b/test/kernel/vecoperation.f90 @@ -46,7 +46,7 @@ program vecoperation ! descriptor type(psb_desc_type) :: desc_a ! vector - type(psb_d_vect_type) :: x,y + type(psb_d_vect_type) :: x,y,z ! blacs parameters integer(psb_ipk_) :: ictxt, iam, np ! other variables @@ -55,7 +55,7 @@ program vecoperation integer(psb_lpk_), allocatable :: myidx(:) real(psb_dpk_) :: zt(1), dotresult, norm2, norm1, norminf character(len=20) :: name,ch_err,readinput - real(psb_dpk_), allocatable :: vx(:), vy(:) + real(psb_dpk_), allocatable :: vx(:), vy(:), vz(:) info=psb_success_ @@ -113,6 +113,7 @@ program vecoperation ! Allocate memory call psb_geall(x,desc_a,info) call psb_geall(y,desc_a,info) + call psb_geall(z,desc_a,info) ! Put entries into the vectors do ii=1,nlr zt(1) = 1.0 @@ -159,6 +160,7 @@ program vecoperation call psb_geaxpby(1.0_psb_dpk_, x, 1.0_psb_dpk_, y, desc_a, info) ! \alpha x + \beta y if (iam == psb_root_) then + write(psb_out_unit,'("axpby : x + y")') vx = x%get_vect() write(psb_out_unit,'("x = ",es12.1)')vx(:) vy = y%get_vect() @@ -168,6 +170,7 @@ program vecoperation call psb_gemlt(x,y,desc_a,info) if (iam == psb_root_) then + write(psb_out_unit,'("mlt : y = x*y ")') vx = x%get_vect() write(psb_out_unit,'("x = ",es12.1)')vx(:) vy = y%get_vect() @@ -177,14 +180,22 @@ program vecoperation call psb_gediv(x,y,desc_a,info) if (iam == psb_root_) then + write(psb_out_unit,'("div : x = x/y")') vx = x%get_vect() write(psb_out_unit,'("x = ",es12.1)')vx(:) vy = y%get_vect() write(psb_out_unit,'("y = ",es12.1)')vy(:) end if + call psb_geinv(x,z,desc_a,info) - + if (iam == psb_root_) then + write(psb_out_unit,'("inv : z = 1/x")') + vx = x%get_vect() + write(psb_out_unit,'("x = ",es12.1)')vx(:) + vz = z%get_vect() + write(psb_out_unit,'("z = ",es12.1)')vy(:) + end if ! ! cleanup storage and exit