psblas3/base/modules/psb_s_base_vect_mod.f90

!!$ 
!!$              Parallel Sparse BLAS  version 3.0
!!$    (C) Copyright 2006, 2007, 2008, 2009, 2010
!!$                       Salvatore Filippone    University of Rome Tor Vergata
!!$                       Alfredo Buttari        CNRS-IRIT, Toulouse
!!$ 
!!$  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.
!!$ 
!!$  
!
! package: psb_s_base_vect_mod
!
! This module contains the definition of the psb_s_base_vect type which
! is a container for dense vectors.
!  This is encapsulated instead of being just a simple array to allow for 
!  more complicated situations, such as GPU programming, where the memory
!  area we are interested in is not easily accessible from the host/Fortran
!  side. It is also meant to be encapsulated in an outer type, to allow
!  runtime switching as per the STATE design pattern, similar to the
!  sparse matrix types.
!
!
module psb_s_base_vect_mod
  
  use psb_const_mod
  use psb_error_mod


  !> \namespace  psb_base_mod  \class psb_s_base_vect_type
  !! The psb_s_base_vect_type 
  !! defines a middle level  real(psb_spk_) encapsulated dense vector.
  !! The encapsulation is needed, in place of a simple array, to allow  
  !! for complicated situations, such as GPU programming, where the memory
  !!  area we are interested in is not easily accessible from the host/Fortran
  !!  side. It is also meant to be encapsulated in an outer type, to allow
  !!  runtime switching as per the STATE design pattern, similar to the
  !!  sparse matrix types.
  !!
  type psb_s_base_vect_type
    real(psb_spk_), allocatable :: v(:)
  contains
    !
    !  Constructors/allocators
    !
    procedure, pass(x) :: bld_x    => s_base_bld_x
    procedure, pass(x) :: bld_n    => s_base_bld_n
    generic, public    :: bld      => bld_x, bld_n
    procedure, pass(x) :: all      => s_base_all
    !
    ! Insert/set. Assembly and free.
    ! Assembly does almost nothing here, but is important
    ! in derived classes. 
    !
    procedure, pass(x) :: ins      => s_base_ins
    procedure, pass(x) :: zero     => s_base_zero
    procedure, pass(x) :: asb      => s_base_asb
    procedure, pass(x) :: free     => s_base_free
    !
    ! Sync: centerpiece of handling of external storage.
    ! Any derived class having extra storage upon sync
    ! will guarantee that both fortran/host side and
    ! external side contain the same data. The base
    ! version is only a placeholder. 
    !
    procedure, pass(x) :: sync     => s_base_sync
    procedure, pass(x) :: is_host  => s_base_is_host
    procedure, pass(x) :: is_dev   => s_base_is_dev
    procedure, pass(x) :: is_sync  => s_base_is_sync
    procedure, pass(x) :: set_host => s_base_set_host
    procedure, pass(x) :: set_dev  => s_base_set_dev
    procedure, pass(x) :: set_sync => s_base_set_sync

    !
    ! Basic info
    procedure, pass(x) :: get_nrows => s_base_get_nrows
    procedure, pass(x) :: sizeof   => s_base_sizeof
    !
    ! Set/get data from/to an external array; also
    ! overload assignment.
    !
    procedure, pass(x) :: get_vect => s_base_get_vect
    procedure, pass(x) :: set_scal => s_base_set_scal
    procedure, pass(x) :: set_vect => s_base_set_vect
    generic, public    :: set      => set_vect, set_scal

    !
    ! Dot product and AXPBY
    !
    procedure, pass(x) :: dot_v    => s_base_dot_v
    procedure, pass(x) :: dot_a    => s_base_dot_a
    generic, public    :: dot      => dot_v, dot_a
    procedure, pass(y) :: axpby_v  => s_base_axpby_v
    procedure, pass(y) :: axpby_a  => s_base_axpby_a
    generic, public    :: axpby    => axpby_v, axpby_a
    !
    ! Vector by vector multiplication. Need all variants
    ! to handle multiple requirements from preconditioners
    !
    procedure, pass(y) :: mlt_v    => s_base_mlt_v
    procedure, pass(y) :: mlt_a    => s_base_mlt_a
    procedure, pass(z) :: mlt_a_2  => s_base_mlt_a_2
    procedure, pass(z) :: mlt_v_2  => s_base_mlt_v_2
    procedure, pass(z) :: mlt_va   => s_base_mlt_va
    procedure, pass(z) :: mlt_av   => s_base_mlt_av
    generic, public    :: mlt      => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
    !
    ! Scaling and norms
    !
    procedure, pass(x) :: scal     => s_base_scal
    procedure, pass(x) :: nrm2     => s_base_nrm2
    procedure, pass(x) :: amax     => s_base_amax
    procedure, pass(x) :: asum     => s_base_asum
    !
    ! Gather/scatter. These are needed for MPI interfacing.
    ! May have to be reworked. 
    !
    procedure, pass(x) :: gthab    => s_base_gthab
    procedure, pass(x) :: gthzv    => s_base_gthzv
    generic, public    :: gth      => gthab, gthzv
    procedure, pass(y) :: sctb     => s_base_sctb
    generic, public    :: sct      => sctb
  end type psb_s_base_vect_type

  public  :: psb_s_base_vect
  private :: constructor, size_const
  interface psb_s_base_vect
    module procedure constructor, size_const
  end interface psb_s_base_vect

contains
  
  !
  ! Constructors. 
  !
  
  function constructor(x) result(this)
    real(psb_spk_)   :: x(:)
    type(psb_s_base_vect_type) :: this
    integer(psb_ipk_) :: info

    this%v = x
    call this%asb(size(x,kind=psb_ipk_),info)
  end function constructor
    
  
  function size_const(n) result(this)
    integer(psb_ipk_), intent(in) :: n
    type(psb_s_base_vect_type) :: this
    integer(psb_ipk_) :: info

    call this%asb(n,info)

  end function size_const
  
  !
  ! Build from a sample
  !

  subroutine s_base_bld_x(x,this)
    use psb_realloc_mod
    real(psb_spk_), intent(in) :: this(:)
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(size(this),x%v,info)
    if (info /= 0) then 
      call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
      return
    end if
    x%v(:)  = this(:)

  end subroutine s_base_bld_x
    
  !
  ! Create with size, but no initialization
  !
  subroutine s_base_bld_n(x,n)
    use psb_realloc_mod
    integer(psb_ipk_), intent(in) :: n
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_) :: info

    call psb_realloc(n,x%v,info)
    call x%asb(n,info)

  end subroutine s_base_bld_n
  
  subroutine s_base_all(n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none 
    integer(psb_ipk_), intent(in)               :: n
    class(psb_s_base_vect_type), intent(out)    :: x
    integer(psb_ipk_), intent(out)              :: info
    
    call psb_realloc(n,x%v,info)
    
  end subroutine s_base_all

  !
  ! Insert a bunch of values at specified positions.
  !
  subroutine s_base_ins(n,irl,val,dupl,x,info)
    use psi_serial_mod
    implicit none 
    class(psb_s_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(in)               :: n, dupl
    integer(psb_ipk_), intent(in)               :: irl(:)
    real(psb_spk_), intent(in)        :: val(:)
    integer(psb_ipk_), intent(out)              :: info

    integer(psb_ipk_) :: i

    info = 0
    if (psb_errstatus_fatal()) return 

    if (.not.allocated(x%v)) then 
      info = psb_err_invalid_vect_state_
    else if (n > min(size(irl),size(val))) then 
      info = psb_err_invalid_input_

    else 
      select case(dupl) 
      case(psb_dupl_ovwrt_) 
        do i = 1, n
          !loop over all val's rows

          ! row actual block row 
          if (irl(i) > 0) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i)) = val(i)
          end if
        enddo

      case(psb_dupl_add_) 

        do i = 1, n
          !loop over all val's rows

          if (irl(i) > 0) then
            ! this row belongs to me
            ! copy i-th row of block val in x
            x%v(irl(i)) = x%v(irl(i)) +  val(i)
          end if
        enddo

      case default
        info = 321
!!$      call psb_errpush(info,name)
!!$      goto 9999
      end select
    end if
    if (info /= 0) then 
      call psb_errpush(info,'base_vect_ins')
      return
    end if

  end subroutine s_base_ins

  !
  subroutine s_base_zero(x)
    use psi_serial_mod
    implicit none 
    class(psb_s_base_vect_type), intent(inout)    :: x
    
    if (allocated(x%v)) x%v=szero

  end subroutine s_base_zero

  
  !
  ! Assembly.
  ! For derived classes: after this the vector
  ! storage is supposed to be in sync.
  !
  
  subroutine s_base_asb(n, x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none 
    integer(psb_ipk_), intent(in)              :: n
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(out)             :: info
    
    if (x%get_nrows() < n) &
         & call psb_realloc(n,x%v,info)
    if (info /= 0) &
         & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')

  end subroutine s_base_asb


  subroutine s_base_free(x, info)
    use psi_serial_mod
    use psb_realloc_mod
    implicit none 
    class(psb_s_base_vect_type), intent(inout)  :: x
    integer(psb_ipk_), intent(out)              :: info
    
    info = 0
    if (allocated(x%v)) deallocate(x%v, stat=info)
    if (info /= 0) call & 
         & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
        
  end subroutine s_base_free

  

  !
  ! The base version of SYNC & friends does nothing, it's just
  ! a placeholder.
  ! 
  subroutine s_base_sync(x)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    
  end subroutine s_base_sync


  subroutine s_base_set_host(x)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    
  end subroutine s_base_set_host

  subroutine s_base_set_dev(x)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    
  end subroutine s_base_set_dev

  subroutine s_base_set_sync(x)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    
  end subroutine s_base_set_sync

  function s_base_is_dev(x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(in) :: x
    logical  :: res
  
    res = .false.
  end function s_base_is_dev
  
  function s_base_is_host(x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function s_base_is_host

  function s_base_is_sync(x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(in) :: x
    logical  :: res

    res = .true.
  end function s_base_is_sync


  !
  ! Size info. 
  !
    
  function s_base_get_nrows(x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(in) :: x
    integer(psb_ipk_) :: res

    res = 0
    if (allocated(x%v)) res = size(x%v)

  end function s_base_get_nrows

  function s_base_sizeof(x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(in) :: x
    integer(psb_long_int_k_) :: res
    
    ! Force 8-byte integers.
    res = (1_psb_long_int_k_ * psb_sizeof_sp) * x%get_nrows()

  end function s_base_sizeof

  
  !
  ! Two versions of extracting an array: one of them
  ! overload the assignment. 
  !
    
  function  s_base_get_vect(x) result(res)
    class(psb_s_base_vect_type), intent(inout) :: x
    real(psb_spk_), allocatable                 :: res(:)
    integer(psb_ipk_) :: info
    
    if (.not.allocated(x%v)) return 
    call x%sync()
    allocate(res(x%get_nrows()),stat=info) 
    if (info /= 0) then 
      call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
      return
    end if
    res(:) = x%v(:)
  end function s_base_get_vect
    
  !
  ! Reset all values 
  !
  subroutine s_base_set_scal(x,val)
    class(psb_s_base_vect_type), intent(inout)  :: x
    real(psb_spk_), intent(in) :: val
        
    integer(psb_ipk_) :: info
    x%v = val
    
  end subroutine s_base_set_scal

  subroutine s_base_set_vect(x,val)
    class(psb_s_base_vect_type), intent(inout)  :: x
    real(psb_spk_), intent(in) :: val(:)
    integer(psb_ipk_) :: nr
    integer(psb_ipk_) :: info

    if (allocated(x%v)) then 
      nr = min(size(x%v),size(val))
      x%v(1:nr) = val(1:nr)
    else
      x%v = val
    end if

  end subroutine s_base_set_vect

  !
  ! Dot products 
  ! 
  function s_base_dot_v(n,x,y) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x, y
    integer(psb_ipk_), intent(in)           :: n
    real(psb_spk_)                :: res
    real(psb_spk_), external      :: sdot
    
    res = szero
    !
    ! Note: this is the base implementation.
    !  When we get here, we are sure that X is of
    !  TYPE psb_s_base_vect.
    !  If Y is not, throw the burden on it, implicitly
    !  calling dot_a
    !
    select type(yy => y)
    type is (psb_s_base_vect_type)
      res = sdot(n,x%v,1,y%v,1)
    class default
      res = y%dot(n,x%v)
    end select

  end function s_base_dot_v

  !
  ! Base workhorse is good old BLAS1
  !
  function s_base_dot_a(n,x,y) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    real(psb_spk_), intent(in)    :: y(:)
    integer(psb_ipk_), intent(in)           :: n
    real(psb_spk_)                :: res
    real(psb_spk_), external      :: sdot
    
    res = sdot(n,y,1,x%v,1)

  end function s_base_dot_a
    
  !
  ! AXPBY is invoked via Y, hence the structure below. 
  !
  subroutine s_base_axpby_v(m,alpha, x, beta, y, info)
    use psi_serial_mod
    implicit none 
    integer(psb_ipk_), intent(in)               :: m
    class(psb_s_base_vect_type), intent(inout)  :: x
    class(psb_s_base_vect_type), intent(inout)  :: y
    real(psb_spk_), intent (in)       :: alpha, beta
    integer(psb_ipk_), intent(out)              :: info
    
    select type(xx => x)
    type is (psb_s_base_vect_type)
      call psb_geaxpby(m,alpha,x%v,beta,y%v,info)
    class default
      call y%axpby(m,alpha,x%v,beta,info)
    end select

  end subroutine s_base_axpby_v

  subroutine s_base_axpby_a(m,alpha, x, beta, y, info)
    use psi_serial_mod
    implicit none 
    integer(psb_ipk_), intent(in)               :: m
    real(psb_spk_), intent(in)        :: x(:)
    class(psb_s_base_vect_type), intent(inout)  :: y
    real(psb_spk_), intent (in)       :: alpha, beta
    integer(psb_ipk_), intent(out)              :: info
    
    call psb_geaxpby(m,alpha,x,beta,y%v,info)
    
  end subroutine s_base_axpby_a

  
  !
  !  Multiple variants of two operations:
  !  Simple multiplication  Y(:) = X(:)*Y(:)
  !  blas-like:   Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
  !
  !  Variants expanded according to the dynamic type
  !  of the involved entities
  !

  subroutine s_base_mlt_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
    select type(xx => x)
    type is (psb_s_base_vect_type)
      n = min(size(y%v), size(xx%v))
      do i=1, n 
        y%v(i) = y%v(i)*xx%v(i)
      end do
    class default
      call y%mlt(x%v,info)
    end select

  end subroutine s_base_mlt_v

  subroutine s_base_mlt_a(x, y, info)
    use psi_serial_mod
    implicit none 
    real(psb_spk_), intent(in)        :: x(:)
    class(psb_s_base_vect_type), intent(inout)  :: y
    integer(psb_ipk_), intent(out)              :: info
    integer(psb_ipk_) :: i, n

    info = 0
    n = min(size(y%v), size(x))
    do i=1, n 
      y%v(i) = y%v(i)*x(i)
    end do
    
  end subroutine s_base_mlt_a


  subroutine s_base_mlt_a_2(alpha,x,y,beta,z,info)
    use psi_serial_mod
    implicit none 
    real(psb_spk_), intent(in)        :: alpha,beta
    real(psb_spk_), intent(in)        :: y(:)
    real(psb_spk_), intent(in)        :: x(:)
    class(psb_s_base_vect_type), intent(inout)  :: z
    integer(psb_ipk_), intent(out)              :: info
    integer(psb_ipk_) :: i, n

    info = 0    
    n = min(size(z%v), size(x), size(y))
!!$    write(0,*) 'Mlt_a_2: ',n
    if (alpha == szero) then 
      if (beta == sone) then 
        return 
      else
        do i=1, n
          z%v(i) = beta*z%v(i)
        end do
      end if
    else
      if (alpha == sone) then 
        if (beta == szero) then 
          do i=1, n 
            z%v(i) = y(i)*x(i)
          end do
        else if (beta == sone) then 
          do i=1, n 
            z%v(i) = z%v(i) + y(i)*x(i)
          end do
        else 
          do i=1, n 
            z%v(i) = beta*z%v(i) + y(i)*x(i)
          end do
        end if
      else if (alpha == -sone) then 
        if (beta == szero) then 
          do i=1, n 
            z%v(i) = -y(i)*x(i)
          end do
        else if (beta == sone) then 
          do i=1, n 
            z%v(i) = z%v(i) - y(i)*x(i)
          end do
        else 
          do i=1, n 
            z%v(i) = beta*z%v(i) - y(i)*x(i)
          end do
        end if
      else
        if (beta == szero) then 
          do i=1, n 
            z%v(i) = alpha*y(i)*x(i)
          end do
        else if (beta == sone) then 
          do i=1, n 
            z%v(i) = z%v(i) + alpha*y(i)*x(i)
          end do
        else 
          do i=1, n 
            z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
          end do
        end if
      end if
    end if
  end subroutine s_base_mlt_a_2

  subroutine s_base_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
    use psi_serial_mod
    use psb_string_mod
    implicit none 
    real(psb_spk_), intent(in)        :: alpha,beta
    class(psb_s_base_vect_type), intent(inout)  :: x
    class(psb_s_base_vect_type), intent(inout)  :: y
    class(psb_s_base_vect_type), intent(inout)  :: z
    integer(psb_ipk_), intent(out)              :: info    
    character(len=1), intent(in), optional     :: conjgx, conjgy
    integer(psb_ipk_) :: i, n
    logical :: conjgx_, conjgy_

    info = 0
    if (.not.psb_s_is_complex_) then
      call z%mlt(alpha,x%v,y%v,beta,info)
    else 
      conjgx_=.false.
      if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
      conjgy_=.false.
      if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
      if (conjgx_) x%v=(x%v)
      if (conjgy_) y%v=(y%v)
      call z%mlt(alpha,x%v,y%v,beta,info)
      if (conjgx_) x%v=(x%v)
      if (conjgy_) y%v=(y%v)
    end if
  end subroutine s_base_mlt_v_2

  subroutine s_base_mlt_av(alpha,x,y,beta,z,info)
    use psi_serial_mod
    implicit none 
    real(psb_spk_), intent(in)        :: alpha,beta
    real(psb_spk_), intent(in)        :: x(:)
    class(psb_s_base_vect_type), intent(inout)  :: y
    class(psb_s_base_vect_type), intent(inout)  :: z
    integer(psb_ipk_), intent(out)              :: info    
    integer(psb_ipk_) :: i, n

    info = 0
    
    call z%mlt(alpha,x,y%v,beta,info)

  end subroutine s_base_mlt_av

  subroutine s_base_mlt_va(alpha,x,y,beta,z,info)
    use psi_serial_mod
    implicit none 
    real(psb_spk_), intent(in)        :: alpha,beta
    real(psb_spk_), intent(in)        :: y(:)
    class(psb_s_base_vect_type), intent(inout)  :: x
    class(psb_s_base_vect_type), intent(inout)  :: z
    integer(psb_ipk_), intent(out)              :: info    
    integer(psb_ipk_) :: i, n

    info = 0
    
    call z%mlt(alpha,y,x,beta,info)

  end subroutine s_base_mlt_va


  !
  ! Simple scaling 
  !
  
  subroutine s_base_scal(alpha, x)
    use psi_serial_mod
    implicit none 
    class(psb_s_base_vect_type), intent(inout)  :: x
    real(psb_spk_), intent (in)       :: alpha
    
    if (allocated(x%v)) x%v = alpha*x%v

  end subroutine s_base_scal
  
  !
  ! Norms 1, 2 and infinity
  !

  function s_base_nrm2(n,x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)           :: n
    real(psb_spk_)                :: res
    real(psb_spk_), external      :: snrm2
    
    res =  snrm2(n,x%v,1)

  end function s_base_nrm2
  
  function s_base_amax(n,x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)           :: n
    real(psb_spk_)                :: res
    
    res =  maxval(abs(x%v(1:n)))

  end function s_base_amax

  function s_base_asum(n,x) result(res)
    implicit none 
    class(psb_s_base_vect_type), intent(inout) :: x
    integer(psb_ipk_), intent(in)           :: n
    real(psb_spk_)                :: res
    
    res =  sum(abs(x%v(1:n)))

  end function s_base_asum
  
  
  !
  ! Gather: Y = beta * Y + alpha * X(IDX(:))
  !

  subroutine s_base_gthab(n,idx,alpha,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    real(psb_spk_) :: alpha, beta, y(:)
    class(psb_s_base_vect_type) :: x
    
    call x%sync()
    call psi_gth(n,idx,alpha,x%v,beta,y)

  end subroutine s_base_gthab
  !
  ! shortcut alpha=1 beta=0
  ! 
  subroutine s_base_gthzv(n,idx,x,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    real(psb_spk_) ::  y(:)
    class(psb_s_base_vect_type) :: x
    
    call x%sync()
    call psi_gth(n,idx,x%v,y)

  end subroutine s_base_gthzv

  !
  ! Scatter: 
  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
  ! 
  
  subroutine s_base_sctb(n,idx,x,beta,y)
    use psi_serial_mod
    integer(psb_ipk_) :: n, idx(:)
    real(psb_spk_) :: beta, x(:)
    class(psb_s_base_vect_type) :: y
    
    call y%sync()
    call psi_sct(n,idx,x,beta,y%v)
    call y%set_host()

  end subroutine s_base_sctb

end module psb_s_base_vect_mod