diff --git a/base/modules/Makefile b/base/modules/Makefile
index 0ccf1133..5356395c 100644
--- a/base/modules/Makefile
+++ b/base/modules/Makefile
@@ -23,7 +23,8 @@ UTIL_MODS = psb_string_mod.o psb_desc_const_mod.o psb_indx_map_mod.o\
 	psb_vect_mod.o\
 	psb_s_psblas_mod.o psb_c_psblas_mod.o \
 	psb_d_psblas_mod.o psb_z_psblas_mod.o psb_psblas_mod.o \
-	psi_serial_mod.o \
+	psi_serial_mod.o psi_i_serial_mod.o \
+	psi_s_serial_mod.o psi_d_serial_mod.o psi_c_serial_mod.o psi_z_serial_mod.o \
 	psi_mod.o psi_i_mod.o psi_s_mod.o psi_d_mod.o psi_c_mod.o psi_z_mod.o\
 	psb_ip_reord_mod.o\
 	psb_i_sort_mod.o psb_s_sort_mod.o psb_d_sort_mod.o \
@@ -140,6 +141,8 @@ psb_c_comm_mod.o: psb_c_vect_mod.o psb_desc_mod.o psb_mat_mod.o
 psb_z_comm_mod.o: psb_z_vect_mod.o psb_desc_mod.o psb_mat_mod.o 
 psb_sort_mod.o: psb_i_sort_mod.o psb_s_sort_mod.o psb_d_sort_mod.o \
 	psb_c_sort_mod.o psb_z_sort_mod.o psb_ip_reord_mod.o psi_serial_mod.o
+psi_serial_mod.o:  psi_i_serial_mod.o \
+	psi_s_serial_mod.o psi_d_serial_mod.o psi_c_serial_mod.o psi_z_serial_mod.o 
 
 psb_base_mod.o: $(MODULES)
 
diff --git a/base/modules/psb_c_base_vect_mod.f90 b/base/modules/psb_c_base_vect_mod.f90
index 0fb91efd..afe28f03 100644
--- a/base/modules/psb_c_base_vect_mod.f90
+++ b/base/modules/psb_c_base_vect_mod.f90
@@ -1439,7 +1439,7 @@ module psb_c_base_multivect_mod
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => c_base_mlv_gthab
+    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
@@ -2406,30 +2406,35 @@ contains
     end if
     
   end subroutine c_base_mlv_absval2
-!!$  
-!!$  
-!!$  !
-!!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
-!!$  !
-!!$  !
-!!$  !> Function  base_mlv_gthab
-!!$  !! \memberof  psb_c_base_multivect_type
-!!$  !! \brief gather into an array
-!!$  !!    Y = beta * Y + alpha * X(IDX(:))
-!!$  !! \param n  how many entries to consider
-!!$  !! \param idx(:) indices
-!!$  !! \param alpha
-!!$  !! \param beta
-!!$  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(:)
-!!$    class(psb_c_base_multivect_type) :: x
-!!$    
-!!$    call x%sync()
-!!$    call psi_gth(n,idx,alpha,x%v,beta,y)
-!!$
-!!$  end subroutine c_base_mlv_gthab
+  
+  
+  !
+  ! Gather: Y = beta * Y + alpha * X(IDX(:))
+  !
+  !
+  !> Function  base_mlv_gthab
+  !! \memberof  psb_c_base_multivect_type
+  !! \brief gather into an array
+  !!    Y = beta * Y + alpha * X(IDX(:))
+  !! \param n  how many entries to consider
+  !! \param idx(:) indices
+  !! \param alpha
+  !! \param beta
+  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(:)
+    class(psb_c_base_multivect_type) :: x
+    integer(psb_ipk_) :: nc
+
+    call x%sync()
+    if (.not.allocated(x%v)) then
+      return
+    end if
+    nc = psb_size(x%v,2)
+    call psi_gth(n,nc,idx,alpha,x%v,beta,y)
+    
+  end subroutine c_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
diff --git a/base/modules/psb_d_base_vect_mod.f90 b/base/modules/psb_d_base_vect_mod.f90
index dc061325..5892892a 100644
--- a/base/modules/psb_d_base_vect_mod.f90
+++ b/base/modules/psb_d_base_vect_mod.f90
@@ -1439,7 +1439,7 @@ module psb_d_base_multivect_mod
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => d_base_mlv_gthab
+    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
@@ -2406,30 +2406,35 @@ contains
     end if
     
   end subroutine d_base_mlv_absval2
-!!$  
-!!$  
-!!$  !
-!!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
-!!$  !
-!!$  !
-!!$  !> Function  base_mlv_gthab
-!!$  !! \memberof  psb_d_base_multivect_type
-!!$  !! \brief gather into an array
-!!$  !!    Y = beta * Y + alpha * X(IDX(:))
-!!$  !! \param n  how many entries to consider
-!!$  !! \param idx(:) indices
-!!$  !! \param alpha
-!!$  !! \param beta
-!!$  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(:)
-!!$    class(psb_d_base_multivect_type) :: x
-!!$    
-!!$    call x%sync()
-!!$    call psi_gth(n,idx,alpha,x%v,beta,y)
-!!$
-!!$  end subroutine d_base_mlv_gthab
+  
+  
+  !
+  ! Gather: Y = beta * Y + alpha * X(IDX(:))
+  !
+  !
+  !> Function  base_mlv_gthab
+  !! \memberof  psb_d_base_multivect_type
+  !! \brief gather into an array
+  !!    Y = beta * Y + alpha * X(IDX(:))
+  !! \param n  how many entries to consider
+  !! \param idx(:) indices
+  !! \param alpha
+  !! \param beta
+  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(:)
+    class(psb_d_base_multivect_type) :: x
+    integer(psb_ipk_) :: nc
+
+    call x%sync()
+    if (.not.allocated(x%v)) then
+      return
+    end if
+    nc = psb_size(x%v,2)
+    call psi_gth(n,nc,idx,alpha,x%v,beta,y)
+    
+  end subroutine d_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
diff --git a/base/modules/psb_s_base_vect_mod.f90 b/base/modules/psb_s_base_vect_mod.f90
index 457c1e7d..a64db92d 100644
--- a/base/modules/psb_s_base_vect_mod.f90
+++ b/base/modules/psb_s_base_vect_mod.f90
@@ -1439,7 +1439,7 @@ module psb_s_base_multivect_mod
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => s_base_mlv_gthab
+    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
@@ -2406,30 +2406,35 @@ contains
     end if
     
   end subroutine s_base_mlv_absval2
-!!$  
-!!$  
-!!$  !
-!!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
-!!$  !
-!!$  !
-!!$  !> Function  base_mlv_gthab
-!!$  !! \memberof  psb_s_base_multivect_type
-!!$  !! \brief gather into an array
-!!$  !!    Y = beta * Y + alpha * X(IDX(:))
-!!$  !! \param n  how many entries to consider
-!!$  !! \param idx(:) indices
-!!$  !! \param alpha
-!!$  !! \param beta
-!!$  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(:)
-!!$    class(psb_s_base_multivect_type) :: x
-!!$    
-!!$    call x%sync()
-!!$    call psi_gth(n,idx,alpha,x%v,beta,y)
-!!$
-!!$  end subroutine s_base_mlv_gthab
+  
+  
+  !
+  ! Gather: Y = beta * Y + alpha * X(IDX(:))
+  !
+  !
+  !> Function  base_mlv_gthab
+  !! \memberof  psb_s_base_multivect_type
+  !! \brief gather into an array
+  !!    Y = beta * Y + alpha * X(IDX(:))
+  !! \param n  how many entries to consider
+  !! \param idx(:) indices
+  !! \param alpha
+  !! \param beta
+  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(:)
+    class(psb_s_base_multivect_type) :: x
+    integer(psb_ipk_) :: nc
+
+    call x%sync()
+    if (.not.allocated(x%v)) then
+      return
+    end if
+    nc = psb_size(x%v,2)
+    call psi_gth(n,nc,idx,alpha,x%v,beta,y)
+    
+  end subroutine s_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
diff --git a/base/modules/psb_z_base_vect_mod.f90 b/base/modules/psb_z_base_vect_mod.f90
index ce1785fc..21bda8ce 100644
--- a/base/modules/psb_z_base_vect_mod.f90
+++ b/base/modules/psb_z_base_vect_mod.f90
@@ -1439,7 +1439,7 @@ module psb_z_base_multivect_mod
 !!$    ! Gather/scatter. These are needed for MPI interfacing.
 !!$    ! May have to be reworked. 
 !!$    !
-!!$    procedure, pass(x) :: gthab    => z_base_mlv_gthab
+    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
@@ -2406,30 +2406,35 @@ contains
     end if
     
   end subroutine z_base_mlv_absval2
-!!$  
-!!$  
-!!$  !
-!!$  ! Gather: Y = beta * Y + alpha * X(IDX(:))
-!!$  !
-!!$  !
-!!$  !> Function  base_mlv_gthab
-!!$  !! \memberof  psb_z_base_multivect_type
-!!$  !! \brief gather into an array
-!!$  !!    Y = beta * Y + alpha * X(IDX(:))
-!!$  !! \param n  how many entries to consider
-!!$  !! \param idx(:) indices
-!!$  !! \param alpha
-!!$  !! \param beta
-!!$  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(:)
-!!$    class(psb_z_base_multivect_type) :: x
-!!$    
-!!$    call x%sync()
-!!$    call psi_gth(n,idx,alpha,x%v,beta,y)
-!!$
-!!$  end subroutine z_base_mlv_gthab
+  
+  
+  !
+  ! Gather: Y = beta * Y + alpha * X(IDX(:))
+  !
+  !
+  !> Function  base_mlv_gthab
+  !! \memberof  psb_z_base_multivect_type
+  !! \brief gather into an array
+  !!    Y = beta * Y + alpha * X(IDX(:))
+  !! \param n  how many entries to consider
+  !! \param idx(:) indices
+  !! \param alpha
+  !! \param beta
+  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(:)
+    class(psb_z_base_multivect_type) :: x
+    integer(psb_ipk_) :: nc
+
+    call x%sync()
+    if (.not.allocated(x%v)) then
+      return
+    end if
+    nc = psb_size(x%v,2)
+    call psi_gth(n,nc,idx,alpha,x%v,beta,y)
+    
+  end subroutine z_base_mlv_gthab
 !!$  !
 !!$  ! shortcut alpha=1 beta=0
 !!$  ! 
diff --git a/base/modules/psi_c_serial_mod.f90 b/base/modules/psi_c_serial_mod.f90
new file mode 100644
index 00000000..ef6f6994
--- /dev/null
+++ b/base/modules/psi_c_serial_mod.f90
@@ -0,0 +1,89 @@
+module psi_c_serial_mod
+  use psb_const_mod, only : psb_ipk_, psb_spk_
+
+  interface psb_gelp
+    ! 2-D version
+    subroutine psb_cgelp(trans,iperm,x,info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      complex(psb_spk_), intent(inout)     ::  x(:,:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_cgelp
+    subroutine psb_cgelpv(trans,iperm,x,info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      complex(psb_spk_), intent(inout)     ::  x(:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_cgelpv
+  end interface psb_gelp
+  
+  interface psb_geaxpby 
+    subroutine psi_caxpby(m,n,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m, n
+      complex(psb_spk_), intent (in)       ::  x(:,:)
+      complex(psb_spk_), intent (inout)    ::  y(:,:)
+      complex(psb_spk_), intent (in)       ::  alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_caxpby
+    subroutine psi_caxpbyv(m,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m
+      complex(psb_spk_), intent (in)       ::  x(:)
+      complex(psb_spk_), intent (inout)    ::  y(:)
+      complex(psb_spk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_caxpbyv
+  end interface psb_geaxpby
+
+  interface psi_gth
+    subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_spk_) :: x(:,:), y(:),alpha,beta
+    end subroutine psi_cgthmv
+    subroutine psi_cgthv(n,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_spk_) :: x(:), y(:),alpha,beta
+    end subroutine psi_cgthv
+    subroutine psi_cgthzmv(n,k,idx,x,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_spk_) :: x(:,:), y(:)
+      
+    end subroutine psi_cgthzmv
+    subroutine psi_cgthzv(n,idx,x,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_spk_) :: x(:), y(:)
+    end subroutine psi_cgthzv
+  end interface psi_gth
+
+  interface psi_sct
+    subroutine psi_csctmv(n,k,idx,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_spk_) :: beta, x(:), y(:,:)
+    end subroutine psi_csctmv
+    subroutine psi_csctv(n,idx,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_spk_) :: beta, x(:), y(:)
+    end subroutine psi_csctv
+  end interface psi_sct
+
+end module psi_c_serial_mod
diff --git a/base/modules/psi_d_serial_mod.f90 b/base/modules/psi_d_serial_mod.f90
new file mode 100644
index 00000000..1eb1bc36
--- /dev/null
+++ b/base/modules/psi_d_serial_mod.f90
@@ -0,0 +1,89 @@
+module psi_d_serial_mod
+  use psb_const_mod, only : psb_ipk_, psb_dpk_
+
+  interface psb_gelp
+    ! 2-D version
+    subroutine psb_dgelp(trans,iperm,x,info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      real(psb_dpk_), intent(inout)     ::  x(:,:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_dgelp
+    subroutine psb_dgelpv(trans,iperm,x,info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      real(psb_dpk_), intent(inout)     ::  x(:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_dgelpv
+  end interface psb_gelp
+  
+  interface psb_geaxpby 
+    subroutine psi_daxpby(m,n,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m, n
+      real(psb_dpk_), intent (in)       ::  x(:,:)
+      real(psb_dpk_), intent (inout)    ::  y(:,:)
+      real(psb_dpk_), intent (in)       ::  alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_daxpby
+    subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m
+      real(psb_dpk_), intent (in)       ::  x(:)
+      real(psb_dpk_), intent (inout)    ::  y(:)
+      real(psb_dpk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_daxpbyv
+  end interface psb_geaxpby
+
+  interface psi_gth
+    subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_dpk_) :: x(:,:), y(:),alpha,beta
+    end subroutine psi_dgthmv
+    subroutine psi_dgthv(n,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_dpk_) :: x(:), y(:),alpha,beta
+    end subroutine psi_dgthv
+    subroutine psi_dgthzmv(n,k,idx,x,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_dpk_) :: x(:,:), y(:)
+      
+    end subroutine psi_dgthzmv
+    subroutine psi_dgthzv(n,idx,x,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_dpk_) :: x(:), y(:)
+    end subroutine psi_dgthzv
+  end interface psi_gth
+
+  interface psi_sct
+    subroutine psi_dsctmv(n,k,idx,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_dpk_) :: beta, x(:), y(:,:)
+    end subroutine psi_dsctmv
+    subroutine psi_dsctv(n,idx,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_dpk_) :: beta, x(:), y(:)
+    end subroutine psi_dsctv
+  end interface psi_sct
+
+end module psi_d_serial_mod
diff --git a/base/modules/psi_i_serial_mod.f90 b/base/modules/psi_i_serial_mod.f90
new file mode 100644
index 00000000..f5274150
--- /dev/null
+++ b/base/modules/psi_i_serial_mod.f90
@@ -0,0 +1,89 @@
+module psi_i_serial_mod
+  use psb_const_mod, only : psb_ipk_
+
+  interface psb_gelp
+    ! 2-D version
+    subroutine psb_igelp(trans,iperm,x,info)
+      import :: psb_ipk_
+      implicit none 
+      integer(psb_ipk_), intent(inout)     ::  x(:,:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_igelp
+    subroutine psb_igelpv(trans,iperm,x,info)
+      import :: psb_ipk_
+      implicit none 
+      integer(psb_ipk_), intent(inout)     ::  x(:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_igelpv
+  end interface psb_gelp
+  
+  interface psb_geaxpby 
+    subroutine psi_iaxpby(m,n,alpha, x, beta, y, info)
+      import :: psb_ipk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m, n
+      integer(psb_ipk_), intent (in)       ::  x(:,:)
+      integer(psb_ipk_), intent (inout)    ::  y(:,:)
+      integer(psb_ipk_), intent (in)       ::  alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_iaxpby
+    subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)
+      import :: psb_ipk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m
+      integer(psb_ipk_), intent (in)       ::  x(:)
+      integer(psb_ipk_), intent (inout)    ::  y(:)
+      integer(psb_ipk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_iaxpbyv
+  end interface psb_geaxpby
+
+  interface psi_gth
+    subroutine psi_igthmv(n,k,idx,alpha,x,beta,y)
+      import :: psb_ipk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      integer(psb_ipk_) :: x(:,:), y(:),alpha,beta
+    end subroutine psi_igthmv
+    subroutine psi_igthv(n,idx,alpha,x,beta,y)
+      import :: psb_ipk_
+      implicit none
+      integer(psb_ipk_) :: n, idx(:)
+      integer(psb_ipk_) :: x(:), y(:),alpha,beta
+    end subroutine psi_igthv
+    subroutine psi_igthzmv(n,k,idx,x,y)
+      import :: psb_ipk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      integer(psb_ipk_) :: x(:,:), y(:)
+      
+    end subroutine psi_igthzmv
+    subroutine psi_igthzv(n,idx,x,y)
+      import :: psb_ipk_
+      implicit none 
+      integer(psb_ipk_) :: n, idx(:)
+      integer(psb_ipk_) :: x(:), y(:)
+    end subroutine psi_igthzv
+  end interface psi_gth
+
+  interface psi_sct
+    subroutine psi_isctmv(n,k,idx,x,beta,y)
+      import :: psb_ipk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      integer(psb_ipk_) :: beta, x(:), y(:,:)
+    end subroutine psi_isctmv
+    subroutine psi_isctv(n,idx,x,beta,y)
+      import :: psb_ipk_
+      implicit none
+      
+      integer(psb_ipk_) :: n, idx(:)
+      integer(psb_ipk_) :: beta, x(:), y(:)
+    end subroutine psi_isctv
+  end interface psi_sct
+
+end module psi_i_serial_mod
diff --git a/base/modules/psi_s_serial_mod.f90 b/base/modules/psi_s_serial_mod.f90
new file mode 100644
index 00000000..64dc604e
--- /dev/null
+++ b/base/modules/psi_s_serial_mod.f90
@@ -0,0 +1,89 @@
+module psi_s_serial_mod
+  use psb_const_mod, only : psb_ipk_, psb_spk_
+
+  interface psb_gelp
+    ! 2-D version
+    subroutine psb_sgelp(trans,iperm,x,info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      real(psb_spk_), intent(inout)     ::  x(:,:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_sgelp
+    subroutine psb_sgelpv(trans,iperm,x,info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      real(psb_spk_), intent(inout)     ::  x(:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_sgelpv
+  end interface psb_gelp
+  
+  interface psb_geaxpby 
+    subroutine psi_saxpby(m,n,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m, n
+      real(psb_spk_), intent (in)       ::  x(:,:)
+      real(psb_spk_), intent (inout)    ::  y(:,:)
+      real(psb_spk_), intent (in)       ::  alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_saxpby
+    subroutine psi_saxpbyv(m,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m
+      real(psb_spk_), intent (in)       ::  x(:)
+      real(psb_spk_), intent (inout)    ::  y(:)
+      real(psb_spk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_saxpbyv
+  end interface psb_geaxpby
+
+  interface psi_gth
+    subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_spk_) :: x(:,:), y(:),alpha,beta
+    end subroutine psi_sgthmv
+    subroutine psi_sgthv(n,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_spk_) :: x(:), y(:),alpha,beta
+    end subroutine psi_sgthv
+    subroutine psi_sgthzmv(n,k,idx,x,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_spk_) :: x(:,:), y(:)
+      
+    end subroutine psi_sgthzmv
+    subroutine psi_sgthzv(n,idx,x,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none 
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_spk_) :: x(:), y(:)
+    end subroutine psi_sgthzv
+  end interface psi_gth
+
+  interface psi_sct
+    subroutine psi_ssctmv(n,k,idx,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      real(psb_spk_) :: beta, x(:), y(:,:)
+    end subroutine psi_ssctmv
+    subroutine psi_ssctv(n,idx,x,beta,y)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      
+      integer(psb_ipk_) :: n, idx(:)
+      real(psb_spk_) :: beta, x(:), y(:)
+    end subroutine psi_ssctv
+  end interface psi_sct
+
+end module psi_s_serial_mod
diff --git a/base/modules/psi_serial_mod.f90 b/base/modules/psi_serial_mod.f90
index ca793999..d985e14d 100644
--- a/base/modules/psi_serial_mod.f90
+++ b/base/modules/psi_serial_mod.f90
@@ -30,290 +30,9 @@
 !!$ 
 !!$  
 module psi_serial_mod
-  use psb_const_mod, only : psb_ipk_, psb_spk_, psb_dpk_
-  interface psb_gelp
-    ! 2-D version
-    subroutine psb_sgelp(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      real(psb_spk_), intent(inout)      ::  x(:,:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)                :: trans
-    end subroutine psb_sgelp
-    ! 1-D version
-    subroutine psb_sgelpv(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      real(psb_spk_), intent(inout)    ::  x(:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)              :: trans
-    end subroutine psb_sgelpv
-    subroutine psb_dgelp(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      real(psb_dpk_), intent(inout)      ::  x(:,:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)                :: trans
-    end subroutine psb_dgelp
-    ! 1-D version
-    subroutine psb_dgelpv(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      real(psb_dpk_), intent(inout)    ::  x(:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)              :: trans
-    end subroutine psb_dgelpv
-    ! 2-D version
-    subroutine psb_cgelp(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      complex(psb_spk_), intent(inout)      ::  x(:,:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)                :: trans
-    end subroutine psb_cgelp
-    ! 1-D version
-    subroutine psb_cgelpv(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      complex(psb_spk_), intent(inout)    ::  x(:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)              :: trans
-    end subroutine psb_cgelpv
-    ! 2-D version
-    subroutine psb_zgelp(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      complex(psb_dpk_), intent(inout)      ::  x(:,:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)                :: trans
-    end subroutine psb_zgelp
-    ! 1-D version
-    subroutine psb_zgelpv(trans,iperm,x,info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      complex(psb_dpk_), intent(inout)    ::  x(:)
-      integer(psb_ipk_), intent(in)                  ::  iperm(:)
-      integer(psb_ipk_), intent(out)                 ::  info
-      character, intent(in)              :: trans
-    end subroutine psb_zgelpv
-  end interface
-
-
-
-  interface psi_gth
-    subroutine psi_igthv(n,idx,alpha,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_ 
-      integer(psb_ipk_) :: n, idx(:)
-      integer(psb_ipk_) :: x(:), y(:), alpha, beta
-    end subroutine psi_igthv
-    subroutine psi_sgthv(n,idx,alpha,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_spk_) :: x(:), y(:), alpha, beta
-    end subroutine psi_sgthv
-    subroutine psi_dgthv(n,idx,alpha,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_dpk_) :: x(:), y(:), alpha, beta
-    end subroutine psi_dgthv
-    subroutine psi_cgthv(n,idx,alpha,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_spk_) :: x(:), y(:),alpha,beta
-    end subroutine psi_cgthv
-    subroutine psi_zgthv(n,idx,alpha,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_dpk_) :: x(:), y(:),alpha,beta
-    end subroutine psi_zgthv
-    subroutine psi_sgthzmv(n,k,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      real(psb_spk_) :: x(:,:), y(:)
-    end subroutine psi_sgthzmv
-    subroutine psi_dgthzmv(n,k,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      real(psb_dpk_) :: x(:,:), y(:)
-    end subroutine psi_dgthzmv
-    subroutine psi_igthzmv(n,k,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      integer(psb_ipk_) :: x(:,:), y(:)
-    end subroutine psi_igthzmv
-    subroutine psi_cgthzmv(n,k,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      complex(psb_spk_) :: x(:,:), y(:)
-    end subroutine psi_cgthzmv
-    subroutine psi_zgthzmv(n,k,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      complex(psb_dpk_) :: x(:,:), y(:)
-    end subroutine psi_zgthzmv
-    subroutine psi_sgthzv(n,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_spk_) :: x(:), y(:)
-    end subroutine psi_sgthzv
-    subroutine psi_dgthzv(n,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_dpk_) :: x(:), y(:)
-    end subroutine psi_dgthzv
-    subroutine psi_igthzv(n,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      integer(psb_ipk_) :: x(:), y(:)
-    end subroutine psi_igthzv
-    subroutine psi_cgthzv(n,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_spk_) :: x(:), y(:)
-    end subroutine psi_cgthzv
-    subroutine psi_zgthzv(n,idx,x,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_dpk_) :: x(:), y(:)
-    end subroutine psi_zgthzv
-  end interface
-
-
-  interface psi_sct
-    subroutine psi_ssctmv(n,k,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      real(psb_spk_) :: beta, x(:), y(:,:)
-    end subroutine psi_ssctmv
-    subroutine psi_ssctv(n,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_spk_) :: beta, x(:), y(:)
-    end subroutine psi_ssctv
-    subroutine psi_dsctmv(n,k,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      real(psb_dpk_) :: beta, x(:), y(:,:)
-    end subroutine psi_dsctmv
-    subroutine psi_dsctv(n,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      real(psb_dpk_) :: beta, x(:), y(:)
-    end subroutine psi_dsctv
-    subroutine psi_isctmv(n,k,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      integer(psb_ipk_) :: beta, x(:), y(:,:)
-    end subroutine psi_isctmv
-    subroutine psi_isctv(n,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      integer(psb_ipk_) :: beta, x(:), y(:)
-    end subroutine psi_isctv
-    subroutine psi_csctmv(n,k,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      complex(psb_spk_) :: beta, x(:), y(:,:)
-    end subroutine psi_csctmv
-    subroutine psi_csctv(n,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_spk_) :: beta, x(:), y(:)
-    end subroutine psi_csctv
-    subroutine psi_zsctmv(n,k,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, k, idx(:)
-      complex(psb_dpk_) :: beta, x(:), y(:,:)
-    end subroutine psi_zsctmv
-    subroutine psi_zsctv(n,idx,x,beta,y)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_) :: n, idx(:)
-      complex(psb_dpk_) :: beta, x(:), y(:)
-    end subroutine psi_zsctv
-  end interface
-
-
-  interface psb_geaxpby
-    subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m
-      integer(psb_ipk_), intent (in)       ::  x(:)
-      integer(psb_ipk_), intent (inout)    ::  y(:)
-      integer(psb_ipk_), intent (in)       :: alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_iaxpbyv
-    subroutine psi_iaxpby(m,n,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m, n
-      integer(psb_ipk_), intent (in)       ::  x(:,:)
-      integer(psb_ipk_), intent (inout)    ::  y(:,:)
-      integer(psb_ipk_), intent (in)       ::  alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_iaxpby
-    subroutine psi_saxpbyv(m,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m
-      real(psb_spk_), intent (in)       ::  x(:)
-      real(psb_spk_), intent (inout)    ::  y(:)
-      real(psb_spk_), intent (in)       :: alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_saxpbyv
-    subroutine psi_saxpby(m,n,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m, n
-      real(psb_spk_), intent (in)       ::  x(:,:)
-      real(psb_spk_), intent (inout)    ::  y(:,:)
-      real(psb_spk_), intent (in)       ::  alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_saxpby
-    subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m
-      real(psb_dpk_), intent (in)       ::  x(:)
-      real(psb_dpk_), intent (inout)    ::  y(:)
-      real(psb_dpk_), intent (in)       :: alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_daxpbyv
-    subroutine psi_daxpby(m,n,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m, n
-      real(psb_dpk_), intent (in)       ::  x(:,:)
-      real(psb_dpk_), intent (inout)    ::  y(:,:)
-      real(psb_dpk_), intent (in)       ::  alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_daxpby
-    subroutine psi_caxpbyv(m,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m
-      complex(psb_spk_), intent (in)    ::  x(:)
-      complex(psb_spk_), intent (inout) ::  y(:)
-      complex(psb_spk_), intent (in)    :: alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_caxpbyv
-    subroutine psi_caxpby(m,n,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      implicit none 
-      integer(psb_ipk_), intent(in)               :: m, n
-      complex(psb_spk_), intent (in)    ::  x(:,:)
-      complex(psb_spk_), intent (inout) ::  y(:,:)
-      complex(psb_spk_), intent (in)    ::  alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_caxpby
-    subroutine psi_zaxpbyv(m,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_
-      integer(psb_ipk_), intent(in)               :: m
-      complex(psb_dpk_), intent (in)    ::  x(:)
-      complex(psb_dpk_), intent (inout) ::  y(:)
-      complex(psb_dpk_), intent (in)    :: alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_zaxpbyv
-    subroutine psi_zaxpby(m,n,alpha, x, beta, y, info)
-      import :: psb_ipk_, psb_spk_, psb_dpk_ 
-      integer(psb_ipk_), intent(in)               :: m, n
-      complex(psb_dpk_), intent (in)    ::  x(:,:)
-      complex(psb_dpk_), intent (inout) ::  y(:,:)
-      complex(psb_dpk_), intent (in)    ::  alpha, beta
-      integer(psb_ipk_), intent(out)              :: info
-    end subroutine psi_zaxpby
-  end interface
-
+  use psi_i_serial_mod
+  use psi_s_serial_mod
+  use psi_d_serial_mod
+  use psi_c_serial_mod
+  use psi_z_serial_mod
 end module psi_serial_mod
diff --git a/base/modules/psi_z_serial_mod.f90 b/base/modules/psi_z_serial_mod.f90
new file mode 100644
index 00000000..f1e37597
--- /dev/null
+++ b/base/modules/psi_z_serial_mod.f90
@@ -0,0 +1,89 @@
+module psi_z_serial_mod
+  use psb_const_mod, only : psb_ipk_, psb_dpk_
+
+  interface psb_gelp
+    ! 2-D version
+    subroutine psb_zgelp(trans,iperm,x,info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      complex(psb_dpk_), intent(inout)     ::  x(:,:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_zgelp
+    subroutine psb_zgelpv(trans,iperm,x,info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      complex(psb_dpk_), intent(inout)     ::  x(:)
+      integer(psb_ipk_), intent(in)      ::  iperm(:)
+      integer(psb_ipk_), intent(out)     ::  info
+      character, intent(in)              :: trans
+    end subroutine psb_zgelpv
+  end interface psb_gelp
+  
+  interface psb_geaxpby 
+    subroutine psi_zaxpby(m,n,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m, n
+      complex(psb_dpk_), intent (in)       ::  x(:,:)
+      complex(psb_dpk_), intent (inout)    ::  y(:,:)
+      complex(psb_dpk_), intent (in)       ::  alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_zaxpby
+    subroutine psi_zaxpbyv(m,alpha, x, beta, y, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_), intent(in)      :: m
+      complex(psb_dpk_), intent (in)       ::  x(:)
+      complex(psb_dpk_), intent (inout)    ::  y(:)
+      complex(psb_dpk_), intent (in)       :: alpha, beta
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_zaxpbyv
+  end interface psb_geaxpby
+
+  interface psi_gth
+    subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_dpk_) :: x(:,:), y(:),alpha,beta
+    end subroutine psi_zgthmv
+    subroutine psi_zgthv(n,idx,alpha,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_dpk_) :: x(:), y(:),alpha,beta
+    end subroutine psi_zgthv
+    subroutine psi_zgthzmv(n,k,idx,x,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_dpk_) :: x(:,:), y(:)
+      
+    end subroutine psi_zgthzmv
+    subroutine psi_zgthzv(n,idx,x,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none 
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_dpk_) :: x(:), y(:)
+    end subroutine psi_zgthzv
+  end interface psi_gth
+
+  interface psi_sct
+    subroutine psi_zsctmv(n,k,idx,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_) :: n, k, idx(:)
+      complex(psb_dpk_) :: beta, x(:), y(:,:)
+    end subroutine psi_zsctmv
+    subroutine psi_zsctv(n,idx,x,beta,y)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      
+      integer(psb_ipk_) :: n, idx(:)
+      complex(psb_dpk_) :: beta, x(:), y(:)
+    end subroutine psi_zsctv
+  end interface psi_sct
+
+end module psi_z_serial_mod
diff --git a/base/serial/Makefile b/base/serial/Makefile
index b0245f23..75a8c757 100644
--- a/base/serial/Makefile
+++ b/base/serial/Makefile
@@ -1,7 +1,9 @@
 include ../../Make.inc
 
 
-FOBJS = psb_lsame.o psi_serial_impl.o \
+FOBJS = psb_lsame.o psi_i_serial_impl.o \
+	psi_s_serial_impl.o psi_d_serial_impl.o \
+	psi_c_serial_impl.o psi_z_serial_impl.o \
 	psb_srwextd.o psb_drwextd.o psb_crwextd.o psb_zrwextd.o \
 	psb_sspspmm.o psb_dspspmm.o psb_cspspmm.o psb_zspspmm.o \
 	psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \
diff --git a/base/serial/psi_c_serial_impl.f90 b/base/serial/psi_c_serial_impl.f90
new file mode 100644
index 00000000..bacc90f2
--- /dev/null
+++ b/base/serial/psi_c_serial_impl.f90
@@ -0,0 +1,265 @@
+subroutine psi_caxpby(m,n,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m, n
+  complex(psb_spk_), intent (in)       ::  x(:,:)
+  complex(psb_spk_), intent (inout)    ::  y(:,:)
+  complex(psb_spk_), intent (in)       ::  alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (n < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 2; ierr(2) = n
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_
+    ierr(1) = 4; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 6; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psi_caxpby
+subroutine psi_caxpbyv(m,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m
+  complex(psb_spk_), intent (in)       ::  x(:)
+  complex(psb_spk_), intent (inout)    ::  y(:)
+  complex(psb_spk_), intent (in)       :: alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 3; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 5; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psi_caxpbyv
+subroutine psi_cgthv(n,idx,alpha,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_spk_) :: x(:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i
+  if (beta == czero) then 
+    if (alpha == czero) then 
+      do i=1,n
+        y(i) = czero
+      end do
+    else if (alpha == cone) then 
+      do i=1,n
+        y(i) = x(idx(i))
+      end do
+    else if (alpha == -cone) then 
+      do i=1,n
+        y(i) = -x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = alpha*x(idx(i))
+      end do
+    end if
+  else 
+    if (beta == cone) then 
+      ! Do nothing
+    else if (beta == -cone) then 
+      y(1:n) = -y(1:n) 
+    else
+      y(1:n) = beta*y(1:n) 
+    end if
+
+    if (alpha == czero) then 
+      ! do nothing
+    else if (alpha == cone) then 
+      do i=1,n
+        y(i) = y(i) + x(idx(i))
+      end do
+    else if (alpha == -cone) then 
+      do i=1,n
+        y(i) = y(i) - x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = y(i) + alpha*x(idx(i))
+      end do
+    end if
+  end if
+
+end subroutine psi_cgthv
+subroutine psi_cgthzmv(n,k,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  complex(psb_spk_) :: x(:,:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  pt=0
+  do j=1,k
+    do i=1,n
+      pt=pt+1
+      y(pt)=x(idx(i),j)
+    end do
+  end do
+
+end subroutine psi_cgthzmv
+subroutine psi_cgthzv(n,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_spk_) :: x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  do i=1,n
+    y(i)=x(idx(i))
+  end do
+
+end subroutine psi_cgthzv
+subroutine psi_csctmv(n,k,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  complex(psb_spk_) :: beta, x(:), y(:,:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == czero) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = x(pt)
+      end do
+    end do
+  else if (beta == cone) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = y(idx(i),j)+x(pt)
+      end do
+    end do
+  else
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
+      end do
+    end do
+  end if
+end subroutine psi_csctmv
+subroutine psi_csctv(n,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_spk_) :: beta, x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  if (beta == czero) then
+    do i=1,n
+      y(idx(i)) = x(i)
+    end do
+  else if (beta == cone) then
+    do i=1,n
+      y(idx(i)) = y(idx(i))+x(i)
+    end do
+  else
+    do i=1,n
+      y(idx(i)) = beta*y(idx(i))+x(i)
+    end do
+  end if
+end subroutine psi_csctv
diff --git a/base/serial/psi_d_serial_impl.f90 b/base/serial/psi_d_serial_impl.f90
new file mode 100644
index 00000000..38031e3a
--- /dev/null
+++ b/base/serial/psi_d_serial_impl.f90
@@ -0,0 +1,265 @@
+subroutine psi_daxpby(m,n,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m, n
+  real(psb_dpk_), intent (in)       ::  x(:,:)
+  real(psb_dpk_), intent (inout)    ::  y(:,:)
+  real(psb_dpk_), intent (in)       ::  alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (n < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 2; ierr(2) = n
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_
+    ierr(1) = 4; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 6; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psi_daxpby
+subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m
+  real(psb_dpk_), intent (in)       ::  x(:)
+  real(psb_dpk_), intent (inout)    ::  y(:)
+  real(psb_dpk_), intent (in)       :: alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 3; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 5; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psi_daxpbyv
+subroutine psi_dgthv(n,idx,alpha,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_dpk_) :: x(:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i
+  if (beta == dzero) then 
+    if (alpha == dzero) then 
+      do i=1,n
+        y(i) = dzero
+      end do
+    else if (alpha == done) then 
+      do i=1,n
+        y(i) = x(idx(i))
+      end do
+    else if (alpha == -done) then 
+      do i=1,n
+        y(i) = -x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = alpha*x(idx(i))
+      end do
+    end if
+  else 
+    if (beta == done) then 
+      ! Do nothing
+    else if (beta == -done) then 
+      y(1:n) = -y(1:n) 
+    else
+      y(1:n) = beta*y(1:n) 
+    end if
+
+    if (alpha == dzero) then 
+      ! do nothing
+    else if (alpha == done) then 
+      do i=1,n
+        y(i) = y(i) + x(idx(i))
+      end do
+    else if (alpha == -done) then 
+      do i=1,n
+        y(i) = y(i) - x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = y(i) + alpha*x(idx(i))
+      end do
+    end if
+  end if
+
+end subroutine psi_dgthv
+subroutine psi_dgthzmv(n,k,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  real(psb_dpk_) :: x(:,:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  pt=0
+  do j=1,k
+    do i=1,n
+      pt=pt+1
+      y(pt)=x(idx(i),j)
+    end do
+  end do
+
+end subroutine psi_dgthzmv
+subroutine psi_dgthzv(n,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_dpk_) :: x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  do i=1,n
+    y(i)=x(idx(i))
+  end do
+
+end subroutine psi_dgthzv
+subroutine psi_dsctmv(n,k,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  real(psb_dpk_) :: beta, x(:), y(:,:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == dzero) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = x(pt)
+      end do
+    end do
+  else if (beta == done) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = y(idx(i),j)+x(pt)
+      end do
+    end do
+  else
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
+      end do
+    end do
+  end if
+end subroutine psi_dsctmv
+subroutine psi_dsctv(n,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_dpk_) :: beta, x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  if (beta == dzero) then
+    do i=1,n
+      y(idx(i)) = x(i)
+    end do
+  else if (beta == done) then
+    do i=1,n
+      y(idx(i)) = y(idx(i))+x(i)
+    end do
+  else
+    do i=1,n
+      y(idx(i)) = beta*y(idx(i))+x(i)
+    end do
+  end if
+end subroutine psi_dsctv
diff --git a/base/serial/psi_i_serial_impl.f90 b/base/serial/psi_i_serial_impl.f90
new file mode 100644
index 00000000..6b227983
--- /dev/null
+++ b/base/serial/psi_i_serial_impl.f90
@@ -0,0 +1,265 @@
+subroutine psi_iaxpby(m,n,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m, n
+  integer(psb_ipk_), intent (in)       ::  x(:,:)
+  integer(psb_ipk_), intent (inout)    ::  y(:,:)
+  integer(psb_ipk_), intent (in)       ::  alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (n < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 2; ierr(2) = n
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_
+    ierr(1) = 4; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 6; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psi_iaxpby
+subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m
+  integer(psb_ipk_), intent (in)       ::  x(:)
+  integer(psb_ipk_), intent (inout)    ::  y(:)
+  integer(psb_ipk_), intent (in)       :: alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 3; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 5; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psi_iaxpbyv
+subroutine psi_igthv(n,idx,alpha,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  integer(psb_ipk_) :: x(:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i
+  if (beta == izero) then 
+    if (alpha == izero) then 
+      do i=1,n
+        y(i) = izero
+      end do
+    else if (alpha == ione) then 
+      do i=1,n
+        y(i) = x(idx(i))
+      end do
+    else if (alpha == -ione) then 
+      do i=1,n
+        y(i) = -x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = alpha*x(idx(i))
+      end do
+    end if
+  else 
+    if (beta == ione) then 
+      ! Do nothing
+    else if (beta == -ione) then 
+      y(1:n) = -y(1:n) 
+    else
+      y(1:n) = beta*y(1:n) 
+    end if
+
+    if (alpha == izero) then 
+      ! do nothing
+    else if (alpha == ione) then 
+      do i=1,n
+        y(i) = y(i) + x(idx(i))
+      end do
+    else if (alpha == -ione) then 
+      do i=1,n
+        y(i) = y(i) - x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = y(i) + alpha*x(idx(i))
+      end do
+    end if
+  end if
+
+end subroutine psi_igthv
+subroutine psi_igthzmv(n,k,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  integer(psb_ipk_) :: x(:,:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  pt=0
+  do j=1,k
+    do i=1,n
+      pt=pt+1
+      y(pt)=x(idx(i),j)
+    end do
+  end do
+
+end subroutine psi_igthzmv
+subroutine psi_igthzv(n,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  integer(psb_ipk_) :: x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  do i=1,n
+    y(i)=x(idx(i))
+  end do
+
+end subroutine psi_igthzv
+subroutine psi_isctmv(n,k,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  integer(psb_ipk_) :: beta, x(:), y(:,:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == izero) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = x(pt)
+      end do
+    end do
+  else if (beta == ione) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = y(idx(i),j)+x(pt)
+      end do
+    end do
+  else
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
+      end do
+    end do
+  end if
+end subroutine psi_isctmv
+subroutine psi_isctv(n,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  integer(psb_ipk_) :: beta, x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  if (beta == izero) then
+    do i=1,n
+      y(idx(i)) = x(i)
+    end do
+  else if (beta == ione) then
+    do i=1,n
+      y(idx(i)) = y(idx(i))+x(i)
+    end do
+  else
+    do i=1,n
+      y(idx(i)) = beta*y(idx(i))+x(i)
+    end do
+  end if
+end subroutine psi_isctv
diff --git a/base/serial/psi_s_serial_impl.f90 b/base/serial/psi_s_serial_impl.f90
new file mode 100644
index 00000000..aaae65f6
--- /dev/null
+++ b/base/serial/psi_s_serial_impl.f90
@@ -0,0 +1,265 @@
+subroutine psi_saxpby(m,n,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m, n
+  real(psb_spk_), intent (in)       ::  x(:,:)
+  real(psb_spk_), intent (inout)    ::  y(:,:)
+  real(psb_spk_), intent (in)       ::  alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (n < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 2; ierr(2) = n
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_
+    ierr(1) = 4; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 6; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psi_saxpby
+subroutine psi_saxpbyv(m,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m
+  real(psb_spk_), intent (in)       ::  x(:)
+  real(psb_spk_), intent (inout)    ::  y(:)
+  real(psb_spk_), intent (in)       :: alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 3; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 5; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psi_saxpbyv
+subroutine psi_sgthv(n,idx,alpha,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_spk_) :: x(:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i
+  if (beta == szero) then 
+    if (alpha == szero) then 
+      do i=1,n
+        y(i) = szero
+      end do
+    else if (alpha == sone) then 
+      do i=1,n
+        y(i) = x(idx(i))
+      end do
+    else if (alpha == -sone) then 
+      do i=1,n
+        y(i) = -x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = alpha*x(idx(i))
+      end do
+    end if
+  else 
+    if (beta == sone) then 
+      ! Do nothing
+    else if (beta == -sone) then 
+      y(1:n) = -y(1:n) 
+    else
+      y(1:n) = beta*y(1:n) 
+    end if
+
+    if (alpha == szero) then 
+      ! do nothing
+    else if (alpha == sone) then 
+      do i=1,n
+        y(i) = y(i) + x(idx(i))
+      end do
+    else if (alpha == -sone) then 
+      do i=1,n
+        y(i) = y(i) - x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = y(i) + alpha*x(idx(i))
+      end do
+    end if
+  end if
+
+end subroutine psi_sgthv
+subroutine psi_sgthzmv(n,k,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  real(psb_spk_) :: x(:,:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  pt=0
+  do j=1,k
+    do i=1,n
+      pt=pt+1
+      y(pt)=x(idx(i),j)
+    end do
+  end do
+
+end subroutine psi_sgthzmv
+subroutine psi_sgthzv(n,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_spk_) :: x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  do i=1,n
+    y(i)=x(idx(i))
+  end do
+
+end subroutine psi_sgthzv
+subroutine psi_ssctmv(n,k,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  real(psb_spk_) :: beta, x(:), y(:,:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == szero) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = x(pt)
+      end do
+    end do
+  else if (beta == sone) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = y(idx(i),j)+x(pt)
+      end do
+    end do
+  else
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
+      end do
+    end do
+  end if
+end subroutine psi_ssctmv
+subroutine psi_ssctv(n,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_spk_) :: beta, x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  if (beta == szero) then
+    do i=1,n
+      y(idx(i)) = x(i)
+    end do
+  else if (beta == sone) then
+    do i=1,n
+      y(idx(i)) = y(idx(i))+x(i)
+    end do
+  else
+    do i=1,n
+      y(idx(i)) = beta*y(idx(i))+x(i)
+    end do
+  end if
+end subroutine psi_ssctv
diff --git a/base/serial/psi_z_serial_impl.f90 b/base/serial/psi_z_serial_impl.f90
new file mode 100644
index 00000000..a3842126
--- /dev/null
+++ b/base/serial/psi_z_serial_impl.f90
@@ -0,0 +1,265 @@
+subroutine psi_zaxpby(m,n,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m, n
+  complex(psb_dpk_), intent (in)       ::  x(:,:)
+  complex(psb_dpk_), intent (inout)    ::  y(:,:)
+  complex(psb_dpk_), intent (in)       ::  alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (n < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 2; ierr(2) = n
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_
+    ierr(1) = 4; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 6; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end subroutine psi_zaxpby
+subroutine psi_zaxpbyv(m,alpha, x, beta, y, info)
+  
+  use psb_const_mod
+  use psb_error_mod
+  implicit none 
+  integer(psb_ipk_), intent(in)      :: m
+  complex(psb_dpk_), intent (in)       ::  x(:)
+  complex(psb_dpk_), intent (inout)    ::  y(:)
+  complex(psb_dpk_), intent (in)       :: alpha, beta
+  integer(psb_ipk_), intent(out)     :: info
+  integer(psb_ipk_) :: err_act
+  integer(psb_ipk_) :: lx, ly
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20)        :: name, ch_err
+
+  name='psb_geaxpby'
+  if(psb_get_errstatus() /= 0) return 
+  info=psb_success_
+  call psb_erractionsave(err_act)
+
+  if (m < 0) then
+    info = psb_err_iarg_neg_
+    ierr(1) = 1; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  lx = size(x,1)
+  ly = size(y,1)
+  if (lx < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 3; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+  if (ly < m) then 
+    info = psb_err_input_asize_small_i_ 
+    ierr(1) = 5; ierr(2) = m
+    call psb_errpush(info,name,i_err=ierr)
+    goto 9999 
+  end if
+
+  if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psi_zaxpbyv
+subroutine psi_zgthv(n,idx,alpha,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_dpk_) :: x(:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i
+  if (beta == zzero) then 
+    if (alpha == zzero) then 
+      do i=1,n
+        y(i) = zzero
+      end do
+    else if (alpha == zone) then 
+      do i=1,n
+        y(i) = x(idx(i))
+      end do
+    else if (alpha == -zone) then 
+      do i=1,n
+        y(i) = -x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = alpha*x(idx(i))
+      end do
+    end if
+  else 
+    if (beta == zone) then 
+      ! Do nothing
+    else if (beta == -zone) then 
+      y(1:n) = -y(1:n) 
+    else
+      y(1:n) = beta*y(1:n) 
+    end if
+
+    if (alpha == zzero) then 
+      ! do nothing
+    else if (alpha == zone) then 
+      do i=1,n
+        y(i) = y(i) + x(idx(i))
+      end do
+    else if (alpha == -zone) then 
+      do i=1,n
+        y(i) = y(i) - x(idx(i))
+      end do
+    else  
+      do i=1,n
+        y(i) = y(i) + alpha*x(idx(i))
+      end do
+    end if
+  end if
+
+end subroutine psi_zgthv
+subroutine psi_zgthzmv(n,k,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  complex(psb_dpk_) :: x(:,:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  pt=0
+  do j=1,k
+    do i=1,n
+      pt=pt+1
+      y(pt)=x(idx(i),j)
+    end do
+  end do
+
+end subroutine psi_zgthzmv
+subroutine psi_zgthzv(n,idx,x,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_dpk_) :: x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  do i=1,n
+    y(i)=x(idx(i))
+  end do
+
+end subroutine psi_zgthzv
+subroutine psi_zsctmv(n,k,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, k, idx(:)
+  complex(psb_dpk_) :: beta, x(:), y(:,:)
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == zzero) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = x(pt)
+      end do
+    end do
+  else if (beta == zone) then
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = y(idx(i),j)+x(pt)
+      end do
+    end do
+  else
+    pt=0
+    do j=1,k
+      do i=1,n
+        pt=pt+1
+        y(idx(i),j) = beta*y(idx(i),j)+x(pt)
+      end do
+    end do
+  end if
+end subroutine psi_zsctmv
+subroutine psi_zsctv(n,idx,x,beta,y)
+  
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_dpk_) :: beta, x(:), y(:)
+
+  ! Locals
+  integer(psb_ipk_) :: i
+
+  if (beta == zzero) then
+    do i=1,n
+      y(idx(i)) = x(i)
+    end do
+  else if (beta == zone) then
+    do i=1,n
+      y(idx(i)) = y(idx(i))+x(i)
+    end do
+  else
+    do i=1,n
+      y(idx(i)) = beta*y(idx(i))+x(i)
+    end do
+  end if
+end subroutine psi_zsctv