diff --git a/base/modules/psb_c_base_vect_mod.f90 b/base/modules/psb_c_base_vect_mod.f90
index 58735b76..afdcf44e 100644
--- a/base/modules/psb_c_base_vect_mod.f90
+++ b/base/modules/psb_c_base_vect_mod.f90
@@ -738,7 +738,7 @@ contains
   !> Function  base_dot_v
   !! \memberof  psb_c_base_vect_type
   !! \brief  Dot product by another base_vector
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y    The other (base_vect) to be multiplied by
   !!
   function c_base_dot_v(n,x,y) result(res)
@@ -772,7 +772,7 @@ contains
   !> Function  base_dot_a
   !! \memberof  psb_c_base_vect_type
   !! \brief  Dot product by a normal array
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y(:) The array to be multiplied by
   !!
   function c_base_dot_a(n,x,y) result(res)
@@ -795,7 +795,7 @@ contains
   !> Function  base_axpby_v
   !! \memberof  psb_c_base_vect_type
   !! \brief AXPBY  by a (base_vect) y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x     The class(base_vect) to be added
   !! \param beta scalar alpha
@@ -823,7 +823,7 @@ contains
   !> Function  base_axpby_a
   !! \memberof  psb_c_base_vect_type
   !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x(:) The array to be added
   !! \param beta scalar alpha
@@ -1339,7 +1339,7 @@ module psb_c_base_multivect_mod
   
   use psb_const_mod
   use psb_error_mod
-
+  use psb_realloc_mod
 
   !> \namespace  psb_base_mod  \class psb_c_base_vect_type
   !! The psb_c_base_vect_type 
@@ -1361,20 +1361,20 @@ module psb_c_base_multivect_mod
     !
     !  Constructors/allocators
     !
-    procedure, pass(x) :: bld_x    => c_base_mv_bld_x
-    procedure, pass(x) :: bld_n    => c_base_mv_bld_n
+    procedure, pass(x) :: bld_x    => c_base_mlv_bld_x
+    procedure, pass(x) :: bld_n    => c_base_mlv_bld_n
     generic, public    :: bld      => bld_x, bld_n
-    procedure, pass(x) :: all      => c_base_mv_all
-    procedure, pass(x) :: mold     => c_base_mv_mold
+    procedure, pass(x) :: all      => c_base_mlv_all
+    procedure, pass(x) :: mold     => c_base_mlv_mold
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
     ! in derived classes. 
     !
-    procedure, pass(x) :: ins      => c_base_mv_ins
-    procedure, pass(x) :: zero     => c_base_mv_zero
-    procedure, pass(x) :: asb      => c_base_mv_asb
-    procedure, pass(x) :: free     => c_base_mv_free
+    procedure, pass(x) :: ins      => c_base_mlv_ins
+    procedure, pass(x) :: zero     => c_base_mlv_zero
+    procedure, pass(x) :: asb      => c_base_mlv_asb
+    procedure, pass(x) :: free     => c_base_mlv_free
     !
     ! Sync: centerpiece of handling of external storage.
     ! Any derived class having extra storage upon sync
@@ -1382,66 +1382,66 @@ module psb_c_base_multivect_mod
     ! external side contain the same data. The base
     ! version is only a placeholder. 
     !
-    procedure, pass(x) :: sync     => c_base_mv_sync
-    procedure, pass(x) :: is_host  => c_base_mv_is_host
-    procedure, pass(x) :: is_dev   => c_base_mv_is_dev
-    procedure, pass(x) :: is_sync  => c_base_mv_is_sync
-    procedure, pass(x) :: set_host => c_base_mv_set_host
-    procedure, pass(x) :: set_dev  => c_base_mv_set_dev
-    procedure, pass(x) :: set_sync => c_base_mv_set_sync
+    procedure, pass(x) :: sync     => c_base_mlv_sync
+    procedure, pass(x) :: is_host  => c_base_mlv_is_host
+    procedure, pass(x) :: is_dev   => c_base_mlv_is_dev
+    procedure, pass(x) :: is_sync  => c_base_mlv_is_sync
+    procedure, pass(x) :: set_host => c_base_mlv_set_host
+    procedure, pass(x) :: set_dev  => c_base_mlv_set_dev
+    procedure, pass(x) :: set_sync => c_base_mlv_set_sync
 
     !
     ! Basic info
-    procedure, pass(x) :: get_nrows => c_base_mv_get_nrows
-    procedure, pass(x) :: get_ncols => c_base_mv_get_ncols
-    procedure, pass(x) :: sizeof    => c_base_mv_sizeof
-    procedure, nopass  :: get_fmt   => c_base_mv_get_fmt
+    procedure, pass(x) :: get_nrows => c_base_mlv_get_nrows
+    procedure, pass(x) :: get_ncols => c_base_mlv_get_ncols
+    procedure, pass(x) :: sizeof    => c_base_mlv_sizeof
+    procedure, nopass  :: get_fmt   => c_base_mlv_get_fmt
     !
     ! Set/get data from/to an external array; also
     ! overload assignment.
     !
-    procedure, pass(x) :: get_vect => c_base_mv_get_vect
-    procedure, pass(x) :: set_scal => c_base_mv_set_scal
-    procedure, pass(x) :: set_vect => c_base_mv_set_vect
+    procedure, pass(x) :: get_vect => c_base_mlv_get_vect
+    procedure, pass(x) :: set_scal => c_base_mlv_set_scal
+    procedure, pass(x) :: set_vect => c_base_mlv_set_vect
     generic, public    :: set      => set_vect, set_scal
 
     !
     ! Dot product and AXPBY
     !
-!!$    procedure, pass(x) :: dot_v    => c_base_mv_dot_v
-!!$    procedure, pass(x) :: dot_a    => c_base_mv_dot_a
-!!$    generic, public    :: dot      => dot_v, dot_a
-!!$    procedure, pass(y) :: axpby_v  => c_base_mv_axpby_v
-!!$    procedure, pass(y) :: axpby_a  => c_base_mv_axpby_a
-!!$    generic, public    :: axpby    => axpby_v, axpby_a
+    procedure, pass(x) :: dot_v    => c_base_mlv_dot_v
+    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_a  => c_base_mlv_axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a
 !!$    !
 !!$    ! Vector by vector multiplication. Need all variants
 !!$    ! to handle multiple requirements from preconditioners
 !!$    !
-!!$    procedure, pass(y) :: mlt_v    => c_base_mv_mlt_v
-!!$    procedure, pass(y) :: mlt_a    => c_base_mv_mlt_a
-!!$    procedure, pass(z) :: mlt_a_2  => c_base_mv_mlt_a_2
-!!$    procedure, pass(z) :: mlt_v_2  => c_base_mv_mlt_v_2
-!!$    procedure, pass(z) :: mlt_va   => c_base_mv_mlt_va
-!!$    procedure, pass(z) :: mlt_av   => c_base_mv_mlt_av
+!!$    procedure, pass(y) :: mlt_v    => c_base_mlv_mlt_v
+!!$    procedure, pass(y) :: mlt_a    => c_base_mlv_mlt_a
+!!$    procedure, pass(z) :: mlt_a_2  => c_base_mlv_mlt_a_2
+!!$    procedure, pass(z) :: mlt_v_2  => c_base_mlv_mlt_v_2
+!!$    procedure, pass(z) :: mlt_va   => c_base_mlv_mlt_va
+!!$    procedure, pass(z) :: mlt_av   => c_base_mlv_mlt_av
 !!$    generic, public    :: mlt      => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
 !!$    !
 !!$    ! Scaling and norms
 !!$    !
-!!$    procedure, pass(x) :: scal     => c_base_mv_scal
-!!$    procedure, pass(x) :: nrm2     => c_base_mv_nrm2
-!!$    procedure, pass(x) :: amax     => c_base_mv_amax
-!!$    procedure, pass(x) :: asum     => c_base_mv_asum
+!!$    procedure, pass(x) :: scal     => c_base_mlv_scal
+!!$    procedure, pass(x) :: nrm2     => c_base_mlv_nrm2
+!!$    procedure, pass(x) :: amax     => c_base_mlv_amax
+!!$    procedure, pass(x) :: asum     => c_base_mlv_asum
 !!$    !
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => c_base_mv_gthab
-!!$    procedure, pass(x) :: gthzv    => c_base_mv_gthzv
-!!$    procedure, pass(x) :: gthzv_x  => c_base_mv_gthzv_x
+!!$    procedure, pass(x) :: gthab    => c_base_mlv_gthab
+!!$    procedure, pass(x) :: gthzv    => c_base_mlv_gthzv
+!!$    procedure, pass(x) :: gthzv_x  => c_base_mlv_gthzv_x
 !!$    generic, public    :: gth      => gthab, gthzv, gthzv_x
-!!$    procedure, pass(y) :: sctb     => c_base_mv_sctb
-!!$    procedure, pass(y) :: sctb_x   => c_base_mv_sctb_x
+!!$    procedure, pass(y) :: sctb     => c_base_mlv_sctb
+!!$    procedure, pass(y) :: sctb_x   => c_base_mlv_sctb_x
 !!$    generic, public    :: sct      => sctb, sctb_x
   end type psb_c_base_multivect_type
 
@@ -1491,7 +1491,7 @@ contains
   !! \brief     Build method from an array
   !!  \param   x(:)  input array to be copied
   !!
-  subroutine c_base_mv_bld_x(x,this)
+  subroutine c_base_mlv_bld_x(x,this)
     use psb_realloc_mod
     complex(psb_spk_), intent(in) :: this(:,:)
     class(psb_c_base_multivect_type), intent(inout) :: x
@@ -1499,12 +1499,12 @@ contains
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
     x%v(:,:)  = this(:,:)
 
-  end subroutine c_base_mv_bld_x
+  end subroutine c_base_mlv_bld_x
     
   !
   ! Create with size, but no initialization
@@ -1515,7 +1515,7 @@ contains
   !! \brief     Build method with size (uninitialized data)
   !!  \param    n    size to be allocated. 
   !!
-  subroutine c_base_mv_bld_n(x,m,n)
+  subroutine c_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: m,n
     class(psb_c_base_multivect_type), intent(inout) :: x
@@ -1524,16 +1524,16 @@ contains
     call psb_realloc(m,n,x%v,info)
     call x%asb(m,n,info)
 
-  end subroutine c_base_mv_bld_n
+  end subroutine c_base_mlv_bld_n
   
-  !> Function  base_mv_all:
+  !> Function  base_mlv_all:
   !! \memberof  psb_c_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
   !!  \param    n    size to be allocated. 
   !!  \param info  return code
   !!
-  subroutine c_base_mv_all(m,n, x, info)
+  subroutine c_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1543,15 +1543,15 @@ contains
     
     call psb_realloc(m,n,x%v,info)
     
-  end subroutine c_base_mv_all
+  end subroutine c_base_mlv_all
 
-  !> Function  base_mv_mold:
+  !> Function  base_mlv_mold:
   !! \memberof  psb_c_base_multivect_type
   !! \brief     Mold method: return a variable with the same dynamic type
   !!  \param    y returned variable
   !!  \param info  return code
   !!
-  subroutine c_base_mv_mold(x, y, info)
+  subroutine c_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1561,12 +1561,12 @@ contains
     
     allocate(psb_c_base_multivect_type :: y, stat=info)
 
-  end subroutine c_base_mv_mold
+  end subroutine c_base_mlv_mold
 
   !
   ! Insert a bunch of values at specified positions.
   !
-  !> Function  base_mv_ins:
+  !> Function  base_mlv_ins:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Insert coefficients. 
   !!
@@ -1590,7 +1590,7 @@ contains
   !!  \param info  return code
   !!
   !
-  subroutine c_base_mv_ins(n,irl,val,dupl,x,info)
+  subroutine c_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
     implicit none 
     class(psb_c_base_multivect_type), intent(inout)  :: x
@@ -1642,26 +1642,26 @@ contains
       end select
     end if
     if (info /= 0) then 
-      call psb_errpush(info,'base_mv_vect_ins')
+      call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
 
-  end subroutine c_base_mv_ins
+  end subroutine c_base_mlv_ins
 
   !
-  !> Function  base_mv_zero
+  !> Function  base_mlv_zero
   !! \memberof  psb_c_base_multivect_type
   !! \brief Zero out contents
   !!
   !
-  subroutine c_base_mv_zero(x)
+  subroutine c_base_mlv_zero(x)
     use psi_serial_mod
     implicit none 
     class(psb_c_base_multivect_type), intent(inout)    :: x
     
     if (allocated(x%v)) x%v=czero
 
-  end subroutine c_base_mv_zero
+  end subroutine c_base_mlv_zero
 
   
   !
@@ -1669,7 +1669,7 @@ contains
   ! For derived classes: after this the vector
   ! storage is supposed to be in sync.
   !
-  !> Function  base_mv_asb:
+  !> Function  base_mlv_asb:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
   !!           
@@ -1678,7 +1678,7 @@ contains
   !!
   !
  
-  subroutine c_base_mv_asb(m,n, x, info)
+  subroutine c_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1691,18 +1691,18 @@ contains
     if (info /= 0) &
          & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
 
-  end subroutine c_base_mv_asb
+  end subroutine c_base_mlv_asb
 
 
   !
-  !> Function  base_mv_free:
+  !> Function  base_mlv_free:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Free vector
   !!           
   !!  \param info  return code
   !!
   !
-  subroutine c_base_mv_free(x, info)
+  subroutine c_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1714,7 +1714,7 @@ contains
     if (info /= 0) call & 
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
         
-  end subroutine c_base_mv_free
+  end subroutine c_base_mlv_free
 
   
 
@@ -1723,106 +1723,106 @@ contains
   ! a placeholder.
   ! 
   !
-  !> Function  base_mv_sync:
+  !> Function  base_mlv_sync:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Sync: base version is a no-op.
   !!           
   !
-  subroutine c_base_mv_sync(x)
+  subroutine c_base_mlv_sync(x)
     implicit none 
     class(psb_c_base_multivect_type), intent(inout) :: x
     
-  end subroutine c_base_mv_sync
+  end subroutine c_base_mlv_sync
 
   !
-  !> Function  base_mv_set_host:
+  !> Function  base_mlv_set_host:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Set_host: base version is a no-op.
   !!           
   !
-  subroutine c_base_mv_set_host(x)
+  subroutine c_base_mlv_set_host(x)
     implicit none 
     class(psb_c_base_multivect_type), intent(inout) :: x
     
-  end subroutine c_base_mv_set_host
+  end subroutine c_base_mlv_set_host
 
   !
-  !> Function  base_mv_set_dev:
+  !> Function  base_mlv_set_dev:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
   !!           
   !
-  subroutine c_base_mv_set_dev(x)
+  subroutine c_base_mlv_set_dev(x)
     implicit none 
     class(psb_c_base_multivect_type), intent(inout) :: x
     
-  end subroutine c_base_mv_set_dev
+  end subroutine c_base_mlv_set_dev
 
   !
-  !> Function  base_mv_set_sync:
+  !> Function  base_mlv_set_sync:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
   !!           
   !
-  subroutine c_base_mv_set_sync(x)
+  subroutine c_base_mlv_set_sync(x)
     implicit none 
     class(psb_c_base_multivect_type), intent(inout) :: x
     
-  end subroutine c_base_mv_set_sync
+  end subroutine c_base_mlv_set_sync
 
   !
-  !> Function  base_mv_is_dev:
+  !> Function  base_mlv_is_dev:
   !! \memberof  psb_c_base_multivect_type
   !! \brief Is  vector on external device    .
   !!           
   !
-  function c_base_mv_is_dev(x) result(res)
+  function c_base_mlv_is_dev(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     logical  :: res
   
     res = .false.
-  end function c_base_mv_is_dev
+  end function c_base_mlv_is_dev
   
   !
-  !> Function  base_mv_is_host
+  !> Function  base_mlv_is_host
   !! \memberof  psb_c_base_multivect_type
   !! \brief Is  vector on standard memory    .
   !!           
   !
-  function c_base_mv_is_host(x) result(res)
+  function c_base_mlv_is_host(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function c_base_mv_is_host
+  end function c_base_mlv_is_host
 
   !
-  !> Function  base_mv_is_sync
+  !> Function  base_mlv_is_sync
   !! \memberof  psb_c_base_multivect_type
   !! \brief Is  vector on sync               .
   !!           
   !
-  function c_base_mv_is_sync(x) result(res)
+  function c_base_mlv_is_sync(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function c_base_mv_is_sync
+  end function c_base_mlv_is_sync
 
 
   !
   ! Size info. 
   !
   !
-  !> Function  base_mv_get_nrows
+  !> Function  base_mlv_get_nrows
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Number of entries
   !!           
   !
-  function c_base_mv_get_nrows(x) result(res)
+  function c_base_mlv_get_nrows(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1830,9 +1830,9 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,1)
 
-  end function c_base_mv_get_nrows
+  end function c_base_mlv_get_nrows
 
-  function c_base_mv_get_ncols(x) result(res)
+  function c_base_mlv_get_ncols(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1840,15 +1840,15 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,2)
 
-  end function c_base_mv_get_ncols
+  end function c_base_mlv_get_ncols
   
   !
-  !> Function  base_mv_get_sizeof
+  !> Function  base_mlv_get_sizeof
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Size in bytesa
   !!           
   !
-  function c_base_mv_sizeof(x) result(res)
+  function c_base_mlv_sizeof(x) result(res)
     implicit none 
     class(psb_c_base_multivect_type), intent(in) :: x
     integer(psb_long_int_k_) :: res
@@ -1856,30 +1856,30 @@ contains
     ! Force 8-byte integers.
     res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()
 
-  end function c_base_mv_sizeof
+  end function c_base_mlv_sizeof
 
   !
-  !> Function  base_mv_get_fmt
+  !> Function  base_mlv_get_fmt
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Format
   !!           
   !
-  function c_base_mv_get_fmt() result(res)
+  function c_base_mlv_get_fmt() result(res)
     implicit none 
     character(len=5) :: res
     res = 'BASE'
-  end function c_base_mv_get_fmt
+  end function c_base_mlv_get_fmt
   
 
   !
   !
   !
-  !> Function  base_mv_get_vect
+  !> Function  base_mlv_get_vect
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
   !    
-  function  c_base_mv_get_vect(x) result(res)
+  function  c_base_mlv_get_vect(x) result(res)
     class(psb_c_base_multivect_type), intent(inout) :: x
     complex(psb_spk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
@@ -1889,37 +1889,37 @@ contains
     call x%sync()
     allocate(res(m,n),stat=info) 
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
     res(1:m,1:n) = x%v(1:m,1:n)
-  end function c_base_mv_get_vect
+  end function c_base_mlv_get_vect
     
   !
   ! Reset all values 
   !
   !
-  !> Function  base_mv_set_scal
+  !> Function  base_mlv_set_scal
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Set all entries
   !! \param val   The value to set
   !!
-  subroutine c_base_mv_set_scal(x,val)
+  subroutine c_base_mlv_set_scal(x,val)
     class(psb_c_base_multivect_type), intent(inout)  :: x
     complex(psb_spk_), intent(in) :: val
         
     integer(psb_ipk_) :: info
     x%v = val
     
-  end subroutine c_base_mv_set_scal
+  end subroutine c_base_mlv_set_scal
 
   !
-  !> Function  base_mv_set_vect
+  !> Function  base_mlv_set_vect
   !! \memberof  psb_c_base_multivect_type
   !! \brief  Set all entries
   !! \param val(:)  The vector to be copied in 
   !!
-  subroutine c_base_mv_set_vect(x,val)
+  subroutine c_base_mlv_set_vect(x,val)
     class(psb_c_base_multivect_type), intent(inout)  :: x
     complex(psb_spk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr
@@ -1934,121 +1934,148 @@ contains
       x%v = val
     end if
 
-  end subroutine c_base_mv_set_vect
+  end subroutine c_base_mlv_set_vect
 
-!!$  !
-!!$  ! Dot products 
-!!$  ! 
-!!$  !
-!!$  !> Function  base_mv_dot_v
-!!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief  Dot product by another base_mv_vector
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y    The other (base_mv_vect) to be multiplied by
-!!$  !!
-!!$  function c_base_mv_dot_v(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_c_base_multivect_type), intent(inout) :: x, y
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    complex(psb_spk_)                :: res
-!!$    complex(psb_spk_), external      :: ddot
-!!$    
-!!$    res = izero
-!!$    !
-!!$    ! Note: this is the base implementation.
-!!$    !  When we get here, we are sure that X is of
-!!$    !  TYPE psb_c_base_mv_vect.
-!!$    !  If Y is not, throw the burden on it, implicitly
-!!$    !  calling dot_a
-!!$    !
-!!$    select type(yy => y)
-!!$    type is (psb_c_base_multivect_type)
-!!$      res = ddot(n,x%v,1,y%v,1)
-!!$    class default
-!!$      res = y%dot(n,x%v)
-!!$    end select
-!!$
-!!$  end function c_base_mv_dot_v
-!!$
-!!$  !
-!!$  ! Base workhorse is good old BLAS1
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_dot_a
-!!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief  Dot product by a normal array
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y(:) The array to be multiplied by
-!!$  !!
-!!$  function c_base_mv_dot_a(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_c_base_multivect_type), intent(inout) :: x
-!!$    complex(psb_spk_), intent(in)    :: y(:)
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    complex(psb_spk_)                :: res
-!!$    integer(psb_ipk_), external      :: ddot
-!!$    
-!!$    res = ddot(n,y,1,x%v,1)
-!!$
-!!$  end function c_base_mv_dot_a
+  !
+  ! Dot products 
+  ! 
+  !
+  !> Function  base_mlv_dot_v
+  !! \memberof  psb_c_base_multivect_type
+  !! \brief  Dot product by another base_mlv_vector
+  !! \param n    Number of entries to be considered
+  !! \param y    The other (base_mlv_vect) to be multiplied by
+  !!
+  function c_base_mlv_dot_v(n,x,y) result(res)
+    implicit none 
+    class(psb_c_base_multivect_type), intent(inout) :: x, y
+    integer(psb_ipk_), intent(in)           :: n
+    complex(psb_spk_), allocatable   :: res(:)
+    complex(psb_spk_), external      :: cdotc
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    res = czero
+    !
+    ! Note: this is the base implementation.
+    !  When we get here, we are sure that X is of
+    !  TYPE psb_c_base_mlv_vect (or its class does not care). 
+    !  If Y is not, throw the burden on it, implicitly
+    !  calling dot_a
+    !
+    select type(yy => y)
+    type is (psb_c_base_multivect_type)
+      if (y%is_dev()) call y%sync()
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+      allocate(res(nc))
+      do j=1,nc
+        res(j) = cdotc(n,x%v(:,j),1,y%v(:,j),1)
+      end do
+    class default
+      res = y%dot(n,x%v)
+    end select
+
+  end function c_base_mlv_dot_v
+
+  !
+  ! Base workhorse is good old BLAS1
+  !
+  !
+  !> Function  base_mlv_dot_a
+  !! \memberof  psb_c_base_multivect_type
+  !! \brief  Dot product by a normal array
+  !! \param n    Number of entries to be considered
+  !! \param y(:) The array to be multiplied by
+  !!
+  function c_base_mlv_dot_a(n,x,y) result(res)
+    implicit none 
+    class(psb_c_base_multivect_type), intent(inout) :: x
+    complex(psb_spk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    complex(psb_spk_), allocatable     :: res(:)
+    integer(psb_ipk_), external      :: cdotc
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    nc = min(psb_size(x%v,2),size(y,2))
+    allocate(res(nc))
+    do j=1,nc
+      res(j) = cdotc(n,x%v(:,j),1,y(:,j),1)
+    end do
+
+  end function c_base_mlv_dot_a
     
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_v
-!!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief AXPBY  by a (base_mv_vect) y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x     The class(base_mv_vect) to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine c_base_mv_axpby_v(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    class(psb_c_base_multivect_type), intent(inout)  :: x
-!!$    class(psb_c_base_multivect_type), intent(inout)  :: y
-!!$    complex(psb_spk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    select type(xx => x)
-!!$    type is (psb_c_base_multivect_type)
-!!$      call psb_geaxpby(m,alpha,x%v,beta,y%v,info)
-!!$    class default
-!!$      call y%axpby(m,alpha,x%v,beta,info)
-!!$    end select
-!!$
-!!$  end subroutine c_base_mv_axpby_v
-!!$
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_a
-!!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x(:) The array to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine c_base_mv_axpby_a(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    complex(psb_spk_), intent(in)        :: x(:)
-!!$    class(psb_c_base_multivect_type), intent(inout)  :: y
-!!$    complex(psb_spk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    call psb_geaxpby(m,alpha,x,beta,y%v,info)
-!!$    
-!!$  end subroutine c_base_mv_axpby_a
+  !
+  ! 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=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_mlv_vect) to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine c_base_mlv_axpby_v(m,alpha, x, beta, y, info, n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    class(psb_c_base_multivect_type), intent(inout)  :: x
+    class(psb_c_base_multivect_type), intent(inout)  :: y
+    complex(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+    end if
+    select type(xx => x)
+    type is (psb_c_base_multivect_type)
+      call psb_geaxpby(m,nc,alpha,x%v,beta,y%v,info)
+    class default
+      call y%axpby(m,alpha,x%v,beta,info,n=n)
+    end select
+
+  end subroutine c_base_mlv_axpby_v
+
+  !
+  ! AXPBY is invoked via Y, hence the structure below. 
+  !
+  !
+  !> Function  base_mlv_axpby_a
+  !! \memberof  psb_c_base_multivect_type
+  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x(:) The array to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine c_base_mlv_axpby_a(m,alpha, x, beta, y, info,n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    complex(psb_spk_), intent(in)        :: x(:,:)
+    class(psb_c_base_multivect_type), intent(inout)  :: y
+    complex(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(size(x,2),psb_size(y%v,2))
+    end if
+    
+    call psb_geaxpby(m,nc,alpha,x,beta,y%v,info)
+    
+  end subroutine c_base_mlv_axpby_a
 
   
 !!$  !
@@ -2060,13 +2087,13 @@ contains
 !!$  !  of the involved entities
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief Vector entry-by-entry multiply  by a base_mv_vect array y=x*y
-!!$  !! \param x   The class(base_mv_vect) to be multiplied by
+!!$  !! \brief Vector entry-by-entry multiply  by a base_mlv_vect array y=x*y
+!!$  !! \param x   The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine c_base_mv_mlt_v(x, y, info)
+!!$  subroutine c_base_mlv_mlt_v(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_c_base_multivect_type), intent(inout)  :: x
@@ -2085,16 +2112,16 @@ contains
 !!$      call y%mlt(x%v,info)
 !!$    end select
 !!$
-!!$  end subroutine c_base_mv_mlt_v
+!!$  end subroutine c_base_mlv_mlt_v
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief Vector entry-by-entry multiply  by a normal array y=x*y
 !!$  !! \param x(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine c_base_mv_mlt_a(x, y, info)
+!!$  subroutine c_base_mlv_mlt_a(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_spk_), intent(in)        :: x(:)
@@ -2108,11 +2135,11 @@ contains
 !!$      y%v(i) = y%v(i)*x(i)
 !!$    end do
 !!$    
-!!$  end subroutine c_base_mv_mlt_a
+!!$  end subroutine c_base_mlv_mlt_a
 !!$
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a_2
+!!$  !> Function  base_mlv_mlt_a_2
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by normal arrays
 !!$  !!        z=beta*z+alpha*x*y
@@ -2122,7 +2149,7 @@ contains
 !!$  !! \param y(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine c_base_mv_mlt_a_2(alpha,x,y,beta,z,info)
+!!$  subroutine c_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_spk_), intent(in)        :: alpha,beta
@@ -2187,20 +2214,20 @@ contains
 !!$        end if
 !!$      end if
 !!$    end if
-!!$  end subroutine c_base_mv_mlt_a_2
+!!$  end subroutine c_base_mlv_mlt_a_2
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_v_2
+!!$  !> Function  base_mlv_mlt_v_2
 !!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mv_vect)
+!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mlv_vect)
 !!$  !!        z=beta*z+alpha*x*y
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  !! \param x  The class(base_mv_vect) to be multiplied b
-!!$  !! \param y  The class(base_mv_vect) to be multiplied by
+!!$  !! \param x  The class(base_mlv_vect) to be multiplied b
+!!$  !! \param y  The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine c_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
+!!$  subroutine c_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
 !!$    use psi_serial_mod
 !!$    use psb_string_mod
 !!$    implicit none 
@@ -2227,9 +2254,9 @@ contains
 !!$      if (conjgx_) x%v=(x%v)
 !!$      if (conjgy_) y%v=(y%v)
 !!$    end if
-!!$  end subroutine c_base_mv_mlt_v_2
+!!$  end subroutine c_base_mlv_mlt_v_2
 !!$
-!!$  subroutine c_base_mv_mlt_av(alpha,x,y,beta,z,info)
+!!$  subroutine c_base_mlv_mlt_av(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_spk_), intent(in)        :: alpha,beta
@@ -2243,9 +2270,9 @@ contains
 !!$    
 !!$    call z%mlt(alpha,x,y%v,beta,info)
 !!$
-!!$  end subroutine c_base_mv_mlt_av
+!!$  end subroutine c_base_mlv_mlt_av
 !!$
-!!$  subroutine c_base_mv_mlt_va(alpha,x,y,beta,z,info)
+!!$  subroutine c_base_mlv_mlt_va(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_spk_), intent(in)        :: alpha,beta
@@ -2259,18 +2286,18 @@ contains
 !!$    
 !!$    call z%mlt(alpha,y,x,beta,info)
 !!$
-!!$  end subroutine c_base_mv_mlt_va
+!!$  end subroutine c_base_mlv_mlt_va
 !!$
 !!$
 !!$  !
 !!$  ! Simple scaling 
 !!$  !
-!!$  !> Function  base_mv_scal
+!!$  !> Function  base_mlv_scal
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief Scale all entries  x = alpha*x
 !!$  !! \param alpha   The multiplier
 !!$  !!
-!!$  subroutine c_base_mv_scal(alpha, x)
+!!$  subroutine c_base_mlv_scal(alpha, x)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_c_base_multivect_type), intent(inout)  :: x
@@ -2278,16 +2305,16 @@ contains
 !!$    
 !!$    if (allocated(x%v)) x%v = alpha*x%v
 !!$
-!!$  end subroutine c_base_mv_scal
+!!$  end subroutine c_base_mlv_scal
 !!$  
 !!$  !
 !!$  ! Norms 1, 2 and infinity
 !!$  !
-!!$  !> Function  base_mv_nrm2
+!!$  !> Function  base_mlv_nrm2
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief 2-norm |x(1:n)|_2
 !!$  !! \param n  how many entries to consider
-!!$  function c_base_mv_nrm2(n,x) result(res)
+!!$  function c_base_mlv_nrm2(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_c_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2296,14 +2323,14 @@ contains
 !!$    
 !!$    res =  dnrm2(n,x%v,1)
 !!$
-!!$  end function c_base_mv_nrm2
+!!$  end function c_base_mlv_nrm2
 !!$  
 !!$  !
-!!$  !> Function  base_mv_amax
+!!$  !> Function  base_mlv_amax
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief infinity-norm |x(1:n)|_\infty
 !!$  !! \param n  how many entries to consider
-!!$  function c_base_mv_amax(n,x) result(res)
+!!$  function c_base_mlv_amax(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_c_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2311,14 +2338,14 @@ contains
 !!$    
 !!$    res =  maxval(abs(x%v(1:n)))
 !!$
-!!$  end function c_base_mv_amax
+!!$  end function c_base_mlv_amax
 !!$
 !!$  !
-!!$  !> Function  base_mv_asum
+!!$  !> Function  base_mlv_asum
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief 1-norm |x(1:n)|_1
 !!$  !! \param n  how many entries to consider
-!!$  function c_base_mv_asum(n,x) result(res)
+!!$  function c_base_mlv_asum(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_c_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2326,14 +2353,14 @@ contains
 !!$    
 !!$    res =  sum(abs(x%v(1:n)))
 !!$
-!!$  end function c_base_mv_asum
+!!$  end function c_base_mlv_asum
 !!$  
 !!$  
 !!$  !
 !!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_gthab
+!!$  !> Function  base_mlv_gthab
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief gather into an array
 !!$  !!    Y = beta * Y + alpha * X(IDX(:))
@@ -2341,7 +2368,7 @@ contains
 !!$  !! \param idx(:) indices
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  subroutine c_base_mv_gthab(n,idx,alpha,x,beta,y)
+!!$  subroutine c_base_mlv_gthab(n,idx,alpha,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_spk_) :: alpha, beta, y(:)
@@ -2350,17 +2377,17 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,alpha,x%v,beta,y)
 !!$
-!!$  end subroutine c_base_mv_gthab
+!!$  end subroutine c_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine c_base_mv_gthzv_x(i,n,idx,x,y)
+!!$  subroutine c_base_mlv_gthzv_x(i,n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i,n
 !!$    class(psb_c_base_multivect_type) :: idx
@@ -2369,18 +2396,18 @@ contains
 !!$    
 !!$    call x%gth(n,idx%v(i:),y)
 !!$
-!!$  end subroutine c_base_mv_gthzv_x
+!!$  end subroutine c_base_mlv_gthzv_x
 !!$
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_c_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine c_base_mv_gthzv(n,idx,x,y)
+!!$  subroutine c_base_mlv_gthzv(n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_spk_) ::  y(:)
@@ -2389,22 +2416,22 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,x%v,y)
 !!$
-!!$  end subroutine c_base_mv_gthzv
+!!$  end subroutine c_base_mlv_gthzv
 !!$
 !!$  !
 !!$  ! Scatter: 
 !!$  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
 !!$  ! 
 !!$  !
-!!$  !> Function  base_mv_sctb
+!!$  !> Function  base_mlv_sctb
 !!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief scatter into a class(base_mv_vect)
+!!$  !! \brief scatter into a class(base_mlv_vect)
 !!$  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
 !!$  !! \param beta
 !!$  !! \param x(:) 
-!!$  subroutine c_base_mv_sctb(n,idx,x,beta,y)
+!!$  subroutine c_base_mlv_sctb(n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_spk_) :: beta, x(:)
@@ -2414,9 +2441,9 @@ contains
 !!$    call psi_sct(n,idx,x,beta,y%v)
 !!$    call y%set_host()
 !!$
-!!$  end subroutine c_base_mv_sctb
+!!$  end subroutine c_base_mlv_sctb
 !!$
-!!$  subroutine c_base_mv_sctb_x(i,n,idx,x,beta,y)
+!!$  subroutine c_base_mlv_sctb_x(i,n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i, n
 !!$    class(psb_c_base_multivect_type) :: idx
@@ -2425,6 +2452,6 @@ contains
 !!$    
 !!$    call y%sct(n,idx%v(i:),x,beta)
 !!$
-!!$  end subroutine c_base_mv_sctb_x
+!!$  end subroutine c_base_mlv_sctb_x
 end module psb_c_base_multivect_mod
 
diff --git a/base/modules/psb_d_base_vect_mod.f90 b/base/modules/psb_d_base_vect_mod.f90
index eb098253..657b3fc3 100644
--- a/base/modules/psb_d_base_vect_mod.f90
+++ b/base/modules/psb_d_base_vect_mod.f90
@@ -738,7 +738,7 @@ contains
   !> Function  base_dot_v
   !! \memberof  psb_d_base_vect_type
   !! \brief  Dot product by another base_vector
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y    The other (base_vect) to be multiplied by
   !!
   function d_base_dot_v(n,x,y) result(res)
@@ -772,7 +772,7 @@ contains
   !> Function  base_dot_a
   !! \memberof  psb_d_base_vect_type
   !! \brief  Dot product by a normal array
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y(:) The array to be multiplied by
   !!
   function d_base_dot_a(n,x,y) result(res)
@@ -795,7 +795,7 @@ contains
   !> Function  base_axpby_v
   !! \memberof  psb_d_base_vect_type
   !! \brief AXPBY  by a (base_vect) y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x     The class(base_vect) to be added
   !! \param beta scalar alpha
@@ -823,7 +823,7 @@ contains
   !> Function  base_axpby_a
   !! \memberof  psb_d_base_vect_type
   !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x(:) The array to be added
   !! \param beta scalar alpha
@@ -1339,7 +1339,7 @@ module psb_d_base_multivect_mod
   
   use psb_const_mod
   use psb_error_mod
-
+  use psb_realloc_mod
 
   !> \namespace  psb_base_mod  \class psb_d_base_vect_type
   !! The psb_d_base_vect_type 
@@ -1361,20 +1361,20 @@ module psb_d_base_multivect_mod
     !
     !  Constructors/allocators
     !
-    procedure, pass(x) :: bld_x    => d_base_mv_bld_x
-    procedure, pass(x) :: bld_n    => d_base_mv_bld_n
+    procedure, pass(x) :: bld_x    => d_base_mlv_bld_x
+    procedure, pass(x) :: bld_n    => d_base_mlv_bld_n
     generic, public    :: bld      => bld_x, bld_n
-    procedure, pass(x) :: all      => d_base_mv_all
-    procedure, pass(x) :: mold     => d_base_mv_mold
+    procedure, pass(x) :: all      => d_base_mlv_all
+    procedure, pass(x) :: mold     => d_base_mlv_mold
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
     ! in derived classes. 
     !
-    procedure, pass(x) :: ins      => d_base_mv_ins
-    procedure, pass(x) :: zero     => d_base_mv_zero
-    procedure, pass(x) :: asb      => d_base_mv_asb
-    procedure, pass(x) :: free     => d_base_mv_free
+    procedure, pass(x) :: ins      => d_base_mlv_ins
+    procedure, pass(x) :: zero     => d_base_mlv_zero
+    procedure, pass(x) :: asb      => d_base_mlv_asb
+    procedure, pass(x) :: free     => d_base_mlv_free
     !
     ! Sync: centerpiece of handling of external storage.
     ! Any derived class having extra storage upon sync
@@ -1382,66 +1382,66 @@ module psb_d_base_multivect_mod
     ! external side contain the same data. The base
     ! version is only a placeholder. 
     !
-    procedure, pass(x) :: sync     => d_base_mv_sync
-    procedure, pass(x) :: is_host  => d_base_mv_is_host
-    procedure, pass(x) :: is_dev   => d_base_mv_is_dev
-    procedure, pass(x) :: is_sync  => d_base_mv_is_sync
-    procedure, pass(x) :: set_host => d_base_mv_set_host
-    procedure, pass(x) :: set_dev  => d_base_mv_set_dev
-    procedure, pass(x) :: set_sync => d_base_mv_set_sync
+    procedure, pass(x) :: sync     => d_base_mlv_sync
+    procedure, pass(x) :: is_host  => d_base_mlv_is_host
+    procedure, pass(x) :: is_dev   => d_base_mlv_is_dev
+    procedure, pass(x) :: is_sync  => d_base_mlv_is_sync
+    procedure, pass(x) :: set_host => d_base_mlv_set_host
+    procedure, pass(x) :: set_dev  => d_base_mlv_set_dev
+    procedure, pass(x) :: set_sync => d_base_mlv_set_sync
 
     !
     ! Basic info
-    procedure, pass(x) :: get_nrows => d_base_mv_get_nrows
-    procedure, pass(x) :: get_ncols => d_base_mv_get_ncols
-    procedure, pass(x) :: sizeof    => d_base_mv_sizeof
-    procedure, nopass  :: get_fmt   => d_base_mv_get_fmt
+    procedure, pass(x) :: get_nrows => d_base_mlv_get_nrows
+    procedure, pass(x) :: get_ncols => d_base_mlv_get_ncols
+    procedure, pass(x) :: sizeof    => d_base_mlv_sizeof
+    procedure, nopass  :: get_fmt   => d_base_mlv_get_fmt
     !
     ! Set/get data from/to an external array; also
     ! overload assignment.
     !
-    procedure, pass(x) :: get_vect => d_base_mv_get_vect
-    procedure, pass(x) :: set_scal => d_base_mv_set_scal
-    procedure, pass(x) :: set_vect => d_base_mv_set_vect
+    procedure, pass(x) :: get_vect => d_base_mlv_get_vect
+    procedure, pass(x) :: set_scal => d_base_mlv_set_scal
+    procedure, pass(x) :: set_vect => d_base_mlv_set_vect
     generic, public    :: set      => set_vect, set_scal
 
     !
     ! Dot product and AXPBY
     !
-!!$    procedure, pass(x) :: dot_v    => d_base_mv_dot_v
-!!$    procedure, pass(x) :: dot_a    => d_base_mv_dot_a
-!!$    generic, public    :: dot      => dot_v, dot_a
-!!$    procedure, pass(y) :: axpby_v  => d_base_mv_axpby_v
-!!$    procedure, pass(y) :: axpby_a  => d_base_mv_axpby_a
-!!$    generic, public    :: axpby    => axpby_v, axpby_a
+    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(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
 !!$    !
 !!$    ! Vector by vector multiplication. Need all variants
 !!$    ! to handle multiple requirements from preconditioners
 !!$    !
-!!$    procedure, pass(y) :: mlt_v    => d_base_mv_mlt_v
-!!$    procedure, pass(y) :: mlt_a    => d_base_mv_mlt_a
-!!$    procedure, pass(z) :: mlt_a_2  => d_base_mv_mlt_a_2
-!!$    procedure, pass(z) :: mlt_v_2  => d_base_mv_mlt_v_2
-!!$    procedure, pass(z) :: mlt_va   => d_base_mv_mlt_va
-!!$    procedure, pass(z) :: mlt_av   => d_base_mv_mlt_av
+!!$    procedure, pass(y) :: mlt_v    => d_base_mlv_mlt_v
+!!$    procedure, pass(y) :: mlt_a    => d_base_mlv_mlt_a
+!!$    procedure, pass(z) :: mlt_a_2  => d_base_mlv_mlt_a_2
+!!$    procedure, pass(z) :: mlt_v_2  => d_base_mlv_mlt_v_2
+!!$    procedure, pass(z) :: mlt_va   => d_base_mlv_mlt_va
+!!$    procedure, pass(z) :: mlt_av   => d_base_mlv_mlt_av
 !!$    generic, public    :: mlt      => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
 !!$    !
 !!$    ! Scaling and norms
 !!$    !
-!!$    procedure, pass(x) :: scal     => d_base_mv_scal
-!!$    procedure, pass(x) :: nrm2     => d_base_mv_nrm2
-!!$    procedure, pass(x) :: amax     => d_base_mv_amax
-!!$    procedure, pass(x) :: asum     => d_base_mv_asum
+!!$    procedure, pass(x) :: scal     => d_base_mlv_scal
+!!$    procedure, pass(x) :: nrm2     => d_base_mlv_nrm2
+!!$    procedure, pass(x) :: amax     => d_base_mlv_amax
+!!$    procedure, pass(x) :: asum     => d_base_mlv_asum
 !!$    !
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => d_base_mv_gthab
-!!$    procedure, pass(x) :: gthzv    => d_base_mv_gthzv
-!!$    procedure, pass(x) :: gthzv_x  => d_base_mv_gthzv_x
+!!$    procedure, pass(x) :: gthab    => d_base_mlv_gthab
+!!$    procedure, pass(x) :: gthzv    => d_base_mlv_gthzv
+!!$    procedure, pass(x) :: gthzv_x  => d_base_mlv_gthzv_x
 !!$    generic, public    :: gth      => gthab, gthzv, gthzv_x
-!!$    procedure, pass(y) :: sctb     => d_base_mv_sctb
-!!$    procedure, pass(y) :: sctb_x   => d_base_mv_sctb_x
+!!$    procedure, pass(y) :: sctb     => d_base_mlv_sctb
+!!$    procedure, pass(y) :: sctb_x   => d_base_mlv_sctb_x
 !!$    generic, public    :: sct      => sctb, sctb_x
   end type psb_d_base_multivect_type
 
@@ -1491,7 +1491,7 @@ contains
   !! \brief     Build method from an array
   !!  \param   x(:)  input array to be copied
   !!
-  subroutine d_base_mv_bld_x(x,this)
+  subroutine d_base_mlv_bld_x(x,this)
     use psb_realloc_mod
     real(psb_dpk_), intent(in) :: this(:,:)
     class(psb_d_base_multivect_type), intent(inout) :: x
@@ -1499,12 +1499,12 @@ contains
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
     x%v(:,:)  = this(:,:)
 
-  end subroutine d_base_mv_bld_x
+  end subroutine d_base_mlv_bld_x
     
   !
   ! Create with size, but no initialization
@@ -1515,7 +1515,7 @@ contains
   !! \brief     Build method with size (uninitialized data)
   !!  \param    n    size to be allocated. 
   !!
-  subroutine d_base_mv_bld_n(x,m,n)
+  subroutine d_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: m,n
     class(psb_d_base_multivect_type), intent(inout) :: x
@@ -1524,16 +1524,16 @@ contains
     call psb_realloc(m,n,x%v,info)
     call x%asb(m,n,info)
 
-  end subroutine d_base_mv_bld_n
+  end subroutine d_base_mlv_bld_n
   
-  !> Function  base_mv_all:
+  !> Function  base_mlv_all:
   !! \memberof  psb_d_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
   !!  \param    n    size to be allocated. 
   !!  \param info  return code
   !!
-  subroutine d_base_mv_all(m,n, x, info)
+  subroutine d_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1543,15 +1543,15 @@ contains
     
     call psb_realloc(m,n,x%v,info)
     
-  end subroutine d_base_mv_all
+  end subroutine d_base_mlv_all
 
-  !> Function  base_mv_mold:
+  !> Function  base_mlv_mold:
   !! \memberof  psb_d_base_multivect_type
   !! \brief     Mold method: return a variable with the same dynamic type
   !!  \param    y returned variable
   !!  \param info  return code
   !!
-  subroutine d_base_mv_mold(x, y, info)
+  subroutine d_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1561,12 +1561,12 @@ contains
     
     allocate(psb_d_base_multivect_type :: y, stat=info)
 
-  end subroutine d_base_mv_mold
+  end subroutine d_base_mlv_mold
 
   !
   ! Insert a bunch of values at specified positions.
   !
-  !> Function  base_mv_ins:
+  !> Function  base_mlv_ins:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Insert coefficients. 
   !!
@@ -1590,7 +1590,7 @@ contains
   !!  \param info  return code
   !!
   !
-  subroutine d_base_mv_ins(n,irl,val,dupl,x,info)
+  subroutine d_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
     implicit none 
     class(psb_d_base_multivect_type), intent(inout)  :: x
@@ -1642,26 +1642,26 @@ contains
       end select
     end if
     if (info /= 0) then 
-      call psb_errpush(info,'base_mv_vect_ins')
+      call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
 
-  end subroutine d_base_mv_ins
+  end subroutine d_base_mlv_ins
 
   !
-  !> Function  base_mv_zero
+  !> Function  base_mlv_zero
   !! \memberof  psb_d_base_multivect_type
   !! \brief Zero out contents
   !!
   !
-  subroutine d_base_mv_zero(x)
+  subroutine d_base_mlv_zero(x)
     use psi_serial_mod
     implicit none 
     class(psb_d_base_multivect_type), intent(inout)    :: x
     
     if (allocated(x%v)) x%v=dzero
 
-  end subroutine d_base_mv_zero
+  end subroutine d_base_mlv_zero
 
   
   !
@@ -1669,7 +1669,7 @@ contains
   ! For derived classes: after this the vector
   ! storage is supposed to be in sync.
   !
-  !> Function  base_mv_asb:
+  !> Function  base_mlv_asb:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
   !!           
@@ -1678,7 +1678,7 @@ contains
   !!
   !
  
-  subroutine d_base_mv_asb(m,n, x, info)
+  subroutine d_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1691,18 +1691,18 @@ contains
     if (info /= 0) &
          & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
 
-  end subroutine d_base_mv_asb
+  end subroutine d_base_mlv_asb
 
 
   !
-  !> Function  base_mv_free:
+  !> Function  base_mlv_free:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Free vector
   !!           
   !!  \param info  return code
   !!
   !
-  subroutine d_base_mv_free(x, info)
+  subroutine d_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1714,7 +1714,7 @@ contains
     if (info /= 0) call & 
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
         
-  end subroutine d_base_mv_free
+  end subroutine d_base_mlv_free
 
   
 
@@ -1723,106 +1723,106 @@ contains
   ! a placeholder.
   ! 
   !
-  !> Function  base_mv_sync:
+  !> Function  base_mlv_sync:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Sync: base version is a no-op.
   !!           
   !
-  subroutine d_base_mv_sync(x)
+  subroutine d_base_mlv_sync(x)
     implicit none 
     class(psb_d_base_multivect_type), intent(inout) :: x
     
-  end subroutine d_base_mv_sync
+  end subroutine d_base_mlv_sync
 
   !
-  !> Function  base_mv_set_host:
+  !> Function  base_mlv_set_host:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_host: base version is a no-op.
   !!           
   !
-  subroutine d_base_mv_set_host(x)
+  subroutine d_base_mlv_set_host(x)
     implicit none 
     class(psb_d_base_multivect_type), intent(inout) :: x
     
-  end subroutine d_base_mv_set_host
+  end subroutine d_base_mlv_set_host
 
   !
-  !> Function  base_mv_set_dev:
+  !> Function  base_mlv_set_dev:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
   !!           
   !
-  subroutine d_base_mv_set_dev(x)
+  subroutine d_base_mlv_set_dev(x)
     implicit none 
     class(psb_d_base_multivect_type), intent(inout) :: x
     
-  end subroutine d_base_mv_set_dev
+  end subroutine d_base_mlv_set_dev
 
   !
-  !> Function  base_mv_set_sync:
+  !> Function  base_mlv_set_sync:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
   !!           
   !
-  subroutine d_base_mv_set_sync(x)
+  subroutine d_base_mlv_set_sync(x)
     implicit none 
     class(psb_d_base_multivect_type), intent(inout) :: x
     
-  end subroutine d_base_mv_set_sync
+  end subroutine d_base_mlv_set_sync
 
   !
-  !> Function  base_mv_is_dev:
+  !> Function  base_mlv_is_dev:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on external device    .
   !!           
   !
-  function d_base_mv_is_dev(x) result(res)
+  function d_base_mlv_is_dev(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
   
     res = .false.
-  end function d_base_mv_is_dev
+  end function d_base_mlv_is_dev
   
   !
-  !> Function  base_mv_is_host
+  !> Function  base_mlv_is_host
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on standard memory    .
   !!           
   !
-  function d_base_mv_is_host(x) result(res)
+  function d_base_mlv_is_host(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function d_base_mv_is_host
+  end function d_base_mlv_is_host
 
   !
-  !> Function  base_mv_is_sync
+  !> Function  base_mlv_is_sync
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on sync               .
   !!           
   !
-  function d_base_mv_is_sync(x) result(res)
+  function d_base_mlv_is_sync(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function d_base_mv_is_sync
+  end function d_base_mlv_is_sync
 
 
   !
   ! Size info. 
   !
   !
-  !> Function  base_mv_get_nrows
+  !> Function  base_mlv_get_nrows
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Number of entries
   !!           
   !
-  function d_base_mv_get_nrows(x) result(res)
+  function d_base_mlv_get_nrows(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1830,9 +1830,9 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,1)
 
-  end function d_base_mv_get_nrows
+  end function d_base_mlv_get_nrows
 
-  function d_base_mv_get_ncols(x) result(res)
+  function d_base_mlv_get_ncols(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1840,15 +1840,15 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,2)
 
-  end function d_base_mv_get_ncols
+  end function d_base_mlv_get_ncols
   
   !
-  !> Function  base_mv_get_sizeof
+  !> Function  base_mlv_get_sizeof
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Size in bytesa
   !!           
   !
-  function d_base_mv_sizeof(x) result(res)
+  function d_base_mlv_sizeof(x) result(res)
     implicit none 
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_long_int_k_) :: res
@@ -1856,30 +1856,30 @@ contains
     ! Force 8-byte integers.
     res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()
 
-  end function d_base_mv_sizeof
+  end function d_base_mlv_sizeof
 
   !
-  !> Function  base_mv_get_fmt
+  !> Function  base_mlv_get_fmt
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Format
   !!           
   !
-  function d_base_mv_get_fmt() result(res)
+  function d_base_mlv_get_fmt() result(res)
     implicit none 
     character(len=5) :: res
     res = 'BASE'
-  end function d_base_mv_get_fmt
+  end function d_base_mlv_get_fmt
   
 
   !
   !
   !
-  !> Function  base_mv_get_vect
+  !> Function  base_mlv_get_vect
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
   !    
-  function  d_base_mv_get_vect(x) result(res)
+  function  d_base_mlv_get_vect(x) result(res)
     class(psb_d_base_multivect_type), intent(inout) :: x
     real(psb_dpk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
@@ -1889,37 +1889,37 @@ contains
     call x%sync()
     allocate(res(m,n),stat=info) 
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
     res(1:m,1:n) = x%v(1:m,1:n)
-  end function d_base_mv_get_vect
+  end function d_base_mlv_get_vect
     
   !
   ! Reset all values 
   !
   !
-  !> Function  base_mv_set_scal
+  !> Function  base_mlv_set_scal
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Set all entries
   !! \param val   The value to set
   !!
-  subroutine d_base_mv_set_scal(x,val)
+  subroutine d_base_mlv_set_scal(x,val)
     class(psb_d_base_multivect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val
         
     integer(psb_ipk_) :: info
     x%v = val
     
-  end subroutine d_base_mv_set_scal
+  end subroutine d_base_mlv_set_scal
 
   !
-  !> Function  base_mv_set_vect
+  !> Function  base_mlv_set_vect
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Set all entries
   !! \param val(:)  The vector to be copied in 
   !!
-  subroutine d_base_mv_set_vect(x,val)
+  subroutine d_base_mlv_set_vect(x,val)
     class(psb_d_base_multivect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr
@@ -1934,121 +1934,148 @@ contains
       x%v = val
     end if
 
-  end subroutine d_base_mv_set_vect
+  end subroutine d_base_mlv_set_vect
 
-!!$  !
-!!$  ! Dot products 
-!!$  ! 
-!!$  !
-!!$  !> Function  base_mv_dot_v
-!!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief  Dot product by another base_mv_vector
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y    The other (base_mv_vect) to be multiplied by
-!!$  !!
-!!$  function d_base_mv_dot_v(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_d_base_multivect_type), intent(inout) :: x, y
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_dpk_)                :: res
-!!$    real(psb_dpk_), external      :: ddot
-!!$    
-!!$    res = izero
-!!$    !
-!!$    ! Note: this is the base implementation.
-!!$    !  When we get here, we are sure that X is of
-!!$    !  TYPE psb_d_base_mv_vect.
-!!$    !  If Y is not, throw the burden on it, implicitly
-!!$    !  calling dot_a
-!!$    !
-!!$    select type(yy => y)
-!!$    type is (psb_d_base_multivect_type)
-!!$      res = ddot(n,x%v,1,y%v,1)
-!!$    class default
-!!$      res = y%dot(n,x%v)
-!!$    end select
-!!$
-!!$  end function d_base_mv_dot_v
-!!$
-!!$  !
-!!$  ! Base workhorse is good old BLAS1
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_dot_a
-!!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief  Dot product by a normal array
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y(:) The array to be multiplied by
-!!$  !!
-!!$  function d_base_mv_dot_a(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_d_base_multivect_type), intent(inout) :: x
-!!$    real(psb_dpk_), intent(in)    :: y(:)
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_dpk_)                :: res
-!!$    integer(psb_ipk_), external      :: ddot
-!!$    
-!!$    res = ddot(n,y,1,x%v,1)
-!!$
-!!$  end function d_base_mv_dot_a
+  !
+  ! Dot products 
+  ! 
+  !
+  !> Function  base_mlv_dot_v
+  !! \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_mlv_vect) to be multiplied by
+  !!
+  function d_base_mlv_dot_v(n,x,y) result(res)
+    implicit none 
+    class(psb_d_base_multivect_type), intent(inout) :: x, 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()
+    res = dzero
+    !
+    ! Note: this is the base implementation.
+    !  When we get here, we are sure that X is of
+    !  TYPE psb_d_base_mlv_vect (or its class does not care). 
+    !  If Y is not, throw the burden on it, implicitly
+    !  calling dot_a
+    !
+    select type(yy => y)
+    type is (psb_d_base_multivect_type)
+      if (y%is_dev()) call y%sync()
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+      allocate(res(nc))
+      do j=1,nc
+        res(j) = ddot(n,x%v(:,j),1,y%v(:,j),1)
+      end do
+    class default
+      res = y%dot(n,x%v)
+    end select
+
+  end function d_base_mlv_dot_v
+
+  !
+  ! Base workhorse is good old BLAS1
+  !
+  !
+  !> Function  base_mlv_dot_a
+  !! \memberof  psb_d_base_multivect_type
+  !! \brief  Dot product by a normal array
+  !! \param n    Number of entries to be considered
+  !! \param y(:) The array to be multiplied by
+  !!
+  function d_base_mlv_dot_a(n,x,y) result(res)
+    implicit none 
+    class(psb_d_base_multivect_type), intent(inout) :: x
+    real(psb_dpk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_dpk_), allocatable     :: res(:)
+    integer(psb_ipk_), external      :: ddot
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    nc = min(psb_size(x%v,2),size(y,2))
+    allocate(res(nc))
+    do j=1,nc
+      res(j) = ddot(n,x%v(:,j),1,y(:,j),1)
+    end do
+
+  end function d_base_mlv_dot_a
     
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_v
-!!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief AXPBY  by a (base_mv_vect) y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x     The class(base_mv_vect) to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine d_base_mv_axpby_v(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    class(psb_d_base_multivect_type), intent(inout)  :: x
-!!$    class(psb_d_base_multivect_type), intent(inout)  :: y
-!!$    real(psb_dpk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    select type(xx => x)
-!!$    type is (psb_d_base_multivect_type)
-!!$      call psb_geaxpby(m,alpha,x%v,beta,y%v,info)
-!!$    class default
-!!$      call y%axpby(m,alpha,x%v,beta,info)
-!!$    end select
-!!$
-!!$  end subroutine d_base_mv_axpby_v
-!!$
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_a
-!!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x(:) The array to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine d_base_mv_axpby_a(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    real(psb_dpk_), intent(in)        :: x(:)
-!!$    class(psb_d_base_multivect_type), intent(inout)  :: y
-!!$    real(psb_dpk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    call psb_geaxpby(m,alpha,x,beta,y%v,info)
-!!$    
-!!$  end subroutine d_base_mv_axpby_a
+  !
+  ! 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=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_mlv_vect) to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine d_base_mlv_axpby_v(m,alpha, x, beta, y, info, n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    class(psb_d_base_multivect_type), intent(inout)  :: x
+    class(psb_d_base_multivect_type), intent(inout)  :: y
+    real(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+    end if
+    select type(xx => x)
+    type is (psb_d_base_multivect_type)
+      call psb_geaxpby(m,nc,alpha,x%v,beta,y%v,info)
+    class default
+      call y%axpby(m,alpha,x%v,beta,info,n=n)
+    end select
+
+  end subroutine d_base_mlv_axpby_v
+
+  !
+  ! AXPBY is invoked via Y, hence the structure below. 
+  !
+  !
+  !> Function  base_mlv_axpby_a
+  !! \memberof  psb_d_base_multivect_type
+  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x(:) The array to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine d_base_mlv_axpby_a(m,alpha, x, beta, y, info,n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    real(psb_dpk_), intent(in)        :: x(:,:)
+    class(psb_d_base_multivect_type), intent(inout)  :: y
+    real(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(size(x,2),psb_size(y%v,2))
+    end if
+    
+    call psb_geaxpby(m,nc,alpha,x,beta,y%v,info)
+    
+  end subroutine d_base_mlv_axpby_a
 
   
 !!$  !
@@ -2060,13 +2087,13 @@ contains
 !!$  !  of the involved entities
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief Vector entry-by-entry multiply  by a base_mv_vect array y=x*y
-!!$  !! \param x   The class(base_mv_vect) to be multiplied by
+!!$  !! \brief Vector entry-by-entry multiply  by a base_mlv_vect array y=x*y
+!!$  !! \param x   The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine d_base_mv_mlt_v(x, y, info)
+!!$  subroutine d_base_mlv_mlt_v(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: x
@@ -2085,16 +2112,16 @@ contains
 !!$      call y%mlt(x%v,info)
 !!$    end select
 !!$
-!!$  end subroutine d_base_mv_mlt_v
+!!$  end subroutine d_base_mlv_mlt_v
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief Vector entry-by-entry multiply  by a normal array y=x*y
 !!$  !! \param x(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine d_base_mv_mlt_a(x, y, info)
+!!$  subroutine d_base_mlv_mlt_a(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_dpk_), intent(in)        :: x(:)
@@ -2108,11 +2135,11 @@ contains
 !!$      y%v(i) = y%v(i)*x(i)
 !!$    end do
 !!$    
-!!$  end subroutine d_base_mv_mlt_a
+!!$  end subroutine d_base_mlv_mlt_a
 !!$
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a_2
+!!$  !> Function  base_mlv_mlt_a_2
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by normal arrays
 !!$  !!        z=beta*z+alpha*x*y
@@ -2122,7 +2149,7 @@ contains
 !!$  !! \param y(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine d_base_mv_mlt_a_2(alpha,x,y,beta,z,info)
+!!$  subroutine d_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_dpk_), intent(in)        :: alpha,beta
@@ -2187,20 +2214,20 @@ contains
 !!$        end if
 !!$      end if
 !!$    end if
-!!$  end subroutine d_base_mv_mlt_a_2
+!!$  end subroutine d_base_mlv_mlt_a_2
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_v_2
+!!$  !> Function  base_mlv_mlt_v_2
 !!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mv_vect)
+!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mlv_vect)
 !!$  !!        z=beta*z+alpha*x*y
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  !! \param x  The class(base_mv_vect) to be multiplied b
-!!$  !! \param y  The class(base_mv_vect) to be multiplied by
+!!$  !! \param x  The class(base_mlv_vect) to be multiplied b
+!!$  !! \param y  The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine d_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
+!!$  subroutine d_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
 !!$    use psi_serial_mod
 !!$    use psb_string_mod
 !!$    implicit none 
@@ -2227,9 +2254,9 @@ contains
 !!$      if (conjgx_) x%v=(x%v)
 !!$      if (conjgy_) y%v=(y%v)
 !!$    end if
-!!$  end subroutine d_base_mv_mlt_v_2
+!!$  end subroutine d_base_mlv_mlt_v_2
 !!$
-!!$  subroutine d_base_mv_mlt_av(alpha,x,y,beta,z,info)
+!!$  subroutine d_base_mlv_mlt_av(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_dpk_), intent(in)        :: alpha,beta
@@ -2243,9 +2270,9 @@ contains
 !!$    
 !!$    call z%mlt(alpha,x,y%v,beta,info)
 !!$
-!!$  end subroutine d_base_mv_mlt_av
+!!$  end subroutine d_base_mlv_mlt_av
 !!$
-!!$  subroutine d_base_mv_mlt_va(alpha,x,y,beta,z,info)
+!!$  subroutine d_base_mlv_mlt_va(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_dpk_), intent(in)        :: alpha,beta
@@ -2259,18 +2286,18 @@ contains
 !!$    
 !!$    call z%mlt(alpha,y,x,beta,info)
 !!$
-!!$  end subroutine d_base_mv_mlt_va
+!!$  end subroutine d_base_mlv_mlt_va
 !!$
 !!$
 !!$  !
 !!$  ! Simple scaling 
 !!$  !
-!!$  !> Function  base_mv_scal
+!!$  !> Function  base_mlv_scal
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief Scale all entries  x = alpha*x
 !!$  !! \param alpha   The multiplier
 !!$  !!
-!!$  subroutine d_base_mv_scal(alpha, x)
+!!$  subroutine d_base_mlv_scal(alpha, x)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: x
@@ -2278,16 +2305,16 @@ contains
 !!$    
 !!$    if (allocated(x%v)) x%v = alpha*x%v
 !!$
-!!$  end subroutine d_base_mv_scal
+!!$  end subroutine d_base_mlv_scal
 !!$  
 !!$  !
 !!$  ! Norms 1, 2 and infinity
 !!$  !
-!!$  !> Function  base_mv_nrm2
+!!$  !> Function  base_mlv_nrm2
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief 2-norm |x(1:n)|_2
 !!$  !! \param n  how many entries to consider
-!!$  function d_base_mv_nrm2(n,x) result(res)
+!!$  function d_base_mlv_nrm2(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_d_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2296,14 +2323,14 @@ contains
 !!$    
 !!$    res =  dnrm2(n,x%v,1)
 !!$
-!!$  end function d_base_mv_nrm2
+!!$  end function d_base_mlv_nrm2
 !!$  
 !!$  !
-!!$  !> Function  base_mv_amax
+!!$  !> Function  base_mlv_amax
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief infinity-norm |x(1:n)|_\infty
 !!$  !! \param n  how many entries to consider
-!!$  function d_base_mv_amax(n,x) result(res)
+!!$  function d_base_mlv_amax(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_d_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2311,14 +2338,14 @@ contains
 !!$    
 !!$    res =  maxval(abs(x%v(1:n)))
 !!$
-!!$  end function d_base_mv_amax
+!!$  end function d_base_mlv_amax
 !!$
 !!$  !
-!!$  !> Function  base_mv_asum
+!!$  !> Function  base_mlv_asum
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief 1-norm |x(1:n)|_1
 !!$  !! \param n  how many entries to consider
-!!$  function d_base_mv_asum(n,x) result(res)
+!!$  function d_base_mlv_asum(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_d_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2326,14 +2353,14 @@ contains
 !!$    
 !!$    res =  sum(abs(x%v(1:n)))
 !!$
-!!$  end function d_base_mv_asum
+!!$  end function d_base_mlv_asum
 !!$  
 !!$  
 !!$  !
 !!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_gthab
+!!$  !> Function  base_mlv_gthab
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief gather into an array
 !!$  !!    Y = beta * Y + alpha * X(IDX(:))
@@ -2341,7 +2368,7 @@ contains
 !!$  !! \param idx(:) indices
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  subroutine d_base_mv_gthab(n,idx,alpha,x,beta,y)
+!!$  subroutine d_base_mlv_gthab(n,idx,alpha,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_dpk_) :: alpha, beta, y(:)
@@ -2350,17 +2377,17 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,alpha,x%v,beta,y)
 !!$
-!!$  end subroutine d_base_mv_gthab
+!!$  end subroutine d_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine d_base_mv_gthzv_x(i,n,idx,x,y)
+!!$  subroutine d_base_mlv_gthzv_x(i,n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i,n
 !!$    class(psb_d_base_multivect_type) :: idx
@@ -2369,18 +2396,18 @@ contains
 !!$    
 !!$    call x%gth(n,idx%v(i:),y)
 !!$
-!!$  end subroutine d_base_mv_gthzv_x
+!!$  end subroutine d_base_mlv_gthzv_x
 !!$
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_d_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine d_base_mv_gthzv(n,idx,x,y)
+!!$  subroutine d_base_mlv_gthzv(n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_dpk_) ::  y(:)
@@ -2389,22 +2416,22 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,x%v,y)
 !!$
-!!$  end subroutine d_base_mv_gthzv
+!!$  end subroutine d_base_mlv_gthzv
 !!$
 !!$  !
 !!$  ! Scatter: 
 !!$  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
 !!$  ! 
 !!$  !
-!!$  !> Function  base_mv_sctb
+!!$  !> Function  base_mlv_sctb
 !!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief scatter into a class(base_mv_vect)
+!!$  !! \brief scatter into a class(base_mlv_vect)
 !!$  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
 !!$  !! \param beta
 !!$  !! \param x(:) 
-!!$  subroutine d_base_mv_sctb(n,idx,x,beta,y)
+!!$  subroutine d_base_mlv_sctb(n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_dpk_) :: beta, x(:)
@@ -2414,9 +2441,9 @@ contains
 !!$    call psi_sct(n,idx,x,beta,y%v)
 !!$    call y%set_host()
 !!$
-!!$  end subroutine d_base_mv_sctb
+!!$  end subroutine d_base_mlv_sctb
 !!$
-!!$  subroutine d_base_mv_sctb_x(i,n,idx,x,beta,y)
+!!$  subroutine d_base_mlv_sctb_x(i,n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i, n
 !!$    class(psb_d_base_multivect_type) :: idx
@@ -2425,6 +2452,6 @@ contains
 !!$    
 !!$    call y%sct(n,idx%v(i:),x,beta)
 !!$
-!!$  end subroutine d_base_mv_sctb_x
+!!$  end subroutine d_base_mlv_sctb_x
 end module psb_d_base_multivect_mod
 
diff --git a/base/modules/psb_i_base_vect_mod.f90 b/base/modules/psb_i_base_vect_mod.f90
index a5814f3d..a7ba7747 100644
--- a/base/modules/psb_i_base_vect_mod.f90
+++ b/base/modules/psb_i_base_vect_mod.f90
@@ -879,7 +879,7 @@ module psb_i_base_multivect_mod
   
   use psb_const_mod
   use psb_error_mod
-
+  use psb_realloc_mod
 
   !> \namespace  psb_base_mod  \class psb_i_base_vect_type
   !! The psb_i_base_vect_type 
@@ -901,20 +901,20 @@ module psb_i_base_multivect_mod
     !
     !  Constructors/allocators
     !
-    procedure, pass(x) :: bld_x    => i_base_mv_bld_x
-    procedure, pass(x) :: bld_n    => i_base_mv_bld_n
+    procedure, pass(x) :: bld_x    => i_base_mlv_bld_x
+    procedure, pass(x) :: bld_n    => i_base_mlv_bld_n
     generic, public    :: bld      => bld_x, bld_n
-    procedure, pass(x) :: all      => i_base_mv_all
-    procedure, pass(x) :: mold     => i_base_mv_mold
+    procedure, pass(x) :: all      => i_base_mlv_all
+    procedure, pass(x) :: mold     => i_base_mlv_mold
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
     ! in derived classes. 
     !
-    procedure, pass(x) :: ins      => i_base_mv_ins
-    procedure, pass(x) :: zero     => i_base_mv_zero
-    procedure, pass(x) :: asb      => i_base_mv_asb
-    procedure, pass(x) :: free     => i_base_mv_free
+    procedure, pass(x) :: ins      => i_base_mlv_ins
+    procedure, pass(x) :: zero     => i_base_mlv_zero
+    procedure, pass(x) :: asb      => i_base_mlv_asb
+    procedure, pass(x) :: free     => i_base_mlv_free
     !
     ! Sync: centerpiece of handling of external storage.
     ! Any derived class having extra storage upon sync
@@ -922,27 +922,27 @@ module psb_i_base_multivect_mod
     ! external side contain the same data. The base
     ! version is only a placeholder. 
     !
-    procedure, pass(x) :: sync     => i_base_mv_sync
-    procedure, pass(x) :: is_host  => i_base_mv_is_host
-    procedure, pass(x) :: is_dev   => i_base_mv_is_dev
-    procedure, pass(x) :: is_sync  => i_base_mv_is_sync
-    procedure, pass(x) :: set_host => i_base_mv_set_host
-    procedure, pass(x) :: set_dev  => i_base_mv_set_dev
-    procedure, pass(x) :: set_sync => i_base_mv_set_sync
+    procedure, pass(x) :: sync     => i_base_mlv_sync
+    procedure, pass(x) :: is_host  => i_base_mlv_is_host
+    procedure, pass(x) :: is_dev   => i_base_mlv_is_dev
+    procedure, pass(x) :: is_sync  => i_base_mlv_is_sync
+    procedure, pass(x) :: set_host => i_base_mlv_set_host
+    procedure, pass(x) :: set_dev  => i_base_mlv_set_dev
+    procedure, pass(x) :: set_sync => i_base_mlv_set_sync
 
     !
     ! Basic info
-    procedure, pass(x) :: get_nrows => i_base_mv_get_nrows
-    procedure, pass(x) :: get_ncols => i_base_mv_get_ncols
-    procedure, pass(x) :: sizeof    => i_base_mv_sizeof
-    procedure, nopass  :: get_fmt   => i_base_mv_get_fmt
+    procedure, pass(x) :: get_nrows => i_base_mlv_get_nrows
+    procedure, pass(x) :: get_ncols => i_base_mlv_get_ncols
+    procedure, pass(x) :: sizeof    => i_base_mlv_sizeof
+    procedure, nopass  :: get_fmt   => i_base_mlv_get_fmt
     !
     ! Set/get data from/to an external array; also
     ! overload assignment.
     !
-    procedure, pass(x) :: get_vect => i_base_mv_get_vect
-    procedure, pass(x) :: set_scal => i_base_mv_set_scal
-    procedure, pass(x) :: set_vect => i_base_mv_set_vect
+    procedure, pass(x) :: get_vect => i_base_mlv_get_vect
+    procedure, pass(x) :: set_scal => i_base_mlv_set_scal
+    procedure, pass(x) :: set_vect => i_base_mlv_set_vect
     generic, public    :: set      => set_vect, set_scal
 
   end type psb_i_base_multivect_type
@@ -993,7 +993,7 @@ contains
   !! \brief     Build method from an array
   !!  \param   x(:)  input array to be copied
   !!
-  subroutine i_base_mv_bld_x(x,this)
+  subroutine i_base_mlv_bld_x(x,this)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: this(:,:)
     class(psb_i_base_multivect_type), intent(inout) :: x
@@ -1001,12 +1001,12 @@ contains
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
     x%v(:,:)  = this(:,:)
 
-  end subroutine i_base_mv_bld_x
+  end subroutine i_base_mlv_bld_x
     
   !
   ! Create with size, but no initialization
@@ -1017,7 +1017,7 @@ contains
   !! \brief     Build method with size (uninitialized data)
   !!  \param    n    size to be allocated. 
   !!
-  subroutine i_base_mv_bld_n(x,m,n)
+  subroutine i_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: m,n
     class(psb_i_base_multivect_type), intent(inout) :: x
@@ -1026,16 +1026,16 @@ contains
     call psb_realloc(m,n,x%v,info)
     call x%asb(m,n,info)
 
-  end subroutine i_base_mv_bld_n
+  end subroutine i_base_mlv_bld_n
   
-  !> Function  base_mv_all:
+  !> Function  base_mlv_all:
   !! \memberof  psb_i_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
   !!  \param    n    size to be allocated. 
   !!  \param info  return code
   !!
-  subroutine i_base_mv_all(m,n, x, info)
+  subroutine i_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1045,15 +1045,15 @@ contains
     
     call psb_realloc(m,n,x%v,info)
     
-  end subroutine i_base_mv_all
+  end subroutine i_base_mlv_all
 
-  !> Function  base_mv_mold:
+  !> Function  base_mlv_mold:
   !! \memberof  psb_i_base_multivect_type
   !! \brief     Mold method: return a variable with the same dynamic type
   !!  \param    y returned variable
   !!  \param info  return code
   !!
-  subroutine i_base_mv_mold(x, y, info)
+  subroutine i_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1063,12 +1063,12 @@ contains
     
     allocate(psb_i_base_multivect_type :: y, stat=info)
 
-  end subroutine i_base_mv_mold
+  end subroutine i_base_mlv_mold
 
   !
   ! Insert a bunch of values at specified positions.
   !
-  !> Function  base_mv_ins:
+  !> Function  base_mlv_ins:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Insert coefficients. 
   !!
@@ -1092,7 +1092,7 @@ contains
   !!  \param info  return code
   !!
   !
-  subroutine i_base_mv_ins(n,irl,val,dupl,x,info)
+  subroutine i_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
     implicit none 
     class(psb_i_base_multivect_type), intent(inout)  :: x
@@ -1144,26 +1144,26 @@ contains
       end select
     end if
     if (info /= 0) then 
-      call psb_errpush(info,'base_mv_vect_ins')
+      call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
 
-  end subroutine i_base_mv_ins
+  end subroutine i_base_mlv_ins
 
   !
-  !> Function  base_mv_zero
+  !> Function  base_mlv_zero
   !! \memberof  psb_i_base_multivect_type
   !! \brief Zero out contents
   !!
   !
-  subroutine i_base_mv_zero(x)
+  subroutine i_base_mlv_zero(x)
     use psi_serial_mod
     implicit none 
     class(psb_i_base_multivect_type), intent(inout)    :: x
     
     if (allocated(x%v)) x%v=izero
 
-  end subroutine i_base_mv_zero
+  end subroutine i_base_mlv_zero
 
   
   !
@@ -1171,7 +1171,7 @@ contains
   ! For derived classes: after this the vector
   ! storage is supposed to be in sync.
   !
-  !> Function  base_mv_asb:
+  !> Function  base_mlv_asb:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
   !!           
@@ -1180,7 +1180,7 @@ contains
   !!
   !
  
-  subroutine i_base_mv_asb(m,n, x, info)
+  subroutine i_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1193,18 +1193,18 @@ contains
     if (info /= 0) &
          & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
 
-  end subroutine i_base_mv_asb
+  end subroutine i_base_mlv_asb
 
 
   !
-  !> Function  base_mv_free:
+  !> Function  base_mlv_free:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Free vector
   !!           
   !!  \param info  return code
   !!
   !
-  subroutine i_base_mv_free(x, info)
+  subroutine i_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1216,7 +1216,7 @@ contains
     if (info /= 0) call & 
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
         
-  end subroutine i_base_mv_free
+  end subroutine i_base_mlv_free
 
   
 
@@ -1225,106 +1225,106 @@ contains
   ! a placeholder.
   ! 
   !
-  !> Function  base_mv_sync:
+  !> Function  base_mlv_sync:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Sync: base version is a no-op.
   !!           
   !
-  subroutine i_base_mv_sync(x)
+  subroutine i_base_mlv_sync(x)
     implicit none 
     class(psb_i_base_multivect_type), intent(inout) :: x
     
-  end subroutine i_base_mv_sync
+  end subroutine i_base_mlv_sync
 
   !
-  !> Function  base_mv_set_host:
+  !> Function  base_mlv_set_host:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Set_host: base version is a no-op.
   !!           
   !
-  subroutine i_base_mv_set_host(x)
+  subroutine i_base_mlv_set_host(x)
     implicit none 
     class(psb_i_base_multivect_type), intent(inout) :: x
     
-  end subroutine i_base_mv_set_host
+  end subroutine i_base_mlv_set_host
 
   !
-  !> Function  base_mv_set_dev:
+  !> Function  base_mlv_set_dev:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
   !!           
   !
-  subroutine i_base_mv_set_dev(x)
+  subroutine i_base_mlv_set_dev(x)
     implicit none 
     class(psb_i_base_multivect_type), intent(inout) :: x
     
-  end subroutine i_base_mv_set_dev
+  end subroutine i_base_mlv_set_dev
 
   !
-  !> Function  base_mv_set_sync:
+  !> Function  base_mlv_set_sync:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
   !!           
   !
-  subroutine i_base_mv_set_sync(x)
+  subroutine i_base_mlv_set_sync(x)
     implicit none 
     class(psb_i_base_multivect_type), intent(inout) :: x
     
-  end subroutine i_base_mv_set_sync
+  end subroutine i_base_mlv_set_sync
 
   !
-  !> Function  base_mv_is_dev:
+  !> Function  base_mlv_is_dev:
   !! \memberof  psb_i_base_multivect_type
   !! \brief Is  vector on external device    .
   !!           
   !
-  function i_base_mv_is_dev(x) result(res)
+  function i_base_mlv_is_dev(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     logical  :: res
   
     res = .false.
-  end function i_base_mv_is_dev
+  end function i_base_mlv_is_dev
   
   !
-  !> Function  base_mv_is_host
+  !> Function  base_mlv_is_host
   !! \memberof  psb_i_base_multivect_type
   !! \brief Is  vector on standard memory    .
   !!           
   !
-  function i_base_mv_is_host(x) result(res)
+  function i_base_mlv_is_host(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function i_base_mv_is_host
+  end function i_base_mlv_is_host
 
   !
-  !> Function  base_mv_is_sync
+  !> Function  base_mlv_is_sync
   !! \memberof  psb_i_base_multivect_type
   !! \brief Is  vector on sync               .
   !!           
   !
-  function i_base_mv_is_sync(x) result(res)
+  function i_base_mlv_is_sync(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function i_base_mv_is_sync
+  end function i_base_mlv_is_sync
 
 
   !
   ! Size info. 
   !
   !
-  !> Function  base_mv_get_nrows
+  !> Function  base_mlv_get_nrows
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Number of entries
   !!           
   !
-  function i_base_mv_get_nrows(x) result(res)
+  function i_base_mlv_get_nrows(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1332,9 +1332,9 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,1)
 
-  end function i_base_mv_get_nrows
+  end function i_base_mlv_get_nrows
 
-  function i_base_mv_get_ncols(x) result(res)
+  function i_base_mlv_get_ncols(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1342,15 +1342,15 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,2)
 
-  end function i_base_mv_get_ncols
+  end function i_base_mlv_get_ncols
   
   !
-  !> Function  base_mv_get_sizeof
+  !> Function  base_mlv_get_sizeof
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Size in bytesa
   !!           
   !
-  function i_base_mv_sizeof(x) result(res)
+  function i_base_mlv_sizeof(x) result(res)
     implicit none 
     class(psb_i_base_multivect_type), intent(in) :: x
     integer(psb_long_int_k_) :: res
@@ -1358,30 +1358,30 @@ contains
     ! Force 8-byte integers.
     res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()
 
-  end function i_base_mv_sizeof
+  end function i_base_mlv_sizeof
 
   !
-  !> Function  base_mv_get_fmt
+  !> Function  base_mlv_get_fmt
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Format
   !!           
   !
-  function i_base_mv_get_fmt() result(res)
+  function i_base_mlv_get_fmt() result(res)
     implicit none 
     character(len=5) :: res
     res = 'BASE'
-  end function i_base_mv_get_fmt
+  end function i_base_mlv_get_fmt
   
 
   !
   !
   !
-  !> Function  base_mv_get_vect
+  !> Function  base_mlv_get_vect
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
   !    
-  function  i_base_mv_get_vect(x) result(res)
+  function  i_base_mlv_get_vect(x) result(res)
     class(psb_i_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
@@ -1391,37 +1391,37 @@ contains
     call x%sync()
     allocate(res(m,n),stat=info) 
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
     res(1:m,1:n) = x%v(1:m,1:n)
-  end function i_base_mv_get_vect
+  end function i_base_mlv_get_vect
     
   !
   ! Reset all values 
   !
   !
-  !> Function  base_mv_set_scal
+  !> Function  base_mlv_set_scal
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Set all entries
   !! \param val   The value to set
   !!
-  subroutine i_base_mv_set_scal(x,val)
+  subroutine i_base_mlv_set_scal(x,val)
     class(psb_i_base_multivect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(in) :: val
         
     integer(psb_ipk_) :: info
     x%v = val
     
-  end subroutine i_base_mv_set_scal
+  end subroutine i_base_mlv_set_scal
 
   !
-  !> Function  base_mv_set_vect
+  !> Function  base_mlv_set_vect
   !! \memberof  psb_i_base_multivect_type
   !! \brief  Set all entries
   !! \param val(:)  The vector to be copied in 
   !!
-  subroutine i_base_mv_set_vect(x,val)
+  subroutine i_base_mlv_set_vect(x,val)
     class(psb_i_base_multivect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr
@@ -1436,7 +1436,7 @@ contains
       x%v = val
     end if
 
-  end subroutine i_base_mv_set_vect
+  end subroutine i_base_mlv_set_vect
 
 end module psb_i_base_multivect_mod
 
diff --git a/base/modules/psb_s_base_vect_mod.f90 b/base/modules/psb_s_base_vect_mod.f90
index 27940dcb..d82f5765 100644
--- a/base/modules/psb_s_base_vect_mod.f90
+++ b/base/modules/psb_s_base_vect_mod.f90
@@ -738,7 +738,7 @@ contains
   !> Function  base_dot_v
   !! \memberof  psb_s_base_vect_type
   !! \brief  Dot product by another base_vector
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y    The other (base_vect) to be multiplied by
   !!
   function s_base_dot_v(n,x,y) result(res)
@@ -772,7 +772,7 @@ contains
   !> Function  base_dot_a
   !! \memberof  psb_s_base_vect_type
   !! \brief  Dot product by a normal array
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y(:) The array to be multiplied by
   !!
   function s_base_dot_a(n,x,y) result(res)
@@ -795,7 +795,7 @@ contains
   !> Function  base_axpby_v
   !! \memberof  psb_s_base_vect_type
   !! \brief AXPBY  by a (base_vect) y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x     The class(base_vect) to be added
   !! \param beta scalar alpha
@@ -823,7 +823,7 @@ contains
   !> Function  base_axpby_a
   !! \memberof  psb_s_base_vect_type
   !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x(:) The array to be added
   !! \param beta scalar alpha
@@ -1339,7 +1339,7 @@ module psb_s_base_multivect_mod
   
   use psb_const_mod
   use psb_error_mod
-
+  use psb_realloc_mod
 
   !> \namespace  psb_base_mod  \class psb_s_base_vect_type
   !! The psb_s_base_vect_type 
@@ -1361,20 +1361,20 @@ module psb_s_base_multivect_mod
     !
     !  Constructors/allocators
     !
-    procedure, pass(x) :: bld_x    => s_base_mv_bld_x
-    procedure, pass(x) :: bld_n    => s_base_mv_bld_n
+    procedure, pass(x) :: bld_x    => s_base_mlv_bld_x
+    procedure, pass(x) :: bld_n    => s_base_mlv_bld_n
     generic, public    :: bld      => bld_x, bld_n
-    procedure, pass(x) :: all      => s_base_mv_all
-    procedure, pass(x) :: mold     => s_base_mv_mold
+    procedure, pass(x) :: all      => s_base_mlv_all
+    procedure, pass(x) :: mold     => s_base_mlv_mold
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
     ! in derived classes. 
     !
-    procedure, pass(x) :: ins      => s_base_mv_ins
-    procedure, pass(x) :: zero     => s_base_mv_zero
-    procedure, pass(x) :: asb      => s_base_mv_asb
-    procedure, pass(x) :: free     => s_base_mv_free
+    procedure, pass(x) :: ins      => s_base_mlv_ins
+    procedure, pass(x) :: zero     => s_base_mlv_zero
+    procedure, pass(x) :: asb      => s_base_mlv_asb
+    procedure, pass(x) :: free     => s_base_mlv_free
     !
     ! Sync: centerpiece of handling of external storage.
     ! Any derived class having extra storage upon sync
@@ -1382,66 +1382,66 @@ module psb_s_base_multivect_mod
     ! external side contain the same data. The base
     ! version is only a placeholder. 
     !
-    procedure, pass(x) :: sync     => s_base_mv_sync
-    procedure, pass(x) :: is_host  => s_base_mv_is_host
-    procedure, pass(x) :: is_dev   => s_base_mv_is_dev
-    procedure, pass(x) :: is_sync  => s_base_mv_is_sync
-    procedure, pass(x) :: set_host => s_base_mv_set_host
-    procedure, pass(x) :: set_dev  => s_base_mv_set_dev
-    procedure, pass(x) :: set_sync => s_base_mv_set_sync
+    procedure, pass(x) :: sync     => s_base_mlv_sync
+    procedure, pass(x) :: is_host  => s_base_mlv_is_host
+    procedure, pass(x) :: is_dev   => s_base_mlv_is_dev
+    procedure, pass(x) :: is_sync  => s_base_mlv_is_sync
+    procedure, pass(x) :: set_host => s_base_mlv_set_host
+    procedure, pass(x) :: set_dev  => s_base_mlv_set_dev
+    procedure, pass(x) :: set_sync => s_base_mlv_set_sync
 
     !
     ! Basic info
-    procedure, pass(x) :: get_nrows => s_base_mv_get_nrows
-    procedure, pass(x) :: get_ncols => s_base_mv_get_ncols
-    procedure, pass(x) :: sizeof    => s_base_mv_sizeof
-    procedure, nopass  :: get_fmt   => s_base_mv_get_fmt
+    procedure, pass(x) :: get_nrows => s_base_mlv_get_nrows
+    procedure, pass(x) :: get_ncols => s_base_mlv_get_ncols
+    procedure, pass(x) :: sizeof    => s_base_mlv_sizeof
+    procedure, nopass  :: get_fmt   => s_base_mlv_get_fmt
     !
     ! Set/get data from/to an external array; also
     ! overload assignment.
     !
-    procedure, pass(x) :: get_vect => s_base_mv_get_vect
-    procedure, pass(x) :: set_scal => s_base_mv_set_scal
-    procedure, pass(x) :: set_vect => s_base_mv_set_vect
+    procedure, pass(x) :: get_vect => s_base_mlv_get_vect
+    procedure, pass(x) :: set_scal => s_base_mlv_set_scal
+    procedure, pass(x) :: set_vect => s_base_mlv_set_vect
     generic, public    :: set      => set_vect, set_scal
 
     !
     ! Dot product and AXPBY
     !
-!!$    procedure, pass(x) :: dot_v    => s_base_mv_dot_v
-!!$    procedure, pass(x) :: dot_a    => s_base_mv_dot_a
-!!$    generic, public    :: dot      => dot_v, dot_a
-!!$    procedure, pass(y) :: axpby_v  => s_base_mv_axpby_v
-!!$    procedure, pass(y) :: axpby_a  => s_base_mv_axpby_a
-!!$    generic, public    :: axpby    => axpby_v, axpby_a
+    procedure, pass(x) :: dot_v    => s_base_mlv_dot_v
+    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_a  => s_base_mlv_axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a
 !!$    !
 !!$    ! Vector by vector multiplication. Need all variants
 !!$    ! to handle multiple requirements from preconditioners
 !!$    !
-!!$    procedure, pass(y) :: mlt_v    => s_base_mv_mlt_v
-!!$    procedure, pass(y) :: mlt_a    => s_base_mv_mlt_a
-!!$    procedure, pass(z) :: mlt_a_2  => s_base_mv_mlt_a_2
-!!$    procedure, pass(z) :: mlt_v_2  => s_base_mv_mlt_v_2
-!!$    procedure, pass(z) :: mlt_va   => s_base_mv_mlt_va
-!!$    procedure, pass(z) :: mlt_av   => s_base_mv_mlt_av
+!!$    procedure, pass(y) :: mlt_v    => s_base_mlv_mlt_v
+!!$    procedure, pass(y) :: mlt_a    => s_base_mlv_mlt_a
+!!$    procedure, pass(z) :: mlt_a_2  => s_base_mlv_mlt_a_2
+!!$    procedure, pass(z) :: mlt_v_2  => s_base_mlv_mlt_v_2
+!!$    procedure, pass(z) :: mlt_va   => s_base_mlv_mlt_va
+!!$    procedure, pass(z) :: mlt_av   => s_base_mlv_mlt_av
 !!$    generic, public    :: mlt      => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
 !!$    !
 !!$    ! Scaling and norms
 !!$    !
-!!$    procedure, pass(x) :: scal     => s_base_mv_scal
-!!$    procedure, pass(x) :: nrm2     => s_base_mv_nrm2
-!!$    procedure, pass(x) :: amax     => s_base_mv_amax
-!!$    procedure, pass(x) :: asum     => s_base_mv_asum
+!!$    procedure, pass(x) :: scal     => s_base_mlv_scal
+!!$    procedure, pass(x) :: nrm2     => s_base_mlv_nrm2
+!!$    procedure, pass(x) :: amax     => s_base_mlv_amax
+!!$    procedure, pass(x) :: asum     => s_base_mlv_asum
 !!$    !
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => s_base_mv_gthab
-!!$    procedure, pass(x) :: gthzv    => s_base_mv_gthzv
-!!$    procedure, pass(x) :: gthzv_x  => s_base_mv_gthzv_x
+!!$    procedure, pass(x) :: gthab    => s_base_mlv_gthab
+!!$    procedure, pass(x) :: gthzv    => s_base_mlv_gthzv
+!!$    procedure, pass(x) :: gthzv_x  => s_base_mlv_gthzv_x
 !!$    generic, public    :: gth      => gthab, gthzv, gthzv_x
-!!$    procedure, pass(y) :: sctb     => s_base_mv_sctb
-!!$    procedure, pass(y) :: sctb_x   => s_base_mv_sctb_x
+!!$    procedure, pass(y) :: sctb     => s_base_mlv_sctb
+!!$    procedure, pass(y) :: sctb_x   => s_base_mlv_sctb_x
 !!$    generic, public    :: sct      => sctb, sctb_x
   end type psb_s_base_multivect_type
 
@@ -1491,7 +1491,7 @@ contains
   !! \brief     Build method from an array
   !!  \param   x(:)  input array to be copied
   !!
-  subroutine s_base_mv_bld_x(x,this)
+  subroutine s_base_mlv_bld_x(x,this)
     use psb_realloc_mod
     real(psb_spk_), intent(in) :: this(:,:)
     class(psb_s_base_multivect_type), intent(inout) :: x
@@ -1499,12 +1499,12 @@ contains
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
     x%v(:,:)  = this(:,:)
 
-  end subroutine s_base_mv_bld_x
+  end subroutine s_base_mlv_bld_x
     
   !
   ! Create with size, but no initialization
@@ -1515,7 +1515,7 @@ contains
   !! \brief     Build method with size (uninitialized data)
   !!  \param    n    size to be allocated. 
   !!
-  subroutine s_base_mv_bld_n(x,m,n)
+  subroutine s_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: m,n
     class(psb_s_base_multivect_type), intent(inout) :: x
@@ -1524,16 +1524,16 @@ contains
     call psb_realloc(m,n,x%v,info)
     call x%asb(m,n,info)
 
-  end subroutine s_base_mv_bld_n
+  end subroutine s_base_mlv_bld_n
   
-  !> Function  base_mv_all:
+  !> Function  base_mlv_all:
   !! \memberof  psb_s_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
   !!  \param    n    size to be allocated. 
   !!  \param info  return code
   !!
-  subroutine s_base_mv_all(m,n, x, info)
+  subroutine s_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1543,15 +1543,15 @@ contains
     
     call psb_realloc(m,n,x%v,info)
     
-  end subroutine s_base_mv_all
+  end subroutine s_base_mlv_all
 
-  !> Function  base_mv_mold:
+  !> Function  base_mlv_mold:
   !! \memberof  psb_s_base_multivect_type
   !! \brief     Mold method: return a variable with the same dynamic type
   !!  \param    y returned variable
   !!  \param info  return code
   !!
-  subroutine s_base_mv_mold(x, y, info)
+  subroutine s_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1561,12 +1561,12 @@ contains
     
     allocate(psb_s_base_multivect_type :: y, stat=info)
 
-  end subroutine s_base_mv_mold
+  end subroutine s_base_mlv_mold
 
   !
   ! Insert a bunch of values at specified positions.
   !
-  !> Function  base_mv_ins:
+  !> Function  base_mlv_ins:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Insert coefficients. 
   !!
@@ -1590,7 +1590,7 @@ contains
   !!  \param info  return code
   !!
   !
-  subroutine s_base_mv_ins(n,irl,val,dupl,x,info)
+  subroutine s_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
     implicit none 
     class(psb_s_base_multivect_type), intent(inout)  :: x
@@ -1642,26 +1642,26 @@ contains
       end select
     end if
     if (info /= 0) then 
-      call psb_errpush(info,'base_mv_vect_ins')
+      call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
 
-  end subroutine s_base_mv_ins
+  end subroutine s_base_mlv_ins
 
   !
-  !> Function  base_mv_zero
+  !> Function  base_mlv_zero
   !! \memberof  psb_s_base_multivect_type
   !! \brief Zero out contents
   !!
   !
-  subroutine s_base_mv_zero(x)
+  subroutine s_base_mlv_zero(x)
     use psi_serial_mod
     implicit none 
     class(psb_s_base_multivect_type), intent(inout)    :: x
     
     if (allocated(x%v)) x%v=szero
 
-  end subroutine s_base_mv_zero
+  end subroutine s_base_mlv_zero
 
   
   !
@@ -1669,7 +1669,7 @@ contains
   ! For derived classes: after this the vector
   ! storage is supposed to be in sync.
   !
-  !> Function  base_mv_asb:
+  !> Function  base_mlv_asb:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
   !!           
@@ -1678,7 +1678,7 @@ contains
   !!
   !
  
-  subroutine s_base_mv_asb(m,n, x, info)
+  subroutine s_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1691,18 +1691,18 @@ contains
     if (info /= 0) &
          & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
 
-  end subroutine s_base_mv_asb
+  end subroutine s_base_mlv_asb
 
 
   !
-  !> Function  base_mv_free:
+  !> Function  base_mlv_free:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Free vector
   !!           
   !!  \param info  return code
   !!
   !
-  subroutine s_base_mv_free(x, info)
+  subroutine s_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1714,7 +1714,7 @@ contains
     if (info /= 0) call & 
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
         
-  end subroutine s_base_mv_free
+  end subroutine s_base_mlv_free
 
   
 
@@ -1723,106 +1723,106 @@ contains
   ! a placeholder.
   ! 
   !
-  !> Function  base_mv_sync:
+  !> Function  base_mlv_sync:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Sync: base version is a no-op.
   !!           
   !
-  subroutine s_base_mv_sync(x)
+  subroutine s_base_mlv_sync(x)
     implicit none 
     class(psb_s_base_multivect_type), intent(inout) :: x
     
-  end subroutine s_base_mv_sync
+  end subroutine s_base_mlv_sync
 
   !
-  !> Function  base_mv_set_host:
+  !> Function  base_mlv_set_host:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Set_host: base version is a no-op.
   !!           
   !
-  subroutine s_base_mv_set_host(x)
+  subroutine s_base_mlv_set_host(x)
     implicit none 
     class(psb_s_base_multivect_type), intent(inout) :: x
     
-  end subroutine s_base_mv_set_host
+  end subroutine s_base_mlv_set_host
 
   !
-  !> Function  base_mv_set_dev:
+  !> Function  base_mlv_set_dev:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
   !!           
   !
-  subroutine s_base_mv_set_dev(x)
+  subroutine s_base_mlv_set_dev(x)
     implicit none 
     class(psb_s_base_multivect_type), intent(inout) :: x
     
-  end subroutine s_base_mv_set_dev
+  end subroutine s_base_mlv_set_dev
 
   !
-  !> Function  base_mv_set_sync:
+  !> Function  base_mlv_set_sync:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
   !!           
   !
-  subroutine s_base_mv_set_sync(x)
+  subroutine s_base_mlv_set_sync(x)
     implicit none 
     class(psb_s_base_multivect_type), intent(inout) :: x
     
-  end subroutine s_base_mv_set_sync
+  end subroutine s_base_mlv_set_sync
 
   !
-  !> Function  base_mv_is_dev:
+  !> Function  base_mlv_is_dev:
   !! \memberof  psb_s_base_multivect_type
   !! \brief Is  vector on external device    .
   !!           
   !
-  function s_base_mv_is_dev(x) result(res)
+  function s_base_mlv_is_dev(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     logical  :: res
   
     res = .false.
-  end function s_base_mv_is_dev
+  end function s_base_mlv_is_dev
   
   !
-  !> Function  base_mv_is_host
+  !> Function  base_mlv_is_host
   !! \memberof  psb_s_base_multivect_type
   !! \brief Is  vector on standard memory    .
   !!           
   !
-  function s_base_mv_is_host(x) result(res)
+  function s_base_mlv_is_host(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function s_base_mv_is_host
+  end function s_base_mlv_is_host
 
   !
-  !> Function  base_mv_is_sync
+  !> Function  base_mlv_is_sync
   !! \memberof  psb_s_base_multivect_type
   !! \brief Is  vector on sync               .
   !!           
   !
-  function s_base_mv_is_sync(x) result(res)
+  function s_base_mlv_is_sync(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function s_base_mv_is_sync
+  end function s_base_mlv_is_sync
 
 
   !
   ! Size info. 
   !
   !
-  !> Function  base_mv_get_nrows
+  !> Function  base_mlv_get_nrows
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Number of entries
   !!           
   !
-  function s_base_mv_get_nrows(x) result(res)
+  function s_base_mlv_get_nrows(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1830,9 +1830,9 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,1)
 
-  end function s_base_mv_get_nrows
+  end function s_base_mlv_get_nrows
 
-  function s_base_mv_get_ncols(x) result(res)
+  function s_base_mlv_get_ncols(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1840,15 +1840,15 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,2)
 
-  end function s_base_mv_get_ncols
+  end function s_base_mlv_get_ncols
   
   !
-  !> Function  base_mv_get_sizeof
+  !> Function  base_mlv_get_sizeof
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Size in bytesa
   !!           
   !
-  function s_base_mv_sizeof(x) result(res)
+  function s_base_mlv_sizeof(x) result(res)
     implicit none 
     class(psb_s_base_multivect_type), intent(in) :: x
     integer(psb_long_int_k_) :: res
@@ -1856,30 +1856,30 @@ contains
     ! Force 8-byte integers.
     res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()
 
-  end function s_base_mv_sizeof
+  end function s_base_mlv_sizeof
 
   !
-  !> Function  base_mv_get_fmt
+  !> Function  base_mlv_get_fmt
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Format
   !!           
   !
-  function s_base_mv_get_fmt() result(res)
+  function s_base_mlv_get_fmt() result(res)
     implicit none 
     character(len=5) :: res
     res = 'BASE'
-  end function s_base_mv_get_fmt
+  end function s_base_mlv_get_fmt
   
 
   !
   !
   !
-  !> Function  base_mv_get_vect
+  !> Function  base_mlv_get_vect
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
   !    
-  function  s_base_mv_get_vect(x) result(res)
+  function  s_base_mlv_get_vect(x) result(res)
     class(psb_s_base_multivect_type), intent(inout) :: x
     real(psb_spk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
@@ -1889,37 +1889,37 @@ contains
     call x%sync()
     allocate(res(m,n),stat=info) 
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
     res(1:m,1:n) = x%v(1:m,1:n)
-  end function s_base_mv_get_vect
+  end function s_base_mlv_get_vect
     
   !
   ! Reset all values 
   !
   !
-  !> Function  base_mv_set_scal
+  !> Function  base_mlv_set_scal
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Set all entries
   !! \param val   The value to set
   !!
-  subroutine s_base_mv_set_scal(x,val)
+  subroutine s_base_mlv_set_scal(x,val)
     class(psb_s_base_multivect_type), intent(inout)  :: x
     real(psb_spk_), intent(in) :: val
         
     integer(psb_ipk_) :: info
     x%v = val
     
-  end subroutine s_base_mv_set_scal
+  end subroutine s_base_mlv_set_scal
 
   !
-  !> Function  base_mv_set_vect
+  !> Function  base_mlv_set_vect
   !! \memberof  psb_s_base_multivect_type
   !! \brief  Set all entries
   !! \param val(:)  The vector to be copied in 
   !!
-  subroutine s_base_mv_set_vect(x,val)
+  subroutine s_base_mlv_set_vect(x,val)
     class(psb_s_base_multivect_type), intent(inout)  :: x
     real(psb_spk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr
@@ -1934,121 +1934,148 @@ contains
       x%v = val
     end if
 
-  end subroutine s_base_mv_set_vect
+  end subroutine s_base_mlv_set_vect
 
-!!$  !
-!!$  ! Dot products 
-!!$  ! 
-!!$  !
-!!$  !> Function  base_mv_dot_v
-!!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief  Dot product by another base_mv_vector
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y    The other (base_mv_vect) to be multiplied by
-!!$  !!
-!!$  function s_base_mv_dot_v(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_s_base_multivect_type), intent(inout) :: x, y
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_spk_)                :: res
-!!$    real(psb_spk_), external      :: ddot
-!!$    
-!!$    res = izero
-!!$    !
-!!$    ! Note: this is the base implementation.
-!!$    !  When we get here, we are sure that X is of
-!!$    !  TYPE psb_s_base_mv_vect.
-!!$    !  If Y is not, throw the burden on it, implicitly
-!!$    !  calling dot_a
-!!$    !
-!!$    select type(yy => y)
-!!$    type is (psb_s_base_multivect_type)
-!!$      res = ddot(n,x%v,1,y%v,1)
-!!$    class default
-!!$      res = y%dot(n,x%v)
-!!$    end select
-!!$
-!!$  end function s_base_mv_dot_v
-!!$
-!!$  !
-!!$  ! Base workhorse is good old BLAS1
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_dot_a
-!!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief  Dot product by a normal array
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y(:) The array to be multiplied by
-!!$  !!
-!!$  function s_base_mv_dot_a(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_s_base_multivect_type), intent(inout) :: x
-!!$    real(psb_spk_), intent(in)    :: y(:)
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    real(psb_spk_)                :: res
-!!$    integer(psb_ipk_), external      :: ddot
-!!$    
-!!$    res = ddot(n,y,1,x%v,1)
-!!$
-!!$  end function s_base_mv_dot_a
+  !
+  ! Dot products 
+  ! 
+  !
+  !> Function  base_mlv_dot_v
+  !! \memberof  psb_s_base_multivect_type
+  !! \brief  Dot product by another base_mlv_vector
+  !! \param n    Number of entries to be considered
+  !! \param y    The other (base_mlv_vect) to be multiplied by
+  !!
+  function s_base_mlv_dot_v(n,x,y) result(res)
+    implicit none 
+    class(psb_s_base_multivect_type), intent(inout) :: x, y
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_spk_), allocatable   :: res(:)
+    real(psb_spk_), external      :: sdot
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    res = szero
+    !
+    ! Note: this is the base implementation.
+    !  When we get here, we are sure that X is of
+    !  TYPE psb_s_base_mlv_vect (or its class does not care). 
+    !  If Y is not, throw the burden on it, implicitly
+    !  calling dot_a
+    !
+    select type(yy => y)
+    type is (psb_s_base_multivect_type)
+      if (y%is_dev()) call y%sync()
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+      allocate(res(nc))
+      do j=1,nc
+        res(j) = sdot(n,x%v(:,j),1,y%v(:,j),1)
+      end do
+    class default
+      res = y%dot(n,x%v)
+    end select
+
+  end function s_base_mlv_dot_v
+
+  !
+  ! Base workhorse is good old BLAS1
+  !
+  !
+  !> Function  base_mlv_dot_a
+  !! \memberof  psb_s_base_multivect_type
+  !! \brief  Dot product by a normal array
+  !! \param n    Number of entries to be considered
+  !! \param y(:) The array to be multiplied by
+  !!
+  function s_base_mlv_dot_a(n,x,y) result(res)
+    implicit none 
+    class(psb_s_base_multivect_type), intent(inout) :: x
+    real(psb_spk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    real(psb_spk_), allocatable     :: res(:)
+    integer(psb_ipk_), external      :: sdot
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    nc = min(psb_size(x%v,2),size(y,2))
+    allocate(res(nc))
+    do j=1,nc
+      res(j) = sdot(n,x%v(:,j),1,y(:,j),1)
+    end do
+
+  end function s_base_mlv_dot_a
     
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_v
-!!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief AXPBY  by a (base_mv_vect) y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x     The class(base_mv_vect) to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine s_base_mv_axpby_v(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    class(psb_s_base_multivect_type), intent(inout)  :: x
-!!$    class(psb_s_base_multivect_type), intent(inout)  :: y
-!!$    real(psb_spk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    select type(xx => x)
-!!$    type is (psb_s_base_multivect_type)
-!!$      call psb_geaxpby(m,alpha,x%v,beta,y%v,info)
-!!$    class default
-!!$      call y%axpby(m,alpha,x%v,beta,info)
-!!$    end select
-!!$
-!!$  end subroutine s_base_mv_axpby_v
-!!$
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_a
-!!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x(:) The array to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine s_base_mv_axpby_a(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    real(psb_spk_), intent(in)        :: x(:)
-!!$    class(psb_s_base_multivect_type), intent(inout)  :: y
-!!$    real(psb_spk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    call psb_geaxpby(m,alpha,x,beta,y%v,info)
-!!$    
-!!$  end subroutine s_base_mv_axpby_a
+  !
+  ! 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=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_mlv_vect) to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine s_base_mlv_axpby_v(m,alpha, x, beta, y, info, n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    class(psb_s_base_multivect_type), intent(inout)  :: x
+    class(psb_s_base_multivect_type), intent(inout)  :: y
+    real(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+    end if
+    select type(xx => x)
+    type is (psb_s_base_multivect_type)
+      call psb_geaxpby(m,nc,alpha,x%v,beta,y%v,info)
+    class default
+      call y%axpby(m,alpha,x%v,beta,info,n=n)
+    end select
+
+  end subroutine s_base_mlv_axpby_v
+
+  !
+  ! AXPBY is invoked via Y, hence the structure below. 
+  !
+  !
+  !> Function  base_mlv_axpby_a
+  !! \memberof  psb_s_base_multivect_type
+  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x(:) The array to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine s_base_mlv_axpby_a(m,alpha, x, beta, y, info,n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    real(psb_spk_), intent(in)        :: x(:,:)
+    class(psb_s_base_multivect_type), intent(inout)  :: y
+    real(psb_spk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(size(x,2),psb_size(y%v,2))
+    end if
+    
+    call psb_geaxpby(m,nc,alpha,x,beta,y%v,info)
+    
+  end subroutine s_base_mlv_axpby_a
 
   
 !!$  !
@@ -2060,13 +2087,13 @@ contains
 !!$  !  of the involved entities
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief Vector entry-by-entry multiply  by a base_mv_vect array y=x*y
-!!$  !! \param x   The class(base_mv_vect) to be multiplied by
+!!$  !! \brief Vector entry-by-entry multiply  by a base_mlv_vect array y=x*y
+!!$  !! \param x   The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine s_base_mv_mlt_v(x, y, info)
+!!$  subroutine s_base_mlv_mlt_v(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_s_base_multivect_type), intent(inout)  :: x
@@ -2085,16 +2112,16 @@ contains
 !!$      call y%mlt(x%v,info)
 !!$    end select
 !!$
-!!$  end subroutine s_base_mv_mlt_v
+!!$  end subroutine s_base_mlv_mlt_v
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief Vector entry-by-entry multiply  by a normal array y=x*y
 !!$  !! \param x(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine s_base_mv_mlt_a(x, y, info)
+!!$  subroutine s_base_mlv_mlt_a(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_spk_), intent(in)        :: x(:)
@@ -2108,11 +2135,11 @@ contains
 !!$      y%v(i) = y%v(i)*x(i)
 !!$    end do
 !!$    
-!!$  end subroutine s_base_mv_mlt_a
+!!$  end subroutine s_base_mlv_mlt_a
 !!$
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a_2
+!!$  !> Function  base_mlv_mlt_a_2
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by normal arrays
 !!$  !!        z=beta*z+alpha*x*y
@@ -2122,7 +2149,7 @@ contains
 !!$  !! \param y(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine s_base_mv_mlt_a_2(alpha,x,y,beta,z,info)
+!!$  subroutine s_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_spk_), intent(in)        :: alpha,beta
@@ -2187,20 +2214,20 @@ contains
 !!$        end if
 !!$      end if
 !!$    end if
-!!$  end subroutine s_base_mv_mlt_a_2
+!!$  end subroutine s_base_mlv_mlt_a_2
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_v_2
+!!$  !> Function  base_mlv_mlt_v_2
 !!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mv_vect)
+!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mlv_vect)
 !!$  !!        z=beta*z+alpha*x*y
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  !! \param x  The class(base_mv_vect) to be multiplied b
-!!$  !! \param y  The class(base_mv_vect) to be multiplied by
+!!$  !! \param x  The class(base_mlv_vect) to be multiplied b
+!!$  !! \param y  The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine s_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
+!!$  subroutine s_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
 !!$    use psi_serial_mod
 !!$    use psb_string_mod
 !!$    implicit none 
@@ -2227,9 +2254,9 @@ contains
 !!$      if (conjgx_) x%v=(x%v)
 !!$      if (conjgy_) y%v=(y%v)
 !!$    end if
-!!$  end subroutine s_base_mv_mlt_v_2
+!!$  end subroutine s_base_mlv_mlt_v_2
 !!$
-!!$  subroutine s_base_mv_mlt_av(alpha,x,y,beta,z,info)
+!!$  subroutine s_base_mlv_mlt_av(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_spk_), intent(in)        :: alpha,beta
@@ -2243,9 +2270,9 @@ contains
 !!$    
 !!$    call z%mlt(alpha,x,y%v,beta,info)
 !!$
-!!$  end subroutine s_base_mv_mlt_av
+!!$  end subroutine s_base_mlv_mlt_av
 !!$
-!!$  subroutine s_base_mv_mlt_va(alpha,x,y,beta,z,info)
+!!$  subroutine s_base_mlv_mlt_va(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    real(psb_spk_), intent(in)        :: alpha,beta
@@ -2259,18 +2286,18 @@ contains
 !!$    
 !!$    call z%mlt(alpha,y,x,beta,info)
 !!$
-!!$  end subroutine s_base_mv_mlt_va
+!!$  end subroutine s_base_mlv_mlt_va
 !!$
 !!$
 !!$  !
 !!$  ! Simple scaling 
 !!$  !
-!!$  !> Function  base_mv_scal
+!!$  !> Function  base_mlv_scal
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief Scale all entries  x = alpha*x
 !!$  !! \param alpha   The multiplier
 !!$  !!
-!!$  subroutine s_base_mv_scal(alpha, x)
+!!$  subroutine s_base_mlv_scal(alpha, x)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_s_base_multivect_type), intent(inout)  :: x
@@ -2278,16 +2305,16 @@ contains
 !!$    
 !!$    if (allocated(x%v)) x%v = alpha*x%v
 !!$
-!!$  end subroutine s_base_mv_scal
+!!$  end subroutine s_base_mlv_scal
 !!$  
 !!$  !
 !!$  ! Norms 1, 2 and infinity
 !!$  !
-!!$  !> Function  base_mv_nrm2
+!!$  !> Function  base_mlv_nrm2
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief 2-norm |x(1:n)|_2
 !!$  !! \param n  how many entries to consider
-!!$  function s_base_mv_nrm2(n,x) result(res)
+!!$  function s_base_mlv_nrm2(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_s_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2296,14 +2323,14 @@ contains
 !!$    
 !!$    res =  dnrm2(n,x%v,1)
 !!$
-!!$  end function s_base_mv_nrm2
+!!$  end function s_base_mlv_nrm2
 !!$  
 !!$  !
-!!$  !> Function  base_mv_amax
+!!$  !> Function  base_mlv_amax
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief infinity-norm |x(1:n)|_\infty
 !!$  !! \param n  how many entries to consider
-!!$  function s_base_mv_amax(n,x) result(res)
+!!$  function s_base_mlv_amax(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_s_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2311,14 +2338,14 @@ contains
 !!$    
 !!$    res =  maxval(abs(x%v(1:n)))
 !!$
-!!$  end function s_base_mv_amax
+!!$  end function s_base_mlv_amax
 !!$
 !!$  !
-!!$  !> Function  base_mv_asum
+!!$  !> Function  base_mlv_asum
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief 1-norm |x(1:n)|_1
 !!$  !! \param n  how many entries to consider
-!!$  function s_base_mv_asum(n,x) result(res)
+!!$  function s_base_mlv_asum(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_s_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2326,14 +2353,14 @@ contains
 !!$    
 !!$    res =  sum(abs(x%v(1:n)))
 !!$
-!!$  end function s_base_mv_asum
+!!$  end function s_base_mlv_asum
 !!$  
 !!$  
 !!$  !
 !!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_gthab
+!!$  !> Function  base_mlv_gthab
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief gather into an array
 !!$  !!    Y = beta * Y + alpha * X(IDX(:))
@@ -2341,7 +2368,7 @@ contains
 !!$  !! \param idx(:) indices
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  subroutine s_base_mv_gthab(n,idx,alpha,x,beta,y)
+!!$  subroutine s_base_mlv_gthab(n,idx,alpha,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_spk_) :: alpha, beta, y(:)
@@ -2350,17 +2377,17 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,alpha,x%v,beta,y)
 !!$
-!!$  end subroutine s_base_mv_gthab
+!!$  end subroutine s_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine s_base_mv_gthzv_x(i,n,idx,x,y)
+!!$  subroutine s_base_mlv_gthzv_x(i,n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i,n
 !!$    class(psb_s_base_multivect_type) :: idx
@@ -2369,18 +2396,18 @@ contains
 !!$    
 !!$    call x%gth(n,idx%v(i:),y)
 !!$
-!!$  end subroutine s_base_mv_gthzv_x
+!!$  end subroutine s_base_mlv_gthzv_x
 !!$
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_s_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine s_base_mv_gthzv(n,idx,x,y)
+!!$  subroutine s_base_mlv_gthzv(n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_spk_) ::  y(:)
@@ -2389,22 +2416,22 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,x%v,y)
 !!$
-!!$  end subroutine s_base_mv_gthzv
+!!$  end subroutine s_base_mlv_gthzv
 !!$
 !!$  !
 !!$  ! Scatter: 
 !!$  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
 !!$  ! 
 !!$  !
-!!$  !> Function  base_mv_sctb
+!!$  !> Function  base_mlv_sctb
 !!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief scatter into a class(base_mv_vect)
+!!$  !! \brief scatter into a class(base_mlv_vect)
 !!$  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
 !!$  !! \param beta
 !!$  !! \param x(:) 
-!!$  subroutine s_base_mv_sctb(n,idx,x,beta,y)
+!!$  subroutine s_base_mlv_sctb(n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    real(psb_spk_) :: beta, x(:)
@@ -2414,9 +2441,9 @@ contains
 !!$    call psi_sct(n,idx,x,beta,y%v)
 !!$    call y%set_host()
 !!$
-!!$  end subroutine s_base_mv_sctb
+!!$  end subroutine s_base_mlv_sctb
 !!$
-!!$  subroutine s_base_mv_sctb_x(i,n,idx,x,beta,y)
+!!$  subroutine s_base_mlv_sctb_x(i,n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i, n
 !!$    class(psb_s_base_multivect_type) :: idx
@@ -2425,6 +2452,6 @@ contains
 !!$    
 !!$    call y%sct(n,idx%v(i:),x,beta)
 !!$
-!!$  end subroutine s_base_mv_sctb_x
+!!$  end subroutine s_base_mlv_sctb_x
 end module psb_s_base_multivect_mod
 
diff --git a/base/modules/psb_z_base_vect_mod.f90 b/base/modules/psb_z_base_vect_mod.f90
index 19313199..3717135b 100644
--- a/base/modules/psb_z_base_vect_mod.f90
+++ b/base/modules/psb_z_base_vect_mod.f90
@@ -738,7 +738,7 @@ contains
   !> Function  base_dot_v
   !! \memberof  psb_z_base_vect_type
   !! \brief  Dot product by another base_vector
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y    The other (base_vect) to be multiplied by
   !!
   function z_base_dot_v(n,x,y) result(res)
@@ -772,7 +772,7 @@ contains
   !> Function  base_dot_a
   !! \memberof  psb_z_base_vect_type
   !! \brief  Dot product by a normal array
-  !! \param n    Number of entries to be considere
+  !! \param n    Number of entries to be considered
   !! \param y(:) The array to be multiplied by
   !!
   function z_base_dot_a(n,x,y) result(res)
@@ -795,7 +795,7 @@ contains
   !> Function  base_axpby_v
   !! \memberof  psb_z_base_vect_type
   !! \brief AXPBY  by a (base_vect) y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x     The class(base_vect) to be added
   !! \param beta scalar alpha
@@ -823,7 +823,7 @@ contains
   !> Function  base_axpby_a
   !! \memberof  psb_z_base_vect_type
   !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-  !! \param m    Number of entries to be considere
+  !! \param m    Number of entries to be considered
   !! \param alpha scalar alpha
   !! \param x(:) The array to be added
   !! \param beta scalar alpha
@@ -1339,7 +1339,7 @@ module psb_z_base_multivect_mod
   
   use psb_const_mod
   use psb_error_mod
-
+  use psb_realloc_mod
 
   !> \namespace  psb_base_mod  \class psb_z_base_vect_type
   !! The psb_z_base_vect_type 
@@ -1361,20 +1361,20 @@ module psb_z_base_multivect_mod
     !
     !  Constructors/allocators
     !
-    procedure, pass(x) :: bld_x    => z_base_mv_bld_x
-    procedure, pass(x) :: bld_n    => z_base_mv_bld_n
+    procedure, pass(x) :: bld_x    => z_base_mlv_bld_x
+    procedure, pass(x) :: bld_n    => z_base_mlv_bld_n
     generic, public    :: bld      => bld_x, bld_n
-    procedure, pass(x) :: all      => z_base_mv_all
-    procedure, pass(x) :: mold     => z_base_mv_mold
+    procedure, pass(x) :: all      => z_base_mlv_all
+    procedure, pass(x) :: mold     => z_base_mlv_mold
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
     ! in derived classes. 
     !
-    procedure, pass(x) :: ins      => z_base_mv_ins
-    procedure, pass(x) :: zero     => z_base_mv_zero
-    procedure, pass(x) :: asb      => z_base_mv_asb
-    procedure, pass(x) :: free     => z_base_mv_free
+    procedure, pass(x) :: ins      => z_base_mlv_ins
+    procedure, pass(x) :: zero     => z_base_mlv_zero
+    procedure, pass(x) :: asb      => z_base_mlv_asb
+    procedure, pass(x) :: free     => z_base_mlv_free
     !
     ! Sync: centerpiece of handling of external storage.
     ! Any derived class having extra storage upon sync
@@ -1382,66 +1382,66 @@ module psb_z_base_multivect_mod
     ! external side contain the same data. The base
     ! version is only a placeholder. 
     !
-    procedure, pass(x) :: sync     => z_base_mv_sync
-    procedure, pass(x) :: is_host  => z_base_mv_is_host
-    procedure, pass(x) :: is_dev   => z_base_mv_is_dev
-    procedure, pass(x) :: is_sync  => z_base_mv_is_sync
-    procedure, pass(x) :: set_host => z_base_mv_set_host
-    procedure, pass(x) :: set_dev  => z_base_mv_set_dev
-    procedure, pass(x) :: set_sync => z_base_mv_set_sync
+    procedure, pass(x) :: sync     => z_base_mlv_sync
+    procedure, pass(x) :: is_host  => z_base_mlv_is_host
+    procedure, pass(x) :: is_dev   => z_base_mlv_is_dev
+    procedure, pass(x) :: is_sync  => z_base_mlv_is_sync
+    procedure, pass(x) :: set_host => z_base_mlv_set_host
+    procedure, pass(x) :: set_dev  => z_base_mlv_set_dev
+    procedure, pass(x) :: set_sync => z_base_mlv_set_sync
 
     !
     ! Basic info
-    procedure, pass(x) :: get_nrows => z_base_mv_get_nrows
-    procedure, pass(x) :: get_ncols => z_base_mv_get_ncols
-    procedure, pass(x) :: sizeof    => z_base_mv_sizeof
-    procedure, nopass  :: get_fmt   => z_base_mv_get_fmt
+    procedure, pass(x) :: get_nrows => z_base_mlv_get_nrows
+    procedure, pass(x) :: get_ncols => z_base_mlv_get_ncols
+    procedure, pass(x) :: sizeof    => z_base_mlv_sizeof
+    procedure, nopass  :: get_fmt   => z_base_mlv_get_fmt
     !
     ! Set/get data from/to an external array; also
     ! overload assignment.
     !
-    procedure, pass(x) :: get_vect => z_base_mv_get_vect
-    procedure, pass(x) :: set_scal => z_base_mv_set_scal
-    procedure, pass(x) :: set_vect => z_base_mv_set_vect
+    procedure, pass(x) :: get_vect => z_base_mlv_get_vect
+    procedure, pass(x) :: set_scal => z_base_mlv_set_scal
+    procedure, pass(x) :: set_vect => z_base_mlv_set_vect
     generic, public    :: set      => set_vect, set_scal
 
     !
     ! Dot product and AXPBY
     !
-!!$    procedure, pass(x) :: dot_v    => z_base_mv_dot_v
-!!$    procedure, pass(x) :: dot_a    => z_base_mv_dot_a
-!!$    generic, public    :: dot      => dot_v, dot_a
-!!$    procedure, pass(y) :: axpby_v  => z_base_mv_axpby_v
-!!$    procedure, pass(y) :: axpby_a  => z_base_mv_axpby_a
-!!$    generic, public    :: axpby    => axpby_v, axpby_a
+    procedure, pass(x) :: dot_v    => z_base_mlv_dot_v
+    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_a  => z_base_mlv_axpby_a
+    generic, public    :: axpby    => axpby_v, axpby_a
 !!$    !
 !!$    ! Vector by vector multiplication. Need all variants
 !!$    ! to handle multiple requirements from preconditioners
 !!$    !
-!!$    procedure, pass(y) :: mlt_v    => z_base_mv_mlt_v
-!!$    procedure, pass(y) :: mlt_a    => z_base_mv_mlt_a
-!!$    procedure, pass(z) :: mlt_a_2  => z_base_mv_mlt_a_2
-!!$    procedure, pass(z) :: mlt_v_2  => z_base_mv_mlt_v_2
-!!$    procedure, pass(z) :: mlt_va   => z_base_mv_mlt_va
-!!$    procedure, pass(z) :: mlt_av   => z_base_mv_mlt_av
+!!$    procedure, pass(y) :: mlt_v    => z_base_mlv_mlt_v
+!!$    procedure, pass(y) :: mlt_a    => z_base_mlv_mlt_a
+!!$    procedure, pass(z) :: mlt_a_2  => z_base_mlv_mlt_a_2
+!!$    procedure, pass(z) :: mlt_v_2  => z_base_mlv_mlt_v_2
+!!$    procedure, pass(z) :: mlt_va   => z_base_mlv_mlt_va
+!!$    procedure, pass(z) :: mlt_av   => z_base_mlv_mlt_av
 !!$    generic, public    :: mlt      => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
 !!$    !
 !!$    ! Scaling and norms
 !!$    !
-!!$    procedure, pass(x) :: scal     => z_base_mv_scal
-!!$    procedure, pass(x) :: nrm2     => z_base_mv_nrm2
-!!$    procedure, pass(x) :: amax     => z_base_mv_amax
-!!$    procedure, pass(x) :: asum     => z_base_mv_asum
+!!$    procedure, pass(x) :: scal     => z_base_mlv_scal
+!!$    procedure, pass(x) :: nrm2     => z_base_mlv_nrm2
+!!$    procedure, pass(x) :: amax     => z_base_mlv_amax
+!!$    procedure, pass(x) :: asum     => z_base_mlv_asum
 !!$    !
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => z_base_mv_gthab
-!!$    procedure, pass(x) :: gthzv    => z_base_mv_gthzv
-!!$    procedure, pass(x) :: gthzv_x  => z_base_mv_gthzv_x
+!!$    procedure, pass(x) :: gthab    => z_base_mlv_gthab
+!!$    procedure, pass(x) :: gthzv    => z_base_mlv_gthzv
+!!$    procedure, pass(x) :: gthzv_x  => z_base_mlv_gthzv_x
 !!$    generic, public    :: gth      => gthab, gthzv, gthzv_x
-!!$    procedure, pass(y) :: sctb     => z_base_mv_sctb
-!!$    procedure, pass(y) :: sctb_x   => z_base_mv_sctb_x
+!!$    procedure, pass(y) :: sctb     => z_base_mlv_sctb
+!!$    procedure, pass(y) :: sctb_x   => z_base_mlv_sctb_x
 !!$    generic, public    :: sct      => sctb, sctb_x
   end type psb_z_base_multivect_type
 
@@ -1491,7 +1491,7 @@ contains
   !! \brief     Build method from an array
   !!  \param   x(:)  input array to be copied
   !!
-  subroutine z_base_mv_bld_x(x,this)
+  subroutine z_base_mlv_bld_x(x,this)
     use psb_realloc_mod
     complex(psb_dpk_), intent(in) :: this(:,:)
     class(psb_z_base_multivect_type), intent(inout) :: x
@@ -1499,12 +1499,12 @@ contains
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
     x%v(:,:)  = this(:,:)
 
-  end subroutine z_base_mv_bld_x
+  end subroutine z_base_mlv_bld_x
     
   !
   ! Create with size, but no initialization
@@ -1515,7 +1515,7 @@ contains
   !! \brief     Build method with size (uninitialized data)
   !!  \param    n    size to be allocated. 
   !!
-  subroutine z_base_mv_bld_n(x,m,n)
+  subroutine z_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
     integer(psb_ipk_), intent(in) :: m,n
     class(psb_z_base_multivect_type), intent(inout) :: x
@@ -1524,16 +1524,16 @@ contains
     call psb_realloc(m,n,x%v,info)
     call x%asb(m,n,info)
 
-  end subroutine z_base_mv_bld_n
+  end subroutine z_base_mlv_bld_n
   
-  !> Function  base_mv_all:
+  !> Function  base_mlv_all:
   !! \memberof  psb_z_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
   !!  \param    n    size to be allocated. 
   !!  \param info  return code
   !!
-  subroutine z_base_mv_all(m,n, x, info)
+  subroutine z_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1543,15 +1543,15 @@ contains
     
     call psb_realloc(m,n,x%v,info)
     
-  end subroutine z_base_mv_all
+  end subroutine z_base_mlv_all
 
-  !> Function  base_mv_mold:
+  !> Function  base_mlv_mold:
   !! \memberof  psb_z_base_multivect_type
   !! \brief     Mold method: return a variable with the same dynamic type
   !!  \param    y returned variable
   !!  \param info  return code
   !!
-  subroutine z_base_mv_mold(x, y, info)
+  subroutine z_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1561,12 +1561,12 @@ contains
     
     allocate(psb_z_base_multivect_type :: y, stat=info)
 
-  end subroutine z_base_mv_mold
+  end subroutine z_base_mlv_mold
 
   !
   ! Insert a bunch of values at specified positions.
   !
-  !> Function  base_mv_ins:
+  !> Function  base_mlv_ins:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Insert coefficients. 
   !!
@@ -1590,7 +1590,7 @@ contains
   !!  \param info  return code
   !!
   !
-  subroutine z_base_mv_ins(n,irl,val,dupl,x,info)
+  subroutine z_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
     implicit none 
     class(psb_z_base_multivect_type), intent(inout)  :: x
@@ -1642,26 +1642,26 @@ contains
       end select
     end if
     if (info /= 0) then 
-      call psb_errpush(info,'base_mv_vect_ins')
+      call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
 
-  end subroutine z_base_mv_ins
+  end subroutine z_base_mlv_ins
 
   !
-  !> Function  base_mv_zero
+  !> Function  base_mlv_zero
   !! \memberof  psb_z_base_multivect_type
   !! \brief Zero out contents
   !!
   !
-  subroutine z_base_mv_zero(x)
+  subroutine z_base_mlv_zero(x)
     use psi_serial_mod
     implicit none 
     class(psb_z_base_multivect_type), intent(inout)    :: x
     
     if (allocated(x%v)) x%v=zzero
 
-  end subroutine z_base_mv_zero
+  end subroutine z_base_mlv_zero
 
   
   !
@@ -1669,7 +1669,7 @@ contains
   ! For derived classes: after this the vector
   ! storage is supposed to be in sync.
   !
-  !> Function  base_mv_asb:
+  !> Function  base_mlv_asb:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
   !!           
@@ -1678,7 +1678,7 @@ contains
   !!
   !
  
-  subroutine z_base_mv_asb(m,n, x, info)
+  subroutine z_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1691,18 +1691,18 @@ contains
     if (info /= 0) &
          & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb')
 
-  end subroutine z_base_mv_asb
+  end subroutine z_base_mlv_asb
 
 
   !
-  !> Function  base_mv_free:
+  !> Function  base_mlv_free:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Free vector
   !!           
   !!  \param info  return code
   !!
   !
-  subroutine z_base_mv_free(x, info)
+  subroutine z_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
     implicit none 
@@ -1714,7 +1714,7 @@ contains
     if (info /= 0) call & 
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
         
-  end subroutine z_base_mv_free
+  end subroutine z_base_mlv_free
 
   
 
@@ -1723,106 +1723,106 @@ contains
   ! a placeholder.
   ! 
   !
-  !> Function  base_mv_sync:
+  !> Function  base_mlv_sync:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Sync: base version is a no-op.
   !!           
   !
-  subroutine z_base_mv_sync(x)
+  subroutine z_base_mlv_sync(x)
     implicit none 
     class(psb_z_base_multivect_type), intent(inout) :: x
     
-  end subroutine z_base_mv_sync
+  end subroutine z_base_mlv_sync
 
   !
-  !> Function  base_mv_set_host:
+  !> Function  base_mlv_set_host:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Set_host: base version is a no-op.
   !!           
   !
-  subroutine z_base_mv_set_host(x)
+  subroutine z_base_mlv_set_host(x)
     implicit none 
     class(psb_z_base_multivect_type), intent(inout) :: x
     
-  end subroutine z_base_mv_set_host
+  end subroutine z_base_mlv_set_host
 
   !
-  !> Function  base_mv_set_dev:
+  !> Function  base_mlv_set_dev:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
   !!           
   !
-  subroutine z_base_mv_set_dev(x)
+  subroutine z_base_mlv_set_dev(x)
     implicit none 
     class(psb_z_base_multivect_type), intent(inout) :: x
     
-  end subroutine z_base_mv_set_dev
+  end subroutine z_base_mlv_set_dev
 
   !
-  !> Function  base_mv_set_sync:
+  !> Function  base_mlv_set_sync:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
   !!           
   !
-  subroutine z_base_mv_set_sync(x)
+  subroutine z_base_mlv_set_sync(x)
     implicit none 
     class(psb_z_base_multivect_type), intent(inout) :: x
     
-  end subroutine z_base_mv_set_sync
+  end subroutine z_base_mlv_set_sync
 
   !
-  !> Function  base_mv_is_dev:
+  !> Function  base_mlv_is_dev:
   !! \memberof  psb_z_base_multivect_type
   !! \brief Is  vector on external device    .
   !!           
   !
-  function z_base_mv_is_dev(x) result(res)
+  function z_base_mlv_is_dev(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     logical  :: res
   
     res = .false.
-  end function z_base_mv_is_dev
+  end function z_base_mlv_is_dev
   
   !
-  !> Function  base_mv_is_host
+  !> Function  base_mlv_is_host
   !! \memberof  psb_z_base_multivect_type
   !! \brief Is  vector on standard memory    .
   !!           
   !
-  function z_base_mv_is_host(x) result(res)
+  function z_base_mlv_is_host(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function z_base_mv_is_host
+  end function z_base_mlv_is_host
 
   !
-  !> Function  base_mv_is_sync
+  !> Function  base_mlv_is_sync
   !! \memberof  psb_z_base_multivect_type
   !! \brief Is  vector on sync               .
   !!           
   !
-  function z_base_mv_is_sync(x) result(res)
+  function z_base_mlv_is_sync(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     logical  :: res
 
     res = .true.
-  end function z_base_mv_is_sync
+  end function z_base_mlv_is_sync
 
 
   !
   ! Size info. 
   !
   !
-  !> Function  base_mv_get_nrows
+  !> Function  base_mlv_get_nrows
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Number of entries
   !!           
   !
-  function z_base_mv_get_nrows(x) result(res)
+  function z_base_mlv_get_nrows(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1830,9 +1830,9 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,1)
 
-  end function z_base_mv_get_nrows
+  end function z_base_mlv_get_nrows
 
-  function z_base_mv_get_ncols(x) result(res)
+  function z_base_mlv_get_ncols(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
@@ -1840,15 +1840,15 @@ contains
     res = 0
     if (allocated(x%v)) res = size(x%v,2)
 
-  end function z_base_mv_get_ncols
+  end function z_base_mlv_get_ncols
   
   !
-  !> Function  base_mv_get_sizeof
+  !> Function  base_mlv_get_sizeof
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Size in bytesa
   !!           
   !
-  function z_base_mv_sizeof(x) result(res)
+  function z_base_mlv_sizeof(x) result(res)
     implicit none 
     class(psb_z_base_multivect_type), intent(in) :: x
     integer(psb_long_int_k_) :: res
@@ -1856,30 +1856,30 @@ contains
     ! Force 8-byte integers.
     res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols()
 
-  end function z_base_mv_sizeof
+  end function z_base_mlv_sizeof
 
   !
-  !> Function  base_mv_get_fmt
+  !> Function  base_mlv_get_fmt
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Format
   !!           
   !
-  function z_base_mv_get_fmt() result(res)
+  function z_base_mlv_get_fmt() result(res)
     implicit none 
     character(len=5) :: res
     res = 'BASE'
-  end function z_base_mv_get_fmt
+  end function z_base_mlv_get_fmt
   
 
   !
   !
   !
-  !> Function  base_mv_get_vect
+  !> Function  base_mlv_get_vect
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
   !    
-  function  z_base_mv_get_vect(x) result(res)
+  function  z_base_mlv_get_vect(x) result(res)
     class(psb_z_base_multivect_type), intent(inout) :: x
     complex(psb_dpk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
@@ -1889,37 +1889,37 @@ contains
     call x%sync()
     allocate(res(m,n),stat=info) 
     if (info /= 0) then 
-      call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect')
+      call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
     res(1:m,1:n) = x%v(1:m,1:n)
-  end function z_base_mv_get_vect
+  end function z_base_mlv_get_vect
     
   !
   ! Reset all values 
   !
   !
-  !> Function  base_mv_set_scal
+  !> Function  base_mlv_set_scal
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Set all entries
   !! \param val   The value to set
   !!
-  subroutine z_base_mv_set_scal(x,val)
+  subroutine z_base_mlv_set_scal(x,val)
     class(psb_z_base_multivect_type), intent(inout)  :: x
     complex(psb_dpk_), intent(in) :: val
         
     integer(psb_ipk_) :: info
     x%v = val
     
-  end subroutine z_base_mv_set_scal
+  end subroutine z_base_mlv_set_scal
 
   !
-  !> Function  base_mv_set_vect
+  !> Function  base_mlv_set_vect
   !! \memberof  psb_z_base_multivect_type
   !! \brief  Set all entries
   !! \param val(:)  The vector to be copied in 
   !!
-  subroutine z_base_mv_set_vect(x,val)
+  subroutine z_base_mlv_set_vect(x,val)
     class(psb_z_base_multivect_type), intent(inout)  :: x
     complex(psb_dpk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr
@@ -1934,121 +1934,148 @@ contains
       x%v = val
     end if
 
-  end subroutine z_base_mv_set_vect
+  end subroutine z_base_mlv_set_vect
 
-!!$  !
-!!$  ! Dot products 
-!!$  ! 
-!!$  !
-!!$  !> Function  base_mv_dot_v
-!!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief  Dot product by another base_mv_vector
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y    The other (base_mv_vect) to be multiplied by
-!!$  !!
-!!$  function z_base_mv_dot_v(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_z_base_multivect_type), intent(inout) :: x, y
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    complex(psb_dpk_)                :: res
-!!$    complex(psb_dpk_), external      :: ddot
-!!$    
-!!$    res = izero
-!!$    !
-!!$    ! Note: this is the base implementation.
-!!$    !  When we get here, we are sure that X is of
-!!$    !  TYPE psb_z_base_mv_vect.
-!!$    !  If Y is not, throw the burden on it, implicitly
-!!$    !  calling dot_a
-!!$    !
-!!$    select type(yy => y)
-!!$    type is (psb_z_base_multivect_type)
-!!$      res = ddot(n,x%v,1,y%v,1)
-!!$    class default
-!!$      res = y%dot(n,x%v)
-!!$    end select
-!!$
-!!$  end function z_base_mv_dot_v
-!!$
-!!$  !
-!!$  ! Base workhorse is good old BLAS1
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_dot_a
-!!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief  Dot product by a normal array
-!!$  !! \param n    Number of entries to be considere
-!!$  !! \param y(:) The array to be multiplied by
-!!$  !!
-!!$  function z_base_mv_dot_a(n,x,y) result(res)
-!!$    implicit none 
-!!$    class(psb_z_base_multivect_type), intent(inout) :: x
-!!$    complex(psb_dpk_), intent(in)    :: y(:)
-!!$    integer(psb_ipk_), intent(in)           :: n
-!!$    complex(psb_dpk_)                :: res
-!!$    integer(psb_ipk_), external      :: ddot
-!!$    
-!!$    res = ddot(n,y,1,x%v,1)
-!!$
-!!$  end function z_base_mv_dot_a
+  !
+  ! Dot products 
+  ! 
+  !
+  !> Function  base_mlv_dot_v
+  !! \memberof  psb_z_base_multivect_type
+  !! \brief  Dot product by another base_mlv_vector
+  !! \param n    Number of entries to be considered
+  !! \param y    The other (base_mlv_vect) to be multiplied by
+  !!
+  function z_base_mlv_dot_v(n,x,y) result(res)
+    implicit none 
+    class(psb_z_base_multivect_type), intent(inout) :: x, y
+    integer(psb_ipk_), intent(in)           :: n
+    complex(psb_dpk_), allocatable   :: res(:)
+    complex(psb_dpk_), external      :: zdotc
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    res = zzero
+    !
+    ! Note: this is the base implementation.
+    !  When we get here, we are sure that X is of
+    !  TYPE psb_z_base_mlv_vect (or its class does not care). 
+    !  If Y is not, throw the burden on it, implicitly
+    !  calling dot_a
+    !
+    select type(yy => y)
+    type is (psb_z_base_multivect_type)
+      if (y%is_dev()) call y%sync()
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+      allocate(res(nc))
+      do j=1,nc
+        res(j) = zdotc(n,x%v(:,j),1,y%v(:,j),1)
+      end do
+    class default
+      res = y%dot(n,x%v)
+    end select
+
+  end function z_base_mlv_dot_v
+
+  !
+  ! Base workhorse is good old BLAS1
+  !
+  !
+  !> Function  base_mlv_dot_a
+  !! \memberof  psb_z_base_multivect_type
+  !! \brief  Dot product by a normal array
+  !! \param n    Number of entries to be considered
+  !! \param y(:) The array to be multiplied by
+  !!
+  function z_base_mlv_dot_a(n,x,y) result(res)
+    implicit none 
+    class(psb_z_base_multivect_type), intent(inout) :: x
+    complex(psb_dpk_), intent(in)    :: y(:,:)
+    integer(psb_ipk_), intent(in)           :: n
+    complex(psb_dpk_), allocatable     :: res(:)
+    integer(psb_ipk_), external      :: zdotc
+    integer(psb_ipk_) :: j,nc
+
+    if (x%is_dev()) call x%sync()
+    nc = min(psb_size(x%v,2),size(y,2))
+    allocate(res(nc))
+    do j=1,nc
+      res(j) = zdotc(n,x%v(:,j),1,y(:,j),1)
+    end do
+
+  end function z_base_mlv_dot_a
     
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_v
-!!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief AXPBY  by a (base_mv_vect) y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x     The class(base_mv_vect) to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine z_base_mv_axpby_v(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    class(psb_z_base_multivect_type), intent(inout)  :: x
-!!$    class(psb_z_base_multivect_type), intent(inout)  :: y
-!!$    complex(psb_dpk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    select type(xx => x)
-!!$    type is (psb_z_base_multivect_type)
-!!$      call psb_geaxpby(m,alpha,x%v,beta,y%v,info)
-!!$    class default
-!!$      call y%axpby(m,alpha,x%v,beta,info)
-!!$    end select
-!!$
-!!$  end subroutine z_base_mv_axpby_v
-!!$
-!!$  !
-!!$  ! AXPBY is invoked via Y, hence the structure below. 
-!!$  !
-!!$  !
-!!$  !> Function  base_mv_axpby_a
-!!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
-!!$  !! \param m    Number of entries to be considere
-!!$  !! \param alpha scalar alpha
-!!$  !! \param x(:) The array to be added
-!!$  !! \param beta scalar alpha
-!!$  !! \param info   return code
-!!$  !!
-!!$  subroutine z_base_mv_axpby_a(m,alpha, x, beta, y, info)
-!!$    use psi_serial_mod
-!!$    implicit none 
-!!$    integer(psb_ipk_), intent(in)               :: m
-!!$    complex(psb_dpk_), intent(in)        :: x(:)
-!!$    class(psb_z_base_multivect_type), intent(inout)  :: y
-!!$    complex(psb_dpk_), intent (in)       :: alpha, beta
-!!$    integer(psb_ipk_), intent(out)              :: info
-!!$    
-!!$    call psb_geaxpby(m,alpha,x,beta,y%v,info)
-!!$    
-!!$  end subroutine z_base_mv_axpby_a
+  !
+  ! 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=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x     The class(base_mlv_vect) to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine z_base_mlv_axpby_v(m,alpha, x, beta, y, info, n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    class(psb_z_base_multivect_type), intent(inout)  :: x
+    class(psb_z_base_multivect_type), intent(inout)  :: y
+    complex(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(psb_size(x%v,2),psb_size(y%v,2))
+    end if
+    select type(xx => x)
+    type is (psb_z_base_multivect_type)
+      call psb_geaxpby(m,nc,alpha,x%v,beta,y%v,info)
+    class default
+      call y%axpby(m,alpha,x%v,beta,info,n=n)
+    end select
+
+  end subroutine z_base_mlv_axpby_v
+
+  !
+  ! AXPBY is invoked via Y, hence the structure below. 
+  !
+  !
+  !> Function  base_mlv_axpby_a
+  !! \memberof  psb_z_base_multivect_type
+  !! \brief AXPBY  by a normal array y=alpha*x+beta*y
+  !! \param m    Number of entries to be considered
+  !! \param alpha scalar alpha
+  !! \param x(:) The array to be added
+  !! \param beta scalar alpha
+  !! \param info   return code
+  !!
+  subroutine z_base_mlv_axpby_a(m,alpha, x, beta, y, info,n)
+    use psi_serial_mod
+    implicit none 
+    integer(psb_ipk_), intent(in)               :: m
+    complex(psb_dpk_), intent(in)        :: x(:,:)
+    class(psb_z_base_multivect_type), intent(inout)  :: y
+    complex(psb_dpk_), intent (in)       :: alpha, beta
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_), intent(in), optional     :: n
+    integer(psb_ipk_)  :: nc
+    if (present(n)) then
+      nc = n
+    else
+      nc = min(size(x,2),psb_size(y%v,2))
+    end if
+    
+    call psb_geaxpby(m,nc,alpha,x,beta,y%v,info)
+    
+  end subroutine z_base_mlv_axpby_a
 
   
 !!$  !
@@ -2060,13 +2087,13 @@ contains
 !!$  !  of the involved entities
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief Vector entry-by-entry multiply  by a base_mv_vect array y=x*y
-!!$  !! \param x   The class(base_mv_vect) to be multiplied by
+!!$  !! \brief Vector entry-by-entry multiply  by a base_mlv_vect array y=x*y
+!!$  !! \param x   The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine z_base_mv_mlt_v(x, y, info)
+!!$  subroutine z_base_mlv_mlt_v(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_z_base_multivect_type), intent(inout)  :: x
@@ -2085,16 +2112,16 @@ contains
 !!$      call y%mlt(x%v,info)
 !!$    end select
 !!$
-!!$  end subroutine z_base_mv_mlt_v
+!!$  end subroutine z_base_mlv_mlt_v
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a
+!!$  !> Function  base_mlv_mlt_a
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief Vector entry-by-entry multiply  by a normal array y=x*y
 !!$  !! \param x(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine z_base_mv_mlt_a(x, y, info)
+!!$  subroutine z_base_mlv_mlt_a(x, y, info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_dpk_), intent(in)        :: x(:)
@@ -2108,11 +2135,11 @@ contains
 !!$      y%v(i) = y%v(i)*x(i)
 !!$    end do
 !!$    
-!!$  end subroutine z_base_mv_mlt_a
+!!$  end subroutine z_base_mlv_mlt_a
 !!$
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_a_2
+!!$  !> Function  base_mlv_mlt_a_2
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by normal arrays
 !!$  !!        z=beta*z+alpha*x*y
@@ -2122,7 +2149,7 @@ contains
 !!$  !! \param y(:) The array to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine z_base_mv_mlt_a_2(alpha,x,y,beta,z,info)
+!!$  subroutine z_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_dpk_), intent(in)        :: alpha,beta
@@ -2187,20 +2214,20 @@ contains
 !!$        end if
 !!$      end if
 !!$    end if
-!!$  end subroutine z_base_mv_mlt_a_2
+!!$  end subroutine z_base_mlv_mlt_a_2
 !!$
 !!$  !
-!!$  !> Function  base_mv_mlt_v_2
+!!$  !> Function  base_mlv_mlt_v_2
 !!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mv_vect)
+!!$  !! \brief AXPBY-like Vector entry-by-entry multiply  by class(base_mlv_vect)
 !!$  !!        z=beta*z+alpha*x*y
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  !! \param x  The class(base_mv_vect) to be multiplied b
-!!$  !! \param y  The class(base_mv_vect) to be multiplied by
+!!$  !! \param x  The class(base_mlv_vect) to be multiplied b
+!!$  !! \param y  The class(base_mlv_vect) to be multiplied by
 !!$  !! \param info   return code
 !!$  !!
-!!$  subroutine z_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
+!!$  subroutine z_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
 !!$    use psi_serial_mod
 !!$    use psb_string_mod
 !!$    implicit none 
@@ -2227,9 +2254,9 @@ contains
 !!$      if (conjgx_) x%v=(x%v)
 !!$      if (conjgy_) y%v=(y%v)
 !!$    end if
-!!$  end subroutine z_base_mv_mlt_v_2
+!!$  end subroutine z_base_mlv_mlt_v_2
 !!$
-!!$  subroutine z_base_mv_mlt_av(alpha,x,y,beta,z,info)
+!!$  subroutine z_base_mlv_mlt_av(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_dpk_), intent(in)        :: alpha,beta
@@ -2243,9 +2270,9 @@ contains
 !!$    
 !!$    call z%mlt(alpha,x,y%v,beta,info)
 !!$
-!!$  end subroutine z_base_mv_mlt_av
+!!$  end subroutine z_base_mlv_mlt_av
 !!$
-!!$  subroutine z_base_mv_mlt_va(alpha,x,y,beta,z,info)
+!!$  subroutine z_base_mlv_mlt_va(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    complex(psb_dpk_), intent(in)        :: alpha,beta
@@ -2259,18 +2286,18 @@ contains
 !!$    
 !!$    call z%mlt(alpha,y,x,beta,info)
 !!$
-!!$  end subroutine z_base_mv_mlt_va
+!!$  end subroutine z_base_mlv_mlt_va
 !!$
 !!$
 !!$  !
 !!$  ! Simple scaling 
 !!$  !
-!!$  !> Function  base_mv_scal
+!!$  !> Function  base_mlv_scal
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief Scale all entries  x = alpha*x
 !!$  !! \param alpha   The multiplier
 !!$  !!
-!!$  subroutine z_base_mv_scal(alpha, x)
+!!$  subroutine z_base_mlv_scal(alpha, x)
 !!$    use psi_serial_mod
 !!$    implicit none 
 !!$    class(psb_z_base_multivect_type), intent(inout)  :: x
@@ -2278,16 +2305,16 @@ contains
 !!$    
 !!$    if (allocated(x%v)) x%v = alpha*x%v
 !!$
-!!$  end subroutine z_base_mv_scal
+!!$  end subroutine z_base_mlv_scal
 !!$  
 !!$  !
 !!$  ! Norms 1, 2 and infinity
 !!$  !
-!!$  !> Function  base_mv_nrm2
+!!$  !> Function  base_mlv_nrm2
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief 2-norm |x(1:n)|_2
 !!$  !! \param n  how many entries to consider
-!!$  function z_base_mv_nrm2(n,x) result(res)
+!!$  function z_base_mlv_nrm2(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_z_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2296,14 +2323,14 @@ contains
 !!$    
 !!$    res =  dnrm2(n,x%v,1)
 !!$
-!!$  end function z_base_mv_nrm2
+!!$  end function z_base_mlv_nrm2
 !!$  
 !!$  !
-!!$  !> Function  base_mv_amax
+!!$  !> Function  base_mlv_amax
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief infinity-norm |x(1:n)|_\infty
 !!$  !! \param n  how many entries to consider
-!!$  function z_base_mv_amax(n,x) result(res)
+!!$  function z_base_mlv_amax(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_z_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2311,14 +2338,14 @@ contains
 !!$    
 !!$    res =  maxval(abs(x%v(1:n)))
 !!$
-!!$  end function z_base_mv_amax
+!!$  end function z_base_mlv_amax
 !!$
 !!$  !
-!!$  !> Function  base_mv_asum
+!!$  !> Function  base_mlv_asum
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief 1-norm |x(1:n)|_1
 !!$  !! \param n  how many entries to consider
-!!$  function z_base_mv_asum(n,x) result(res)
+!!$  function z_base_mlv_asum(n,x) result(res)
 !!$    implicit none 
 !!$    class(psb_z_base_multivect_type), intent(inout) :: x
 !!$    integer(psb_ipk_), intent(in)           :: n
@@ -2326,14 +2353,14 @@ contains
 !!$    
 !!$    res =  sum(abs(x%v(1:n)))
 !!$
-!!$  end function z_base_mv_asum
+!!$  end function z_base_mlv_asum
 !!$  
 !!$  
 !!$  !
 !!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
 !!$  !
 !!$  !
-!!$  !> Function  base_mv_gthab
+!!$  !> Function  base_mlv_gthab
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief gather into an array
 !!$  !!    Y = beta * Y + alpha * X(IDX(:))
@@ -2341,7 +2368,7 @@ contains
 !!$  !! \param idx(:) indices
 !!$  !! \param alpha
 !!$  !! \param beta
-!!$  subroutine z_base_mv_gthab(n,idx,alpha,x,beta,y)
+!!$  subroutine z_base_mlv_gthab(n,idx,alpha,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_dpk_) :: alpha, beta, y(:)
@@ -2350,17 +2377,17 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,alpha,x%v,beta,y)
 !!$
-!!$  end subroutine z_base_mv_gthab
+!!$  end subroutine z_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine z_base_mv_gthzv_x(i,n,idx,x,y)
+!!$  subroutine z_base_mlv_gthzv_x(i,n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i,n
 !!$    class(psb_z_base_multivect_type) :: idx
@@ -2369,18 +2396,18 @@ contains
 !!$    
 !!$    call x%gth(n,idx%v(i:),y)
 !!$
-!!$  end subroutine z_base_mv_gthzv_x
+!!$  end subroutine z_base_mlv_gthzv_x
 !!$
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
-!!$  !> Function  base_mv_gthzv
+!!$  !> Function  base_mlv_gthzv
 !!$  !! \memberof  psb_z_base_multivect_type
 !!$  !! \brief gather into an array special alpha=1 beta=0
 !!$  !!    Y = X(IDX(:))
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
-!!$  subroutine z_base_mv_gthzv(n,idx,x,y)
+!!$  subroutine z_base_mlv_gthzv(n,idx,x,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_dpk_) ::  y(:)
@@ -2389,22 +2416,22 @@ contains
 !!$    call x%sync()
 !!$    call psi_gth(n,idx,x%v,y)
 !!$
-!!$  end subroutine z_base_mv_gthzv
+!!$  end subroutine z_base_mlv_gthzv
 !!$
 !!$  !
 !!$  ! Scatter: 
 !!$  ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
 !!$  ! 
 !!$  !
-!!$  !> Function  base_mv_sctb
+!!$  !> Function  base_mlv_sctb
 !!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief scatter into a class(base_mv_vect)
+!!$  !! \brief scatter into a class(base_mlv_vect)
 !!$  !!    Y(IDX(:)) = beta * Y(IDX(:)) +  X(:)
 !!$  !! \param n  how many entries to consider
 !!$  !! \param idx(:) indices
 !!$  !! \param beta
 !!$  !! \param x(:) 
-!!$  subroutine z_base_mv_sctb(n,idx,x,beta,y)
+!!$  subroutine z_base_mlv_sctb(n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: n, idx(:)
 !!$    complex(psb_dpk_) :: beta, x(:)
@@ -2414,9 +2441,9 @@ contains
 !!$    call psi_sct(n,idx,x,beta,y%v)
 !!$    call y%set_host()
 !!$
-!!$  end subroutine z_base_mv_sctb
+!!$  end subroutine z_base_mlv_sctb
 !!$
-!!$  subroutine z_base_mv_sctb_x(i,n,idx,x,beta,y)
+!!$  subroutine z_base_mlv_sctb_x(i,n,idx,x,beta,y)
 !!$    use psi_serial_mod
 !!$    integer(psb_ipk_) :: i, n
 !!$    class(psb_z_base_multivect_type) :: idx
@@ -2425,6 +2452,6 @@ contains
 !!$    
 !!$    call y%sct(n,idx%v(i:),x,beta)
 !!$
-!!$  end subroutine z_base_mv_sctb_x
+!!$  end subroutine z_base_mlv_sctb_x
 end module psb_z_base_multivect_mod