diff --git a/base/modules/psblas/psb_c_psblas_mod.F90 b/base/modules/psblas/psb_c_psblas_mod.F90
index 98deebd8..a8c244d8 100644
--- a/base/modules/psblas/psb_c_psblas_mod.F90
+++ b/base/modules/psblas/psb_c_psblas_mod.F90
@@ -33,6 +33,7 @@ module psb_c_psblas_mod
   use psb_desc_mod, only : psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_
   use psb_c_vect_mod, only : psb_c_vect_type
   use psb_c_mat_mod, only : psb_cspmat_type
+  use psb_c_multivect_mod, only: psb_c_multivect_type
 
   interface psb_gedot
     function psb_cdot_vect(x, y, desc_a,info,global) result(res)
@@ -141,6 +142,17 @@ module psb_c_psblas_mod
       integer(psb_ipk_), optional, intent(in) :: n, jx, jy
       integer(psb_ipk_), intent(out)      :: info
     end subroutine psb_caxpby
+    subroutine psb_caxpby_multivect_vect(alpha, x, beta, y,&
+         & j, desc_a, info)
+      import :: psb_desc_type, psb_spk_, psb_ipk_, &
+              & psb_c_vect_type, psb_cspmat_type, psb_c_multivect_type
+      type(psb_c_vect_type), intent (inout) ::  x
+      type(psb_c_multivect_type), intent (inout) ::  y
+      complex(psb_spk_), intent (in)        :: alpha, beta
+      integer(psb_ipk_), intent(in)        :: j
+      type(psb_desc_type), intent (in)      :: desc_a
+      integer(psb_ipk_), intent(out)                  :: info
+    end subroutine psb_caxpby_multivect_vect
   end interface
 
   interface psb_geamax
diff --git a/base/modules/psblas/psb_d_psblas_mod.F90 b/base/modules/psblas/psb_d_psblas_mod.F90
index e4988387..93c63c79 100644
--- a/base/modules/psblas/psb_d_psblas_mod.F90
+++ b/base/modules/psblas/psb_d_psblas_mod.F90
@@ -33,6 +33,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_mat_mod, only : psb_dspmat_type
+  use psb_d_multivect_mod, only: psb_d_multivect_type
 
   interface psb_gedot
     function psb_ddot_vect(x, y, desc_a,info,global) result(res)
@@ -141,6 +142,16 @@ module psb_d_psblas_mod
       integer(psb_ipk_), optional, intent(in) :: n, jx, jy
       integer(psb_ipk_), intent(out)      :: info
     end subroutine psb_daxpby
+    subroutine psb_daxpby_multivect_vect(alpha, x, beta, y, j, desc_a, info)
+      import :: psb_desc_type, psb_dpk_, psb_ipk_, &
+              & psb_d_vect_type, psb_dspmat_type, psb_d_multivect_type
+      type(psb_d_vect_type), intent (inout) ::  x
+      type(psb_d_multivect_type), intent (inout) ::  y
+      real(psb_dpk_), intent (in)        :: alpha, beta
+      integer(psb_ipk_), intent(in)        :: j
+      type(psb_desc_type), intent (in)      :: desc_a
+      integer(psb_ipk_), intent(out)                  :: info
+    end subroutine psb_daxpby_multivect_vect
   end interface
 
   interface psb_geamax
diff --git a/base/modules/psblas/psb_s_psblas_mod.F90 b/base/modules/psblas/psb_s_psblas_mod.F90
index 93fe74b9..ebecee29 100644
--- a/base/modules/psblas/psb_s_psblas_mod.F90
+++ b/base/modules/psblas/psb_s_psblas_mod.F90
@@ -33,6 +33,7 @@ module psb_s_psblas_mod
   use psb_desc_mod, only : psb_desc_type, psb_spk_, psb_ipk_, psb_lpk_
   use psb_s_vect_mod, only : psb_s_vect_type
   use psb_s_mat_mod, only : psb_sspmat_type
+  use psb_s_multivect_mod, only: psb_s_multivect_type
 
   interface psb_gedot
     function psb_sdot_vect(x, y, desc_a,info,global) result(res)
@@ -141,6 +142,17 @@ module psb_s_psblas_mod
       integer(psb_ipk_), optional, intent(in) :: n, jx, jy
       integer(psb_ipk_), intent(out)      :: info
     end subroutine psb_saxpby
+    subroutine psb_saxpby_multivect_vect(alpha, x, beta, y,&
+         & j, desc_a, info)
+      import :: psb_desc_type, psb_spk_, psb_ipk_, &
+              & psb_s_vect_type, psb_sspmat_type, psb_s_multivect_type
+      type(psb_s_vect_type), intent (inout) ::  x
+      type(psb_s_multivect_type), intent (inout) ::  y
+      real(psb_spk_), intent (in)        :: alpha, beta
+      integer(psb_ipk_), intent(in)        :: j
+      type(psb_desc_type), intent (in)      :: desc_a
+      integer(psb_ipk_), intent(out)                  :: info
+    end subroutine psb_saxpby_multivect_vect
   end interface
 
   interface psb_geamax
diff --git a/base/modules/psblas/psb_z_psblas_mod.F90 b/base/modules/psblas/psb_z_psblas_mod.F90
index 06be1b82..55e1698e 100644
--- a/base/modules/psblas/psb_z_psblas_mod.F90
+++ b/base/modules/psblas/psb_z_psblas_mod.F90
@@ -33,6 +33,7 @@ module psb_z_psblas_mod
   use psb_desc_mod, only : psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_
   use psb_z_vect_mod, only : psb_z_vect_type
   use psb_z_mat_mod, only : psb_zspmat_type
+  use psb_z_multivect_mod, only: psb_z_multivect_type
 
   interface psb_gedot
     function psb_zdot_vect(x, y, desc_a,info,global) result(res)
@@ -141,6 +142,17 @@ module psb_z_psblas_mod
       integer(psb_ipk_), optional, intent(in) :: n, jx, jy
       integer(psb_ipk_), intent(out)      :: info
     end subroutine psb_zaxpby
+    subroutine psb_zaxpby_multivect_vect(alpha, x, beta, y,&
+         & j, desc_a, info)
+      import :: psb_desc_type, psb_dpk_, psb_ipk_, &
+              & psb_z_vect_type, psb_zspmat_type, psb_z_multivect_type
+      type(psb_z_vect_type), intent (inout) ::  x
+      type(psb_z_multivect_type), intent (inout) ::  y
+      complex(psb_dpk_), intent (in)        :: alpha, beta
+      integer(psb_ipk_), intent(in)        :: j
+      type(psb_desc_type), intent (in)      :: desc_a
+      integer(psb_ipk_), intent(out)                  :: info
+    end subroutine psb_zaxpby_multivect_vect
   end interface
 
   interface psb_geamax
diff --git a/base/modules/serial/psb_c_base_vect_mod.F90 b/base/modules/serial/psb_c_base_vect_mod.F90
index 44044771..373a7d10 100644
--- a/base/modules/serial/psb_c_base_vect_mod.F90
+++ b/base/modules/serial/psb_c_base_vect_mod.F90
@@ -267,7 +267,7 @@ contains
     class(psb_c_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -800,7 +800,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -808,7 +808,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function c_base_get_vect
 
   !
@@ -836,7 +836,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -864,7 +864,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
@@ -918,7 +918,7 @@ contains
     class(psb_c_base_vect_type), intent(inout)  :: x
 
     integer(psb_ipk_) :: i
-    
+
     if (allocated(x%v)) then
       if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
@@ -1170,7 +1170,7 @@ contains
     info = 0
     if (y%is_dev()) call y%sync()
     n = min(size(y%v), size(x))
-    !$omp parallel do private(i)    
+    !$omp parallel do private(i)
     do i=1, n
       y%v(i) = y%v(i)*x(i)
     end do
@@ -1216,7 +1216,7 @@ contains
     else
       if (alpha == cone) then
         if (beta == czero) then
-          !$omp parallel do private(i)    
+          !$omp parallel do private(i)
           do i=1, n
             z%v(i) = y(i)*x(i)
           end do
@@ -1681,7 +1681,7 @@ contains
     end if
     call x%set_host()
   end subroutine c_base_scal
-  
+
   !
   ! Norms 1, 2 and infinity
   !
@@ -1737,7 +1737,7 @@ contains
     integer(psb_ipk_), intent(in)           :: n
     real(psb_spk_)                :: res
     integer(psb_ipk_) :: i
-    
+
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     res=szero
@@ -1964,7 +1964,7 @@ contains
     z%v = x + b
 #endif
     info = 0
-    
+
   end subroutine c_base_addconst_a2
   !
   !> Function  _base_addconst_v2
@@ -2070,8 +2070,9 @@ module psb_c_base_multivect_mod
     procedure, pass(x) :: dot_a    => c_base_mlv_dot_a
     generic, public    :: dot      => dot_v, dot_a
     procedure, pass(y) :: axpby_v  => c_base_mlv_axpby_v
+    procedure, pass(y) :: axpby_vv => c_base_mlv_axpby_vv
     procedure, pass(y) :: axpby_a  => c_base_mlv_axpby_a
-    generic, public    :: axpby    => axpby_v, axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a, axpby_vv
     !
     ! MultiVector by vector/multivector multiplication. Need all variants
     ! to handle multiple requirements from preconditioners
@@ -2728,6 +2729,35 @@ contains
 
   end subroutine c_base_mlv_axpby_v
 
+  !
+  ! AXPBY is invoked via Y, hence the structure below.
+  !
+  !
+  !
+  !> Function  base_mlv_axpby_v
+  !! \memberof  psb_c_base_multivect_type
+  !! \brief AXPBY  by a (base_mlv_vect) y(j)=alpha*x+beta*y(j)
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_v_vect) to be added
+  !! \param beta scalar alpha
+  !! \param j to what column of the multivector should we add
+  !! \param info   return code
+  !!
+  subroutine c_base_mlv_axpby_vv(m,alpha, x, beta, y, j, info)
+    use psi_serial_mod
+    implicit none
+    integer(psb_ipk_), intent(in)                 :: m
+    class(psb_c_base_vect_type), intent(inout)  :: x
+    class(psb_c_base_multivect_type), intent(inout)  :: y
+    complex(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(in)      :: j
+    integer(psb_ipk_), intent(out)              :: info
+
+    call psb_geaxpby(m,alpha,x%v,beta,y%v(:,j),info)
+
+  end subroutine c_base_mlv_axpby_vv
+
   !
   ! AXPBY is invoked via Y, hence the structure below.
   !
diff --git a/base/modules/serial/psb_c_vect_mod.F90 b/base/modules/serial/psb_c_vect_mod.F90
index 1a336d11..25d69911 100644
--- a/base/modules/serial/psb_c_vect_mod.F90
+++ b/base/modules/serial/psb_c_vect_mod.F90
@@ -57,7 +57,7 @@ module psb_c_vect_mod
     procedure, pass(x) :: is_remote_build => c_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => c_vect_set_remote_build
     procedure, pass(x) :: get_dupl => c_vect_get_dupl
-    procedure, pass(x) :: set_dupl => c_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => c_vect_set_dupl
     procedure, pass(x) :: get_nrmv => c_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => c_vect_set_nrmv
     procedure, pass(x) :: all      => c_vect_all
@@ -214,7 +214,7 @@ contains
 
     x%nrmv = val
   end subroutine c_vect_set_nrmv
-        
+
 
   function c_vect_is_remote_build(x) result(res)
     implicit none
@@ -234,7 +234,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine c_vect_set_remote_build
-        
+
   subroutine  psb_c_set_vect_default(v)
     implicit none
     class(psb_c_base_vect_type), intent(in) :: v
@@ -592,7 +592,7 @@ contains
       allocate(tmp,stat=info,mold=psb_c_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -1149,7 +1149,7 @@ contains
     ! Temp vectors
     type(psb_c_vect_type) :: wtemp
 
-    info = 0 
+    info = 0
     if( allocated(w%v) ) then
       if (.not.present(aux)) then
         allocate(wtemp%v, mold=w%v)
@@ -1316,6 +1316,8 @@ module psb_c_multivect_mod
 !!$    procedure, pass(x) :: nrm2     => c_vect_nrm2
 !!$    procedure, pass(x) :: amax     => c_vect_amax
 !!$    procedure, pass(x) :: asum     => c_vect_asum
+       procedure, pass(y) :: axpby_vv => c_vect_axpby_vv
+       generic, public :: axpby => axpby_vv
   end type psb_c_multivect_type
 
   public  :: psb_c_multivect, psb_c_multivect_type,&
@@ -1341,7 +1343,7 @@ module psb_c_multivect_mod
 
 contains
 
-  
+
   function c_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_c_multivect_type), intent(in) :: x
@@ -1360,7 +1362,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine c_mvect_set_dupl
-        
+
 
   function c_mvect_is_remote_build(x) result(res)
     implicit none
@@ -1380,7 +1382,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine c_mvect_set_remote_build
-        
+
 
   subroutine  psb_c_set_multivect_default(v)
     implicit none
@@ -1769,6 +1771,24 @@ contains
 !!$
 !!$  end subroutine c_vect_axpby_v
 !!$
+    subroutine c_vect_axpby_vv(m,alpha, x, beta, y, j, info)
+      use psi_serial_mod
+      use psb_c_vect_mod
+      implicit none
+      integer(psb_ipk_), intent(in)            :: m
+      type(psb_c_vect_type), intent(inout)  :: x
+      class(psb_c_multivect_type), intent(inout)  :: y
+      complex(psb_spk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(in)      :: j
+      integer(psb_ipk_), intent(out)           :: info
+
+      if (allocated(x%v).and.allocated(y%v)) then
+         call y%v%axpby(m,alpha,x%v,beta,j,info)
+      else
+         info = psb_err_invalid_vect_state_
+      end if
+
+    end subroutine c_vect_axpby_vv
 !!$  subroutine c_vect_axpby_a(m,alpha, x, beta, y, info)
 !!$    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 a28d12f6..b1d7a392 100644
--- a/base/modules/serial/psb_d_base_vect_mod.F90
+++ b/base/modules/serial/psb_d_base_vect_mod.F90
@@ -274,7 +274,7 @@ contains
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -807,7 +807,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -815,7 +815,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function d_base_get_vect
 
   !
@@ -843,7 +843,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -871,7 +871,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
@@ -925,7 +925,7 @@ contains
     class(psb_d_base_vect_type), intent(inout)  :: x
 
     integer(psb_ipk_) :: i
-    
+
     if (allocated(x%v)) then
       if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
@@ -1177,7 +1177,7 @@ contains
     info = 0
     if (y%is_dev()) call y%sync()
     n = min(size(y%v), size(x))
-    !$omp parallel do private(i)    
+    !$omp parallel do private(i)
     do i=1, n
       y%v(i) = y%v(i)*x(i)
     end do
@@ -1223,7 +1223,7 @@ contains
     else
       if (alpha == done) then
         if (beta == dzero) then
-          !$omp parallel do private(i)    
+          !$omp parallel do private(i)
           do i=1, n
             z%v(i) = y(i)*x(i)
           end do
@@ -1688,7 +1688,7 @@ contains
     end if
     call x%set_host()
   end subroutine d_base_scal
-  
+
   !
   ! Norms 1, 2 and infinity
   !
@@ -1824,7 +1824,7 @@ contains
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
     integer(psb_ipk_) :: i
-    
+
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     res=dzero
@@ -2143,7 +2143,7 @@ contains
     z%v = x + b
 #endif
     info = 0
-    
+
   end subroutine d_base_addconst_a2
   !
   !> Function  _base_addconst_v2
@@ -2249,8 +2249,9 @@ module psb_d_base_multivect_mod
     procedure, pass(x) :: dot_a    => d_base_mlv_dot_a
     generic, public    :: dot      => dot_v, dot_a
     procedure, pass(y) :: axpby_v  => d_base_mlv_axpby_v
+    procedure, pass(y) :: axpby_vv => d_base_mlv_axpby_vv
     procedure, pass(y) :: axpby_a  => d_base_mlv_axpby_a
-    generic, public    :: axpby    => axpby_v, axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a, axpby_vv
     !
     ! MultiVector by vector/multivector multiplication. Need all variants
     ! to handle multiple requirements from preconditioners
@@ -2907,6 +2908,35 @@ contains
 
   end subroutine d_base_mlv_axpby_v
 
+  !
+  ! AXPBY is invoked via Y, hence the structure below.
+  !
+  !
+  !
+  !> Function  base_mlv_axpby_v
+  !! \memberof  psb_d_base_multivect_type
+  !! \brief AXPBY  by a (base_mlv_vect) y(j)=alpha*x+beta*y(j)
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_v_vect) to be added
+  !! \param beta scalar alpha
+  !! \param j to what column of the multivector should we add
+  !! \param info   return code
+  !!
+  subroutine d_base_mlv_axpby_vv(m,alpha, x, beta, y, j, info)
+    use psi_serial_mod
+    implicit none
+    integer(psb_ipk_), intent(in)                 :: m
+    class(psb_d_base_vect_type), intent(inout)  :: x
+    class(psb_d_base_multivect_type), intent(inout)  :: y
+    real(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(in)      :: j
+    integer(psb_ipk_), intent(out)              :: info
+
+    call psb_geaxpby(m,alpha,x%v,beta,y%v(:,j),info)
+
+  end subroutine d_base_mlv_axpby_vv
+
   !
   ! AXPBY is invoked via Y, hence the structure below.
   !
diff --git a/base/modules/serial/psb_d_vect_mod.F90 b/base/modules/serial/psb_d_vect_mod.F90
index 88fa3262..64a71368 100644
--- a/base/modules/serial/psb_d_vect_mod.F90
+++ b/base/modules/serial/psb_d_vect_mod.F90
@@ -57,7 +57,7 @@ module psb_d_vect_mod
     procedure, pass(x) :: is_remote_build => d_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => d_vect_set_remote_build
     procedure, pass(x) :: get_dupl => d_vect_get_dupl
-    procedure, pass(x) :: set_dupl => d_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => d_vect_set_dupl
     procedure, pass(x) :: get_nrmv => d_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => d_vect_set_nrmv
     procedure, pass(x) :: all      => d_vect_all
@@ -221,7 +221,7 @@ contains
 
     x%nrmv = val
   end subroutine d_vect_set_nrmv
-        
+
 
   function d_vect_is_remote_build(x) result(res)
     implicit none
@@ -241,7 +241,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine d_vect_set_remote_build
-        
+
   subroutine  psb_d_set_vect_default(v)
     implicit none
     class(psb_d_base_vect_type), intent(in) :: v
@@ -599,7 +599,7 @@ contains
       allocate(tmp,stat=info,mold=psb_d_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -1156,7 +1156,7 @@ contains
     ! Temp vectors
     type(psb_d_vect_type) :: wtemp
 
-    info = 0 
+    info = 0
     if( allocated(w%v) ) then
       if (.not.present(aux)) then
         allocate(wtemp%v, mold=w%v)
@@ -1395,6 +1395,8 @@ module psb_d_multivect_mod
 !!$    procedure, pass(x) :: nrm2     => d_vect_nrm2
 !!$    procedure, pass(x) :: amax     => d_vect_amax
 !!$    procedure, pass(x) :: asum     => d_vect_asum
+       procedure, pass(y) :: axpby_vv => d_vect_axpby_vv
+       generic, public :: axpby => axpby_vv
   end type psb_d_multivect_type
 
   public  :: psb_d_multivect, psb_d_multivect_type,&
@@ -1420,7 +1422,7 @@ module psb_d_multivect_mod
 
 contains
 
-  
+
   function d_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_d_multivect_type), intent(in) :: x
@@ -1439,7 +1441,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine d_mvect_set_dupl
-        
+
 
   function d_mvect_is_remote_build(x) result(res)
     implicit none
@@ -1459,7 +1461,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine d_mvect_set_remote_build
-        
+
 
   subroutine  psb_d_set_multivect_default(v)
     implicit none
@@ -1848,6 +1850,24 @@ contains
 !!$
 !!$  end subroutine d_vect_axpby_v
 !!$
+    subroutine d_vect_axpby_vv(m,alpha, x, beta, y, j, info)
+      use psi_serial_mod
+      use psb_d_vect_mod
+      implicit none
+      integer(psb_ipk_), intent(in)            :: m
+      type(psb_d_vect_type), intent(inout)  :: x
+      class(psb_d_multivect_type), intent(inout)  :: y
+      real(psb_dpk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(in)      :: j
+      integer(psb_ipk_), intent(out)           :: info
+
+      if (allocated(x%v).and.allocated(y%v)) then
+         call y%v%axpby(m,alpha,x%v,beta,j,info)
+      else
+         info = psb_err_invalid_vect_state_
+      end if
+
+    end subroutine d_vect_axpby_vv
 !!$  subroutine d_vect_axpby_a(m,alpha, x, beta, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none
diff --git a/base/modules/serial/psb_i_base_vect_mod.F90 b/base/modules/serial/psb_i_base_vect_mod.F90
index 0289ecd0..78f655a4 100644
--- a/base/modules/serial/psb_i_base_vect_mod.F90
+++ b/base/modules/serial/psb_i_base_vect_mod.F90
@@ -203,7 +203,7 @@ contains
     class(psb_i_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -736,7 +736,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -744,7 +744,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function i_base_get_vect
 
   !
@@ -772,7 +772,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -800,7 +800,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
diff --git a/base/modules/serial/psb_i_vect_mod.F90 b/base/modules/serial/psb_i_vect_mod.F90
index ab371bd5..13f65a6c 100644
--- a/base/modules/serial/psb_i_vect_mod.F90
+++ b/base/modules/serial/psb_i_vect_mod.F90
@@ -56,7 +56,7 @@ module psb_i_vect_mod
     procedure, pass(x) :: is_remote_build => i_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => i_vect_set_remote_build
     procedure, pass(x) :: get_dupl => i_vect_get_dupl
-    procedure, pass(x) :: set_dupl => i_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => i_vect_set_dupl
     procedure, pass(x) :: get_nrmv => i_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => i_vect_set_nrmv
     procedure, pass(x) :: all      => i_vect_all
@@ -161,7 +161,7 @@ contains
 
     x%nrmv = val
   end subroutine i_vect_set_nrmv
-        
+
 
   function i_vect_is_remote_build(x) result(res)
     implicit none
@@ -181,7 +181,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine i_vect_set_remote_build
-        
+
   subroutine  psb_i_set_vect_default(v)
     implicit none
     class(psb_i_base_vect_type), intent(in) :: v
@@ -539,7 +539,7 @@ contains
       allocate(tmp,stat=info,mold=psb_i_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -698,7 +698,7 @@ module psb_i_multivect_mod
 
 contains
 
-  
+
   function i_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_i_multivect_type), intent(in) :: x
@@ -717,7 +717,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine i_mvect_set_dupl
-        
+
 
   function i_mvect_is_remote_build(x) result(res)
     implicit none
@@ -737,7 +737,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine i_mvect_set_remote_build
-        
+
 
   subroutine  psb_i_set_multivect_default(v)
     implicit none
diff --git a/base/modules/serial/psb_l_base_vect_mod.F90 b/base/modules/serial/psb_l_base_vect_mod.F90
index d8654f63..5181a25e 100644
--- a/base/modules/serial/psb_l_base_vect_mod.F90
+++ b/base/modules/serial/psb_l_base_vect_mod.F90
@@ -204,7 +204,7 @@ contains
     class(psb_l_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -737,7 +737,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -745,7 +745,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function l_base_get_vect
 
   !
@@ -773,7 +773,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -801,7 +801,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
diff --git a/base/modules/serial/psb_l_vect_mod.F90 b/base/modules/serial/psb_l_vect_mod.F90
index 779d4723..bcf2b938 100644
--- a/base/modules/serial/psb_l_vect_mod.F90
+++ b/base/modules/serial/psb_l_vect_mod.F90
@@ -57,7 +57,7 @@ module psb_l_vect_mod
     procedure, pass(x) :: is_remote_build => l_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => l_vect_set_remote_build
     procedure, pass(x) :: get_dupl => l_vect_get_dupl
-    procedure, pass(x) :: set_dupl => l_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => l_vect_set_dupl
     procedure, pass(x) :: get_nrmv => l_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => l_vect_set_nrmv
     procedure, pass(x) :: all      => l_vect_all
@@ -162,7 +162,7 @@ contains
 
     x%nrmv = val
   end subroutine l_vect_set_nrmv
-        
+
 
   function l_vect_is_remote_build(x) result(res)
     implicit none
@@ -182,7 +182,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine l_vect_set_remote_build
-        
+
   subroutine  psb_l_set_vect_default(v)
     implicit none
     class(psb_l_base_vect_type), intent(in) :: v
@@ -540,7 +540,7 @@ contains
       allocate(tmp,stat=info,mold=psb_l_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -699,7 +699,7 @@ module psb_l_multivect_mod
 
 contains
 
-  
+
   function l_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_l_multivect_type), intent(in) :: x
@@ -718,7 +718,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine l_mvect_set_dupl
-        
+
 
   function l_mvect_is_remote_build(x) result(res)
     implicit none
@@ -738,7 +738,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine l_mvect_set_remote_build
-        
+
 
   subroutine  psb_l_set_multivect_default(v)
     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 4bd6bbfb..4e89e420 100644
--- a/base/modules/serial/psb_s_base_vect_mod.F90
+++ b/base/modules/serial/psb_s_base_vect_mod.F90
@@ -274,7 +274,7 @@ contains
     class(psb_s_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -807,7 +807,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -815,7 +815,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function s_base_get_vect
 
   !
@@ -843,7 +843,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -871,7 +871,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
@@ -925,7 +925,7 @@ contains
     class(psb_s_base_vect_type), intent(inout)  :: x
 
     integer(psb_ipk_) :: i
-    
+
     if (allocated(x%v)) then
       if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
@@ -1177,7 +1177,7 @@ contains
     info = 0
     if (y%is_dev()) call y%sync()
     n = min(size(y%v), size(x))
-    !$omp parallel do private(i)    
+    !$omp parallel do private(i)
     do i=1, n
       y%v(i) = y%v(i)*x(i)
     end do
@@ -1223,7 +1223,7 @@ contains
     else
       if (alpha == sone) then
         if (beta == szero) then
-          !$omp parallel do private(i)    
+          !$omp parallel do private(i)
           do i=1, n
             z%v(i) = y(i)*x(i)
           end do
@@ -1688,7 +1688,7 @@ contains
     end if
     call x%set_host()
   end subroutine s_base_scal
-  
+
   !
   ! Norms 1, 2 and infinity
   !
@@ -1824,7 +1824,7 @@ contains
     integer(psb_ipk_), intent(in)           :: n
     real(psb_spk_)                :: res
     integer(psb_ipk_) :: i
-    
+
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     res=szero
@@ -2143,7 +2143,7 @@ contains
     z%v = x + b
 #endif
     info = 0
-    
+
   end subroutine s_base_addconst_a2
   !
   !> Function  _base_addconst_v2
@@ -2249,8 +2249,9 @@ module psb_s_base_multivect_mod
     procedure, pass(x) :: dot_a    => s_base_mlv_dot_a
     generic, public    :: dot      => dot_v, dot_a
     procedure, pass(y) :: axpby_v  => s_base_mlv_axpby_v
+    procedure, pass(y) :: axpby_vv => s_base_mlv_axpby_vv
     procedure, pass(y) :: axpby_a  => s_base_mlv_axpby_a
-    generic, public    :: axpby    => axpby_v, axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a, axpby_vv
     !
     ! MultiVector by vector/multivector multiplication. Need all variants
     ! to handle multiple requirements from preconditioners
@@ -2907,6 +2908,35 @@ contains
 
   end subroutine s_base_mlv_axpby_v
 
+  !
+  ! AXPBY is invoked via Y, hence the structure below.
+  !
+  !
+  !
+  !> Function  base_mlv_axpby_v
+  !! \memberof  psb_s_base_multivect_type
+  !! \brief AXPBY  by a (base_mlv_vect) y(j)=alpha*x+beta*y(j)
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_v_vect) to be added
+  !! \param beta scalar alpha
+  !! \param j to what column of the multivector should we add
+  !! \param info   return code
+  !!
+  subroutine s_base_mlv_axpby_vv(m,alpha, x, beta, y, j, 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_multivect_type), intent(inout)  :: y
+    real(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(in)      :: j
+    integer(psb_ipk_), intent(out)              :: info
+
+    call psb_geaxpby(m,alpha,x%v,beta,y%v(:,j),info)
+
+  end subroutine s_base_mlv_axpby_vv
+
   !
   ! AXPBY is invoked via Y, hence the structure below.
   !
diff --git a/base/modules/serial/psb_s_vect_mod.F90 b/base/modules/serial/psb_s_vect_mod.F90
index 7a54ecf0..e69554e3 100644
--- a/base/modules/serial/psb_s_vect_mod.F90
+++ b/base/modules/serial/psb_s_vect_mod.F90
@@ -57,7 +57,7 @@ module psb_s_vect_mod
     procedure, pass(x) :: is_remote_build => s_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => s_vect_set_remote_build
     procedure, pass(x) :: get_dupl => s_vect_get_dupl
-    procedure, pass(x) :: set_dupl => s_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => s_vect_set_dupl
     procedure, pass(x) :: get_nrmv => s_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => s_vect_set_nrmv
     procedure, pass(x) :: all      => s_vect_all
@@ -221,7 +221,7 @@ contains
 
     x%nrmv = val
   end subroutine s_vect_set_nrmv
-        
+
 
   function s_vect_is_remote_build(x) result(res)
     implicit none
@@ -241,7 +241,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine s_vect_set_remote_build
-        
+
   subroutine  psb_s_set_vect_default(v)
     implicit none
     class(psb_s_base_vect_type), intent(in) :: v
@@ -599,7 +599,7 @@ contains
       allocate(tmp,stat=info,mold=psb_s_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -1156,7 +1156,7 @@ contains
     ! Temp vectors
     type(psb_s_vect_type) :: wtemp
 
-    info = 0 
+    info = 0
     if( allocated(w%v) ) then
       if (.not.present(aux)) then
         allocate(wtemp%v, mold=w%v)
@@ -1395,6 +1395,8 @@ module psb_s_multivect_mod
 !!$    procedure, pass(x) :: nrm2     => s_vect_nrm2
 !!$    procedure, pass(x) :: amax     => s_vect_amax
 !!$    procedure, pass(x) :: asum     => s_vect_asum
+       procedure, pass(y) :: axpby_vv => s_vect_axpby_vv
+       generic, public :: axpby => axpby_vv
   end type psb_s_multivect_type
 
   public  :: psb_s_multivect, psb_s_multivect_type,&
@@ -1420,7 +1422,7 @@ module psb_s_multivect_mod
 
 contains
 
-  
+
   function s_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_s_multivect_type), intent(in) :: x
@@ -1439,7 +1441,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine s_mvect_set_dupl
-        
+
 
   function s_mvect_is_remote_build(x) result(res)
     implicit none
@@ -1459,7 +1461,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine s_mvect_set_remote_build
-        
+
 
   subroutine  psb_s_set_multivect_default(v)
     implicit none
@@ -1848,6 +1850,24 @@ contains
 !!$
 !!$  end subroutine s_vect_axpby_v
 !!$
+    subroutine s_vect_axpby_vv(m,alpha, x, beta, y, j, info)
+      use psi_serial_mod
+      use psb_s_vect_mod
+      implicit none
+      integer(psb_ipk_), intent(in)            :: m
+      type(psb_s_vect_type), intent(inout)  :: x
+      class(psb_s_multivect_type), intent(inout)  :: y
+      real(psb_spk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(in)      :: j
+      integer(psb_ipk_), intent(out)           :: info
+
+      if (allocated(x%v).and.allocated(y%v)) then
+         call y%v%axpby(m,alpha,x%v,beta,j,info)
+      else
+         info = psb_err_invalid_vect_state_
+      end if
+
+    end subroutine s_vect_axpby_vv
 !!$  subroutine s_vect_axpby_a(m,alpha, x, beta, y, info)
 !!$    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 c52dcd59..e10c773b 100644
--- a/base/modules/serial/psb_z_base_vect_mod.F90
+++ b/base/modules/serial/psb_z_base_vect_mod.F90
@@ -267,7 +267,7 @@ contains
     class(psb_z_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
     integer(psb_ipk_) :: i
-    
+
     call psb_realloc(size(this),x%v,info)
     if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
@@ -800,7 +800,7 @@ contains
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
-    if (.false.) then 
+    if (.false.) then
       res(1:isz) = x%v(1:isz)
     else
       !$omp parallel do private(i)
@@ -808,7 +808,7 @@ contains
         res(i) = x%v(i)
       end do
     end if
-    
+
   end function z_base_get_vect
 
   !
@@ -836,7 +836,7 @@ contains
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     !$omp parallel do private(i)
-    do i = first_, last_        
+    do i = first_, last_
       x%v(i) = val
     end do
 #else
@@ -864,7 +864,7 @@ contains
     if (.not.allocated(x%v)) then
       call psb_realloc(size(val),x%v,info)
     end if
-    
+
     first_                     = 1
     if (present(first)) first_ = max(1,first)
     last_                      = min(psb_size(x%v),first_+size(val)-1)
@@ -918,7 +918,7 @@ contains
     class(psb_z_base_vect_type), intent(inout)  :: x
 
     integer(psb_ipk_) :: i
-    
+
     if (allocated(x%v)) then
       if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
@@ -1170,7 +1170,7 @@ contains
     info = 0
     if (y%is_dev()) call y%sync()
     n = min(size(y%v), size(x))
-    !$omp parallel do private(i)    
+    !$omp parallel do private(i)
     do i=1, n
       y%v(i) = y%v(i)*x(i)
     end do
@@ -1216,7 +1216,7 @@ contains
     else
       if (alpha == zone) then
         if (beta == zzero) then
-          !$omp parallel do private(i)    
+          !$omp parallel do private(i)
           do i=1, n
             z%v(i) = y(i)*x(i)
           end do
@@ -1681,7 +1681,7 @@ contains
     end if
     call x%set_host()
   end subroutine z_base_scal
-  
+
   !
   ! Norms 1, 2 and infinity
   !
@@ -1737,7 +1737,7 @@ contains
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
     integer(psb_ipk_) :: i
-    
+
     if (x%is_dev()) call x%sync()
 #if defined(OPENMP)
     res=dzero
@@ -1964,7 +1964,7 @@ contains
     z%v = x + b
 #endif
     info = 0
-    
+
   end subroutine z_base_addconst_a2
   !
   !> Function  _base_addconst_v2
@@ -2070,8 +2070,9 @@ module psb_z_base_multivect_mod
     procedure, pass(x) :: dot_a    => z_base_mlv_dot_a
     generic, public    :: dot      => dot_v, dot_a
     procedure, pass(y) :: axpby_v  => z_base_mlv_axpby_v
+    procedure, pass(y) :: axpby_vv => z_base_mlv_axpby_vv
     procedure, pass(y) :: axpby_a  => z_base_mlv_axpby_a
-    generic, public    :: axpby    => axpby_v, axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a, axpby_vv
     !
     ! MultiVector by vector/multivector multiplication. Need all variants
     ! to handle multiple requirements from preconditioners
@@ -2728,6 +2729,35 @@ contains
 
   end subroutine z_base_mlv_axpby_v
 
+  !
+  ! AXPBY is invoked via Y, hence the structure below.
+  !
+  !
+  !
+  !> Function  base_mlv_axpby_v
+  !! \memberof  psb_z_base_multivect_type
+  !! \brief AXPBY  by a (base_mlv_vect) y(j)=alpha*x+beta*y(j)
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_v_vect) to be added
+  !! \param beta scalar alpha
+  !! \param j to what column of the multivector should we add
+  !! \param info   return code
+  !!
+  subroutine z_base_mlv_axpby_vv(m,alpha, x, beta, y, j, info)
+    use psi_serial_mod
+    implicit none
+    integer(psb_ipk_), intent(in)                 :: m
+    class(psb_z_base_vect_type), intent(inout)  :: x
+    class(psb_z_base_multivect_type), intent(inout)  :: y
+    complex(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(in)      :: j
+    integer(psb_ipk_), intent(out)              :: info
+
+    call psb_geaxpby(m,alpha,x%v,beta,y%v(:,j),info)
+
+  end subroutine z_base_mlv_axpby_vv
+
   !
   ! AXPBY is invoked via Y, hence the structure below.
   !
diff --git a/base/modules/serial/psb_z_vect_mod.F90 b/base/modules/serial/psb_z_vect_mod.F90
index e8a34859..87c89e90 100644
--- a/base/modules/serial/psb_z_vect_mod.F90
+++ b/base/modules/serial/psb_z_vect_mod.F90
@@ -57,7 +57,7 @@ module psb_z_vect_mod
     procedure, pass(x) :: is_remote_build => z_vect_is_remote_build
     procedure, pass(x) :: set_remote_build => z_vect_set_remote_build
     procedure, pass(x) :: get_dupl => z_vect_get_dupl
-    procedure, pass(x) :: set_dupl => z_vect_set_dupl 
+    procedure, pass(x) :: set_dupl => z_vect_set_dupl
     procedure, pass(x) :: get_nrmv => z_vect_get_nrmv
     procedure, pass(x) :: set_nrmv => z_vect_set_nrmv
     procedure, pass(x) :: all      => z_vect_all
@@ -214,7 +214,7 @@ contains
 
     x%nrmv = val
   end subroutine z_vect_set_nrmv
-        
+
 
   function z_vect_is_remote_build(x) result(res)
     implicit none
@@ -234,7 +234,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine z_vect_set_remote_build
-        
+
   subroutine  psb_z_set_vect_default(v)
     implicit none
     class(psb_z_base_vect_type), intent(in) :: v
@@ -592,7 +592,7 @@ contains
       allocate(tmp,stat=info,mold=psb_z_get_base_vect_default())
     end if
     if (allocated(x%v)) then
-      if (allocated(x%v%v)) then 
+      if (allocated(x%v%v)) then
         call x%v%sync()
         if (info == psb_success_) call tmp%bld(x%v%v)
         call x%v%free(info)
@@ -1149,7 +1149,7 @@ contains
     ! Temp vectors
     type(psb_z_vect_type) :: wtemp
 
-    info = 0 
+    info = 0
     if( allocated(w%v) ) then
       if (.not.present(aux)) then
         allocate(wtemp%v, mold=w%v)
@@ -1316,6 +1316,8 @@ module psb_z_multivect_mod
 !!$    procedure, pass(x) :: nrm2     => z_vect_nrm2
 !!$    procedure, pass(x) :: amax     => z_vect_amax
 !!$    procedure, pass(x) :: asum     => z_vect_asum
+       procedure, pass(y) :: axpby_vv => z_vect_axpby_vv
+       generic, public :: axpby => axpby_vv
   end type psb_z_multivect_type
 
   public  :: psb_z_multivect, psb_z_multivect_type,&
@@ -1341,7 +1343,7 @@ module psb_z_multivect_mod
 
 contains
 
-  
+
   function z_mvect_get_dupl(x) result(res)
     implicit none
     class(psb_z_multivect_type), intent(in) :: x
@@ -1360,7 +1362,7 @@ contains
       x%dupl = psb_dupl_def_
     end if
   end subroutine z_mvect_set_dupl
-        
+
 
   function z_mvect_is_remote_build(x) result(res)
     implicit none
@@ -1380,7 +1382,7 @@ contains
       x%remote_build = psb_matbld_remote_
     end if
   end subroutine z_mvect_set_remote_build
-        
+
 
   subroutine  psb_z_set_multivect_default(v)
     implicit none
@@ -1769,6 +1771,24 @@ contains
 !!$
 !!$  end subroutine z_vect_axpby_v
 !!$
+    subroutine z_vect_axpby_vv(m,alpha, x, beta, y, j, info)
+      use psi_serial_mod
+      use psb_z_vect_mod
+      implicit none
+      integer(psb_ipk_), intent(in)            :: m
+      type(psb_z_vect_type), intent(inout)  :: x
+      class(psb_z_multivect_type), intent(inout)  :: y
+      complex(psb_dpk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(in)      :: j
+      integer(psb_ipk_), intent(out)           :: info
+
+      if (allocated(x%v).and.allocated(y%v)) then
+         call y%v%axpby(m,alpha,x%v,beta,j,info)
+      else
+         info = psb_err_invalid_vect_state_
+      end if
+
+    end subroutine z_vect_axpby_vv
 !!$  subroutine z_vect_axpby_a(m,alpha, x, beta, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none
diff --git a/base/psblas/psb_caxpby.f90 b/base/psblas/psb_caxpby.f90
index da3dd93b..d76bef74 100644
--- a/base/psblas/psb_caxpby.f90
+++ b/base/psblas/psb_caxpby.f90
@@ -741,3 +741,138 @@ subroutine psb_caddconst_vect(x,b,z,desc_a,info)
   return
 
 end subroutine psb_caddconst_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_caxpby.f90
+
+!
+! Subroutine: psb_caxpby_vect
+!    Adds one distributed vector to another,
+!
+!    Y := beta * Y + alpha * X
+!
+! Arguments:
+!    alpha  -  complex,input        The scalar used to multiply each component of X
+!    x      - type(psb_c_vect_type) The input vector containing the entries of X
+!    beta   -  complex,input        The scalar used to multiply each component of Y
+!    y      - type(psb_c_multivect_type)  The input/output vector Y
+!    j      - integer(psb_ipk) The column component of the multivect to which we should add
+!    desc_a -  type(psb_desc_type)  The communication descriptor.
+!    info   -  integer              Return code
+!
+!  Note: from a functional point of view, X is input, but here
+!        it's declared INOUT because of the sync() methods.
+!
+subroutine psb_caxpby_multivect_vect(alpha, x, beta, y,&
+     & j, desc_a, info)
+  use psb_base_mod, psb_protect_name => psb_caxpby_multivect_vect
+  implicit none
+  type(psb_c_vect_type), intent (inout) ::  x
+  type(psb_c_multivect_type), intent (inout) ::  y
+  complex(psb_spk_), intent (in)        :: alpha, beta
+  integer(psb_ipk_), intent(in)        :: j
+  type(psb_desc_type), intent (in)      :: desc_a
+  integer(psb_ipk_), intent(out)                  :: info
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me,&
+       & err_act, iix, jjx, iiy, jjy
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m
+  character(len=20)        :: name, ch_err
+
+  name='psb_cgeaxpby'
+  if (psb_errstatus_fatal()) return
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  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
+
+
+  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 ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+  end if
+
+  if(desc_a%get_local_rows() > 0) then
+    call y%axpby(desc_a%get_local_rows(),&
+         & alpha,x,beta,j,info)
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end subroutine psb_caxpby_multivect_vect
diff --git a/base/psblas/psb_daxpby.f90 b/base/psblas/psb_daxpby.f90
index c386f8f2..d7573e32 100644
--- a/base/psblas/psb_daxpby.f90
+++ b/base/psblas/psb_daxpby.f90
@@ -741,3 +741,138 @@ subroutine psb_daddconst_vect(x,b,z,desc_a,info)
   return
 
 end subroutine psb_daddconst_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_daxpby.f90
+
+!
+! Subroutine: psb_daxpby_vect
+!    Adds one distributed vector to another,
+!
+!    Y := beta * Y + alpha * X
+!
+! Arguments:
+!    alpha  -  real,input        The scalar used to multiply each component of X
+!    x      - type(psb_d_vect_type) The input vector containing the entries of X
+!    beta   -  real,input        The scalar used to multiply each component of Y
+!    y      - type(psb_d_multivect_type)  The input/output vector Y
+!    j      - integer(psb_ipk) The column component of the multivect to which we should add
+!    desc_a -  type(psb_desc_type)  The communication descriptor.
+!    info   -  integer              Return code
+!
+!  Note: from a functional point of view, X is input, but here
+!        it's declared INOUT because of the sync() methods.
+!
+subroutine psb_daxpby_multivect_vect(alpha, x, beta, y,&
+     & j, desc_a, info)
+  use psb_base_mod, psb_protect_name => psb_daxpby_multivect_vect
+  implicit none
+  type(psb_d_vect_type), intent (inout) ::  x
+  type(psb_d_multivect_type), intent (inout) ::  y
+  real(psb_dpk_), intent (in)        :: alpha, beta
+  integer(psb_ipk_), intent(in)        :: j
+  type(psb_desc_type), intent (in)      :: desc_a
+  integer(psb_ipk_), intent(out)                  :: info
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me,&
+       & err_act, iix, jjx, iiy, jjy
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m
+  character(len=20)        :: name, ch_err
+
+  name='psb_dgeaxpby'
+  if (psb_errstatus_fatal()) return
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  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
+
+
+  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 ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+  end if
+
+  if(desc_a%get_local_rows() > 0) then
+    call y%axpby(desc_a%get_local_rows(),&
+         & alpha,x,beta,j,info)
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end subroutine psb_daxpby_multivect_vect
diff --git a/base/psblas/psb_saxpby.f90 b/base/psblas/psb_saxpby.f90
index 78f4d01a..e50528ca 100644
--- a/base/psblas/psb_saxpby.f90
+++ b/base/psblas/psb_saxpby.f90
@@ -741,3 +741,138 @@ subroutine psb_saddconst_vect(x,b,z,desc_a,info)
   return
 
 end subroutine psb_saddconst_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_saxpby.f90
+
+!
+! Subroutine: psb_saxpby_vect
+!    Adds one distributed vector to another,
+!
+!    Y := beta * Y + alpha * X
+!
+! Arguments:
+!    alpha  -  real,input        The scalar used to multiply each component of X
+!    x      - type(psb_s_vect_type) The input vector containing the entries of X
+!    beta   -  real,input        The scalar used to multiply each component of Y
+!    y      - type(psb_s_multivect_type)  The input/output vector Y
+!    j      - integer(psb_ipk) The column component of the multivect to which we should add
+!    desc_a -  type(psb_desc_type)  The communication descriptor.
+!    info   -  integer              Return code
+!
+!  Note: from a functional point of view, X is input, but here
+!        it's declared INOUT because of the sync() methods.
+!
+subroutine psb_saxpby_multivect_vect(alpha, x, beta, y,&
+     & j, desc_a, info)
+  use psb_base_mod, psb_protect_name => psb_saxpby_multivect_vect
+  implicit none
+  type(psb_s_vect_type), intent (inout) ::  x
+  type(psb_s_multivect_type), intent (inout) ::  y
+  real(psb_spk_), intent (in)        :: alpha, beta
+  integer(psb_ipk_), intent(in)        :: j
+  type(psb_desc_type), intent (in)      :: desc_a
+  integer(psb_ipk_), intent(out)                  :: info
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me,&
+       & err_act, iix, jjx, iiy, jjy
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m
+  character(len=20)        :: name, ch_err
+
+  name='psb_sgeaxpby'
+  if (psb_errstatus_fatal()) return
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  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
+
+
+  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 ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+  end if
+
+  if(desc_a%get_local_rows() > 0) then
+    call y%axpby(desc_a%get_local_rows(),&
+         & alpha,x,beta,j,info)
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end subroutine psb_saxpby_multivect_vect
diff --git a/base/psblas/psb_zaxpby.f90 b/base/psblas/psb_zaxpby.f90
index 2258f38f..444ce0ec 100644
--- a/base/psblas/psb_zaxpby.f90
+++ b/base/psblas/psb_zaxpby.f90
@@ -741,3 +741,138 @@ subroutine psb_zaddconst_vect(x,b,z,desc_a,info)
   return
 
 end subroutine psb_zaddconst_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_zaxpby.f90
+
+!
+! Subroutine: psb_zaxpby_vect
+!    Adds one distributed vector to another,
+!
+!    Y := beta * Y + alpha * X
+!
+! Arguments:
+!    alpha  -  complex,input        The scalar used to multiply each component of X
+!    x      - type(psb_z_vect_type) The input vector containing the entries of X
+!    beta   -  complex,input        The scalar used to multiply each component of Y
+!    y      - type(psb_z_multivect_type)  The input/output vector Y
+!    j      - integer(psb_ipk) The column component of the multivect to which we should add
+!    desc_a -  type(psb_desc_type)  The communication descriptor.
+!    info   -  integer              Return code
+!
+!  Note: from a functional point of view, X is input, but here
+!        it's declared INOUT because of the sync() methods.
+!
+subroutine psb_zaxpby_multivect_vect(alpha, x, beta, y,&
+     & j, desc_a, info)
+  use psb_base_mod, psb_protect_name => psb_zaxpby_multivect_vect
+  implicit none
+  type(psb_z_vect_type), intent (inout) ::  x
+  type(psb_z_multivect_type), intent (inout) ::  y
+  complex(psb_dpk_), intent (in)        :: alpha, beta
+  integer(psb_ipk_), intent(in)        :: j
+  type(psb_desc_type), intent (in)      :: desc_a
+  integer(psb_ipk_), intent(out)                  :: info
+
+  ! locals
+  type(psb_ctxt_type) :: ctxt
+  integer(psb_ipk_) :: np, me,&
+       & err_act, iix, jjx, iiy, jjy
+  integer(psb_lpk_) :: ix, ijx, iy, ijy, m
+  character(len=20)        :: name, ch_err
+
+  name='psb_zgeaxpby'
+  if (psb_errstatus_fatal()) return
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  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
+
+
+  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 ((iix /= ione).or.(iiy /= ione)) then
+    info=psb_err_ix_n1_iy_n1_unsupported_
+    call psb_errpush(info,name)
+  end if
+
+  if(desc_a%get_local_rows() > 0) then
+    call y%axpby(desc_a%get_local_rows(),&
+         & alpha,x,beta,j,info)
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(ctxt,err_act)
+
+  return
+
+end subroutine psb_zaxpby_multivect_vect