diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90
index 729a8745..8aa8a8a0 100644
--- a/base/modules/psb_serial_mod.f90
+++ b/base/modules/psb_serial_mod.f90
@@ -326,6 +326,131 @@ module psb_serial_mod
     end subroutine psb_zgeprt1
   end interface
   
+
+  
+  interface psb_spdot_srtd
+    function psb_s_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1,nv2
+      integer, intent(in) :: iv1(*), iv2(*)
+      real(psb_spk_), intent(in) :: v1(*),v2(*)
+      real(psb_spk_)      :: dot
+    end function psb_s_spdot_srtd
+
+    function psb_d_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1,nv2
+      integer, intent(in) :: iv1(*), iv2(*)
+      real(psb_dpk_), intent(in) :: v1(*),v2(*)
+      real(psb_dpk_)      :: dot
+    end function psb_d_spdot_srtd
+
+    function psb_c_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1,nv2
+      integer, intent(in) :: iv1(*), iv2(*)
+      complex(psb_spk_), intent(in) :: v1(*),v2(*)
+      complex(psb_spk_)      :: dot
+    end function psb_c_spdot_srtd
+
+    function psb_z_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1,nv2
+      integer, intent(in) :: iv1(*), iv2(*)
+      complex(psb_dpk_), intent(in) :: v1(*),v2(*)
+      complex(psb_dpk_)      :: dot
+    end function psb_z_spdot_srtd
+  end interface
+
+
+  interface psb_spge_dot
+    function psb_s_spge_dot(nv1,iv1,v1,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1
+      integer, intent(in) :: iv1(*)
+      real(psb_spk_), intent(in) :: v1(*),v2(*)
+      real(psb_spk_)      :: dot
+    end function psb_s_spge_dot
+
+    function psb_d_spge_dot(nv1,iv1,v1,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1
+      integer, intent(in) :: iv1(*)
+      real(psb_dpk_), intent(in) :: v1(*),v2(*)
+      real(psb_dpk_)      :: dot
+    end function psb_d_spge_dot
+
+    function psb_c_spge_dot(nv1,iv1,v1,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1
+      integer, intent(in) :: iv1(*)
+      complex(psb_spk_), intent(in) :: v1(*),v2(*)
+      complex(psb_spk_)      :: dot
+    end function psb_c_spge_dot
+
+    function psb_z_spge_dot(nv1,iv1,v1,v2) result(dot) 
+      use psb_const_mod
+      integer, intent(in) :: nv1
+      integer, intent(in) :: iv1(*)
+      complex(psb_dpk_), intent(in) :: v1(*),v2(*)
+      complex(psb_dpk_)      :: dot
+    end function psb_z_spge_dot
+  end interface
+
+
+  interface psb_nspaxpby
+    subroutine psb_d_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
+      use psb_const_mod
+      integer, intent(out)              :: nz
+      integer, intent(out)              :: iz(:)
+      real(psb_dpk_), intent (out)      :: z(:)
+      integer, intent(in)               :: nx, ny
+      integer, intent(in)               :: ix(:), iy(:)
+      real(psb_dpk_), intent (in)       :: x(:), y(:)
+      real(psb_dpk_), intent (in)       :: alpha, beta
+      integer, intent(out)              :: info
+    end subroutine psb_d_nspaxpby
+  end interface
+
+  interface psb_aspxpby
+    subroutine psb_s_aspxpby(alpha, nx, ix, x, beta, y, info)
+      use psb_const_mod
+      integer, intent(in)               :: nx
+      integer, intent(in)               :: ix(:)
+      real(psb_spk_), intent (in)       :: x(:)
+      real(psb_spk_), intent (inout)    :: y(:)
+      real(psb_spk_), intent (in)       :: alpha, beta
+      integer, intent(out)              :: info
+    end subroutine psb_s_aspxpby
+    subroutine psb_d_aspxpby(alpha, nx, ix, x, beta, y, info)
+      use psb_const_mod
+      integer, intent(in)               :: nx
+      integer, intent(in)               :: ix(:)
+      real(psb_dpk_), intent (in)       :: x(:)
+      real(psb_dpk_), intent (inout)    :: y(:)
+      real(psb_dpk_), intent (in)       :: alpha, beta
+      integer, intent(out)              :: info
+    end subroutine psb_d_aspxpby
+    subroutine psb_c_aspxpby(alpha, nx, ix, x, beta, y, info)
+      use psb_const_mod
+      integer, intent(in)               :: nx
+      integer, intent(in)               :: ix(:)
+      complex(psb_spk_), intent (in)    :: x(:)
+      complex(psb_spk_), intent (inout) :: y(:)
+      complex(psb_spk_), intent (in)    :: alpha, beta
+      integer, intent(out)              :: info
+    end subroutine psb_c_aspxpby
+    subroutine psb_z_aspxpby(alpha, nx, ix, x, beta, y, info)
+      use psb_const_mod
+      integer, intent(in)               :: nx
+      integer, intent(in)               :: ix(:)
+      complex(psb_dpk_), intent (in)    :: x(:)
+      complex(psb_dpk_), intent (inout) :: y(:)
+      complex(psb_dpk_), intent (in)    :: alpha, beta
+      integer, intent(out)              :: info
+    end subroutine psb_z_aspxpby
+  end interface
+
   
 end module psb_serial_mod
 
diff --git a/base/serial/Makefile b/base/serial/Makefile
index 1e226e52..684d7baa 100644
--- a/base/serial/Makefile
+++ b/base/serial/Makefile
@@ -5,7 +5,8 @@ FOBJS = psb_lsame.o psi_serial_impl.o psb_sort_impl.o \
 	psb_srwextd.o psb_drwextd.o psb_crwextd.o psb_zrwextd.o \
 	psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \
 	psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \
-	psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o
+	psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o\
+	psb_spdot_srtd.o psb_aspxpby.o psb_spge_dot.o 
 
 LIBDIR=..
 MODDIR=../modules
diff --git a/base/serial/psb_aspxpby.f90 b/base/serial/psb_aspxpby.f90
new file mode 100644
index 00000000..edfa77a4
--- /dev/null
+++ b/base/serial/psb_aspxpby.f90
@@ -0,0 +1,193 @@
+subroutine psb_s_aspxpby(alpha, nx, ix, x, beta, y, info)
+  use psb_const_mod
+  integer, intent(in)               :: nx
+  integer, intent(in)               :: ix(:)
+  real(psb_spk_), intent (in)       :: x(:)
+  real(psb_spk_), intent (inout)    :: y(:)
+  real(psb_spk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+  integer :: i, ip 
+
+  info=psb_success_
+
+  if (nx > max(size(ix),size(x))) then 
+    info = -2
+    return
+  end if
+
+  if (beta /= sone) then 
+    if (beta == -sone) then 
+      y(:) = -y(:)
+    else if (beta == szero) then 
+      y(:) = szero
+    else
+      y(:) = beta * y(:)
+    end if
+  end if
+
+  if (alpha == szero) return 
+
+  if (alpha == sone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + x(i)
+    end do
+  else   if (alpha == -sone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) - x(i)
+    end do
+  else
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + alpha*x(i)
+    end do
+  end if
+
+
+end subroutine psb_s_aspxpby
+
+subroutine psb_d_aspxpby(alpha, nx, ix, x, beta, y, info)
+  use psb_const_mod
+  integer, intent(in)               :: nx
+  integer, intent(in)               :: ix(:)
+  real(psb_dpk_), intent (in)       :: x(:)
+  real(psb_dpk_), intent (inout)    :: y(:)
+  real(psb_dpk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+  integer :: i, ip 
+
+  info=psb_success_
+  
+  if (nx > max(size(ix),size(x))) then 
+    info = -2
+    return
+  end if
+
+  if (beta /= done) then 
+    if (beta == -done) then 
+      y(:) = -y(:)
+    else if (beta == dzero) then 
+      y(:) = dzero
+    else
+      y(:) = beta * y(:)
+    end if
+  end if
+
+  if (alpha == dzero) return 
+
+  if (alpha == done) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + x(i)
+    end do
+  else   if (alpha == -done) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) - x(i)
+    end do
+  else
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + alpha*x(i)
+    end do
+  end if
+
+
+end subroutine psb_d_aspxpby
+subroutine psb_c_aspxpby(alpha, nx, ix, x, beta, y, info)
+  use psb_const_mod
+  integer, intent(in)               :: nx
+  integer, intent(in)               :: ix(:)
+  complex(psb_spk_), intent (in)    :: x(:)
+  complex(psb_spk_), intent (inout) :: y(:)
+  complex(psb_spk_), intent (in)    :: alpha, beta
+  integer, intent(out)              :: info
+  integer :: i, ip 
+
+  info=psb_success_
+  
+  if (nx > max(size(ix),size(x))) then 
+    info = -2
+    return
+  end if
+  
+  if (beta /= cone) then 
+    if (beta == -cone) then 
+      y(:) = -y(:)
+    else if (beta == czero) then 
+      y(:) = czero
+    else
+      y(:) = beta * y(:)
+    end if
+  end if
+
+  if (alpha == czero) return 
+
+  if (alpha == cone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + x(i)
+    end do
+  else   if (alpha == -cone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) - x(i)
+    end do
+  else
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + alpha*x(i)
+    end do
+  end if
+
+
+end subroutine psb_c_aspxpby
+
+subroutine psb_z_aspxpby(alpha, nx, ix, x, beta, y, info)
+  use psb_const_mod
+  integer, intent(in)               :: nx
+  integer, intent(in)               :: ix(:)
+  complex(psb_dpk_), intent (in)    :: x(:)
+  complex(psb_dpk_), intent (inout) :: y(:)
+  complex(psb_dpk_), intent (in)    :: alpha, beta
+  integer, intent(out)              :: info
+  integer :: i, ip 
+
+  info=psb_success_
+  
+  if (nx > max(size(ix),size(x))) then 
+    info = -2
+    return
+  end if
+  
+  if (beta /= zone) then 
+    if (beta == -zone) then 
+      y(:) = -y(:)
+    else if (beta == zzero) then 
+      y(:) = zzero
+    else
+      y(:) = beta * y(:)
+    end if
+  end if
+
+  if (alpha == zzero) return 
+
+  if (alpha == zone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + x(i)
+    end do
+  else if (alpha == -zone) then 
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) - x(i)
+    end do
+  else
+    do i=1, nx
+      ip = ix(i) 
+      y(ip) = y(ip) + alpha*x(i)
+    end do
+  end if
+
+end subroutine psb_z_aspxpby
diff --git a/base/serial/psb_spdot_srtd.f90 b/base/serial/psb_spdot_srtd.f90
new file mode 100644
index 00000000..e45c7666
--- /dev/null
+++ b/base/serial/psb_spdot_srtd.f90
@@ -0,0 +1,521 @@
+subroutine psb_s_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
+  use psb_const_mod
+  integer, intent(out)              :: nz
+  integer, intent(out)              :: iz(:)
+  real(psb_spk_), intent (out)      :: z(:)
+  integer, intent(in)               :: nx, ny
+  integer, intent(in)               :: ix(:), iy(:)
+  real(psb_spk_), intent (in)       :: x(:), y(:)
+  real(psb_spk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+
+  integer        :: i,j,k, ipx, ipy, isz
+  real(psb_spk_) :: acc
+
+  info=psb_success_
+  nz   = 0
+  ipx  = 1
+  ipy  = 1 
+  isz  = min(size(iz),size(z))
+
+  if (beta == szero) then 
+    if (alpha == szero) return 
+    nz = nx
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = ix(1:nx)
+    z(1:nz)  = alpha*x(1:nx) 
+
+  else if (alpha == szero) then
+
+    nz = ny
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = iy(1:ny)
+    z(1:nz)  = beta*y(1:ny) 
+
+  else
+    
+    do 
+      if (ipx > nx) exit
+      if (ipy > ny) exit
+      if (ix(ipx) == iy(ipy)) then 
+        acc = beta*y(ipy) + alpha*x(ipx) 
+        if (acc /= szero) then 
+          nz = nz + 1 
+          if (nz > isz) then 
+            info = -1 
+            return
+          endif
+          iz(nz) = ix(ipx) 
+          z(nz)  = acc
+          ipx    = ipx + 1
+          ipy    = ipy + 1
+        end if
+      else 
+        nz = nz + 1 
+        if (nz > isz) then 
+          info = -1 
+          return
+        endif
+        if (ix(ipx) < iy(ipy)) then 
+          iz(nz) = ix(ipx) 
+          z(nz)  = alpha*x(ipx) 
+          ipx    = ipx + 1
+        else
+          iz(nz) = iy(ipy) 
+          z(nz)  = beta*y(ipy) 
+          ipy    = ipy + 1
+        end if
+      end if
+    end do
+    do 
+      if (ipx > nx) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = ix(ipx) 
+      z(nz)  = alpha*x(ipx) 
+      ipx    = ipx + 1
+    end do
+    do 
+      if (ipy > ny) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = iy(ipy) 
+      z(nz)  = beta*y(ipy) 
+      ipy    = ipy + 1
+    end do
+  end if
+
+end subroutine psb_s_nspaxpby
+
+function psb_s_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1,nv2
+  integer, intent(in) :: iv1(*), iv2(*)
+  real(psb_spk_), intent(in) :: v1(*),v2(*)
+  real(psb_spk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = szero 
+  ip1 = 1
+  ip2 = 1
+
+  do 
+    if (ip1 > nv1) exit
+    if (ip2 > nv2) exit
+    if (iv1(ip1) == iv2(ip2)) then 
+      dot = dot + v1(ip1)*v2(ip2)
+      ip1 = ip1 + 1
+      ip2 = ip2 + 1
+    else if (iv1(ip1) < iv2(ip2)) then 
+      ip1 = ip1 + 1 
+    else
+      ip2 = ip2 + 1 
+    end if
+  end do
+
+end function psb_s_spdot_srtd
+
+subroutine psb_d_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
+  use psb_const_mod
+  integer, intent(out)              :: nz
+  integer, intent(out)              :: iz(:)
+  real(psb_dpk_), intent (out)      :: z(:)
+  integer, intent(in)               :: nx, ny
+  integer, intent(in)               :: ix(:), iy(:)
+  real(psb_dpk_), intent (in)       :: x(:), y(:)
+  real(psb_dpk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+
+  integer        :: i,j,k, ipx, ipy, isz
+  real(psb_dpk_) :: acc
+
+  info=psb_success_
+  nz   = 0
+  ipx  = 1
+  ipy  = 1 
+  isz  = min(size(iz),size(z))
+
+  if (beta == dzero) then 
+    if (alpha == dzero) return 
+    nz = nx
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = ix(1:nx)
+    z(1:nz)  = alpha*x(1:nx) 
+
+  else if (alpha == dzero) then
+
+    nz = ny
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = iy(1:ny)
+    z(1:nz)  = beta*y(1:ny) 
+
+  else
+    
+    do 
+      if (ipx > nx) exit
+      if (ipy > ny) exit
+      if (ix(ipx) == iy(ipy)) then 
+        acc = beta*y(ipy) + alpha*x(ipx) 
+        if (acc /= dzero) then 
+          nz = nz + 1 
+          if (nz > isz) then 
+            info = -1 
+            return
+          endif
+          iz(nz) = ix(ipx) 
+          z(nz)  = acc
+          ipx    = ipx + 1
+          ipy    = ipy + 1
+        end if
+      else 
+        nz = nz + 1 
+        if (nz > isz) then 
+          info = -1 
+          return
+        endif
+        if (ix(ipx) < iy(ipy)) then 
+          iz(nz) = ix(ipx) 
+          z(nz)  = alpha*x(ipx) 
+          ipx    = ipx + 1
+        else
+          iz(nz) = iy(ipy) 
+          z(nz)  = beta*y(ipy) 
+          ipy    = ipy + 1
+        end if
+      end if
+    end do
+    do 
+      if (ipx > nx) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = ix(ipx) 
+      z(nz)  = alpha*x(ipx) 
+      ipx    = ipx + 1
+    end do
+    do 
+      if (ipy > ny) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = iy(ipy) 
+      z(nz)  = beta*y(ipy) 
+      ipy    = ipy + 1
+    end do
+  end if
+
+end subroutine psb_d_nspaxpby
+
+
+
+function psb_d_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1,nv2
+  integer, intent(in) :: iv1(*), iv2(*)
+  real(psb_dpk_), intent(in) :: v1(*),v2(*)
+  real(psb_dpk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = dzero 
+  ip1 = 1
+  ip2 = 1
+  if (nv1 == 0) return
+  if (nv2 == 0) return
+  do 
+    if (iv1(ip1) == iv2(ip2)) then 
+      dot = dot + v1(ip1)*v2(ip2)
+      ip1 = ip1 + 1
+      ip2 = ip2 + 1
+    else if (iv1(ip1) < iv2(ip2)) then 
+      ip1 = ip1 + 1 
+    else
+      ip2 = ip2 + 1 
+    end if
+    if ((ip1 > nv1) .or. (ip2 > nv2)) exit
+  end do
+
+end function psb_d_spdot_srtd
+
+subroutine psb_c_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
+  use psb_const_mod
+  integer, intent(out)              :: nz
+  integer, intent(out)              :: iz(:)
+  complex(psb_spk_), intent (out)      :: z(:)
+  integer, intent(in)               :: nx, ny
+  integer, intent(in)               :: ix(:), iy(:)
+  complex(psb_spk_), intent (in)       :: x(:), y(:)
+  complex(psb_spk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+
+  integer        :: i,j,k, ipx, ipy, isz
+  complex(psb_spk_) :: acc
+
+  info=psb_success_
+  nz   = 0
+  ipx  = 1
+  ipy  = 1 
+  isz  = min(size(iz),size(z))
+
+  if (beta == czero) then 
+    if (alpha == czero) return 
+    nz = nx
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = ix(1:nx)
+    z(1:nz)  = alpha*x(1:nx) 
+
+  else if (alpha == czero) then
+
+    nz = ny
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = iy(1:ny)
+    z(1:nz)  = beta*y(1:ny) 
+
+  else
+    
+    do 
+      if (ipx > nx) exit
+      if (ipy > ny) exit
+      if (ix(ipx) == iy(ipy)) then 
+        acc = beta*y(ipy) + alpha*x(ipx) 
+        if (acc /= czero) then 
+          nz = nz + 1 
+          if (nz > isz) then 
+            info = -1 
+            return
+          endif
+          iz(nz) = ix(ipx) 
+          z(nz)  = acc
+          ipx    = ipx + 1
+          ipy    = ipy + 1
+        end if
+      else 
+        nz = nz + 1 
+        if (nz > isz) then 
+          info = -1 
+          return
+        endif
+        if (ix(ipx) < iy(ipy)) then 
+          iz(nz) = ix(ipx) 
+          z(nz)  = alpha*x(ipx) 
+          ipx    = ipx + 1
+        else
+          iz(nz) = iy(ipy) 
+          z(nz)  = beta*y(ipy) 
+          ipy    = ipy + 1
+        end if
+      end if
+    end do
+    do 
+      if (ipx > nx) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = ix(ipx) 
+      z(nz)  = alpha*x(ipx) 
+      ipx    = ipx + 1
+    end do
+    do 
+      if (ipy > ny) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = iy(ipy) 
+      z(nz)  = beta*y(ipy) 
+      ipy    = ipy + 1
+    end do
+  end if
+
+end subroutine psb_c_nspaxpby
+
+function psb_c_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1,nv2
+  integer, intent(in) :: iv1(*), iv2(*)
+  complex(psb_spk_), intent(in) :: v1(*),v2(*)
+  complex(psb_spk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = czero 
+  ip1 = 1
+  ip2 = 1
+
+  do 
+    if (ip1 > nv1) exit
+    if (ip2 > nv2) exit
+    if (iv1(ip1) == iv2(ip2)) then 
+      dot = dot + conjg(v1(ip1))*v2(ip2)
+      ip1 = ip1 + 1
+      ip2 = ip2 + 1
+    else if (iv1(ip1) < iv2(ip2)) then 
+      ip1 = ip1 + 1 
+    else
+      ip2 = ip2 + 1 
+    end if
+  end do
+
+end function psb_c_spdot_srtd
+
+subroutine psb_z_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
+  use psb_const_mod
+  integer, intent(out)              :: nz
+  integer, intent(out)              :: iz(:)
+  complex(psb_dpk_), intent (out)      :: z(:)
+  integer, intent(in)               :: nx, ny
+  integer, intent(in)               :: ix(:), iy(:)
+  complex(psb_dpk_), intent (in)       :: x(:), y(:)
+  complex(psb_dpk_), intent (in)       :: alpha, beta
+  integer, intent(out)              :: info
+
+  integer        :: i,j,k, ipx, ipy, isz
+  complex(psb_dpk_) :: acc
+
+  info=psb_success_
+  nz   = 0
+  ipx  = 1
+  ipy  = 1 
+  isz  = min(size(iz),size(z))
+
+  if (beta == zzero) then 
+    if (alpha == zzero) return 
+    nz = nx
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = ix(1:nx)
+    z(1:nz)  = alpha*x(1:nx) 
+
+  else if (alpha == zzero) then
+
+    nz = ny
+    if (nz > isz) then 
+      info = -1 
+      return
+    endif
+    iz(1:nz) = iy(1:ny)
+    z(1:nz)  = beta*y(1:ny) 
+
+  else
+    
+    do 
+      if (ipx > nx) exit
+      if (ipy > ny) exit
+      if (ix(ipx) == iy(ipy)) then 
+        acc = beta*y(ipy) + alpha*x(ipx) 
+        if (acc /= zzero) then 
+          nz = nz + 1 
+          if (nz > isz) then 
+            info = -1 
+            return
+          endif
+          iz(nz) = ix(ipx) 
+          z(nz)  = acc
+          ipx    = ipx + 1
+          ipy    = ipy + 1
+        end if
+      else 
+        nz = nz + 1 
+        if (nz > isz) then 
+          info = -1 
+          return
+        endif
+        if (ix(ipx) < iy(ipy)) then 
+          iz(nz) = ix(ipx) 
+          z(nz)  = alpha*x(ipx) 
+          ipx    = ipx + 1
+        else
+          iz(nz) = iy(ipy) 
+          z(nz)  = beta*y(ipy) 
+          ipy    = ipy + 1
+        end if
+      end if
+    end do
+    do 
+      if (ipx > nx) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = ix(ipx) 
+      z(nz)  = alpha*x(ipx) 
+      ipx    = ipx + 1
+    end do
+    do 
+      if (ipy > ny) exit
+      nz = nz + 1 
+      if (nz > isz) then 
+        info = -1 
+        return
+      endif
+      iz(nz) = iy(ipy) 
+      z(nz)  = beta*y(ipy) 
+      ipy    = ipy + 1
+    end do
+  end if
+
+end subroutine psb_z_nspaxpby
+
+function psb_z_spdot_srtd(nv1,iv1,v1,nv2,iv2,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1,nv2
+  integer, intent(in) :: iv1(*), iv2(*)
+  complex(psb_dpk_), intent(in) :: v1(*),v2(*)
+  complex(psb_dpk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = zzero 
+  ip1 = 1
+  ip2 = 1
+
+  do 
+    if (ip1 > nv1) exit
+    if (ip2 > nv2) exit
+    if (iv1(ip1) == iv2(ip2)) then 
+      dot = dot + conjg(v1(ip1))*v2(ip2)
+      ip1 = ip1 + 1
+      ip2 = ip2 + 1
+    else if (iv1(ip1) < iv2(ip2)) then 
+      ip1 = ip1 + 1 
+    else
+      ip2 = ip2 + 1 
+    end if
+  end do
+
+end function psb_z_spdot_srtd
diff --git a/base/serial/psb_spge_dot.f90 b/base/serial/psb_spge_dot.f90
new file mode 100644
index 00000000..dbd88955
--- /dev/null
+++ b/base/serial/psb_spge_dot.f90
@@ -0,0 +1,68 @@
+function psb_s_spge_dot(nv1,iv1,v1,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1
+  integer, intent(in) :: iv1(*)
+  real(psb_spk_), intent(in) :: v1(*),v2(*)
+  real(psb_spk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = szero 
+  ip1 = 1
+
+  do i=1, nv1
+    dot = dot + v1(i)*v2(iv1(i))
+  end do
+
+end function psb_s_spge_dot
+
+function psb_d_spge_dot(nv1,iv1,v1,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1
+  integer, intent(in) :: iv1(*)
+  real(psb_dpk_), intent(in) :: v1(*),v2(*)
+  real(psb_dpk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = dzero 
+
+  do i=1, nv1
+    dot = dot + v1(i)*v2(iv1(i))
+  end do
+
+end function psb_d_spge_dot
+
+function psb_c_spge_dot(nv1,iv1,v1,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1
+  integer, intent(in) :: iv1(*)
+  complex(psb_spk_), intent(in) :: v1(*),v2(*)
+  complex(psb_spk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = czero 
+
+  do i=1, nv1
+    dot = dot + v1(i)*v2(iv1(i))
+  end do
+
+end function psb_c_spge_dot
+
+function psb_z_spge_dot(nv1,iv1,v1,v2) result(dot) 
+  use psb_const_mod
+  integer, intent(in) :: nv1
+  integer, intent(in) :: iv1(*)
+  complex(psb_dpk_), intent(in) :: v1(*),v2(*)
+  complex(psb_dpk_)      :: dot
+
+  integer :: i,j,k, ip1, ip2
+
+  dot = zzero 
+
+  do i=1, nv1
+    dot = dot + v1(i)*v2(iv1(i))
+  end do
+
+end function psb_z_spge_dot