diff --git a/base/comm/Makefile b/base/comm/Makefile
index c6c97f25..ee99c2d7 100644
--- a/base/comm/Makefile
+++ b/base/comm/Makefile
@@ -6,7 +6,9 @@ OBJS =  psb_dgather.o psb_dhalo.o psb_dovrl.o \
 	psb_cgather.o psb_chalo.o psb_covrl.o \
 	psb_zgather.o psb_zhalo.o psb_zovrl.o 
 
-MPFOBJS=psb_dscatter.o psb_zscatter.o psb_iscatter.o psb_cscatter.o psb_sscatter.o
+
+MPFOBJS=psb_dscatter.o psb_zscatter.o psb_iscatter.o psb_cscatter.o psb_sscatter.o\
+	psb_dspgather.o
 LIBDIR=..
 MODDIR=../modules
 FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG)$(MODDIR) $(FMFLAG).
diff --git a/base/comm/psb_dspgather.F90 b/base/comm/psb_dspgather.F90
new file mode 100644
index 00000000..aed45035
--- /dev/null
+++ b/base/comm/psb_dspgather.F90
@@ -0,0 +1,113 @@
+subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_descriptor_type
+  use psb_error_mod
+  use psb_mat_mod
+
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_d_sparse_mat), intent(inout) :: loca
+  type(psb_d_sparse_mat), intent(out)   :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer, intent(out)            :: info
+  integer, intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_d_coo_sparse_mat)      :: loc_coo, glob_coo
+  integer :: ictxt,np,me, err_act, icomm, dupl_, nrg, ncg, nzg
+  integer :: ip, ndx,naggrm1,naggrp1, i, j, k
+  logical :: keepnum_, keeploc_
+  integer, allocatable :: nzbr(:), idisp(:)
+  character(len=20) :: name
+  integer            :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=0
+
+  call psb_erractionsave(err_act)
+  ictxt = psb_cd_get_context(desc_a)
+  icomm = psb_cd_get_mpic(desc_a)
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+
+  if (keepnum_) then 
+    nrg = psb_cd_get_global_rows(desc_a)
+    ncg = psb_cd_get_global_rows(desc_a)
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= 0) then 
+      info=4025
+      call psb_errpush(info,name,i_err=(/2*np,0,0,0,0/),&
+           & a_err='integer')
+      goto 9999      
+    end if
+    call loca%mv_to(loc_coo)
+    nzbr(:) = 0
+    nzbr(me+1) = loc_coo%get_nzeros()
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == 0) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= 0) goto 9999
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,mpi_double_precision,&
+         & glob_coo%val,nzbr,idisp,&
+         & mpi_double_precision,icomm,info)
+    if (info == 0) call mpi_allgatherv(loc_coo%ia,ndx,mpi_integer,&
+         & glob_coo%ia,nzbr,idisp,&
+         & mpi_integer,icomm,info)
+    if (info == 0) call mpi_allgatherv(loc_coo%ja,ndx,mpi_integer,&
+         & glob_coo%ja,nzbr,idisp,&
+         & mpi_integer,icomm,info)
+    
+    if (info /= 0) then 
+      call psb_errpush(4001,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if
+    
+    if (keeploc_) then
+      call loca%mv_from(loc_coo)
+    else
+      call loc_coo%free()
+    end if
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(0,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_erractionrestore(err_act)
+  if (err_act.eq.psb_act_abort_) then
+    call psb_error()
+    return
+  end if
+  return
+
+ 
+end subroutine psb_dsp_allgather
diff --git a/base/modules/Makefile b/base/modules/Makefile
index 9bc82c91..89a75b5b 100644
--- a/base/modules/Makefile
+++ b/base/modules/Makefile
@@ -53,6 +53,7 @@ psi_mod.o: psb_penv_mod.o psb_error_mod.o psb_desc_type.o psb_const_mod.o psi_se
 psb_desc_type.o: psb_const_mod.o psb_error_mod.o psb_penv_mod.o psb_realloc_mod.o  psb_hash_mod.o
 psb_linmap_mod.o: psb_linmap_type_mod.o psb_mat_mod.o 
 psb_linmap_type_mod.o: psb_desc_type.o psb_error_mod.o psb_serial_mod.o psb_comm_mod.o psb_mat_mod.o
+psb_comm_mod.o: psb_desc_type.o psb_mat_mod.o
 psb_check_mod.o: psb_desc_type.o
 psb_serial_mod.o: psb_mat_mod.o psb_string_mod.o psb_sort_mod.o psi_serial_mod.o
 psb_sort_mod.o: psb_error_mod.o psb_realloc_mod.o psb_const_mod.o
diff --git a/base/modules/psb_comm_mod.f90 b/base/modules/psb_comm_mod.f90
index 6f79a774..f8738b82 100644
--- a/base/modules/psb_comm_mod.f90
+++ b/base/modules/psb_comm_mod.f90
@@ -302,6 +302,17 @@ module psb_comm_mod
   end interface
 
   interface psb_gather
+    subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+      use psb_descriptor_type
+      use psb_mat_mod
+      implicit none
+      type(psb_d_sparse_mat), intent(inout) :: loca
+      type(psb_d_sparse_mat), intent(out)   :: globa
+      type(psb_desc_type), intent(in) :: desc_a
+      integer, intent(out)            :: info
+      integer, intent(in), optional   :: root,dupl
+      logical, intent(in), optional   :: keepnum,keeploc
+    end subroutine psb_dsp_allgather
     subroutine  psb_igatherm(globx, locx, desc_a, info, root)
       use psb_descriptor_type
       integer, intent(in)             :: locx(:,:)
diff --git a/base/modules/psb_d_mat_mod.f03 b/base/modules/psb_d_mat_mod.f03
index ac75facd..c50b94ef 100644
--- a/base/modules/psb_d_mat_mod.f03
+++ b/base/modules/psb_d_mat_mod.f03
@@ -52,18 +52,25 @@ module psb_d_mat_mod
     procedure, pass(a) :: d_csgetrow
     procedure, pass(a) :: d_csgetblk
     generic, public    :: csget => d_csgetptn, d_csgetrow, d_csgetblk 
-    procedure, pass(a) :: csclip
+    procedure, pass(a) :: d_csclip
+    procedure, pass(a) :: d_b_csclip
+    generic, public    :: csclip => d_b_csclip, d_csclip
     procedure, pass(a) :: reall => reallocate_nz
     procedure, pass(a) :: get_neigh
     procedure, pass(a) :: d_cscnv
     procedure, pass(a) :: d_cscnv_ip
-    generic, public    :: cscnv => d_cscnv, d_cscnv_ip
+    procedure, pass(a) :: d_cscnv_base
+    generic, public    :: cscnv => d_cscnv, d_cscnv_ip, d_cscnv_base
     procedure, pass(a) :: reinit
     procedure, pass(a) :: print => sparse_print
     procedure, pass(a) :: d_mv_from
     generic, public    :: mv_from => d_mv_from
+    procedure, pass(a) :: d_mv_to
+    generic, public    :: mv_to => d_mv_to
     procedure, pass(a) :: d_cp_from
     generic, public    :: cp_from => d_cp_from
+    procedure, pass(a) :: d_cp_to
+    generic, public    :: cp_to => d_cp_to
     procedure, pass(a) :: d_transp_1mat
     procedure, pass(a) :: d_transp_2mat
     generic, public    :: transp => d_transp_1mat, d_transp_2mat
@@ -92,7 +99,7 @@ module psb_d_mat_mod
        & get_state, get_dupl, is_null, is_bld, is_upd, &
        & is_asb, is_sorted, is_upper, is_lower, is_triangle, &
        & is_unit, get_neigh, csall, csput, d_csgetrow,&
-       & d_csgetblk, csclip, d_cscnv, d_cscnv_ip, &
+       & d_csgetblk, d_csclip, d_b_csclip, d_cscnv, d_cscnv_ip, &
        & reallocate_nz, free, trim, &
        & sparse_print, reinit, &
        & set_nrows, set_ncols, set_dupl, &
@@ -100,7 +107,7 @@ module psb_d_mat_mod
        & set_upd, set_asb, set_sorted, &
        & set_upper, set_lower, set_triangle, &
        & set_unit, get_diag, get_nz_row, d_csgetptn, &
-       & d_mv_from, d_cp_from, &
+       & d_mv_from, d_mv_to, d_cp_from, d_cp_to,&
        & d_transp_1mat, d_transp_2mat, &
        & d_transc_1mat, d_transc_2mat
 
@@ -1262,7 +1269,7 @@ contains
 
 
 
-  subroutine csclip(a,b,info,&
+  subroutine d_csclip(a,b,info,&
        & imin,imax,jmin,jmax,rscale,cscale)
     ! Output is always in  COO format 
     use psb_error_mod
@@ -1306,7 +1313,52 @@ contains
       return
     end if
 
-  end subroutine csclip
+  end subroutine d_csclip
+
+  subroutine d_b_csclip(a,b,info,&
+       & imin,imax,jmin,jmax,rscale,cscale)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    use psb_d_base_mat_mod
+    implicit none
+    
+    class(psb_d_sparse_mat), intent(in) :: a
+    type(psb_d_coo_sparse_mat), intent(out) :: b
+    integer,intent(out)                  :: info
+    integer, intent(in), optional        :: imin,imax,jmin,jmax
+    logical, intent(in), optional        :: rscale,cscale
+
+    Integer :: err_act
+    character(len=20)  :: name='csclip'
+    logical, parameter :: debug=.false.
+    type(psb_d_coo_sparse_mat), allocatable  :: acoo
+
+    info = 0
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    allocate(acoo,stat=info)    
+    if (info == 0) call a%a%csclip(b,info,&
+       & imin,imax,jmin,jmax,rscale,cscale)
+    if (info /= 0) goto 9999 
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine d_b_csclip
 
 
 
@@ -1494,6 +1546,61 @@ contains
 
   end subroutine d_cscnv_ip
 
+
+  subroutine d_cscnv_base(a,b,info,dupl)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_d_sparse_mat), intent(in)       :: a
+    class(psb_d_base_sparse_mat), intent(out) :: b
+    integer, intent(out)                   :: info
+    integer,optional, intent(in)           :: dupl
+
+
+    type(psb_d_coo_sparse_mat)  :: altmp
+    Integer :: err_act
+    character(len=20)  :: name='cscnv'
+    logical, parameter :: debug=.false.
+
+    info = 0
+    call psb_erractionsave(err_act)
+
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    call a%a%cp_to_coo(altmp,info )
+    if ((info == 0).and.present(dupl)) then 
+      call altmp%set_dupl(dupl)
+    end if
+    call altmp%fix(info)
+    if (info == 0) call altmp%trim()
+    if (info == 0) call altmp%set_asb() 
+    if (info == 0) call b%mv_from_coo(altmp,info)
+
+    if (info /= 0) then
+      info = 4010
+      call psb_errpush(info,name,a_err="mv_from")
+      goto 9999
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine d_cscnv_base
+
+
+
   subroutine d_mv_from(a,b)
     use psb_error_mod
     use psb_string_mod
@@ -1538,6 +1645,33 @@ contains
     end if
   end subroutine d_cp_from
 
+  subroutine d_mv_to(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_d_sparse_mat), intent(inout) :: a
+    class(psb_d_base_sparse_mat), intent(out) :: b
+    integer :: info
+
+    call b%mv_from_fmt(a%a,info)
+    
+    return
+  end subroutine d_mv_to
+
+  subroutine d_cp_to(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_d_sparse_mat), intent(in) :: a
+    class(psb_d_base_sparse_mat), intent(out) :: b
+    integer :: info
+
+    call b%cp_from_fmt(a%a,info)
+    
+    return
+  end subroutine d_cp_to
+
+
   subroutine d_sparse_mat_move(a,b,info)
     use psb_error_mod
     use psb_string_mod
@@ -1979,7 +2113,7 @@ contains
     integer, intent(out)                 :: info
 
     Integer :: err_act
-    character(len=20)  :: name='csnmi'
+    character(len=20)  :: name='get_diag'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -2015,7 +2149,7 @@ contains
     integer, intent(out)                    :: info
 
     Integer :: err_act
-    character(len=20)  :: name='csnmi'
+    character(len=20)  :: name='scal'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
@@ -2051,7 +2185,7 @@ contains
     integer, intent(out)                    :: info
 
     Integer :: err_act
-    character(len=20)  :: name='csnmi'
+    character(len=20)  :: name='scal'
     logical, parameter :: debug=.false.
 
     call psb_erractionsave(err_act)
diff --git a/base/modules/psb_serial_mod.f90 b/base/modules/psb_serial_mod.f90
index fef55e9e..88acca98 100644
--- a/base/modules/psb_serial_mod.f90
+++ b/base/modules/psb_serial_mod.f90
@@ -41,5 +41,59 @@ module psb_serial_mod
        & psb_sct => psi_sct
 
   use psb_mat_mod
+
+  interface psb_symbmm
+    subroutine psb_dsymbmm(a,b,c,info)
+      use psb_mat_mod
+      implicit none 
+      type(psb_d_sparse_mat), intent(in)  :: a,b
+      type(psb_d_sparse_mat), intent(out) :: c
+      integer, intent(out)                :: info
+    end subroutine psb_dsymbmm
+    subroutine psb_dbase_symbmm(a,b,c,info)
+      use psb_mat_mod
+      implicit none 
+      class(psb_d_base_sparse_mat), intent(in) :: a,b
+      type(psb_d_csr_sparse_mat), intent(out)  :: c
+      integer, intent(out)                     :: info
+    end subroutine psb_dbase_symbmm
+  end interface
+  interface psb_numbmm
+    subroutine psb_dnumbmm(a,b,c)
+      use psb_mat_mod
+      implicit none 
+      type(psb_d_sparse_mat), intent(in) :: a,b
+      type(psb_d_sparse_mat), intent(inout)  :: c
+    end subroutine psb_dnumbmm
+    subroutine psb_dbase_numbmm(a,b,c)
+      use psb_mat_mod
+      implicit none 
+      class(psb_d_base_sparse_mat), intent(in) :: a,b
+      type(psb_d_csr_sparse_mat), intent(inout)  :: c
+    end subroutine psb_dbase_numbmm
+  end interface
+
+  interface psb_rwextd
+    subroutine psb_drwextd(nr,a,info,b,rowscale)
+      use psb_mat_mod
+      implicit none
+      integer, intent(in)                          :: nr
+      type(psb_d_sparse_mat), intent(inout)        :: a
+      integer,intent(out)                          :: info
+      type(psb_d_sparse_mat), intent(in), optional :: b
+      logical,intent(in), optional                 :: rowscale
+    end subroutine psb_drwextd
+    subroutine psb_dbase_rwextd(nr,a,info,b,rowscale)
+      use psb_mat_mod
+      implicit none
+      integer, intent(in)                                :: nr
+      class(psb_d_base_sparse_mat), intent(inout)        :: a
+      integer,intent(out)                                :: info
+      class(psb_d_base_sparse_mat), intent(in), optional :: b
+      logical,intent(in), optional                       :: rowscale
+    end subroutine psb_dbase_rwextd
+  end interface
+  
+  
 end module psb_serial_mod
 
diff --git a/base/serial/Makefile b/base/serial/Makefile
index cbe5f71a..872287c9 100644
--- a/base/serial/Makefile
+++ b/base/serial/Makefile
@@ -1,7 +1,8 @@
 include ../../Make.inc
 
 
-FOBJS = psb_lsame.o    
+FOBJS = psb_lsame.o   psb_dsymbmm.o psb_dnumbmm.o \
+	psb_drwextd.o
 #	psb_regen_mod.o   psb_lsame.o psb_zspgetrow.o\
 #	psb_zcsmm.o psb_zcsmv.o  psb_zspgtdiag.o psb_zspgtblk.o\
 #	psb_zcsnmi.o psb_zcsrws.o psb_zcssm.o psb_zcssv.o psb_zspcnv.o\
diff --git a/base/serial/f77/Makefile b/base/serial/f77/Makefile
index cc3c4bad..5cc8e333 100644
--- a/base/serial/f77/Makefile
+++ b/base/serial/f77/Makefile
@@ -4,7 +4,7 @@ include ../../../Make.inc
 # The object files
 #
 FOBJS =  daxpby.o 	saxpby.o \
-	caxpby.o  zaxpby.o 
+	caxpby.o  zaxpby.o smmp.o
 
 
 #	clpupd.o ccsmm.o cswmm.o ccsnmi.o ccsrws.o\
diff --git a/base/serial/psb_dnumbmm.f90 b/base/serial/psb_dnumbmm.f90
index 428ebdb0..ea943f1e 100644
--- a/base/serial/psb_dnumbmm.f90
+++ b/base/serial/psb_dnumbmm.f90
@@ -41,45 +41,146 @@
 !
 
 subroutine psb_dnumbmm(a,b,c)
-  use psb_spmat_type
+  use psb_mat_mod
+  use psb_string_mod
   use psb_serial_mod, psb_protect_name => psb_dnumbmm
-  implicit none
+  implicit none 
 
-  type(psb_dspmat_type)         :: a,b,c
+  type(psb_d_sparse_mat), intent(in) :: a,b
+  type(psb_d_sparse_mat), intent(inout)  :: c
+  integer               :: info
+  integer               :: err_act
+  character(len=*), parameter ::  name='psb_numbmm'
+
+  call psb_erractionsave(err_act)
+  info = 0
+
+  if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
+    info = 1121
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  select type(aa=>c%a)
+  type is (psb_d_csr_sparse_mat)
+    call psb_numbmm(a%a,b%a,aa)
+  class default
+    info = 1121
+    call psb_errpush(info,name)
+    goto 9999
+  end select
+
+  call c%set_asb()
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act == psb_act_abort_) then
+     call psb_error()
+     return
+  end if
+  return
+
+end subroutine psb_dnumbmm
+
+subroutine psb_dbase_numbmm(a,b,c)
+  use psb_mat_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_dbase_numbmm
+  implicit none 
+
+  class(psb_d_base_sparse_mat), intent(in) :: a,b
+  type(psb_d_csr_sparse_mat), intent(inout)  :: c
+  integer, allocatable  :: itemp(:)
+  integer               :: nze, ma,na,mb,nb  
+  character(len=20)     :: name
   real(psb_dpk_), allocatable :: temp(:)
-  integer                       :: info
-  logical                   :: csra, csrb
-  
-  allocate(temp(max(a%m,a%k,b%m,b%k)),stat=info)
-  if (info /= 0) then
-    return
+  integer               :: info
+  integer               :: err_act
+  name='psb_numbmm'
+  call psb_erractionsave(err_act)
+  info = 0
+
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
   endif
-  call psb_realloc(size(c%ia1),c%aspk,info)
+  allocate(temp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= 0) then 
+    info = 4000 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+
   !
   ! Note: we still have to test about possible performance hits. 
   !
   !
-  csra = (psb_toupper(a%fida(1:3))=='CSR')
-  csrb = (psb_toupper(b%fida(1:3))=='CSR')
-
-  if (csra.and.csrb) then 
-    call numbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,a%aspk,&
-         & b%ia2,b%ia1,0,b%aspk,&
-         & c%ia2,c%ia1,0,c%aspk,temp)
-  else
-    call inner_numbmm(a,b,c,temp,info)
-    if (info /= 0) then 
-      write(0,*) 'Error ',info,' from inner numbmm'
-    end if
+  call psb_ensure_size(size(c%ja),c%val,info)
+  select type(a)
+  type is (psb_d_csr_sparse_mat) 
+    select type(b)
+    type is (psb_d_csr_sparse_mat) 
+      call csr_numbmm(a,b,c,temp,info)
+    class default
+      call gen_numbmm(a,b,c,temp,info)
+    end select
+  class default
+    call gen_numbmm(a,b,c,temp,info)
+  end select
+  
+  if (info /= 0) then 
+    call psb_errpush(info,name)
+    goto 9999
   end if
-  call psb_sp_setifld(psb_spmat_asb_,psb_state_,c,info)
+
+  call c%set_asb()
   deallocate(temp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act == psb_act_abort_) then
+     call psb_error()
+     return
+  end if
   return
 
 contains 
 
-  subroutine inner_numbmm(a,b,c,temp,info)
-    type(psb_dspmat_type) :: a,b,c
+  subroutine csr_numbmm(a,b,c,temp,info)
+    type(psb_d_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_d_csr_sparse_mat), intent(inout) :: c
+    real(psb_dpk_)                          :: temp(:)
+    integer, intent(out)                    :: info
+    integer               :: nze, ma,na,mb,nb
+
+    info = 0 
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    call numbmm(ma,na,nb,a%irp,a%ja,0,a%val,&
+         & b%irp,b%ja,0,b%val,&
+         & c%irp,c%ja,0,c%val,temp)
+
+    
+  end subroutine csr_numbmm
+
+  subroutine gen_numbmm(a,b,c,temp,info)
+    class(psb_d_base_sparse_mat), intent(in)  :: a,b
+    type(psb_d_csr_sparse_mat), intent(inout) :: c
     integer               :: info
     real(psb_dpk_)      :: temp(:)
     integer, allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
@@ -87,13 +188,14 @@ contains
     integer  :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
     real(psb_dpk_)      :: ajj
 
-    n = a%m
-    m = a%k 
-    l = b%k 
+    n = a%get_nrows()
+    m = a%get_ncols() 
+    l = b%get_ncols()
     maxlmn = max(l,m,n)
     allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
          & aval(maxlmn),bval(maxlmn), stat=info)
-    if (info /= 0) then 
+    if (info /= 0) then
+      info = 4000
       return
     endif
 
@@ -105,8 +207,7 @@ contains
     minmn = min(m,n)
     do  i = 1,n
 
-      call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
-
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
       do jj=1, nazr
         j=iacl(jj)
         ajj = aval(jj)
@@ -116,7 +217,7 @@ contains
             return
           
         endif
-        call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
         do k=1,nbzr
           if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
             write(0,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
@@ -127,19 +228,19 @@ contains
           endif
         enddo
       end do
-      do  j = c%ia2(i),c%ia2(i+1)-1
-        if((c%ia1(j)<1).or. (c%ia1(j) > maxlmn))  then 
-          write(0,*) ' NUMBMM: output problem',i,j,c%ia1(j),maxlmn
+      do  j = c%irp(i),c%irp(i+1)-1
+        if((c%ja(j)<1).or. (c%ja(j) > maxlmn))  then 
+          write(0,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
             info = 3
             return
         else
-          c%aspk(j) = temp(c%ia1(j))
-          temp(c%ia1(j)) = dzero
+          c%val(j) = temp(c%ja(j))
+          temp(c%ja(j)) = dzero
         endif
       end do
     end do
 
+    
+  end subroutine gen_numbmm
 
-  end subroutine inner_numbmm
-
-end subroutine psb_dnumbmm
+end subroutine psb_dbase_numbmm
diff --git a/base/serial/psb_drwextd.f90 b/base/serial/psb_drwextd.f90
index 3d4cc214..bbf0ea54 100644
--- a/base/serial/psb_drwextd.f90
+++ b/base/serial/psb_drwextd.f90
@@ -39,127 +39,209 @@
 !
 !
 subroutine psb_drwextd(nr,a,info,b,rowscale)
-  use psb_spmat_type
   use psb_error_mod
   use psb_string_mod
   use psb_serial_mod, psb_protect_name => psb_drwextd
   implicit none
 
   ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
-  integer, intent(in)                            :: nr
-  type(psb_dspmat_type), intent(inout)           :: a
-  integer,intent(out)                            :: info
-  type(psb_dspmat_type), intent(in), optional    :: b
-  logical,intent(in), optional                   :: rowscale
+  integer, intent(in)                          :: nr
+  type(psb_d_sparse_mat), intent(inout)        :: a
+  integer,intent(out)                          :: info
+  type(psb_d_sparse_mat), intent(in), optional :: b
+  logical,intent(in), optional                 :: rowscale
 
   integer :: i,j,ja,jb,err_act,nza,nzb
   character(len=20)                 :: name, ch_err
+  type(psb_d_coo_sparse_mat)        :: actmp
   logical  rowscale_ 
 
   name='psb_drwextd'
   info  = 0
   call psb_erractionsave(err_act)
 
+  if (nr > a%get_nrows()) then 
+    select type(aa=> a%a) 
+    type is (psb_d_csr_sparse_mat)
+      if (present(b)) then 
+        call psb_rwextd(nr,aa,info,b%a,rowscale)
+      else
+        call psb_rwextd(nr,aa,info,rowscale=rowscale)
+      end if
+    type is (psb_d_coo_sparse_mat) 
+      if (present(b)) then 
+        call psb_rwextd(nr,aa,info,b%a,rowscale=rowscale)
+      else
+        call psb_rwextd(nr,aa,info,rowscale=rowscale)
+      end if
+    class default
+      call aa%mv_to_coo(actmp,info)
+      if (info == 0) then 
+        if (present(b)) then 
+          call psb_rwextd(nr,actmp,info,b%a,rowscale=rowscale)
+        else
+          call psb_rwextd(nr,actmp,info,rowscale=rowscale)
+        end if
+      end if
+      if (info == 0) call aa%mv_from_coo(actmp,info)
+    end select
+  end if
+  if (info /= 0) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act == psb_act_abort_) then
+     call psb_error()
+     return
+  end if
+  return
+
+end subroutine psb_drwextd
+subroutine psb_dbase_rwextd(nr,a,info,b,rowscale)
+  use psb_error_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_dbase_rwextd
+  implicit none
+
+  ! Extend matrix A up to NR rows with empty ones (i.e.: all zeroes)
+  integer, intent(in)                                :: nr
+  class(psb_d_base_sparse_mat), intent(inout)        :: a
+  integer,intent(out)                                :: info
+  class(psb_d_base_sparse_mat), intent(in), optional :: b
+  logical,intent(in), optional                       :: rowscale
+
+  integer :: i,j,ja,jb,err_act,nza,nzb, ma, mb, na, nb
+  character(len=20)                 :: name, ch_err
+  logical  rowscale_ 
+
+  name='psb_dbase_rwextd'
+  info  = 0
+  call psb_erractionsave(err_act)
+
   if (present(rowscale)) then 
     rowscale_ = rowscale
   else
     rowscale_ = .true.
   end if
 
-  if (nr > a%m) then 
-    if (psb_toupper(a%fida) == 'CSR') then 
-      call psb_ensure_size(nr+1,a%ia2,info)
-      if (present(b)) then 
-        nzb = psb_sp_get_nnzeros(b)
-        call psb_ensure_size(size(a%ia1)+nzb,a%ia1,info)
-        call psb_ensure_size(size(a%aspk)+nzb,a%aspk,info)
-        if (psb_toupper(b%fida)=='CSR') then 
-
-          do i=1, min(nr-a%m,b%m)
-            a%ia2(a%m+i+1) =  a%ia2(a%m+i) + b%ia2(i+1) - b%ia2(i)
-            ja = a%ia2(a%m+i)
-            jb = b%ia2(i)
-            do 
-              if (jb >=  b%ia2(i+1)) exit
-              a%aspk(ja) = b%aspk(jb)
-              a%ia1(ja) = b%ia1(jb)
-              ja = ja + 1
-              jb = jb + 1
-            end do
-          end do
-          do j=i,nr-a%m
-            a%ia2(a%m+i+1) = a%ia2(a%m+i)
+  ma = a%get_nrows()
+  na = a%get_ncols()
+
+
+  select type(a) 
+  type is (psb_d_csr_sparse_mat)
+
+    call psb_ensure_size(nr+1,a%irp,info)
+
+    if (present(b)) then 
+      mb = b%get_nrows()
+      nb = b%get_ncols()
+      nzb = b%get_nzeros()
+
+      select type (b) 
+      type is (psb_d_csr_sparse_mat)
+        call psb_ensure_size(size(a%ja)+nzb,a%ja,info)
+        call psb_ensure_size(size(a%val)+nzb,a%val,info)
+        do i=1, min(nr-ma,mb)
+          a%irp(ma+i+1) =  a%irp(ma+i) + b%irp(i+1) - b%irp(i)
+          ja = a%irp(ma+i)
+          jb = b%irp(i)
+          do 
+            if (jb >=  b%irp(i+1)) exit
+            a%val(ja) = b%val(jb)
+            a%ja(ja)  = b%ja(jb)
+            ja = ja + 1
+            jb = jb + 1
           end do
-        else 
-          write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
-        endif
-      else
-        do i=a%m+2,nr+1
-          a%ia2(i) = a%ia2(i-1)
         end do
-      end if
-      a%m = nr
-      a%k = max(a%k,b%k)
+        do j=i,nr-ma
+          a%irp(ma+i+1) = a%irp(ma+i)
+        end do
+        class default 
 
-    else if (psb_toupper(a%fida) == 'COO') then 
+        write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
+      end select
+      call a%set_ncols(max(na,nb))
 
-      if (present(b)) then 
-        nza = psb_sp_get_nnzeros(a)
-        nzb = psb_sp_get_nnzeros(b)
-        call psb_sp_reall(a,nza+nzb,info)
-        if (psb_toupper(b%fida)=='COO') then 
-          if (rowscale_) then 
-            do j=1,nzb
-              if ((a%m + b%ia1(j)) <= nr) then 
-                nza = nza + 1
-                a%ia1(nza)  = a%m + b%ia1(j)
-                a%ia2(nza)  = b%ia2(j)
-                a%aspk(nza) = b%aspk(j)
-              end if
-            enddo
-          else
-            do j=1,nzb
-              if ((b%ia1(j)) <= nr) then 
-                nza = nza + 1
-                a%ia1(nza)  = b%ia1(j)
-                a%ia2(nza)  = b%ia2(j)
-                a%aspk(nza) = b%aspk(j)                
-              endif
-            enddo
-          endif
-          a%infoa(psb_nnz_) = nza
-        else if(psb_toupper(b%fida)=='CSR') then 
-          do i=1, min(nr-a%m,b%m)
-            do 
-              jb = b%ia2(i)
-              if (jb >=  b%ia2(i+1)) exit
-              nza = nza + 1 
-              a%aspk(nza) = b%aspk(jb)
-              a%ia1(nza)  = a%m + i
-              a%ia2(nza)  = b%ia1(jb)
-              jb = jb + 1
-            end do
-          end do
-          a%infoa(psb_nnz_) = nza
-        else
-          write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
-        endif
-      endif
-      a%m = nr
-      a%k = max(a%k,b%k)
-    else if (a%fida == 'JAD') then 
-       info=135
-       ch_err=a%fida(1:3)
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
     else
-       info=136
-       ch_err=a%fida(1:3)
-       call psb_errpush(info,name,a_err=ch_err)
-       goto 9999
+
+      do i=ma+2,nr+1
+        a%irp(i) = a%irp(i-1)
+      end do
+
     end if
 
-  end if
+    call a%set_nrows(nr)
+
+
+  type is (psb_d_coo_sparse_mat) 
+    nza = a%get_nzeros()
+
+    if (present(b)) then 
+      mb = b%get_nrows()
+      nb = b%get_ncols()
+      nzb = b%get_nzeros()
+      call a%reallocate(nza+nzb)
+
+      select type(b)
+      type is (psb_d_coo_sparse_mat) 
+
+        if (rowscale_) then 
+          do j=1,nzb
+            if ((ma + b%ia(j)) <= nr) then 
+              nza = nza + 1
+              a%ia(nza)  = ma + b%ia(j)
+              a%ja(nza)  = b%ja(j)
+              a%val(nza) = b%val(j)
+            end if
+          enddo
+        else
+          do j=1,nzb
+            if ((ma + b%ia(j)) <= nr) then 
+              nza = nza + 1
+              a%ia(nza)  = b%ia(j)
+              a%ja(nza)  = b%ja(j)
+              a%val(nza) = b%val(j)
+            end if
+          enddo
+        endif
+        call a%set_nzeros(nza)
+
+      type is (psb_d_csr_sparse_mat) 
+        
+        do i=1, min(nr-ma,mb)
+          do 
+            jb = b%irp(i)
+            if (jb >=  b%irp(i+1)) exit
+            nza = nza + 1 
+            a%val(nza) = b%val(jb)
+            a%ia(nza)  = ma + i
+            a%ja(nza)  = b%ja(jb)
+            jb = jb + 1
+          end do
+        end do
+        call a%set_nzeros(nza)
+
+      class default 
+        write(0,*) 'Implement SPGETBLK in RWEXTD!!!!!!!'
+
+      end select
+
+      call a%set_ncols(max(na,nb))
+    endif
+
+    call a%set_nrows(nr)
+
+  class default 
+    info = 135
+    ch_err=a%get_fmt()
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+
+  end select
 
   call psb_erractionrestore(err_act)
   return
@@ -167,9 +249,9 @@ subroutine psb_drwextd(nr,a,info,b,rowscale)
 9999 continue
   call psb_erractionrestore(err_act)
   if (err_act == psb_act_abort_) then
-     call psb_error()
-     return
+    call psb_error()
+    return
   end if
   return
 
-end subroutine psb_drwextd
+end subroutine psb_dbase_rwextd
diff --git a/base/serial/psb_dsymbmm.f90 b/base/serial/psb_dsymbmm.f90
index 251e03da..211bed40 100644
--- a/base/serial/psb_dsymbmm.f90
+++ b/base/serial/psb_dsymbmm.f90
@@ -40,62 +40,107 @@
 !
 
 subroutine psb_dsymbmm(a,b,c,info)
-  use psb_spmat_type
+  use psb_mat_mod
   use psb_string_mod
   use psb_serial_mod, psb_protect_name => psb_dsymbmm
   implicit none 
 
-  type(psb_dspmat_type) :: a,b,c
+  type(psb_d_sparse_mat), intent(in)    :: a,b
+  type(psb_d_sparse_mat), intent(out)   :: c
+  integer, intent(out)                  :: info
+  type(psb_d_csr_sparse_mat), allocatable :: ccsr
+  integer               :: err_act
+  character(len=*), parameter ::  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = 0
+
+  if ((a%is_null()) .or.(b%is_null())) then
+    info = 1121
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  allocate(ccsr, stat=info)
+  if (info /= 0)  then
+    info = 4000 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call psb_symbmm(a%a,b%a,ccsr,info)
+  if (info /= 0) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call move_alloc(ccsr,c%a)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+  if (err_act == psb_act_abort_) then
+     call psb_error()
+     return
+  end if
+  return
+
+end subroutine psb_dsymbmm
+
+subroutine psb_dbase_symbmm(a,b,c,info)
+  use psb_mat_mod
+  use psb_serial_mod, psb_protect_name => psb_dbase_symbmm
+  implicit none 
+
+  class(psb_d_base_sparse_mat), intent(in) :: a,b
+  type(psb_d_csr_sparse_mat), intent(out)  :: c
+  integer, intent(out)                     :: info
   integer, allocatable  :: itemp(:)
-  integer               :: nze,info
-
-  interface 
-    subroutine symbmm (n, m, l, ia, ja, diaga, &
-         & ib, jb, diagb, ic, jc, diagc, index)
-      integer  n,m,l,  ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
-           & diagc,  index(*)
-      integer, allocatable :: ic(:),jc(:)
-    end subroutine symbmm 
-  end interface
+  integer               :: nze, ma,na,mb,nb
   character(len=20)     :: name
   integer               :: err_act
-  logical               :: csra, csrb
   name='psb_symbmm'
   call psb_erractionsave(err_act)
+  info = 0
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
 
-  csra = (psb_toupper(a%fida(1:3))=='CSR')
-  csrb = (psb_toupper(b%fida(1:3))=='CSR')
 
-  if (b%m /= a%k) then 
-    write(0,*) 'Mismatch in SYMBMM: ',a%m,a%k,b%m,b%k
+  if ( mb /= na ) then 
+    write(0,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
   endif
-  allocate(itemp(max(a%m,a%k,b%m,b%k)),stat=info)    
+  allocate(itemp(max(ma,na,mb,nb)),stat=info)    
   if (info /= 0) then 
-    return
+    info = 4000 
+    call psb_Errpush(info,name)
+    goto 9999
   endif
-  nze = max(a%m+1,2*a%m)
-  call psb_sp_reall(c,nze,info)
   !
   ! Note: we need to test whether there is a performance impact 
   !       in not using the original Douglas & Bank code. 
   !
-  if (csra.and.csrb) then 
-    call symbmm(a%m,a%k,b%k,a%ia2,a%ia1,0,&
-         & b%ia2,b%ia1,0,&
-         & c%ia2,c%ia1,0,itemp)
-  else 
-    call inner_symbmm(a,b,c,itemp,info)
-  endif
-  call psb_realloc(size(c%ia1),c%aspk,info)
-
-  c%pl(1) = 0
-  c%pr(1) = 0
-  c%m=a%m
-  c%k=b%k
-  c%fida='CSR'
-  c%descra='GUN'
-  
+  select type(a)
+  type is (psb_d_csr_sparse_mat) 
+    select type(b)
+    type is (psb_d_csr_sparse_mat) 
+      call csr_symbmm(a,b,c,itemp,info)
+    class default
+      call gen_symbmm(a,b,c,itemp,info)
+    end select
+  class default
+    call gen_symbmm(a,b,c,itemp,info)
+  end select
+
+  if (info /= 0) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+
+  call psb_realloc(size(c%ja),c%val,info)
   deallocate(itemp) 
+
   call psb_erractionrestore(err_act)
   return
 
@@ -109,86 +154,119 @@ subroutine psb_dsymbmm(a,b,c,info)
 
 contains
 
-  subroutine inner_symbmm(a,b,c,index,info)
-    type(psb_dspmat_type) :: a,b,c
+  subroutine csr_symbmm(a,b,c,itemp,info)
+    type(psb_d_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_d_csr_sparse_mat), intent(out) :: c
+    integer                                 :: itemp(:) 
+    integer, intent(out)                    :: info
+    interface 
+      subroutine symbmm (n, m, l, ia, ja, diaga, &
+           & ib, jb, diagb, ic, jc, diagc, index)
+        integer  n,m,l,  ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
+             & diagc,  index(*)
+        integer, allocatable :: ic(:),jc(:)
+      end subroutine symbmm
+    end interface
+    integer               :: nze, ma,na,mb,nb
+
+    info = 0 
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+    call symbmm(ma,na,nb,a%irp,a%ja,0,&
+         & b%irp,b%ja,0,&
+         & c%irp,c%ja,0,itemp)
+    
+  end subroutine csr_symbmm
+  subroutine gen_symbmm(a,b,c,index,info)
+    class(psb_d_base_sparse_mat), intent(in)  :: a,b
+    type(psb_d_csr_sparse_mat), intent(out) :: c
     integer               :: index(:),info
     integer, allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
     real(psb_dpk_), allocatable :: aval(:),bval(:)
     integer  :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
+    integer               :: nze, ma,na,mb,nb
 
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
 
-    n = a%m
-    m = a%k 
-    l = b%k 
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+
+    n = ma
+    m = na 
+    l = nb 
     maxlmn = max(l,m,n)
 
     allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
          & aval(maxlmn),bval(maxlmn), stat=info)
     if (info /= 0) then 
+      info = 4000
       return
     endif
 
-
-    if (size(c%ia2) < n+1) then 
-
-      call psb_realloc(n+1,c%ia2,info)
-    endif
     do i=1,maxlmn
       index(i)=0
     end do
 
-      c%ia2(1)=1
-      minlm = min(l,m)
-      minmn = min(m,n)
-
-      main: do  i=1,n
-        istart=-1
-        length=0
-        call psb_sp_getrow(i,a,nazr,iarw,iacl,aval,info)
-        do jj=1, nazr
-          
-          j=iacl(jj)
-          
-          if ((j<1).or.(j>m)) then 
-            write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
-            info = 1
+    c%irp(1)=1
+    minlm = min(l,m)
+    minmn = min(m,n)
+
+    main: do  i=1,n
+      istart=-1
+      length=0
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
+      do jj=1, nazr
+
+        j=iacl(jj)
+
+        if ((j<1).or.(j>m)) then 
+          write(0,*) ' SymbMM: Problem with A ',i,jj,j,m
+          info = 1
+          return
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
+            info=2
             return
-          endif
-          call psb_sp_getrow(j,b,nbzr,ibrw,ibcl,bval,info)
-          do k=1,nbzr
-            if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
-                write(0,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
-                info=2
-                return
-            else
-              if(index(ibcl(k)) == 0) then
-                index(ibcl(k))=istart
-                istart=ibcl(k)
-                length=length+1
-              endif
+          else
+            if(index(ibcl(k)) == 0) then
+              index(ibcl(k))=istart
+              istart=ibcl(k)
+              length=length+1
             endif
-          end do
+          endif
         end do
+      end do
 
-        c%ia2(i+1)=c%ia2(i)+length
-        
-        if (c%ia2(i+1) > size(c%ia1)) then 
-          if (n > (2*i)) then 
-            nze = max(c%ia2(i+1), c%ia2(i)*((n+i-1)/i))
-          else
-            nze = max(c%ia2(i+1), nint((dble(c%ia2(i))*(dble(n)/i)))   )
-          endif 
-          call psb_realloc(nze,c%ia1,info)
-        end if 
-        do j= c%ia2(i),c%ia2(i+1)-1
-          c%ia1(j)=istart
-          istart=index(istart)
-          index(c%ia1(j))=0
-        end do
-        call psb_msort(c%ia1(c%ia2(i):c%ia2(i)+length-1))
-        index(i) = 0
-      end do main
+      c%irp(i+1)=c%irp(i)+length
 
-  end subroutine inner_symbmm
+      if (c%irp(i+1) > size(c%ja)) then 
+        if (n > (2*i)) then 
+          nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
+        else
+          nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i)))   )
+        endif
+        call psb_realloc(nze,c%ja,info)
+      end if
+      do j= c%irp(i),c%irp(i+1)-1
+        c%ja(j)=istart
+        istart=index(istart)
+        index(c%ja(j))=0
+      end do
+      call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
+      index(i) = 0
+    end do main
 
-end subroutine psb_dsymbmm
+  end subroutine gen_symbmm
+
+end subroutine psb_dbase_symbmm
diff --git a/configure b/configure
index 38140f9f..08e1699c 100755
--- a/configure
+++ b/configure
@@ -771,7 +771,7 @@ SHELL'
 ac_subst_files=''
 ac_user_opts='
 enable_option_checking
-with_serial_mpi
+enable_serial
 with_ccopt
 with_fcopt
 with_f90copt
@@ -1429,14 +1429,14 @@ Optional Features:
   --disable-option-checking  ignore unrecognized --enable/--with options
   --disable-FEATURE       do not include FEATURE (same as --enable-FEATURE=no)
   --enable-FEATURE[=ARG]  include FEATURE [ARG=yes]
+  --enable-serial         Specify whether to enable a fake mpi library to run
+                          in serial mode.
   --disable-dependency-tracking  speeds up one-time build
   --enable-dependency-tracking   do not reject slow dependency extractors
 
 Optional Packages:
   --with-PACKAGE[=ARG]    use PACKAGE [ARG=yes]
   --without-PACKAGE       do not use PACKAGE (same as --with-PACKAGE=no)
-  --with-serial-mpi       Specify whether to enable a fake mpi library to run
-                          in serial mode. --with-serial-mpi={yes|no}
   --with-ccopt            additional CCOPT flags to be added: will prepend
                           to CCOPT
   --with-fcopt            additional FCOPT flags to be added: will prepend
@@ -3394,14 +3394,11 @@ fi
 
 { $as_echo "$as_me:$LINENO: checking whether we want serial (fake) mpi" >&5
 $as_echo_n "checking whether we want serial (fake) mpi... " >&6; }
+# Check whether --enable-serial was given.
+if test "${enable_serial+set}" = set; then
+  enableval=$enable_serial;
+pac_cv_serial_mpi="yes";
 
-# Check whether --with-serial-mpi was given.
-if test "${with_serial_mpi+set}" = set; then
-  withval=$with_serial_mpi;
-pac_cv_serial_mpi="${withval}";
-
-else
-  pac_cv_serial_mpi="no";
 
 fi
 
@@ -3409,7 +3406,8 @@ if test x"$pac_cv_serial_mpi" == x"yes" ; then
    { $as_echo "$as_me:$LINENO: result: yes." >&5
 $as_echo "yes." >&6; }
 else
-   { $as_echo "$as_me:$LINENO: result: no." >&5
+ pac_cv_serial_mpi="no";
+ { $as_echo "$as_me:$LINENO: result: no." >&5
 $as_echo "no." >&6; }
 fi
 
@@ -7119,7 +7117,7 @@ if test "X$FCOPT" == "X" ; then
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options
-   FCOPT="$FCOPT  -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
+   FCOPT="$FCOPT  -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
 fi
 FFLAGS="${FCOPT}"
 
@@ -7152,7 +7150,7 @@ else
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options
-     F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
+     F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
      EXTRA_OPT="-mismatch_all"
      F03COPT="${F90COPT}"
 else
diff --git a/configure.ac b/configure.ac
index a97731b4..e1be73aa 100755
--- a/configure.ac
+++ b/configure.ac
@@ -399,7 +399,7 @@ if test "X$FCOPT" == "X" ; then
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options 
-   FCOPT="$FCOPT  -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
+   FCOPT="$FCOPT  -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
 fi
 FFLAGS="${FCOPT}"
 
@@ -432,7 +432,7 @@ else
 fi
 if test "X$psblas_cv_fc" == X"nag" ; then
    # Add needed options 
-     F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv"
+     F90COPT="$F90COPT -dcfuns -f2003 -wmismatch=mpi_scatterv,mpi_alltoallv,mpi_gatherv,mpi_allgatherv"
      EXTRA_OPT="-mismatch_all"
      F03COPT="${F90COPT}"
 else