Implementation of d type multivec-vect dot products

base/modules/psblas/psb_d_psblas_mod.F90

@@ -32,6 +32,7 @@
 module psb_d_psblas_mod
   use psb_desc_mod, only : psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_
   use psb_d_vect_mod, only : psb_d_vect_type
+  use psb_d_multivect_mod, only : psb_d_multivect_type
   use psb_d_mat_mod, only : psb_dspmat_type
 
   interface psb_gedot
@@ -63,6 +64,26 @@ module psb_d_psblas_mod
       integer(psb_ipk_), intent(out)     :: info
       logical, intent(in), optional      :: global
     end function psb_ddot
+    function psb_ddot_multivect(x, y, desc_a,info,global) result(res)
+      import :: psb_desc_type, psb_dpk_, psb_ipk_, &
+           & psb_d_multivect_type, psb_dspmat_type
+      real(psb_dpk_), dimension(:), allocatable :: res
+      type(psb_d_multivect_type), intent(inout) :: x, y
+      type(psb_desc_type), intent(in)      :: desc_a
+      integer(psb_ipk_), intent(out)       :: info
+      logical, intent(in), optional        :: global
+    end function psb_ddot_multivect
+    function psb_ddot_mvect_vect(x, y, desc_a,info,global) result(res)
+      import :: psb_desc_type, psb_dpk_, psb_ipk_, &
+        & psb_d_multivect_type, psb_dspmat_type, &
+        & psb_d_vect_type
+      real(psb_dpk_), dimension(:), allocatable :: res
+      type(psb_d_multivect_type), intent(inout) :: x
+      type(psb_d_vect_type), intent(inout) :: y
+      type(psb_desc_type), intent(in)      :: desc_a
+      integer(psb_ipk_), intent(out)       :: info
+      logical, intent(in), optional        :: global
+    end function psb_ddot_mvect_vect
   end interface
 
 

base/modules/serial/psb_d_base_vect_mod.F90

@@ -149,7 +149,8 @@ module psb_d_base_vect_mod
     !
     procedure, pass(x) :: dot_v    => d_base_dot_v
     procedure, pass(x) :: dot_a    => d_base_dot_a
-    generic, public    :: dot      => dot_v, dot_a
+    procedure, pass(x) :: dot_a2   => d_base_dot_a2
+    generic, public    :: dot      => dot_v, dot_a, dot_a2
     procedure, pass(y) :: axpby_v  => d_base_axpby_v
     procedure, pass(y) :: axpby_a  => d_base_axpby_a
     procedure, pass(z) :: axpby_v2  => d_base_axpby_v2
@@ -1017,6 +1018,34 @@ contains
 
   end function d_base_dot_a
 
+    !
+  ! Base workhorse is good old BLAS1
+  !
+  !
+  !> Function  base_dot_a
+  !! \memberof  psb_d_base_vect_type
+  !! \brief  Dot product by a normal array
+  !! \param n    Number of entries to be considered
+  !! \param y(:,:) The matrix to be multiplied by
+  !!
+  function d_base_dot_a2(n,x,y) result(res)
+    implicit none
+    class(psb_d_base_vect_type), intent(inout) :: x
+    real(psb_dpk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), allocatable, dimension(:) :: res
+    real(psb_dpk_), external      :: ddot
+
+    ! local
+    integer(psb_ipk_) :: ncol
+
+    ncol = size(y,2)
+    allocate(res(ncol))
+
+    call dgemv('T',n,ncol,done,y,n,x%v,1,dzero,res,1)
+
+  end function d_base_dot_a2
+
   !
   ! AXPBY is invoked via Y, hence the structure below.
   !
@@ -2312,7 +2341,8 @@ module psb_d_base_multivect_mod
     !
     procedure, pass(x) :: dot_v    => d_base_mlv_dot_v
     procedure, pass(x) :: dot_a    => d_base_mlv_dot_a
-    generic, public    :: dot      => dot_v, dot_a
+    procedure, pass(x) :: dot_vect => d_base_mlv_dot_vect
+    generic, public    :: dot      => dot_v, dot_a, dot_vect
     procedure, pass(y) :: axpby_v  => d_base_mlv_axpby_v
     procedure, pass(y) :: axpby_a  => d_base_mlv_axpby_a
     generic, public    :: axpby    => axpby_v, axpby_a
@@ -2883,7 +2913,6 @@ contains
     integer(psb_ipk_) :: j,nc
 
     if (x%is_dev()) call x%sync()
-    res = dzero
     !
     ! Note: this is the base implementation.
     !  When we get here, we are sure that X is of
@@ -2905,6 +2934,36 @@ contains
 
   end function d_base_mlv_dot_v
 
+  !> Function  d_base_mlv_dot_vect
+  !! \memberof  psb_d_base_multivect_type
+  !! \brief  Dot product by another base_mlv_vector
+  !! \param n    Number of entries to be considered
+  !! \param y    The other (base_vect) to be multiplied by
+  !!
+  function d_base_mlv_dot_vect(n,x,y) result(res)
+    implicit none
+    class(psb_d_base_multivect_type), intent(inout) :: x
+    class(psb_d_base_vect_type), intent(inout) :: y
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), allocatable   :: res(:)
+    real(psb_dpk_), external      :: ddot
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+
+    select type(yy => y)
+    type is (psb_d_base_vect_type)
+      nc = psb_size(x%v,2_psb_ipk_)
+      allocate(res(nc))
+      do j=1,nc
+        res(j) = ddot(n,x%v(:,j),1,y%v,1)
+      end do
+    class default
+      res = y%dot(n,x%v)
+    end select
+
+  end function d_base_mlv_dot_vect
+
   !
   ! Base workhorse is good old BLAS1
   !

base/modules/serial/psb_d_vect_mod.F90

@@ -1366,6 +1366,7 @@ end module psb_d_vect_mod
 module psb_d_multivect_mod
 
   use psb_d_base_multivect_mod
+  use psb_d_vect_mod
   use psb_const_mod
   use psb_i_vect_mod
 
@@ -1388,11 +1389,18 @@ module psb_d_multivect_mod
     procedure, pass(x) :: get_dupl => d_mvect_get_dupl
     procedure, pass(x) :: set_dupl => d_mvect_set_dupl
 
+    procedure, pass(x) :: sync     => d_mvect_sync
+    procedure, pass(x) :: is_host  => d_mvect_is_host
+    procedure, pass(x) :: is_dev   => d_mvect_is_dev
+    procedure, pass(x) :: is_sync  => d_mvect_is_sync
+    procedure, pass(x) :: set_host => d_mvect_set_host
+    procedure, pass(x) :: set_dev  => d_mvect_set_dev
+    procedure, pass(x) :: set_sync => d_mvect_set_sync
+
     procedure, pass(x) :: all      => d_mvect_all
     procedure, pass(x) :: reall    => d_mvect_reall
     procedure, pass(x) :: zero     => d_mvect_zero
     procedure, pass(x) :: asb      => d_mvect_asb
-    procedure, pass(x) :: sync     => d_mvect_sync
     procedure, pass(x) :: free     => d_mvect_free
     procedure, pass(x) :: ins      => d_mvect_ins
     procedure, pass(x) :: bld_x    => d_mvect_bld_x
@@ -1411,9 +1419,10 @@ module psb_d_multivect_mod
     procedure, pass(y) :: sctb     => d_mvect_sctb
     procedure, pass(y) :: sctb_x   => d_mvect_sctb_x
     generic, public    :: sct      => sctb, sctb_x
-!!$    procedure, pass(x) :: dot_v    => d_mvect_dot_v
-!!$    procedure, pass(x) :: dot_a    => d_mvect_dot_a
-!!$    generic, public    :: dot      => dot_v, dot_a
+    procedure, pass(x) :: dot_v    => d_mvect_dot_v
+    procedure, pass(x) :: dot_a    => d_mvect_dot_a
+    procedure, pass(x) :: dot_vect => d_mvect_dot_vect
+    generic, public    :: dot      => dot_v, dot_a, dot_vect
 !!$    procedure, pass(y) :: axpby_v  => d_mvect_axpby_v
 !!$    procedure, pass(y) :: axpby_a  => d_mvect_axpby_a
 !!$    generic, public    :: axpby    => axpby_v, axpby_a
@@ -1473,6 +1482,75 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine d_mvect_set_dupl
+
+  subroutine d_mvect_sync(x)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    if (allocated(x%v)) &
+         & call x%v%sync()
+
+  end subroutine d_mvect_sync
+
+  subroutine d_mvect_set_sync(x)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    if (allocated(x%v)) &
+         & call x%v%set_sync()
+
+  end subroutine d_mvect_set_sync
+
+  subroutine d_mvect_set_host(x)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    if (allocated(x%v)) &
+         & call x%v%set_host()
+
+  end subroutine d_mvect_set_host
+
+  subroutine d_mvect_set_dev(x)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    if (allocated(x%v)) &
+         & call x%v%set_dev()
+
+  end subroutine d_mvect_set_dev
+
+  function d_mvect_is_sync(x) result(res)
+    implicit none
+    logical :: res
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    res = .true.
+    if (allocated(x%v)) &
+         & res = x%v%is_sync()
+
+  end function d_mvect_is_sync
+
+  function d_mvect_is_host(x) result(res)
+    implicit none
+    logical :: res
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    res = .true.
+    if (allocated(x%v)) &
+         & res = x%v%is_host()
+
+  end function d_mvect_is_host
+
+  function d_mvect_is_dev(x) result(res)
+    implicit none
+    logical :: res
+    class(psb_d_multivect_type), intent(inout) :: x
+
+    res = .false.
+    if (allocated(x%v)) &
+         & res =  x%v%is_dev()
+
+  end function d_mvect_is_dev
         
 
   function d_mvect_is_remote_build(x) result(res)
@@ -1717,15 +1795,6 @@ contains
 
   end subroutine d_mvect_asb
 
-  subroutine d_mvect_sync(x)
-    implicit none
-    class(psb_d_multivect_type), intent(inout) :: x
-
-    if (allocated(x%v)) &
-         & call x%v%sync()
-
-  end subroutine d_mvect_sync
-
   subroutine d_mvect_gthab(n,idx,alpha,x,beta,y)
     use psi_serial_mod
     integer(psb_ipk_) :: n, idx(:)
@@ -1840,30 +1909,50 @@ contains
   end subroutine d_mvect_cnv
 
 
-!!$  function d_mvect_dot_v(n,x,y) result(res)
-!!$    implicit none
-!!$    class(psb_d_multivect_type), intent(inout) :: x, y
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_dpk_)                :: res
-!!$
-!!$    res = dzero
-!!$    if (allocated(x%v).and.allocated(y%v)) &
-!!$         & res = x%v%dot(n,y%v)
-!!$
-!!$  end function d_mvect_dot_v
-!!$
-!!$  function d_mvect_dot_a(n,x,y) result(res)
-!!$    implicit none
-!!$    class(psb_d_multivect_type), intent(inout) :: x
-!!$    real(psb_dpk_), intent(in)    :: y(:)
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_dpk_)                :: res
-!!$
-!!$    res = dzero
-!!$    if (allocated(x%v)) &
-!!$         & res = x%v%dot(n,y)
-!!$
-!!$  end function d_mvect_dot_a
+  function d_mvect_dot_v(n,x,y) result(res)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x, y
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), dimension(:), allocatable :: res
+
+    if (allocated(x%v).and.allocated(y%v)) then 
+      res = x%v%dot(n,y%v)
+    else
+      allocate(res(1))
+      res(1) = psb_err_invalid_vect_state_
+    end if
+
+  end function d_mvect_dot_v
+
+  function d_mvect_dot_vect(n,x,y) result(res)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+    class(psb_d_vect_type), intent(inout) :: y
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), dimension(:), allocatable :: res
+
+    if (allocated(x%v).and.allocated(y%v)) then 
+      res = x%v%dot(n,y%v)
+    else
+      allocate(res(1))
+      res(1) = psb_err_invalid_vect_state_
+    end if
+
+  end function d_mvect_dot_vect
+
+  function d_mvect_dot_a(n,x,y) result(res)
+    implicit none
+    class(psb_d_multivect_type), intent(inout) :: x
+    real(psb_dpk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), dimension(:), allocatable :: res
+
+    res = dzero
+    if (allocated(x%v)) then
+      res = x%v%dot(n,y)
+    end if
+
+  end function d_mvect_dot_a
 !!$
 !!$  subroutine d_mvect_axpby_v(m,alpha, x, beta, y, info)
 !!$    use psi_serial_mod

base/psblas/psb_ddot.f90

@@ -157,6 +157,347 @@ function psb_ddot_vect(x, y, desc_a,info,global) result(res)
   return
 
 end function psb_ddot_vect
+!   
+!                Parallel Sparse BLAS  version 3.5
+!      (C) Copyright 2006-2018
+!        Salvatore Filippone    
+!        Alfredo Buttari      
+!   
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific written permission.
+!   
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!   
+!    
+! File: psb_ddot.f90
+!
+! Function: psb_ddot_vect
+!    psb_ddot computes the dot product of two distributed vectors,
+!
+!    dot := ( X )**C * ( Y )
+!
+!
+! Arguments:
+!    x      -  type(psb_d_vect_type) The input vector containing the entries of sub( X ).
+!    y      -  type(psb_d_vect_type) The input vector containing the entries of sub( Y ).
+!    desc_a -  type(psb_desc_type).  The communication descriptor.
+!    info   -  integer.              Return code
+!    global -  logical(optional)     Whether to perform the global sum, default: .true.
+!
+!  Note: from a functional point of view, X and Y are input, but here
+!        they are declared INOUT because of the sync() methods. 
+!
+!
+function psb_ddot_multivect(x, y, desc_a,info,global) result(res)
+  use psb_desc_mod
+  use psb_d_base_mat_mod
+  use psb_check_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_d_multivect_mod
+  use psb_d_psblas_mod, psb_protect_name => psb_ddot_multivect
+  implicit none 
+  real(psb_dpk_), dimension(:), allocatable :: res
+  type(psb_d_multivect_type), intent(inout) :: x, y
+  type(psb_desc_type), intent(in)      :: desc_a
+  integer(psb_ipk_), intent(out)       :: info
+  logical, intent(in), optional        :: global
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me, idx, ndm,&
+       & err_act, iix, jjx, iiy, jjy, i, nr
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n
+  logical :: global_
+  character(len=20)      :: name, ch_err
+
+  name='psb_ddot_multivect'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+
+  ctxt=desc_a%get_context()
+  call psb_info(ctxt, me, np)
+  if (np == -ione) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  if (.not.allocated(x%v)) then 
+    info = psb_err_invalid_vect_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  if (.not.allocated(y%v)) then 
+    info = psb_err_invalid_vect_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  if (present(global)) then
+    global_ = global
+  else
+    global_ = .true.
+  end if
+
+  ix = ione
+  ijx = ione
+
+  iy = ione
+  ijy = ione
+
+  m = desc_a%get_global_rows()
+
+  ! check vector correctness
+  call psb_chkvect(m,x%get_ncols(),x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  if (info == psb_success_) &
+       & call psb_chkvect(m,y%get_ncols(),y%get_nrows(),iy,ijy,desc_a,info,iiy,jjy)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkvect'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  if ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  if (x%get_ncols() /= y%get_ncols()) then
+    info=psb_err_invalid_vect_state_
+    call psb_errpush(info,name)
+    goto 9999
+  else
+    allocate(res(x%get_ncols()),stat=info)
+    if (info /= 0) then
+      info=psb_err_alloc_dealloc_
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+  end if
+
+  nr = desc_a%get_local_rows() 
+  if(nr > 0) then
+    res = x%dot(nr,y)
+    ! FIXME
+    ! adjust dot_local because overlapped elements are computed more than once
+    if (size(desc_a%ovrlap_elem,1)>0) then
+      if (x%is_dev()) call x%sync()
+      if (y%is_dev()) call y%sync()
+      do i=1,size(desc_a%ovrlap_elem,1)
+        idx = desc_a%ovrlap_elem(i,1)
+        ndm = desc_a%ovrlap_elem(i,2)
+        res(:) = res(:) - (real(ndm-1)/real(ndm))*(x%v%v(idx,:)*y%v%v(idx,:))
+      end do
+    end if
+  else
+    res = dzero
+  end if
+
+  ! compute global sum
+  if (global_) call psb_sum(ctxt, res)
+
+  call psb_erractionrestore(err_act)
+  return  
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end function psb_ddot_multivect
+!   
+!                Parallel Sparse BLAS  version 3.5
+!      (C) Copyright 2006-2018
+!        Salvatore Filippone    
+!        Alfredo Buttari      
+!   
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific written permission.
+!   
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!   
+!    
+! File: psb_ddot.f90
+!
+! Function: psb_ddot_vect
+!    psb_ddot computes the dot product of two distributed vectors,
+!
+!    dot := ( X )**C * ( Y )
+!
+!
+! Arguments:
+!    x      -  type(psb_d_vect_type) The input vector containing the entries of sub( X ).
+!    y      -  type(psb_d_vect_type) The input vector containing the entries of sub( Y ).
+!    desc_a -  type(psb_desc_type).  The communication descriptor.
+!    info   -  integer.              Return code
+!    global -  logical(optional)     Whether to perform the global sum, default: .true.
+!
+!  Note: from a functional point of view, X and Y are input, but here
+!        they are declared INOUT because of the sync() methods. 
+!
+!
+function psb_ddot_mvect_vect(x, y, desc_a,info,global) result(res)
+  use psb_desc_mod
+  use psb_d_base_mat_mod
+  use psb_check_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_d_multivect_mod
+  use psb_d_vect_mod
+  use psb_d_psblas_mod, psb_protect_name => psb_ddot_mvect_vect
+  implicit none 
+  real(psb_dpk_), dimension(:), allocatable :: res
+  type(psb_d_multivect_type), intent(inout) :: x
+  type(psb_d_vect_type), intent(inout) :: y
+  type(psb_desc_type), intent(in)      :: desc_a
+  integer(psb_ipk_), intent(out)       :: info
+  logical, intent(in), optional        :: global
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me, idx, ndm,&
+       & err_act, iix, jjx, iiy, jjy, i, nr
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m, n
+  logical :: global_
+  character(len=20)      :: name, ch_err
+
+  name='psb_ddot_mvect_vect'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+
+  ctxt=desc_a%get_context()
+  call psb_info(ctxt, me, np)
+  if (np == -ione) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  if (.not.allocated(x%v)) then 
+    info = psb_err_invalid_vect_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  if (.not.allocated(y%v)) then 
+    info = psb_err_invalid_vect_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  if (present(global)) then
+    global_ = global
+  else
+    global_ = .true.
+  end if
+
+  ix = ione
+  ijx = ione
+
+  iy = ione
+  ijy = ione
+
+  m = desc_a%get_global_rows()
+
+  ! check vector correctness
+  call psb_chkvect(m,x%get_ncols(),x%get_nrows(),ix,ijx,desc_a,info,iix,jjx)
+  if (info == psb_success_) &
+       & call psb_chkvect(m,lone,y%get_nrows(),iy,ijy,desc_a,info,iiy,jjy)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkvect'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  if ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  allocate(res(x%get_ncols()),stat=info)
+  if (info /= 0) then
+    info=psb_err_alloc_dealloc_
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  nr = desc_a%get_local_rows() 
+  if(nr > 0) then
+    res = x%dot(nr,y)
+    ! FIXME
+    ! adjust dot_local because overlapped elements are computed more than once
+    if (size(desc_a%ovrlap_elem,1)>0) then
+      if (x%is_dev()) call x%sync()
+      if (y%is_dev()) call y%sync()
+      do i=1,size(desc_a%ovrlap_elem,1)
+        idx = desc_a%ovrlap_elem(i,1)
+        ndm = desc_a%ovrlap_elem(i,2)
+        ! Remove the overlapped elements via dgemv calls 
+        ! res = - (real(ndm-1)/real(ndm))* x(idx,:)^T y(idx)  + 1.0 res
+        call dgemv('T',size(x%v%v,1),size(x%v%v,2),-(real(ndm-1)/real(ndm)), &
+          & size(x%v%v,1),y%v%v(idx),ione,done,res,ione)
+      end do
+    end if
+  else
+    res = dzero
+  end if
+
+  ! compute global sum
+  if (global_) call psb_sum(ctxt, res)
+
+  call psb_erractionrestore(err_act)
+  return  
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end function psb_ddot_mvect_vect
 !
 ! Function: psb_ddot
 !    psb_ddot computes the dot product of two distributed vectors,

test/kernel/Makefile

@@ -18,7 +18,7 @@ DPGOBJS=pdgenspmv.o
 DVECOBJS=vecoperation.o
 EXEDIR=./runs
 
-all: runsd pdgenspmv
+all: runsd pdgenspmv vecoperation
 #d_file_spmv s_file_spmv pdgenspmv vecoperation
 
 runsd:

test/kernel/vecoperation.f90

@@ -35,15 +35,20 @@ module unittestvector_mod
 
   use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,&
        &  psb_dspmat_type, psb_d_vect_type, dzero, psb_ctxt_type,&
-       &  psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type
+       &  psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type,&
+       &  psb_d_base_multivect_type
 
   interface psb_gen_const
-   module procedure  psb_d_gen_const
+   module procedure  psb_d_gen_const, psb_d_gen_const_multi
   end interface psb_gen_const
 
+  interface psb_check_ans
+   module procedure  psb_check_ans_v, psb_check_ans_mv
+  end interface psb_check_ans
+
 contains
 
-  function psb_check_ans(v,val,ctxt) result(ans)
+  function psb_check_ans_v(v,val,ctxt) result(ans)
     use psb_base_mod
 
     implicit none
@@ -73,7 +78,39 @@ contains
       ans = .false.
     end if
 
-  end function psb_check_ans
+  end function psb_check_ans_v
+
+  function psb_check_ans_mv(v,val,ctxt) result(ans)
+    use psb_base_mod
+
+    implicit none
+
+    type(psb_d_multivect_type) :: v
+    real(psb_dpk_)        :: val
+    type(psb_ctxt_type) :: ctxt
+    logical               :: ans
+
+    ! Local variables
+    integer(psb_ipk_) :: np, iam, info
+    real(psb_dpk_)    :: check
+    real(psb_dpk_), allocatable    :: va(:,:)
+
+    call psb_info(ctxt,iam,np)
+
+    va = v%get_vect()
+    va = va - val;
+
+    check = maxval(va);
+
+    call psb_sum(ctxt,check)
+
+    if(check == 0.d0) then
+      ans = .true.
+    else
+      ans = .false.
+    end if
+
+  end function psb_check_ans_mv
   !
   !  subroutine to fill a vector with constant entries
   !
@@ -160,6 +197,93 @@ contains
 
     return
   end subroutine psb_d_gen_const
+    !
+  !  subroutine to fill a multivectorvector with constant entries
+  !
+  subroutine psb_d_gen_const_multi(v,val,idim,jdim,ctxt,desc_a,info)
+    use psb_base_mod
+    implicit none
+
+    type(psb_d_multivect_type) :: v
+    type(psb_desc_type)   :: desc_a
+    integer(psb_lpk_)     :: idim
+    integer(psb_ipk_)     :: jdim !! Number of columns of the multivector
+    type(psb_ctxt_type) :: ctxt
+    integer(psb_ipk_)     :: info
+    real(psb_dpk_)        :: val
+
+    ! Local variables
+    integer(psb_ipk_), parameter    :: nb=20
+    real(psb_dpk_)                  :: zt(nb,jdim) ! Temporary array to fill the vector
+    character(len=40)               :: name, ch_err
+    integer(psb_ipk_)               :: np, iam, nr, nt
+    integer(psb_ipk_)               :: n,nlr,ib,ii
+    integer(psb_ipk_)               :: err_act
+    integer(psb_lpk_), allocatable  :: myidx(:)
+
+
+    info = psb_success_
+    name = 'create_constant_multivector'
+    call psb_erractionsave(err_act)
+
+    call psb_info(ctxt, iam, np)
+
+    n = idim*np         ! The global dimension is the number of process times
+                        ! the input size
+
+    ! We use a simple minded block distribution
+    nt = (n+np-1)/np
+    nr = max(0,min(nt,n-(iam*nt)))
+    nt = nr
+
+    call psb_sum(ctxt,nt)
+    if (nt /= n) then
+      write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,n
+      info = -1
+      call psb_barrier(ctxt)
+      call psb_abort(ctxt)
+      return
+    end if
+    ! Allocate the descriptor with simple minded data distribution
+    call psb_cdall(ctxt,desc_a,info,nl=nr)
+    ! Allocate the vector on the recently build descriptor
+    if (info == psb_success_) call psb_geall(v,desc_a,info,n=jdim)
+    ! Check that allocation has gone good
+    if (info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='allocation rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    do ii=1,nlr,nb
+      ib = min(nb,nlr-ii+1)
+      zt(:,:) = val
+      call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib,1:jdim),v,desc_a,info)
+      if(info /= psb_success_) exit
+    end do
+
+    if(info /= psb_success_) then
+      info=psb_err_from_subroutine_
+      ch_err='insert rout.'
+      call psb_errpush(info,name,a_err=ch_err)
+      goto 9999
+    end if
+
+    ! Assembly of communicator and vector
+    call psb_cdasb(desc_a,info)
+    if (info == psb_success_) call psb_geasb(v,desc_a,info)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+    return
+  end subroutine psb_d_gen_const_multi
 
 end module unittestvector_mod
 
@@ -171,7 +295,8 @@ program vecoperation
   implicit none
 
   ! input parameters
-  integer(psb_lpk_) :: idim = 100
+  integer(psb_lpk_) :: idim = 100 ! Local vector size
+  integer(psb_ipk_) :: nmv = 10 ! Number of columns of the multivector
 
   ! miscellaneous
   real(psb_dpk_), parameter :: one = 1.d0
@@ -184,13 +309,17 @@ program vecoperation
   type(psb_desc_type)   :: desc_a
   ! vector
   type(psb_d_vect_type)  :: x,y,z
+  ! multivector
+  type(psb_d_multivect_type) :: mv1, mv2
   ! blacs parameters
   type(psb_ctxt_type) :: ctxt
   integer(psb_ipk_) :: iam, np
   ! auxiliary parameters
+  integer(psb_ipk_) :: ii
   integer(psb_ipk_) :: info
   character(len=20) :: name,ch_err,readinput
   real(psb_dpk_)    :: ans
+  real(psb_dpk_), allocatable :: ansmv(:)
   logical           :: hasitnotfailed
   integer(psb_lpk_), allocatable  :: myidx(:)
   integer(psb_ipk_) :: ib = 1
@@ -363,9 +492,54 @@ program vecoperation
     if(ans /= onehalf) write(psb_out_unit,'("TEST FAILED --- MaxNorm")')
   end if
 
+
+  !
+  ! Test of multivector operation
+  !
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+  if (iam == psb_root_) write(psb_out_unit,'("Multivector Operations")')
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+  ! X = 1
+  call psb_d_gen_const_multi(mv1,one,idim,nmv,ctxt,desc_a,info)
+  hasitnotfailed = psb_check_ans(mv1,one,ctxt)
+  if (iam == psb_root_) then
+    if(hasitnotfailed) write(psb_out_unit,'("TEST PASSED >>> Constant multivector ")')
+    if(.not.hasitnotfailed) write(psb_out_unit,'("TEST FAILED --- Constant multivector ")')
+  end if
+
+  ! 
+  ! Multivector to field operation
+  !
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+  if (iam == psb_root_) write(psb_out_unit,'("Multivector to Field Operations")')
+  if (iam == psb_root_) write(psb_out_unit,'(" ")')
+
+  ! Dot product: multivector vs multivector
+  call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info)
+  call psb_d_gen_const_multi(mv2,onehalf,idim,nmv,ctxt,desc_a,info)
+  ansmv = psb_gedot(mv1,mv2,desc_a,info)
+  if (iam == psb_root_) then
+    ! write ansmv to check
+    if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")')
+    if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")')
+  end if
+  ! Dot product: multivector vs vector
+  call psb_d_gen_const_multi(mv1,two,idim,nmv,ctxt,desc_a,info)
+  call psb_d_gen_const(x,onehalf,idim,ctxt,desc_a,info)
+  ansmv = psb_gedot(mv1,x,desc_a,info)
+  if (iam == psb_root_) then
+    ! write ansmv to check
+    if(all(ansmv(:) == np*idim)) write(psb_out_unit,'("TEST PASSED >>> Dot product")')
+    if(any(ansmv(:) /= np*idim)) write(psb_out_unit,'("TEST FAILED --- Dot product")')
+  end if
+
+
+
   call psb_gefree(x,desc_a,info)
   call psb_gefree(y,desc_a,info)
   call psb_gefree(z,desc_a,info)
+  call psb_gefree(mv1,desc_a,info)
+  call psb_gefree(mv2,desc_a,info)
   call psb_cdfree(desc_a,info)
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_