diff --git a/base/comm/internals/psi_dswapdata.F90 b/base/comm/internals/psi_dswapdata.F90
index 9f342098..9b3f0df3 100644
--- a/base/comm/internals/psi_dswapdata.F90
+++ b/base/comm/internals/psi_dswapdata.F90
@@ -227,17 +227,13 @@ subroutine psi_dswap_vidx_vect(iictxt,iicomm,flag,beta,y,idx, &
   logical :: swap_mpi, swap_sync, swap_send, swap_recv,&
        & albf,do_send,do_recv, swap_persistent, do_persistent
   logical, parameter :: usersend=.false., debug=.false.
-  integer(psb_ipk_), allocatable :: snd_counts(:), rcv_counts(:), &
-      snd_displs(:), rcv_displs(:), snd_to(:), rcv_from(:), snd_ws(:), rcv_ws(:)
   character(len=20)  :: name
 
   !remove
-  integer :: status(MPI_STATUS_SIZE), string_len, num_neighbors, snd_count, rcv_count
-  character(len=2*MPI_MAX_ERROR_STRING) :: mpistring
+  integer :: status(MPI_STATUS_SIZE)
   logical :: weight
   !character :: mpistring(16384)
 
-
   info=psb_success_
   name='psi_swap_datav'
   call psb_erractionsave(err_act)
@@ -270,34 +266,33 @@ subroutine psi_dswap_vidx_vect(iictxt,iicomm,flag,beta,y,idx, &
   ! check do_persistent twice, here and calling psi_dswapdata_vect, redudent but needed for now
   if (do_persistent) then
     ! if not allocated, allocate buffers and create request
-    if (.not. allocated(y%sndbuf)) then
-      allocate(y%sndbuf(totsnd))
-      allocate(y%rcvbuf(totrcv))
-      ! allocate(y%sndbuf(4))
-      ! allocate(y%rcvbuf(4))
+    if (.not. allocated(y%p)) then
+      allocate(y%p)
+    else
+      print *, "ALLLLLLOOOOCCCCAATTEED"
+    end if
+    if (.not. allocated(y%p%sndbuf)) then
+      allocate(y%p%sndbuf(totsnd), y%p%rcvbuf(totrcv))
+      allocate(y%p%rcv_count, y%p%snd_count)
       si = 1 ! sndbuf index
-      ri = 1 ! rcvbuf index
-      ! y%sndbuf = 10 + me ! remove this after working
-      y%rcvbuf = -1 ! remove this after working
 
-      ! call MPI_Graph_neighbors_count(icomm, me, num_neighbors, ierr)
-      ! print *, me, "~~~~~ NUM NEIGHBORS = ", num_neighbors
-      call MPI_Dist_graph_neighbors_count(icomm, rcv_count, snd_count, weight, ierr)
-      allocate(rcv_from(rcv_count), rcv_ws(rcv_count), snd_to(snd_count), snd_ws(snd_count))
-      call MPI_Dist_graph_neighbors(icomm, rcv_count, rcv_from, rcv_ws, snd_count, snd_to, &
-          snd_ws, ierr)
+      call MPI_Dist_graph_neighbors_count(icomm, y%p%rcv_count, &
+          y%p%snd_count, weight, ierr) ! should weight go into psb_d_persis_vect_type??
+
+      allocate(y%p%rcv_from(y%p%rcv_count), y%p%rcv_ws(y%p%rcv_count))
+      allocate(y%p%snd_to(y%p%snd_count), y%p%snd_ws(y%p%snd_count))
 
-      allocate(snd_counts(snd_count), rcv_counts(rcv_count), &
-          snd_displs(snd_count), rcv_displs(rcv_count))
+      call MPI_Dist_graph_neighbors(icomm, y%p%rcv_count, y%p%rcv_from, &
+          y%p%rcv_ws, y%p%snd_count, y%p%snd_to, y%p%snd_ws, ierr)
 
+      ! artless
+      allocate(y%p%snd_counts(y%p%snd_count), y%p%rcv_counts(y%p%rcv_count), &
+          y%p%snd_displs(y%p%snd_count), y%p%rcv_displs(y%p%rcv_count))
 
-      ! old
-      ! allocate(snd_counts(0:np-1), rcv_counts(0:np-1), &
-      !     snd_displs(0:np-1), rcv_displs(0:np-1))
-      snd_counts=0
-      rcv_counts=0
-      snd_displs=0
-      rcv_displs=0
+      y%p%snd_counts=0
+      y%p%rcv_counts=0
+      y%p%snd_displs=0
+      y%p%rcv_displs=0
 
       pnti   = 1
       snd_pt = 1
@@ -311,136 +306,129 @@ subroutine psi_dswap_vidx_vect(iictxt,iicomm,flag,beta,y,idx, &
         snd_pt = 1+pnti+nerv+psb_n_elem_send_
         rcv_pt = 1+pnti+psb_n_elem_recv_
 
-        do ii=1, snd_count
-          if (rcv_from(ii) == proc_to_comm) then
-            rcv_counts(ii) = nerv
-            rcv_counts(ii) = 4
-            rcv_displs(ii) = si - 1
+        do ii=1, y%p%snd_count
+          if (y%p%rcv_from(ii) == proc_to_comm) then
+            y%p%rcv_counts(ii) = nerv
+            y%p%rcv_displs(ii) = si - 1
           end if
-          if (snd_to(ii) == proc_to_comm) then
-            snd_counts(ii) = nesd
-            snd_counts(ii) = 4
-            snd_displs(ii) = si - 1
+          if (y%p%snd_to(ii) == proc_to_comm) then
+            y%p%snd_counts(ii) = nesd
+            y%p%snd_displs(ii) = si - 1
           end if
         end do
 
-        ! prepare sndbuf
-        ! si = si + 1
+        ! pack sndbuf
         do ii=0,nesd-1
-          ! y%sndbuf(si) = y%v(idx%v(ii+snd_pt))
-          y%sndbuf(si) = me + 10! si
+          y%p%sndbuf(si) = y%v(idx%v(ii+snd_pt))
           si = si + 1
         end do
-
-        ! subroutine psi_dgthv(n,idx,alpha,x,beta,y)  ! THIS THIS
-        ! y = beta*y(1:n) + alpha*x(idx(1:n))
-        ! print *, y%combuf  , idx // index of y
-        ! print *, "```idx_pt = ", idx_pt, "nesd = ", nesd
-        ! print *, me, ":::::::pnti=", pnti, "proc_to_comm=", proc_to_comm,"nerv=", &
-        !     nerv , "nesd=", nesd, "idx_pt=",idx_pt , "snd_pt=",snd_pt, &
-        !     "rcv_pt=", rcv_pt
-        ! print *, me, ":y%v", y%v
         pnti   = pnti + nerv + nesd + 3
-
       end do
 
       ! code for this in ~/src/psblas/psblas3/base/internals/psi_desc_impl.f90
       print *, me, ":totxch=", totxch,"totsnd", totsnd, "totrcv", totrcv
-      print *, me, ": rcv_count = ", rcv_count, "snd_count = ", snd_count
-      print *, me, ": rcv_from =", rcv_from, "snd_to =", snd_to
-      print *, me, ":snd_counts", snd_counts, "snd_displs = ",snd_displs, &
-          "rcv_counts", rcv_counts, "rcv_displs",rcv_displs
-
-      allocate(y%init_request)
-      y%init_request = -1
-      call MPIX_Neighbor_alltoallv_init(y%sndbuf, snd_counts, snd_displs, &
-          MPI_DOUBLE_PRECISION, y%rcvbuf, rcv_counts, rcv_displs, &
+      print *, me, ": rcv_count = ", y%p%rcv_count, "snd_count = ", y%p%snd_count
+      print *, me, ": y%p%rcv_from =", y%p%rcv_from, "y%p%snd_to =", y%p%snd_to
+      print *, me, ":y%p%snd_counts", y%p%snd_counts, "y%p%snd_displs = ",y%p%snd_displs, &
+          "y%p%rcv_counts", y%p%rcv_counts, "y%p%rcv_displs",y%p%rcv_displs
+      ! print *, me, ":y%p%sndbuf", y%p%sndbuf
+      allocate(y%p%init_request)
+      y%p%init_request = -1
+      call MPIX_Neighbor_alltoallv_init(y%p%sndbuf, y%p%snd_counts, y%p%snd_displs, &
+          MPI_DOUBLE_PRECISION, y%p%rcvbuf, y%p%rcv_counts, y%p%rcv_displs, &
           MPI_DOUBLE_PRECISION, icomm, MPI_INFO_NULL, &
-          ! MPI_REAL, y%rcvbuf, rcv_counts, rcv_displs, &
+          ! MPI_REAL, y%p%rcvbuf, y%p%rcv_counts, y%p%rcv_displs, &
           ! MPI_REAL, icomm, MPI_INFO_NULL, &
-          y%init_request, ierr)
+          y%p%init_request, ierr)
 
       if (ierr .ne. 0) then
         print *, "ERROR: MPIX_Neighbor_alltoallvinit ierr = ", ierr
         goto 9999
       end if
 
-      ! print *, me, ":y%init_request = ", y%init_request
-      deallocate(snd_counts, rcv_counts, snd_displs, rcv_displs)
-    else ! send and recv buffers exist, need to pack them
-      print *, "!!!!!!!!!!!!!!!ELSE"
+      ! print *, me, ":y%p%init_request = ", y%p%init_request
+      ! deallocate(y%p%snd_counts, y%p%rcv_counts, y%p%snd_displs, y%p%rcv_displs)
+    else ! send and recv buffers exist, need to pack send buffer
+      ! pnti   = 1
+      ! snd_pt = 1
+      ! rcv_pt = 1
+      ! do i=1, totxch
+      !   proc_to_comm = idx%v(pnti+psb_proc_id_)
+      !   nerv = idx%v(pnti+psb_n_elem_recv_)
+      !   nesd = idx%v(pnti+nerv+psb_n_elem_send_)
+      !   idx_pt = 1+pnti+psb_n_elem_recv_          !
+
+      !   snd_pt = 1+pnti+nerv+psb_n_elem_send_
+      !   rcv_pt = 1+pnti+psb_n_elem_recv_
+
+      !   do ii=1, snd_count
+      !     if (y%p%rcv_from(ii) == proc_to_comm) then
+      !       y%p%rcv_counts(ii) = nerv
+      !       y%p%rcv_displs(ii) = si - 1
+      !     end if
+      !     if (y%p%snd_to(ii) == proc_to_comm) then
+      !       y%p%snd_counts(ii) = nesd
+      !       y%p%snd_displs(ii) = si - 1
+      !     end if
+      !   end do
+
+      !   ! pack sndbuf
+      !   do ii=0,nesd-1
+      !     y%p%sndbuf(si) = y%v(idx%v(ii+snd_pt))
+      !     si = si + 1
+      !   end do
+      !   pnti   = pnti + nerv + nesd + 3
+      ! end do
+
+      print *, "!!!!!!!!!!!!!!!ELSE!!!!!!!!!!!!!!!!"
     end if
 
-    if (me == 3) then
-      ! print *, me,"PRE: y%sndbuf=", y%sndbuf!, "y%rcvbuf=",y%rcvbuf
-    end if
     ! start communication
-    if (.not. allocated(y%init_request)) then
-      print *, "error: y%init_request should be allocated"
+    if (.not. allocated(y%p%init_request)) then
+      print *, "error: y%p%init_request should be allocated"
       goto 9999
     end if
-    call MPI_Start(y%init_request, ierr)
+
+    call MPI_Start(y%p%init_request, ierr)
     if (ierr .ne. 0) then
-      print *, "ERROR: MPI_Start ierr = ", ierr
+      print *, "ERROR: rank ",me,"has MPI_Start status(MPI_ERROR) = ", &
+          status(MPI_ERROR), "and ierr = ", ierr
       goto 9999
     end if
 
-    call MPI_Wait(y%init_request, status, ierr)
-
-    print *, me,": Y%SNDBUF=", y%sndbuf, "Y%RCVBUF=",y%rcvbuf
-    if (status(MPI_ERROR) .ne. 0) then
-      ! print *, me,"ERROR: Y%SNDBUF=", y%sndbuf, "Y%RCVBUF=",y%rcvbuf
-      print *, me, "---ERROR---"
-      print *, "----"
-      ! call MPI_Error_string(status(MPI_ERROR), mpistring, string_len, ierr)
-      ! print *, "MPI_Error_string = ", mpistring(:string_len)
-      print *, "ERROR: rank ",me,"has MPI_Wait status(MPI_ERROR) = ", status(MPI_ERROR)
-      ! goto 9999
-    end if
-
+    call MPI_Wait(y%p%init_request, status, ierr)
 
-    if (me == 1) then
-      ! print *, me,"POST: y%sndbuf=", y%sndbuf, "y%rcvbuf=",y%rcvbuf
-      ! print *, me,"POST: y%rcvbuf=",y%rcvbuf
+    if (ierr .ne. 0) then
+      print *, "ERROR: rank ",me,"has MPI_Wait status(MPI_ERROR) = ", &
+          status(MPI_ERROR), "and ierr = ", ierr
+      goto 9999
     end if
-    ! print *, me,"POST: y%rcvbuf=",y%rcvbuf(:4)
-  end if
-
-
-      ! Then gather for sending.
-    !
-    ! pnti   = 1
-    ! snd_pt = totrcv_+1
-    ! rcv_pt = 1
-    ! do i=1, totxch
-    !   nerv = idx%v(pnti+psb_n_elem_recv_)
-    !   nesd = idx%v(pnti+nerv+psb_n_elem_send_)
-    !   idx_pt = 1+pnti+nerv+psb_n_elem_send_
-    !   call y%gth(idx_pt,snd_pt,nesd,idx)
-    !   rcv_pt = rcv_pt + n*nerv
-    !   snd_pt = snd_pt + n*nesd
-    !   pnti   = pnti + nerv + nesd + 3
-    ! end do
-
-  ! from later, is this needed??
-  ! pnti   = 1
-  ! snd_pt = 1
-  ! rcv_pt = 1
-  ! do i=1, totxch
-  !   proc_to_comm = idx%v(pnti+psb_proc_id_)
-  !   nerv = idx%v(pnti+psb_n_elem_recv_)
-  !   nesd = idx%v(pnti+nerv+psb_n_elem_send_)
-  !   idx_pt = 1+pnti+psb_n_elem_recv_
-  !   snd_pt = 1+pnti+nerv+psb_n_elem_send_
-  !   rcv_pt = 1+pnti+psb_n_elem_recv_
-
-  !   if (debug) write(0,*)me,' Received from: ',prcid(i),&
-  !        & y%combuf(rcv_pt:rcv_pt+nerv-1)
-  !   call y%sct(rcv_pt,nerv,idx,beta)
-  !   pnti   = pnti + nerv + nesd + 3
-  ! end do
 
+    ! ----scatter----
+    pnti   = 1
+    rcv_pt = 1
+    do i=1, totxch
+      proc_to_comm = idx%v(pnti+psb_proc_id_)
+      nerv = idx%v(pnti+psb_n_elem_recv_)
+      nesd = idx%v(pnti+nerv+psb_n_elem_send_)
+      idx_pt = 1+pnti+psb_n_elem_recv_          !
+      rcv_pt = 1+pnti+psb_n_elem_recv_
 
+      do ii=1, y%p%rcv_count
+        if (proc_to_comm .eq. y%p%rcv_from(ii)) then ! gather from rcvbuf
+          if (nerv .ne. y%p%rcv_counts(ii)) then
+            print *, "Error: dsi_dswapdata.F90:", &
+                "number of persistent elements received differs from nerv"
+            goto 9999
+          end if
+          do ri=1,nerv
+            y%v(idx%v(ri+rcv_pt-1)) = y%p%rcvbuf(ri + y%p%rcv_displs(ii))
+          end do
+        end if
+      end do
+      pnti   = pnti + nerv + nesd + 3
+    end do
+  end if ! persistent communication is complete
 
 
   if (debug) write(*,*) me,'Internal buffer'
diff --git a/base/modules/serial/psb_d_base_vect_mod.f90 b/base/modules/serial/psb_d_base_vect_mod.f90
index 8a59b513..d6dfad9f 100644
--- a/base/modules/serial/psb_d_base_vect_mod.f90
+++ b/base/modules/serial/psb_d_base_vect_mod.f90
@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
-!        Alfredo Buttari      
-!   
+!        Salvatore Filippone
+!        Alfredo Buttari
+!
 !    Redistribution and use in source and binary forms, with or without
 !    modification, are permitted provided that the following conditions
 !    are met:
@@ -15,7 +15,7 @@
 !      3. The name of the PSBLAS group or the names of its contributors may
 !         not be used to endorse or promote products derived from this
 !         software without specific written permission.
-!   
+!
 !    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 !    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 !    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
@@ -27,14 +27,14 @@
 !    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 !    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 !    POSSIBILITY OF SUCH DAMAGE.
-!   
-!    
+!
+!
 !
 ! package: psb_d_base_vect_mod
 !
 ! This module contains the definition of the psb_d_base_vect type which
 ! is a container for dense vectors.
-!  This is encapsulated instead of being just a simple array to allow for 
+!  This is encapsulated instead of being just a simple array to allow for
 !  more complicated situations, such as GPU programming, where the memory
 !  area we are interested in is not easily accessible from the host/Fortran
 !  side. It is also meant to be encapsulated in an outer type, to allow
@@ -43,7 +43,7 @@
 !
 !
 module psb_d_base_vect_mod
-  
+
   use psb_const_mod
   use psb_error_mod
   use psb_realloc_mod
@@ -51,20 +51,38 @@ module psb_d_base_vect_mod
   use psb_l_base_vect_mod
 
   !> \namespace  psb_base_mod  \class psb_d_base_vect_type
-  !! The psb_d_base_vect_type 
+  !! The psb_d_base_vect_type
   !! defines a middle level  real(psb_dpk_) encapsulated dense vector.
-  !! The encapsulation is needed, in place of a simple array, to allow  
+  !! The encapsulation is needed, in place of a simple array, to allow
   !! for complicated situations, such as GPU programming, where the memory
   !!  area we are interested in is not easily accessible from the host/Fortran
   !!  side. It is also meant to be encapsulated in an outer type, to allow
   !!  runtime switching as per the STATE design pattern, similar to the
   !!  sparse matrix types.
   !!
+  type psb_d_persis_vect_type
+    integer, allocatable :: init_request
+    real(psb_dpk_), allocatable :: sndbuf(:)
+    real(psb_dpk_), allocatable :: rcvbuf(:)
+    ! required for packing and unpacking when using persistent communication
+    integer(psb_ipk_), allocatable :: snd_count
+    integer(psb_ipk_), allocatable :: rcv_count
+    integer(psb_ipk_), allocatable :: snd_counts(:)
+    integer(psb_ipk_), allocatable :: rcv_counts(:)
+    integer(psb_ipk_), allocatable :: snd_ws(:)
+    integer(psb_ipk_), allocatable :: rcv_ws(:)
+    integer(psb_ipk_), allocatable :: snd_to(:)
+    integer(psb_ipk_), allocatable :: rcv_from(:)
+    integer(psb_ipk_), allocatable :: snd_displs(:) ! snd array displacements
+    integer(psb_ipk_), allocatable :: rcv_displs(:) ! rcv array displacements
+ end type psb_d_persis_vect_type
+
   type psb_d_base_vect_type
-    !> Values. 
+    !> Values.
     real(psb_dpk_), allocatable :: v(:)
-    real(psb_dpk_), allocatable :: combuf(:) 
+    real(psb_dpk_), allocatable :: combuf(:)
     integer(psb_mpk_), allocatable :: comid(:,:)
+    type(psb_d_persis_vect_type), allocatable  :: p
   contains
     !
     !  Constructors/allocators
@@ -78,7 +96,7 @@ module psb_d_base_vect_mod
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
-    ! in derived classes. 
+    ! in derived classes.
     !
     procedure, pass(x) :: ins_a    => d_base_ins_a
     procedure, pass(x) :: ins_v    => d_base_ins_v
@@ -93,7 +111,7 @@ module psb_d_base_vect_mod
     ! Any derived class having extra storage upon sync
     ! will guarantee that both fortran/host side and
     ! external side contain the same data. The base
-    ! version is only a placeholder. 
+    ! version is only a placeholder.
     !
     procedure, pass(x) :: sync     => d_base_sync
     procedure, pass(x) :: is_host  => d_base_is_host
@@ -130,7 +148,7 @@ module psb_d_base_vect_mod
     generic, public    :: set      => set_vect, set_scal
     !
     ! Gather/scatter. These are needed for MPI interfacing.
-    ! May have to be reworked. 
+    ! May have to be reworked.
     !
     procedure, pass(x) :: gthab    => d_base_gthab
     procedure, pass(x) :: gthzv    => d_base_gthzv
@@ -183,11 +201,11 @@ module psb_d_base_vect_mod
   end interface psb_d_base_vect
 
 contains
-  
+
   !
-  ! Constructors. 
+  ! Constructors.
   !
-  
+
   !> Function  constructor:
   !! \brief     Constructor from an array
   !!  \param   x(:)  input array to be copied
@@ -200,11 +218,11 @@ contains
     this%v = x
     call this%asb(size(x,kind=psb_ipk_),info)
   end function constructor
-    
-  
+
+
   !> Function  constructor:
   !! \brief     Constructor from size
-  !!  \param    n   Size of vector to be built. 
+  !!  \param    n   Size of vector to be built.
   !!
   function size_const(n) result(this)
     integer(psb_ipk_), intent(in) :: n
@@ -214,7 +232,7 @@ contains
     call this%asb(n,info)
 
   end function size_const
-  
+
   !
   ! Build from a sample
   !
@@ -226,20 +244,20 @@ contains
   !!
   subroutine d_base_bld_x(x,this)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in) :: this(:)
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
 
     call psb_realloc(size(this),x%v,info)
-    if (info /= 0) then 
+    if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_vect_bld')
       return
     end if
     x%v(:)  = this(:)
 
   end subroutine d_base_bld_x
-    
+
   !
   ! Create with size, but no initialization
   !
@@ -247,11 +265,11 @@ contains
   !> Function  bld_mn:
   !! \memberof  psb_d_base_vect_type
   !! \brief     Build method with size (uninitialized data)
-  !!  \param    n    size to be allocated. 
+  !!  \param    n    size to be allocated.
   !!
   subroutine d_base_bld_mn(x,n)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_mpk_), intent(in) :: n
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
@@ -260,15 +278,15 @@ contains
     call x%asb(n,info)
 
   end subroutine d_base_bld_mn
-  
+
   !> Function  bld_en:
   !! \memberof  psb_d_base_vect_type
   !! \brief     Build method with size (uninitialized data)
-  !!  \param    n    size to be allocated. 
+  !!  \param    n    size to be allocated.
   !!
   subroutine d_base_bld_en(x,n)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_epk_), intent(in) :: n
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_) :: info
@@ -277,24 +295,24 @@ contains
     call x%asb(n,info)
 
   end subroutine d_base_bld_en
-  
+
   !> Function  base_all:
   !! \memberof  psb_d_base_vect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
-  !!  \param    n    size to be allocated. 
+  !!  \param    n    size to be allocated.
   !!  \param info  return code
   !!
   subroutine d_base_all(n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)               :: n
     class(psb_d_base_vect_type), intent(out)    :: x
     integer(psb_ipk_), intent(out)              :: info
-    
+
     call psb_realloc(n,x%v,info)
-    
+
   end subroutine d_base_all
 
   !> Function  base_mold:
@@ -306,11 +324,11 @@ contains
   subroutine d_base_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in)   :: x
     class(psb_d_base_vect_type), intent(out), allocatable :: y
     integer(psb_ipk_), intent(out)              :: info
-    
+
     allocate(psb_d_base_vect_type :: y, stat=info)
 
   end subroutine d_base_mold
@@ -320,21 +338,21 @@ contains
   !
   !> Function  base_ins:
   !! \memberof  psb_d_base_vect_type
-  !! \brief Insert coefficients. 
+  !! \brief Insert coefficients.
   !!
   !!
   !!         Given  a list of N pairs
   !!           (IRL(i),VAL(i))
   !!         record a new coefficient in X such that
   !!            X(IRL(1:N)) = VAL(1:N).
-  !!            
+  !!
   !!         - the update operation will perform either
   !!               X(IRL(1:n)) = VAL(1:N)
   !!           or
   !!               X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N)
   !!           according to the value of DUPLICATE.
-  !!           
-  !!           
+  !!
+  !!
   !!  \param n     number of pairs in input
   !!  \param irl(:)  the input row indices
   !!  \param val(:)  the input coefficients
@@ -344,7 +362,7 @@ contains
   !
   subroutine d_base_ins_a(n,irl,val,dupl,x,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(in)               :: n, dupl
     integer(psb_ipk_), intent(in)               :: irl(:)
@@ -354,21 +372,21 @@ contains
     integer(psb_ipk_) :: i, isz
 
     info = 0
-    if (psb_errstatus_fatal()) return 
+    if (psb_errstatus_fatal()) return
 
-    if (.not.allocated(x%v)) then 
+    if (.not.allocated(x%v)) then
       info = psb_err_invalid_vect_state_
-    else if (n > min(size(irl),size(val))) then 
+    else if (n > min(size(irl),size(val))) then
       info = psb_err_invalid_input_
 
-    else 
+    else
       isz = size(x%v)
-      select case(dupl) 
-      case(psb_dupl_ovwrt_) 
+      select case(dupl)
+      case(psb_dupl_ovwrt_)
         do i = 1, n
           !loop over all val's rows
 
-          ! row actual block row 
+          ! row actual block row
           if ((1 <= irl(i)).and.(irl(i) <= isz)) then
             ! this row belongs to me
             ! copy i-th row of block val in x
@@ -376,7 +394,7 @@ contains
           end if
         enddo
 
-      case(psb_dupl_add_) 
+      case(psb_dupl_add_)
 
         do i = 1, n
           !loop over all val's rows
@@ -394,7 +412,7 @@ contains
       end select
     end if
     call x%set_host()
-    if (info /= 0) then 
+    if (info /= 0) then
       call psb_errpush(info,'base_vect_ins')
       return
     end if
@@ -403,7 +421,7 @@ contains
 
   subroutine d_base_ins_v(n,irl,val,dupl,x,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(in)               :: n, dupl
     class(psb_i_base_vect_type), intent(inout)  :: irl
@@ -413,14 +431,14 @@ contains
     integer(psb_ipk_) :: i, isz
 
     info = 0
-    if (psb_errstatus_fatal()) return 
+    if (psb_errstatus_fatal()) return
 
     if (irl%is_dev()) call irl%sync()
     if (val%is_dev()) call val%sync()
     if (x%is_dev())   call x%sync()
     call x%ins(n,irl%v,val%v,dupl,info)
 
-    if (info /= 0) then 
+    if (info /= 0) then
       call psb_errpush(info,'base_vect_ins')
       return
     end if
@@ -436,14 +454,14 @@ contains
   !
   subroutine d_base_zero(x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)    :: x
-    
+
     if (allocated(x%v)) x%v=dzero
     call x%set_host()
   end subroutine d_base_zero
 
-  
+
   !
   ! Assembly.
   ! For derived classes: after this the vector
@@ -452,20 +470,20 @@ contains
   !> Function  base_asb:
   !! \memberof  psb_d_base_vect_type
   !! \brief Assemble vector: reallocate as necessary.
-  !!           
+  !!
   !!  \param n     final size
   !!  \param info  return code
   !!
   !
- 
+
   subroutine d_base_asb_m(n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_mpk_), intent(in)              :: n
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
-    
+
     info = 0
     if (x%get_nrows() < n) &
          & call psb_realloc(n,x%v,info)
@@ -482,20 +500,20 @@ contains
   !> Function  base_asb:
   !! \memberof  psb_d_base_vect_type
   !! \brief Assemble vector: reallocate as necessary.
-  !!           
+  !!
   !!  \param n     final size
   !!  \param info  return code
   !!
   !
- 
+
   subroutine d_base_asb_e(n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_epk_), intent(in)              :: n
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
-    
+
     info = 0
     if (x%get_nrows() < n) &
          & call psb_realloc(n,x%v,info)
@@ -508,39 +526,39 @@ contains
   !> Function  base_free:
   !! \memberof  psb_d_base_vect_type
   !! \brief Free vector
-  !!           
+  !!
   !!  \param info  return code
   !!
   !
   subroutine d_base_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(out)              :: info
-    
+
     info = 0
     if (allocated(x%v)) deallocate(x%v, stat=info)
     if (info == 0) call x%free_buffer(info)
     if (info == 0) call x%free_comid(info)
-    if (info /= 0) call & 
+    if (info /= 0) call &
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
-        
+
   end subroutine d_base_free
 
-  
+
 
   !
   !> Function  base_free_buffer:
   !! \memberof  psb_d_base_vect_type
   !! \brief Free aux buffer
-  !!           
+  !!
   !!  \param info  return code
   !!
   !
   subroutine d_base_free_buffer(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
@@ -555,17 +573,17 @@ contains
   !!        In some derived classes, e.g. GPU,
   !!        does not really frees to avoid  runtime
   !!        costs
-  !!           
+  !!
   !!  \param info  return code
   !!
   !
   subroutine d_base_maybe_free_buffer(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
-    info = 0 
+    info = 0
     if (psb_get_maybe_free_buffer())&
          &  call x%free_buffer(info)
 
@@ -575,13 +593,13 @@ contains
   !> Function  base_free_comid:
   !! \memberof  psb_d_base_vect_type
   !! \brief Free aux MPI communication id buffer
-  !!           
+  !!
   !!  \param info  return code
   !!
   !
   subroutine d_base_free_comid(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
@@ -593,77 +611,77 @@ contains
   !
   ! The base version of SYNC & friends does nothing, it's just
   ! a placeholder.
-  ! 
+  !
   !
   !> Function  base_sync:
   !! \memberof  psb_d_base_vect_type
   !! \brief Sync: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_sync(x)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
-    
+
   end subroutine d_base_sync
 
   !
   !> Function  base_set_host:
   !! \memberof  psb_d_base_vect_type
   !! \brief Set_host: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_set_host(x)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
-    
+
   end subroutine d_base_set_host
 
   !
   !> Function  base_set_dev:
   !! \memberof  psb_d_base_vect_type
   !! \brief Set_dev: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_set_dev(x)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
-    
+
   end subroutine d_base_set_dev
 
   !
   !> Function  base_set_sync:
   !! \memberof  psb_d_base_vect_type
   !! \brief Set_sync: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_set_sync(x)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
-    
+
   end subroutine d_base_set_sync
 
   !
   !> Function  base_is_dev:
   !! \memberof  psb_d_base_vect_type
   !! \brief Is  vector on external device    .
-  !!           
+  !!
   !
   function d_base_is_dev(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in) :: x
     logical  :: res
-  
+
     res = .false.
   end function d_base_is_dev
-  
+
   !
   !> Function  base_is_host
   !! \memberof  psb_d_base_vect_type
   !! \brief Is  vector on standard memory    .
-  !!           
+  !!
   !
   function d_base_is_host(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in) :: x
     logical  :: res
 
@@ -674,10 +692,10 @@ contains
   !> Function  base_is_sync
   !! \memberof  psb_d_base_vect_type
   !! \brief Is  vector on sync               .
-  !!           
+  !!
   !
   function d_base_is_sync(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in) :: x
     logical  :: res
 
@@ -686,16 +704,16 @@ contains
 
 
   !
-  ! Size info. 
+  ! Size info.
   !
   !
   !> Function  base_get_nrows
   !! \memberof  psb_d_base_vect_type
   !! \brief  Number of entries
-  !!           
+  !!
   !
   function d_base_get_nrows(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in) :: x
     integer(psb_ipk_) :: res
 
@@ -708,13 +726,13 @@ contains
   !> Function  base_get_sizeof
   !! \memberof  psb_d_base_vect_type
   !! \brief  Size in bytes
-  !!           
+  !!
   !
   function d_base_sizeof(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(in) :: x
     integer(psb_epk_) :: res
-    
+
     ! Force 8-byte integers.
     res = (1_psb_epk_ * psb_sizeof_dp) * x%get_nrows()
 
@@ -724,14 +742,14 @@ contains
   !> Function  base_get_fmt
   !! \memberof  psb_d_base_vect_type
   !! \brief  Format
-  !!           
+  !!
   !
   function d_base_get_fmt() result(res)
-    implicit none 
+    implicit none
     character(len=5) :: res
     res = 'BASE'
   end function d_base_get_fmt
-  
+
 
   !
   !
@@ -740,7 +758,7 @@ contains
   !! \memberof  psb_d_base_vect_type
   !! \brief  Extract a copy of the contents
   !!
-  !    
+  !
   function  d_base_get_vect(x,n) result(res)
     class(psb_d_base_vect_type), intent(inout) :: x
     real(psb_dpk_), allocatable                 :: res(:)
@@ -748,21 +766,21 @@ contains
     integer(psb_ipk_), optional :: n
     ! Local variables
     integer(psb_ipk_) :: isz
-    
-    if (.not.allocated(x%v)) return 
+
+    if (.not.allocated(x%v)) return
     if (.not.x%is_host()) call x%sync()
     isz = x%get_nrows()
     if (present(n)) isz = max(0,min(isz,n))
-    allocate(res(isz),stat=info) 
-    if (info /= 0) then 
+    allocate(res(isz),stat=info)
+    if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_get_vect')
       return
     end if
     res(1:isz) = x%v(1:isz)
   end function d_base_get_vect
-    
+
   !
-  ! Reset all values 
+  ! Reset all values
   !
   !
   !> Function  base_set_scal
@@ -771,18 +789,18 @@ contains
   !! \param val   The value to set
   !!
   subroutine d_base_set_scal(x,val,first,last)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val
     integer(psb_ipk_), optional :: first, last
-        
+
     integer(psb_ipk_) :: info, first_, last_
 
     first_=1
     last_=size(x%v)
     if (present(first)) first_ = max(1,first)
     if (present(last))  last_  = min(last,last_)
-    
+
     if (x%is_dev()) call x%sync()
     x%v(first_:last_) = val
     call x%set_host()
@@ -794,14 +812,14 @@ contains
   !> Function  base_set_vect
   !! \memberof  psb_d_base_vect_type
   !! \brief  Set all entries
-  !! \param val(:)  The vector to be copied in 
+  !! \param val(:)  The vector to be copied in
   !!
   subroutine d_base_set_vect(x,val,first,last)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val(:)
     integer(psb_ipk_), optional :: first, last
-        
+
     integer(psb_ipk_) :: info, first_, last_, nr
 
     first_                     = 1
@@ -809,7 +827,7 @@ contains
     last_                      = min(psb_size(x%v),first_+size(val)-1)
     if (present(last))  last_  = min(last,last_)
 
-    if (allocated(x%v)) then 
+    if (allocated(x%v)) then
       if (x%is_dev()) call x%sync()
       x%v(first_:last_) = val(1:last_-first_+1)
     else
@@ -829,7 +847,7 @@ contains
   !! \brief  Set all entries to their respective absolute values.
   !!
   subroutine d_base_absval1(x)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
 
     if (allocated(x%v)) then
@@ -841,21 +859,21 @@ contains
   end subroutine d_base_absval1
 
   subroutine d_base_absval2(x,y)
-    implicit none 
-    class(psb_d_base_vect_type), intent(inout) :: x    
+    implicit none
+    class(psb_d_base_vect_type), intent(inout) :: x
     class(psb_d_base_vect_type), intent(inout) :: y
     integer(psb_ipk_) :: info
     if (.not.x%is_host()) call x%sync()
-    if (allocated(x%v)) then 
+    if (allocated(x%v)) then
       call y%axpby(ione*min(x%get_nrows(),y%get_nrows()),done,x,dzero,info)
       call y%absval()
     end if
-    
+
   end subroutine d_base_absval2
 
   !
-  ! Dot products 
-  ! 
+  ! Dot products
+  !
   !
   !> Function  base_dot_v
   !! \memberof  psb_d_base_vect_type
@@ -864,12 +882,12 @@ contains
   !! \param y    The other (base_vect) to be multiplied by
   !!
   function d_base_dot_v(n,x,y) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x, y
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
     real(psb_dpk_), external      :: ddot
-    
+
     res = dzero
     !
     ! Note: this is the base implementation.
@@ -898,19 +916,19 @@ contains
   !! \param y(:) The array to be multiplied by
   !!
   function d_base_dot_a(n,x,y) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     real(psb_dpk_), intent(in)    :: y(:)
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
     real(psb_dpk_), external      :: ddot
-    
+
     res = ddot(n,y,1,x%v,1)
 
   end function d_base_dot_a
-    
+
   !
-  ! AXPBY is invoked via Y, hence the structure below. 
+  ! AXPBY is invoked via Y, hence the structure below.
   !
   !
   !
@@ -925,13 +943,13 @@ contains
   !!
   subroutine d_base_axpby_v(m,alpha, x, beta, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)               :: m
     class(psb_d_base_vect_type), intent(inout)  :: x
     class(psb_d_base_vect_type), intent(inout)  :: y
     real(psb_dpk_), intent (in)       :: alpha, beta
     integer(psb_ipk_), intent(out)              :: info
-    
+
     if (x%is_dev()) call x%sync()
 
     call y%axpby(m,alpha,x%v,beta,info)
@@ -939,7 +957,7 @@ contains
   end subroutine d_base_axpby_v
 
   !
-  ! AXPBY is invoked via Y, hence the structure below. 
+  ! AXPBY is invoked via Y, hence the structure below.
   !
   !
   !> Function  base_axpby_a
@@ -953,20 +971,20 @@ contains
   !!
   subroutine d_base_axpby_a(m,alpha, x, beta, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)               :: m
     real(psb_dpk_), intent(in)        :: x(:)
     class(psb_d_base_vect_type), intent(inout)  :: y
     real(psb_dpk_), intent (in)       :: alpha, beta
     integer(psb_ipk_), intent(out)              :: info
-    
+
     if (y%is_dev()) call y%sync()
     call psb_geaxpby(m,alpha,x,beta,y%v,info)
     call y%set_host()
-    
+
   end subroutine d_base_axpby_a
 
-  
+
   !
   !  Multiple variants of two operations:
   !  Simple multiplication  Y(:) = X(:)*Y(:)
@@ -984,10 +1002,10 @@ contains
   !!
   subroutine d_base_mlt_v(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     class(psb_d_base_vect_type), intent(inout)  :: y
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_) :: i, n
 
     info = 0
@@ -1005,7 +1023,7 @@ contains
   !!
   subroutine d_base_mlt_a(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: x(:)
     class(psb_d_base_vect_type), intent(inout)  :: y
     integer(psb_ipk_), intent(out)              :: info
@@ -1014,7 +1032,7 @@ contains
     info = 0
     if (y%is_dev()) call y%sync()
     n = min(size(y%v), size(x))
-    do i=1, n 
+    do i=1, n
       y%v(i) = y%v(i)*x(i)
     end do
     call y%set_host()
@@ -1035,7 +1053,7 @@ contains
   !!
   subroutine d_base_mlt_a_2(alpha,x,y,beta,z,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     real(psb_dpk_), intent(in)        :: y(:)
     real(psb_dpk_), intent(in)        :: x(:)
@@ -1043,58 +1061,58 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_) :: i, n
 
-    info = 0    
+    info = 0
     if (z%is_dev()) call z%sync()
 
     n = min(size(z%v), size(x), size(y))
-    if (alpha == dzero) then 
-      if (beta == done) then 
-        return 
+    if (alpha == dzero) then
+      if (beta == done) then
+        return
       else
         do i=1, n
           z%v(i) = beta*z%v(i)
         end do
       end if
     else
-      if (alpha == done) then 
-        if (beta == dzero) then 
-          do i=1, n 
+      if (alpha == done) then
+        if (beta == dzero) then
+          do i=1, n
             z%v(i) = y(i)*x(i)
           end do
-        else if (beta == done) then 
-          do i=1, n 
+        else if (beta == done) then
+          do i=1, n
             z%v(i) = z%v(i) + y(i)*x(i)
           end do
-        else 
-          do i=1, n 
+        else
+          do i=1, n
             z%v(i) = beta*z%v(i) + y(i)*x(i)
           end do
         end if
-      else if (alpha == -done) then 
-        if (beta == dzero) then 
-          do i=1, n 
+      else if (alpha == -done) then
+        if (beta == dzero) then
+          do i=1, n
             z%v(i) = -y(i)*x(i)
           end do
-        else if (beta == done) then 
-          do i=1, n 
+        else if (beta == done) then
+          do i=1, n
             z%v(i) = z%v(i) - y(i)*x(i)
           end do
-        else 
-          do i=1, n 
+        else
+          do i=1, n
             z%v(i) = beta*z%v(i) - y(i)*x(i)
           end do
         end if
       else
-        if (beta == dzero) then 
-          do i=1, n 
+        if (beta == dzero) then
+          do i=1, n
             z%v(i) = alpha*y(i)*x(i)
           end do
-        else if (beta == done) then 
-          do i=1, n 
+        else if (beta == done) then
+          do i=1, n
             z%v(i) = z%v(i) + alpha*y(i)*x(i)
           end do
-        else 
-          do i=1, n 
+        else
+          do i=1, n
             z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
           end do
         end if
@@ -1118,12 +1136,12 @@ contains
   subroutine d_base_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
     use psi_serial_mod
     use psb_string_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     class(psb_d_base_vect_type), intent(inout)  :: x
     class(psb_d_base_vect_type), intent(inout)  :: y
     class(psb_d_base_vect_type), intent(inout)  :: z
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
     character(len=1), intent(in), optional     :: conjgx, conjgy
     integer(psb_ipk_) :: i, n
     logical :: conjgx_, conjgy_
@@ -1133,7 +1151,7 @@ contains
     if (x%is_dev()) call x%sync()
     if (.not.psb_d_is_complex_) then
       call z%mlt(alpha,x%v,y%v,beta,info)
-    else 
+    else
       conjgx_=.false.
       if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
       conjgy_=.false.
@@ -1148,12 +1166,12 @@ contains
 
   subroutine d_base_mlt_av(alpha,x,y,beta,z,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     real(psb_dpk_), intent(in)        :: x(:)
     class(psb_d_base_vect_type), intent(inout)  :: y
     class(psb_d_base_vect_type), intent(inout)  :: z
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_) :: i, n
 
     info = 0
@@ -1164,12 +1182,12 @@ contains
 
   subroutine d_base_mlt_va(alpha,x,y,beta,z,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     real(psb_dpk_), intent(in)        :: y(:)
     class(psb_d_base_vect_type), intent(inout)  :: x
     class(psb_d_base_vect_type), intent(inout)  :: z
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_) :: i, n
 
     info = 0
@@ -1180,7 +1198,7 @@ contains
 
 
   !
-  ! Simple scaling 
+  ! Simple scaling
   !
   !> Function  base_scal
   !! \memberof  psb_d_base_vect_type
@@ -1189,17 +1207,17 @@ contains
   !!
   subroutine d_base_scal(alpha, x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     real(psb_dpk_), intent (in)       :: alpha
-    
-    if (allocated(x%v)) then 
+
+    if (allocated(x%v)) then
       x%v = alpha*x%v
       call x%set_host()
     end if
 
   end subroutine d_base_scal
-  
+
   !
   ! Norms 1, 2 and infinity
   !
@@ -1208,28 +1226,28 @@ contains
   !! \brief 2-norm |x(1:n)|_2
   !! \param n  how many entries to consider
   function d_base_nrm2(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
     real(psb_dpk_), external      :: dnrm2
-    
+
     if (x%is_dev()) call x%sync()
     res =  dnrm2(n,x%v,1)
 
   end function d_base_nrm2
-  
+
   !
   !> Function  base_amax
   !! \memberof  psb_d_base_vect_type
   !! \brief infinity-norm |x(1:n)|_\infty
   !! \param n  how many entries to consider
   function d_base_amax(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
-    
+
     if (x%is_dev()) call x%sync()
     res =  maxval(abs(x%v(1:n)))
 
@@ -1241,17 +1259,17 @@ contains
   !! \brief 1-norm |x(1:n)|_1
   !! \param n  how many entries to consider
   function d_base_asum(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_)                :: res
-    
+
     if (x%is_dev()) call x%sync()
     res =  sum(abs(x%v(1:n)))
 
   end function d_base_asum
-  
-  
+
+
   !
   ! Gather: Y = beta * Y + alpha * X(IDX(:))
   !
@@ -1266,18 +1284,18 @@ contains
   !! \param beta
   subroutine d_base_gthab(n,idx,alpha,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) :: alpha, beta, y(:)
     class(psb_d_base_vect_type) :: x
-    
+
     if (x%is_dev()) call x%sync()
     call psi_gth(n,idx,alpha,x%v,beta,y)
 
   end subroutine d_base_gthab
   !
   ! shortcut alpha=1 beta=0
-  ! 
+  !
   !> Function  base_gthzv
   !! \memberof  psb_d_base_vect_type
   !! \brief gather into an array special alpha=1 beta=0
@@ -1286,28 +1304,28 @@ contains
   !! \param idx(:) indices
   subroutine d_base_gthzv_x(i,n,idx,x,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i,n
     class(psb_i_base_vect_type) :: idx
     real(psb_dpk_) ::  y(:)
     class(psb_d_base_vect_type) :: x
-    
+
     if (idx%is_dev()) call idx%sync()
     call x%gth(n,idx%v(i:),y)
 
   end subroutine d_base_gthzv_x
 
   !
-  ! New comm internals impl. 
+  ! New comm internals impl.
   !
   subroutine d_base_gthzbuf(i,n,idx,x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i,n
     class(psb_i_base_vect_type) :: idx
     class(psb_d_base_vect_type) :: x
-    
-    if (.not.allocated(x%combuf)) then 
+
+    if (.not.allocated(x%combuf)) then
       call psb_errpush(psb_err_alloc_dealloc_,'gthzbuf')
       return
     end if
@@ -1320,22 +1338,22 @@ contains
   !> Function  base_device_wait:
   !! \memberof  psb_d_base_vect_type
   !! \brief device_wait: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_device_wait()
-    implicit none 
-    
+    implicit none
+
   end subroutine d_base_device_wait
 
   function d_base_use_buffer() result(res)
     logical :: res
-    
+
     res = .true.
   end function d_base_use_buffer
 
   subroutine d_base_new_buffer(n,x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)              :: n
     integer(psb_ipk_), intent(out)             :: info
@@ -1345,7 +1363,7 @@ contains
 
   subroutine d_base_new_comid(n,x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)              :: n
     integer(psb_ipk_), intent(out)             :: info
@@ -1356,7 +1374,7 @@ contains
 
   !
   ! shortcut alpha=1 beta=0
-  ! 
+  !
   !> Function  base_gthzv
   !! \memberof  psb_d_base_vect_type
   !! \brief gather into an array special alpha=1 beta=0
@@ -1365,20 +1383,20 @@ contains
   !! \param idx(:) indices
   subroutine d_base_gthzv(n,idx,x,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) ::  y(:)
     class(psb_d_base_vect_type) :: x
-    
+
     if (x%is_dev()) call x%sync()
     call psi_gth(n,idx,x%v,y)
 
   end subroutine d_base_gthzv
 
   !
-  ! Scatter: 
+  ! Scatter:
   ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
-  ! 
+  !
   !
   !> Function  base_sctb
   !! \memberof  psb_d_base_vect_type
@@ -1387,14 +1405,14 @@ contains
   !! \param n  how many entries to consider
   !! \param idx(:) indices
   !! \param beta
-  !! \param x(:) 
+  !! \param x(:)
   subroutine d_base_sctb(n,idx,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) :: beta, x(:)
     class(psb_d_base_vect_type) :: y
-    
+
     if (y%is_dev()) call y%sync()
     call psi_sct(n,idx,x,beta,y%v)
     call y%set_host()
@@ -1403,12 +1421,12 @@ contains
 
   subroutine d_base_sctb_x(i,n,idx,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i, n
     class(psb_i_base_vect_type) :: idx
     real(psb_dpk_) :: beta, x(:)
     class(psb_d_base_vect_type) :: y
-    
+
     if (idx%is_dev()) call idx%sync()
     call y%sct(n,idx%v(i:),x,beta)
     call y%set_host()
@@ -1417,14 +1435,14 @@ contains
 
   subroutine d_base_sctb_buf(i,n,idx,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i, n
     class(psb_i_base_vect_type) :: idx
     real(psb_dpk_) :: beta
     class(psb_d_base_vect_type) :: y
-    
-    
-    if (.not.allocated(y%combuf)) then 
+
+
+    if (.not.allocated(y%combuf)) then
       call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
       return
     end if
@@ -1449,22 +1467,22 @@ module psb_d_base_multivect_mod
   use psb_d_base_vect_mod
 
   !> \namespace  psb_base_mod  \class psb_d_base_vect_type
-  !! The psb_d_base_vect_type 
+  !! The psb_d_base_vect_type
   !! defines a middle level  integer(psb_ipk_) encapsulated dense vector.
-  !! The encapsulation is needed, in place of a simple array, to allow  
+  !! The encapsulation is needed, in place of a simple array, to allow
   !! for complicated situations, such as GPU programming, where the memory
   !!  area we are interested in is not easily accessible from the host/Fortran
   !!  side. It is also meant to be encapsulated in an outer type, to allow
   !!  runtime switching as per the STATE design pattern, similar to the
   !!  sparse matrix types.
   !!
-  private 
+  private
   public  :: psb_d_base_multivect, psb_d_base_multivect_type
 
   type psb_d_base_multivect_type
-    !> Values. 
+    !> Values.
     real(psb_dpk_), allocatable :: v(:,:)
-    real(psb_dpk_), allocatable :: combuf(:) 
+    real(psb_dpk_), allocatable :: combuf(:)
     integer(psb_mpk_), allocatable :: comid(:,:)
   contains
     !
@@ -1478,7 +1496,7 @@ module psb_d_base_multivect_mod
     !
     ! Insert/set. Assembly and free.
     ! Assembly does almost nothing here, but is important
-    ! in derived classes. 
+    ! in derived classes.
     !
     procedure, pass(x) :: ins      => d_base_mlv_ins
     procedure, pass(x) :: zero     => d_base_mlv_zero
@@ -1489,7 +1507,7 @@ module psb_d_base_multivect_mod
     ! Any derived class having extra storage upon sync
     ! will guarantee that both fortran/host side and
     ! external side contain the same data. The base
-    ! version is only a placeholder. 
+    ! version is only a placeholder.
     !
     procedure, pass(x) :: sync     => d_base_mlv_sync
     procedure, pass(x) :: is_host  => d_base_mlv_is_host
@@ -1562,7 +1580,7 @@ module psb_d_base_multivect_mod
 
     !
     ! Gather/scatter. These are needed for MPI interfacing.
-    ! May have to be reworked. 
+    ! May have to be reworked.
     !
     procedure, pass(x) :: gthab    => d_base_mlv_gthab
     procedure, pass(x) :: gthzv    => d_base_mlv_gthzv
@@ -1584,7 +1602,7 @@ module psb_d_base_multivect_mod
 contains
 
   !
-  ! Constructors. 
+  ! Constructors.
   !
 
   !> Function  constructor:
@@ -1603,7 +1621,7 @@ contains
 
   !> Function  constructor:
   !! \brief     Constructor from size
-  !!  \param    n   Size of vector to be built. 
+  !!  \param    n   Size of vector to be built.
   !!
   function size_const(m,n) result(this)
     integer(psb_ipk_), intent(in) :: m,n
@@ -1630,7 +1648,7 @@ contains
     integer(psb_ipk_) :: info
 
     call psb_realloc(size(this,1),size(this,2),x%v,info)
-    if (info /= 0) then 
+    if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_vect_bld')
       return
     end if
@@ -1645,7 +1663,7 @@ contains
   !> Function  bld_n:
   !! \memberof  psb_d_base_multivect_type
   !! \brief     Build method with size (uninitialized data)
-  !!  \param    n    size to be allocated. 
+  !!  \param    n    size to be allocated.
   !!
   subroutine d_base_mlv_bld_n(x,m,n)
     use psb_realloc_mod
@@ -1662,13 +1680,13 @@ contains
   !! \memberof  psb_d_base_multivect_type
   !! \brief     Build method with size (uninitialized data) and
   !!            allocation return code.
-  !!  \param    n    size to be allocated. 
+  !!  \param    n    size to be allocated.
   !!  \param info  return code
   !!
   subroutine d_base_mlv_all(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)               :: m,n
     class(psb_d_base_multivect_type), intent(out) :: x
     integer(psb_ipk_), intent(out)              :: info
@@ -1686,7 +1704,7 @@ contains
   subroutine d_base_mlv_mold(x, y, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in)   :: x
     class(psb_d_base_multivect_type), intent(out), allocatable :: y
     integer(psb_ipk_), intent(out)              :: info
@@ -1700,21 +1718,21 @@ contains
   !
   !> Function  base_mlv_ins:
   !! \memberof  psb_d_base_multivect_type
-  !! \brief Insert coefficients. 
+  !! \brief Insert coefficients.
   !!
   !!
   !!         Given  a list of N pairs
   !!           (IRL(i),VAL(i))
   !!         record a new coefficient in X such that
   !!            X(IRL(1:N)) = VAL(1:N).
-  !!            
+  !!
   !!         - the update operation will perform either
   !!               X(IRL(1:n)) = VAL(1:N)
   !!           or
   !!               X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N)
   !!           according to the value of DUPLICATE.
-  !!           
-  !!           
+  !!
+  !!
   !!  \param n     number of pairs in input
   !!  \param irl(:)  the input row indices
   !!  \param val(:)  the input coefficients
@@ -1724,7 +1742,7 @@ contains
   !
   subroutine d_base_mlv_ins(n,irl,val,dupl,x,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(in)               :: n, dupl
     integer(psb_ipk_), intent(in)               :: irl(:)
@@ -1734,21 +1752,21 @@ contains
     integer(psb_ipk_) :: i, isz
 
     info = 0
-    if (psb_errstatus_fatal()) return 
+    if (psb_errstatus_fatal()) return
 
-    if (.not.allocated(x%v)) then 
+    if (.not.allocated(x%v)) then
       info = psb_err_invalid_vect_state_
-    else if (n > min(size(irl),size(val))) then 
+    else if (n > min(size(irl),size(val))) then
       info = psb_err_invalid_input_
 
-    else 
+    else
       isz = size(x%v,1)
-      select case(dupl) 
-      case(psb_dupl_ovwrt_) 
+      select case(dupl)
+      case(psb_dupl_ovwrt_)
         do i = 1, n
           !loop over all val's rows
 
-          ! row actual block row 
+          ! row actual block row
           if ((1 <= irl(i)).and.(irl(i) <= isz)) then
             ! this row belongs to me
             ! copy i-th row of block val in x
@@ -1756,7 +1774,7 @@ contains
           end if
         enddo
 
-      case(psb_dupl_add_) 
+      case(psb_dupl_add_)
 
         do i = 1, n
           !loop over all val's rows
@@ -1773,7 +1791,7 @@ contains
         ! !$      goto 9999
       end select
     end if
-    if (info /= 0) then 
+    if (info /= 0) then
       call psb_errpush(info,'base_mlv_vect_ins')
       return
     end if
@@ -1788,7 +1806,7 @@ contains
   !
   subroutine d_base_mlv_zero(x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)    :: x
 
     if (allocated(x%v)) x%v=dzero
@@ -1804,7 +1822,7 @@ contains
   !> Function  base_mlv_asb:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Assemble vector: reallocate as necessary.
-  !!           
+  !!
   !!  \param n     final size
   !!  \param info  return code
   !!
@@ -1813,7 +1831,7 @@ contains
   subroutine d_base_mlv_asb(m,n, x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)              :: m,n
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
@@ -1830,20 +1848,20 @@ contains
   !> Function  base_mlv_free:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Free vector
-  !!           
+  !!
   !!  \param info  return code
   !!
   !
   subroutine d_base_mlv_free(x, info)
     use psi_serial_mod
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     integer(psb_ipk_), intent(out)              :: info
 
     info = 0
     if (allocated(x%v)) deallocate(x%v, stat=info)
-    if (info /= 0) call & 
+    if (info /= 0) call &
          & psb_errpush(psb_err_alloc_dealloc_,'vect_free')
 
   end subroutine d_base_mlv_free
@@ -1853,15 +1871,15 @@ contains
   !
   ! The base version of SYNC & friends does nothing, it's just
   ! a placeholder.
-  ! 
+  !
   !
   !> Function  base_mlv_sync:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Sync: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_mlv_sync(x)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
 
   end subroutine d_base_mlv_sync
@@ -1870,10 +1888,10 @@ contains
   !> Function  base_mlv_set_host:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_host: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_mlv_set_host(x)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
 
   end subroutine d_base_mlv_set_host
@@ -1882,10 +1900,10 @@ contains
   !> Function  base_mlv_set_dev:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_dev: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_mlv_set_dev(x)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
 
   end subroutine d_base_mlv_set_dev
@@ -1894,10 +1912,10 @@ contains
   !> Function  base_mlv_set_sync:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Set_sync: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_mlv_set_sync(x)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
 
   end subroutine d_base_mlv_set_sync
@@ -1906,10 +1924,10 @@ contains
   !> Function  base_mlv_is_dev:
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on external device    .
-  !!           
+  !!
   !
   function d_base_mlv_is_dev(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
 
@@ -1920,10 +1938,10 @@ contains
   !> Function  base_mlv_is_host
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on standard memory    .
-  !!           
+  !!
   !
   function d_base_mlv_is_host(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
 
@@ -1934,10 +1952,10 @@ contains
   !> Function  base_mlv_is_sync
   !! \memberof  psb_d_base_multivect_type
   !! \brief Is  vector on sync               .
-  !!           
+  !!
   !
   function d_base_mlv_is_sync(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     logical  :: res
 
@@ -1946,16 +1964,16 @@ contains
 
 
   !
-  ! Size info. 
+  ! Size info.
   !
   !
   !> Function  base_mlv_get_nrows
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Number of entries
-  !!           
+  !!
   !
   function d_base_mlv_get_nrows(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
 
@@ -1965,7 +1983,7 @@ contains
   end function d_base_mlv_get_nrows
 
   function d_base_mlv_get_ncols(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_ipk_) :: res
 
@@ -1978,10 +1996,10 @@ contains
   !> Function  base_mlv_get_sizeof
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Size in bytesa
-  !!           
+  !!
   !
   function d_base_mlv_sizeof(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(in) :: x
     integer(psb_epk_) :: res
 
@@ -1994,10 +2012,10 @@ contains
   !> Function  base_mlv_get_fmt
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Format
-  !!           
+  !!
   !
   function d_base_mlv_get_fmt() result(res)
-    implicit none 
+    implicit none
     character(len=5) :: res
     res = 'BASE'
   end function d_base_mlv_get_fmt
@@ -2010,18 +2028,18 @@ contains
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Extract a copy of the contents
   !!
-  !    
+  !
   function  d_base_mlv_get_vect(x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     real(psb_dpk_), allocatable                 :: res(:,:)
     integer(psb_ipk_) :: info,m,n
     m = x%get_nrows()
     n = x%get_ncols()
-    if (.not.allocated(x%v)) return 
+    if (.not.allocated(x%v)) return
     call x%sync()
-    allocate(res(m,n),stat=info) 
-    if (info /= 0) then 
+    allocate(res(m,n),stat=info)
+    if (info /= 0) then
       call psb_errpush(psb_err_alloc_dealloc_,'base_mlv_get_vect')
       return
     end if
@@ -2029,7 +2047,7 @@ contains
   end function d_base_mlv_get_vect
 
   !
-  ! Reset all values 
+  ! Reset all values
   !
   !
   !> Function  base_mlv_set_scal
@@ -2038,7 +2056,7 @@ contains
   !! \param val   The value to set
   !!
   subroutine d_base_mlv_set_scal(x,val)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val
 
@@ -2051,16 +2069,16 @@ contains
   !> Function  base_mlv_set_vect
   !! \memberof  psb_d_base_multivect_type
   !! \brief  Set all entries
-  !! \param val(:)  The vector to be copied in 
+  !! \param val(:)  The vector to be copied in
   !!
   subroutine d_base_mlv_set_vect(x,val)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     real(psb_dpk_), intent(in) :: val(:,:)
     integer(psb_ipk_) :: nr, nc
     integer(psb_ipk_) :: info
 
-    if (allocated(x%v)) then 
+    if (allocated(x%v)) then
       nr = min(size(x%v,1),size(val,1))
       nc = min(size(x%v,2),size(val,2))
 
@@ -2072,8 +2090,8 @@ contains
   end subroutine d_base_mlv_set_vect
 
   !
-  ! Dot products 
-  ! 
+  ! Dot products
+  !
   !
   !> Function  base_mlv_dot_v
   !! \memberof  psb_d_base_multivect_type
@@ -2082,7 +2100,7 @@ contains
   !! \param y    The other (base_mlv_vect) to be multiplied by
   !!
   function d_base_mlv_dot_v(n,x,y) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x, y
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_), allocatable   :: res(:)
@@ -2094,7 +2112,7 @@ contains
     !
     ! 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). 
+    !  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
     !
@@ -2123,7 +2141,7 @@ contains
   !! \param y(:) The array to be multiplied by
   !!
   function d_base_mlv_dot_a(n,x,y) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     real(psb_dpk_), intent(in)    :: y(:,:)
     integer(psb_ipk_), intent(in)           :: n
@@ -2141,7 +2159,7 @@ contains
   end function d_base_mlv_dot_a
 
   !
-  ! AXPBY is invoked via Y, hence the structure below. 
+  ! AXPBY is invoked via Y, hence the structure below.
   !
   !
   !
@@ -2156,7 +2174,7 @@ contains
   !!
   subroutine d_base_mlv_axpby_v(m,alpha, x, beta, y, info, n)
     use psi_serial_mod
-    implicit none 
+    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
@@ -2180,7 +2198,7 @@ contains
   end subroutine d_base_mlv_axpby_v
 
   !
-  ! AXPBY is invoked via Y, hence the structure below. 
+  ! AXPBY is invoked via Y, hence the structure below.
   !
   !
   !> Function  base_mlv_axpby_a
@@ -2194,7 +2212,7 @@ contains
   !!
   subroutine d_base_mlv_axpby_a(m,alpha, x, beta, y, info,n)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_), intent(in)               :: m
     real(psb_dpk_), intent(in)        :: x(:,:)
     class(psb_d_base_multivect_type), intent(inout)  :: y
@@ -2230,10 +2248,10 @@ contains
   !!
   subroutine d_base_mlv_mlt_mv(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     class(psb_d_base_multivect_type), intent(inout)  :: y
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
 
     info = 0
     if (x%is_dev()) call x%sync()
@@ -2243,10 +2261,10 @@ contains
 
   subroutine d_base_mlv_mlt_mv_v(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_vect_type), intent(inout)  :: x
     class(psb_d_base_multivect_type), intent(inout)  :: y
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
 
     info = 0
     if (x%is_dev()) call x%sync()
@@ -2263,7 +2281,7 @@ contains
   !!
   subroutine d_base_mlv_mlt_ar1(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: x(:)
     class(psb_d_base_multivect_type), intent(inout)  :: y
     integer(psb_ipk_), intent(out)              :: info
@@ -2271,7 +2289,7 @@ contains
 
     info = 0
     n = min(psb_size(y%v,1_psb_ipk_), size(x))
-    do i=1, n 
+    do i=1, n
       y%v(i,:) = y%v(i,:)*x(i)
     end do
 
@@ -2286,7 +2304,7 @@ contains
   !!
   subroutine d_base_mlv_mlt_ar2(x, y, info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: x(:,:)
     class(psb_d_base_multivect_type), intent(inout)  :: y
     integer(psb_ipk_), intent(out)              :: info
@@ -2313,7 +2331,7 @@ contains
   !!
   subroutine d_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
     use psi_serial_mod
-    implicit none 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     real(psb_dpk_), intent(in)        :: y(:,:)
     real(psb_dpk_), intent(in)        :: x(:,:)
@@ -2321,38 +2339,38 @@ contains
     integer(psb_ipk_), intent(out)              :: info
     integer(psb_ipk_) :: i, nr, nc
 
-    info = 0    
+    info = 0
     nr = min(psb_size(z%v,1_psb_ipk_), size(x,1), size(y,1))
-    nc = min(psb_size(z%v,2_psb_ipk_), size(x,2), size(y,2))    
-    if (alpha == dzero) then 
-      if (beta == done) then 
-        return 
+    nc = min(psb_size(z%v,2_psb_ipk_), size(x,2), size(y,2))
+    if (alpha == dzero) then
+      if (beta == done) then
+        return
       else
         z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc)
       end if
     else
-      if (alpha == done) then 
-        if (beta == dzero) then 
+      if (alpha == done) then
+        if (beta == dzero) then
           z%v(1:nr,1:nc) = y(1:nr,1:nc)*x(1:nr,1:nc)
-        else if (beta == done) then 
+        else if (beta == done) then
           z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
-        else 
+        else
           z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
         end if
-      else if (alpha == -done) then 
-        if (beta == dzero) then 
+      else if (alpha == -done) then
+        if (beta == dzero) then
           z%v(1:nr,1:nc) = -y(1:nr,1:nc)*x(1:nr,1:nc)
-        else if (beta == done) then 
+        else if (beta == done) then
           z%v(1:nr,1:nc) = z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
-        else 
+        else
           z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
         end if
       else
-        if (beta == dzero) then 
+        if (beta == dzero) then
           z%v(1:nr,1:nc) = alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
-        else if (beta == done) then 
+        else if (beta == done) then
           z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
-        else 
+        else
           z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
         end if
       end if
@@ -2373,12 +2391,12 @@ contains
   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 
+    implicit none
     real(psb_dpk_), intent(in)        :: alpha,beta
     class(psb_d_base_multivect_type), intent(inout)  :: x
     class(psb_d_base_multivect_type), intent(inout)  :: y
     class(psb_d_base_multivect_type), intent(inout)  :: z
-    integer(psb_ipk_), intent(out)              :: info    
+    integer(psb_ipk_), intent(out)              :: info
     character(len=1), intent(in), optional     :: conjgx, conjgy
     integer(psb_ipk_) :: i, n
     logical :: conjgx_, conjgy_
@@ -2389,7 +2407,7 @@ contains
     if (z%is_dev()) call z%sync()
     if (.not.psb_d_is_complex_) then
       call z%mlt(alpha,x%v,y%v,beta,info)
-    else 
+    else
       conjgx_=.false.
       if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
       conjgy_=.false.
@@ -2404,39 +2422,39 @@ contains
 !!$
 !!$  subroutine d_base_mlv_mlt_av(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
-!!$    implicit none 
+!!$    implicit none
 !!$    real(psb_dpk_), intent(in)        :: alpha,beta
 !!$    real(psb_dpk_), intent(in)        :: x(:)
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: y
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: z
-!!$    integer(psb_ipk_), intent(out)              :: info    
+!!$    integer(psb_ipk_), intent(out)              :: info
 !!$    integer(psb_ipk_) :: i, n
 !!$
 !!$    info = 0
-!!$    
+!!$
 !!$    call z%mlt(alpha,x,y%v,beta,info)
 !!$
 !!$  end subroutine d_base_mlv_mlt_av
 !!$
 !!$  subroutine d_base_mlv_mlt_va(alpha,x,y,beta,z,info)
 !!$    use psi_serial_mod
-!!$    implicit none 
+!!$    implicit none
 !!$    real(psb_dpk_), intent(in)        :: alpha,beta
 !!$    real(psb_dpk_), intent(in)        :: y(:)
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: x
 !!$    class(psb_d_base_multivect_type), intent(inout)  :: z
-!!$    integer(psb_ipk_), intent(out)              :: info    
+!!$    integer(psb_ipk_), intent(out)              :: info
 !!$    integer(psb_ipk_) :: i, n
 !!$
 !!$    info = 0
-!!$    
+!!$
 !!$    call z%mlt(alpha,y,x,beta,info)
 !!$
 !!$  end subroutine d_base_mlv_mlt_va
 !!$
 !!$
   !
-  ! Simple scaling 
+  ! Simple scaling
   !
   !> Function  base_mlv_scal
   !! \memberof  psb_d_base_multivect_type
@@ -2445,7 +2463,7 @@ contains
   !!
   subroutine d_base_mlv_scal(alpha, x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
     real(psb_dpk_), intent (in)       :: alpha
 
@@ -2462,7 +2480,7 @@ contains
   !! \brief 2-norm |x(1:n)|_2
   !! \param n  how many entries to consider
   function d_base_mlv_nrm2(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_), allocatable    :: res(:)
@@ -2484,7 +2502,7 @@ contains
   !! \brief infinity-norm |x(1:n)|_\infty
   !! \param n  how many entries to consider
   function d_base_mlv_amax(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_), allocatable    :: res(:)
@@ -2505,7 +2523,7 @@ contains
   !! \brief 1-norm |x(1:n)|_1
   !! \param n  how many entries to consider
   function d_base_mlv_asum(n,x) result(res)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)           :: n
     real(psb_dpk_), allocatable    :: res(:)
@@ -2528,7 +2546,7 @@ contains
   !! \brief  Set all entries to their respective absolute values.
   !!
   subroutine d_base_mlv_absval1(x)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout)  :: x
 
     if (allocated(x%v)) then
@@ -2540,13 +2558,13 @@ contains
   end subroutine d_base_mlv_absval1
 
   subroutine d_base_mlv_absval2(x,y)
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     class(psb_d_base_multivect_type), intent(inout) :: y
     integer(psb_ipk_) :: info
-    
+
     if (x%is_dev()) call x%sync()
-    if (allocated(x%v)) then 
+    if (allocated(x%v)) then
       call y%axpby(min(x%get_nrows(),y%get_nrows()),done,x,dzero,info)
       call y%absval()
     end if
@@ -2555,15 +2573,15 @@ contains
 
 
   function d_base_mlv_use_buffer() result(res)
-    implicit none 
+    implicit none
     logical :: res
-    
+
     res = .true.
   end function d_base_mlv_use_buffer
 
   subroutine d_base_mlv_new_buffer(n,x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)              :: n
     integer(psb_ipk_), intent(out)             :: info
@@ -2575,7 +2593,7 @@ contains
 
   subroutine d_base_mlv_new_comid(n,x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(in)              :: n
     integer(psb_ipk_), intent(out)             :: info
@@ -2586,12 +2604,12 @@ contains
 
   subroutine d_base_mlv_maybe_free_buffer(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
 
-    info = 0 
+    info = 0
     if (psb_get_maybe_free_buffer())&
          &  call x%free_buffer(info)
 
@@ -2599,7 +2617,7 @@ contains
 
   subroutine d_base_mlv_free_buffer(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
@@ -2609,7 +2627,7 @@ contains
 
   subroutine d_base_mlv_free_comid(x,info)
     use psb_realloc_mod
-    implicit none 
+    implicit none
     class(psb_d_base_multivect_type), intent(inout) :: x
     integer(psb_ipk_), intent(out)             :: info
 
@@ -2632,7 +2650,7 @@ contains
   !! \param beta
   subroutine d_base_mlv_gthab(n,idx,alpha,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) :: alpha, beta, y(:)
     class(psb_d_base_multivect_type) :: x
@@ -2648,7 +2666,7 @@ contains
   end subroutine d_base_mlv_gthab
   !
   ! shortcut alpha=1 beta=0
-  ! 
+  !
   !> Function  base_mlv_gthzv
   !! \memberof  psb_d_base_multivect_type
   !! \brief gather into an array special alpha=1 beta=0
@@ -2657,7 +2675,7 @@ contains
   !! \param idx(:) indices
   subroutine d_base_mlv_gthzv_x(i,n,idx,x,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i,n
     class(psb_i_base_vect_type) :: idx
     real(psb_dpk_) ::  y(:)
@@ -2670,7 +2688,7 @@ contains
 
   !
   ! shortcut alpha=1 beta=0
-  ! 
+  !
   !> Function  base_mlv_gthzv
   !! \memberof  psb_d_base_multivect_type
   !! \brief gather into an array special alpha=1 beta=0
@@ -2679,7 +2697,7 @@ contains
   !! \param idx(:) indices
   subroutine d_base_mlv_gthzv(n,idx,x,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) ::  y(:)
     class(psb_d_base_multivect_type) :: x
@@ -2696,7 +2714,7 @@ contains
   end subroutine d_base_mlv_gthzv
   !
   ! shortcut alpha=1 beta=0
-  ! 
+  !
   !> Function  base_mlv_gthzv
   !! \memberof  psb_d_base_multivect_type
   !! \brief gather into an array special alpha=1 beta=0
@@ -2705,7 +2723,7 @@ contains
   !! \param idx(:) indices
   subroutine d_base_mlv_gthzm(n,idx,x,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) ::  y(:,:)
     class(psb_d_base_multivect_type) :: x
@@ -2722,17 +2740,17 @@ contains
   end subroutine d_base_mlv_gthzm
 
   !
-  ! New comm internals impl. 
+  ! New comm internals impl.
   !
   subroutine d_base_mlv_gthzbuf(i,ixb,n,idx,x)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i, ixb, n
     class(psb_i_base_vect_type) :: idx
     class(psb_d_base_multivect_type) :: x
     integer(psb_ipk_) :: nc
-    
-    if (.not.allocated(x%combuf)) then 
+
+    if (.not.allocated(x%combuf)) then
       call psb_errpush(psb_err_alloc_dealloc_,'gthzbuf')
       return
     end if
@@ -2744,9 +2762,9 @@ contains
   end subroutine d_base_mlv_gthzbuf
 
   !
-  ! Scatter: 
+  ! Scatter:
   ! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
-  ! 
+  !
   !
   !> Function  base_mlv_sctb
   !! \memberof  psb_d_base_multivect_type
@@ -2755,10 +2773,10 @@ contains
   !! \param n  how many entries to consider
   !! \param idx(:) indices
   !! \param beta
-  !! \param x(:) 
+  !! \param x(:)
   subroutine d_base_mlv_sctb(n,idx,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) :: beta, x(:)
     class(psb_d_base_multivect_type) :: y
@@ -2773,7 +2791,7 @@ contains
 
   subroutine d_base_mlv_sctbr2(n,idx,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: n, idx(:)
     real(psb_dpk_) :: beta, x(:,:)
     class(psb_d_base_multivect_type) :: y
@@ -2788,7 +2806,7 @@ contains
 
   subroutine d_base_mlv_sctb_x(i,n,idx,x,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i, n
     class(psb_i_base_vect_type) :: idx
     real( psb_dpk_) :: beta, x(:)
@@ -2800,14 +2818,14 @@ contains
 
   subroutine d_base_mlv_sctb_buf(i,iyb,n,idx,beta,y)
     use psi_serial_mod
-    implicit none 
+    implicit none
     integer(psb_ipk_) :: i, iyb, n
     class(psb_i_base_vect_type) :: idx
     real(psb_dpk_) :: beta
     class(psb_d_base_multivect_type) :: y
     integer(psb_ipk_) :: nc
-    
-    if (.not.allocated(y%combuf)) then 
+
+    if (.not.allocated(y%combuf)) then
       call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
       return
     end if
@@ -2816,19 +2834,18 @@ contains
     nc = y%get_ncols()
     call y%sct(n,idx%v(i:),y%combuf(iyb:),beta)
     call y%set_host()
-    
+
   end subroutine d_base_mlv_sctb_buf
 
   !
   !> Function  base_device_wait:
   !! \memberof  psb_d_base_vect_type
   !! \brief device_wait: base version is a no-op.
-  !!           
+  !!
   !
   subroutine d_base_mlv_device_wait()
-    implicit none 
-    
+    implicit none
+
   end subroutine d_base_mlv_device_wait
 
 end module psb_d_base_multivect_mod
-