diff --git a/base/modules/psb_c_base_mat_mod.f90 b/base/modules/psb_c_base_mat_mod.f90
index 90db864c9..a9b7ce5e0 100644
--- a/base/modules/psb_c_base_mat_mod.f90
+++ b/base/modules/psb_c_base_mat_mod.f90
@@ -42,7 +42,7 @@ module psb_c_base_mat_mod
   !! The psb_c_base_sparse_mat type, extending psb_base_sparse_mat,
   !! defines a middle level  complex(psb_spk_) sparse matrix object.
   !! This class object itself does not have any additional members
-  !! with respect to those of the base class. No methods can be fully
+  !! with respect to those of the base class. Most methods cannot be fully
   !! implemented at this level, but we can define the interface for the
   !! computational methods requiring the knowledge of the underlying
   !! field, such as the matrix-vector product; this interface is defined,
@@ -78,6 +78,7 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: mold          => psb_c_base_mold 
     procedure, pass(a) :: clone         => psb_c_base_clone
     procedure, pass(a) :: make_nonunit  => psb_c_base_make_nonunit
+    procedure, pass(a) :: clean_zeros   => psb_c_base_clean_zeros
     
     !
     ! Transpose methods: defined here but not implemented. 
@@ -161,6 +162,7 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: get_nz_row   => psb_c_coo_get_nz_row
     procedure, pass(a) :: fix          => psb_c_fix_coo
     procedure, pass(a) :: trim         => psb_c_coo_trim
+    procedure, pass(a) :: clean_zeros  => psb_c_coo_clean_zeros
     procedure, pass(a) :: print        => psb_c_coo_print
     procedure, pass(a) :: free         => c_coo_free
     procedure, pass(a) :: mold         => psb_c_coo_mold
@@ -697,6 +699,18 @@ module psb_c_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_c_base_mv_from_fmt
   end interface
+  !
+  !> 
+  !! \memberof  psb_c_base_sparse_mat
+  !! \see psb_c_base_mat_mod::psb_c_base_clean_zeros
+  !
+  interface
+    subroutine  psb_c_base_clean_zeros(a, info)
+      import :: psb_ipk_, psb_c_base_sparse_mat
+      class(psb_c_base_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)              :: info
+    end subroutine psb_c_base_clean_zeros
+  end interface
   
   !
   !> Function  transp:
@@ -1235,6 +1249,18 @@ module psb_c_base_mat_mod
       class(psb_c_coo_sparse_mat), intent(inout) :: a
     end subroutine psb_c_coo_trim
   end interface
+  !
+  !> 
+  !! \memberof  psb_c_coo_sparse_mat
+  !! \see psb_c_base_mat_mod::psb_c_base_clean_zeros
+  !
+  interface
+    subroutine  psb_c_coo_clean_zeros(a,info)
+      import :: psb_ipk_, psb_c_coo_sparse_mat
+      class(psb_c_coo_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)             :: info
+    end subroutine psb_c_coo_clean_zeros
+  end interface
   
   !
   !> 
diff --git a/base/modules/psb_c_comm_mod.f90 b/base/modules/psb_c_comm_mod.f90
index 713dedce5..5cb8ff7ec 100644
--- a/base/modules/psb_c_comm_mod.f90
+++ b/base/modules/psb_c_comm_mod.f90
@@ -97,7 +97,7 @@ module psb_c_comm_mod
   interface psb_scatter
     subroutine  psb_cscatterm(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      complex(psb_spk_), intent(out) :: locx(:,:)
+      complex(psb_spk_), intent(out), allocatable :: locx(:,:)
       complex(psb_spk_), intent(in)  :: globx(:,:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
@@ -105,7 +105,7 @@ module psb_c_comm_mod
     end subroutine psb_cscatterm
     subroutine  psb_cscatterv(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      complex(psb_spk_), intent(out) :: locx(:)
+      complex(psb_spk_), intent(out), allocatable :: locx(:)
       complex(psb_spk_), intent(in)  :: globx(:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
diff --git a/base/modules/psb_c_mat_mod.f90 b/base/modules/psb_c_mat_mod.f90
index 973997527..48652ca7f 100644
--- a/base/modules/psb_c_mat_mod.f90
+++ b/base/modules/psb_c_mat_mod.f90
@@ -44,7 +44,7 @@
 ! the functionalities to have the encapsulated class change its
 ! type dynamically, and to extract/input an inner object.
 !
-! A sparse matric has a state corresponding to its progression
+! A sparse matrix has a state corresponding to its progression
 ! through the application life.
 ! In particular, computational methods can only be invoked when
 ! the matrix is in the ASSEMBLED state, whereas the other states are
@@ -133,6 +133,7 @@ module psb_c_mat_mod
     procedure, pass(a) :: m_csclip    => psb_c_csclip
     procedure, pass(a) :: b_csclip    => psb_c_b_csclip
     generic, public    :: csclip      => b_csclip, m_csclip
+    procedure, pass(a) :: clean_zeros => psb_c_clean_zeros
     procedure, pass(a) :: reall       => psb_c_reallocate_nz
     procedure, pass(a) :: get_neigh   => psb_c_get_neigh
     procedure, pass(a) :: reinit      => psb_c_reinit
@@ -1177,5 +1178,15 @@ contains
 
   end function psb_c_get_nz_row
 
+  subroutine psb_c_clean_zeros(a,info)
+    implicit none 
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_cspmat_type), intent(inout) :: a
+
+    info = 0 
+    if (allocated(a%a)) call a%a%clean_zeros(info)
+
+  end subroutine psb_c_clean_zeros
+
 
 end module psb_c_mat_mod
diff --git a/base/modules/psb_d_base_mat_mod.f90 b/base/modules/psb_d_base_mat_mod.f90
index 304b02b31..c0d9de7e1 100644
--- a/base/modules/psb_d_base_mat_mod.f90
+++ b/base/modules/psb_d_base_mat_mod.f90
@@ -42,7 +42,7 @@ module psb_d_base_mat_mod
   !! The psb_d_base_sparse_mat type, extending psb_base_sparse_mat,
   !! defines a middle level  real(psb_dpk_) sparse matrix object.
   !! This class object itself does not have any additional members
-  !! with respect to those of the base class. No methods can be fully
+  !! with respect to those of the base class. Most methods cannot be fully
   !! implemented at this level, but we can define the interface for the
   !! computational methods requiring the knowledge of the underlying
   !! field, such as the matrix-vector product; this interface is defined,
@@ -78,6 +78,7 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: mold          => psb_d_base_mold 
     procedure, pass(a) :: clone         => psb_d_base_clone
     procedure, pass(a) :: make_nonunit  => psb_d_base_make_nonunit
+    procedure, pass(a) :: clean_zeros   => psb_d_base_clean_zeros
     
     !
     ! Transpose methods: defined here but not implemented. 
@@ -161,6 +162,7 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: get_nz_row   => psb_d_coo_get_nz_row
     procedure, pass(a) :: fix          => psb_d_fix_coo
     procedure, pass(a) :: trim         => psb_d_coo_trim
+    procedure, pass(a) :: clean_zeros  => psb_d_coo_clean_zeros
     procedure, pass(a) :: print        => psb_d_coo_print
     procedure, pass(a) :: free         => d_coo_free
     procedure, pass(a) :: mold         => psb_d_coo_mold
@@ -697,6 +699,18 @@ module psb_d_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_d_base_mv_from_fmt
   end interface
+  !
+  !> 
+  !! \memberof  psb_d_base_sparse_mat
+  !! \see psb_d_base_mat_mod::psb_d_base_clean_zeros
+  !
+  interface
+    subroutine  psb_d_base_clean_zeros(a, info)
+      import :: psb_ipk_, psb_d_base_sparse_mat
+      class(psb_d_base_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)              :: info
+    end subroutine psb_d_base_clean_zeros
+  end interface
   
   !
   !> Function  transp:
@@ -1235,6 +1249,18 @@ module psb_d_base_mat_mod
       class(psb_d_coo_sparse_mat), intent(inout) :: a
     end subroutine psb_d_coo_trim
   end interface
+  !
+  !> 
+  !! \memberof  psb_d_coo_sparse_mat
+  !! \see psb_d_base_mat_mod::psb_d_base_clean_zeros
+  !
+  interface
+    subroutine  psb_d_coo_clean_zeros(a,info)
+      import :: psb_ipk_, psb_d_coo_sparse_mat
+      class(psb_d_coo_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)             :: info
+    end subroutine psb_d_coo_clean_zeros
+  end interface
   
   !
   !> 
diff --git a/base/modules/psb_d_comm_mod.f90 b/base/modules/psb_d_comm_mod.f90
index e99345bdd..e4036f564 100644
--- a/base/modules/psb_d_comm_mod.f90
+++ b/base/modules/psb_d_comm_mod.f90
@@ -97,7 +97,7 @@ module psb_d_comm_mod
   interface psb_scatter
     subroutine  psb_dscatterm(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      real(psb_dpk_), intent(out) :: locx(:,:)
+      real(psb_dpk_), intent(out), allocatable :: locx(:,:)
       real(psb_dpk_), intent(in)  :: globx(:,:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
@@ -105,7 +105,7 @@ module psb_d_comm_mod
     end subroutine psb_dscatterm
     subroutine  psb_dscatterv(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      real(psb_dpk_), intent(out) :: locx(:)
+      real(psb_dpk_), intent(out), allocatable :: locx(:)
       real(psb_dpk_), intent(in)  :: globx(:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
diff --git a/base/modules/psb_d_mat_mod.f90 b/base/modules/psb_d_mat_mod.f90
index 1ad7d0307..c5b2999db 100644
--- a/base/modules/psb_d_mat_mod.f90
+++ b/base/modules/psb_d_mat_mod.f90
@@ -44,7 +44,7 @@
 ! the functionalities to have the encapsulated class change its
 ! type dynamically, and to extract/input an inner object.
 !
-! A sparse matric has a state corresponding to its progression
+! A sparse matrix has a state corresponding to its progression
 ! through the application life.
 ! In particular, computational methods can only be invoked when
 ! the matrix is in the ASSEMBLED state, whereas the other states are
@@ -133,6 +133,7 @@ module psb_d_mat_mod
     procedure, pass(a) :: m_csclip    => psb_d_csclip
     procedure, pass(a) :: b_csclip    => psb_d_b_csclip
     generic, public    :: csclip      => b_csclip, m_csclip
+    procedure, pass(a) :: clean_zeros => psb_d_clean_zeros
     procedure, pass(a) :: reall       => psb_d_reallocate_nz
     procedure, pass(a) :: get_neigh   => psb_d_get_neigh
     procedure, pass(a) :: reinit      => psb_d_reinit
@@ -1177,5 +1178,15 @@ contains
 
   end function psb_d_get_nz_row
 
+  subroutine psb_d_clean_zeros(a,info)
+    implicit none 
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_dspmat_type), intent(inout) :: a
+
+    info = 0 
+    if (allocated(a%a)) call a%a%clean_zeros(info)
+
+  end subroutine psb_d_clean_zeros
+
 
 end module psb_d_mat_mod
diff --git a/base/modules/psb_s_base_mat_mod.f90 b/base/modules/psb_s_base_mat_mod.f90
index 10d4e40b1..24f18a276 100644
--- a/base/modules/psb_s_base_mat_mod.f90
+++ b/base/modules/psb_s_base_mat_mod.f90
@@ -42,7 +42,7 @@ module psb_s_base_mat_mod
   !! The psb_s_base_sparse_mat type, extending psb_base_sparse_mat,
   !! defines a middle level  real(psb_spk_) sparse matrix object.
   !! This class object itself does not have any additional members
-  !! with respect to those of the base class. No methods can be fully
+  !! with respect to those of the base class. Most methods cannot be fully
   !! implemented at this level, but we can define the interface for the
   !! computational methods requiring the knowledge of the underlying
   !! field, such as the matrix-vector product; this interface is defined,
@@ -78,6 +78,7 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: mold          => psb_s_base_mold 
     procedure, pass(a) :: clone         => psb_s_base_clone
     procedure, pass(a) :: make_nonunit  => psb_s_base_make_nonunit
+    procedure, pass(a) :: clean_zeros   => psb_s_base_clean_zeros
     
     !
     ! Transpose methods: defined here but not implemented. 
@@ -161,6 +162,7 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: get_nz_row   => psb_s_coo_get_nz_row
     procedure, pass(a) :: fix          => psb_s_fix_coo
     procedure, pass(a) :: trim         => psb_s_coo_trim
+    procedure, pass(a) :: clean_zeros  => psb_s_coo_clean_zeros
     procedure, pass(a) :: print        => psb_s_coo_print
     procedure, pass(a) :: free         => s_coo_free
     procedure, pass(a) :: mold         => psb_s_coo_mold
@@ -697,6 +699,18 @@ module psb_s_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_s_base_mv_from_fmt
   end interface
+  !
+  !> 
+  !! \memberof  psb_s_base_sparse_mat
+  !! \see psb_s_base_mat_mod::psb_s_base_clean_zeros
+  !
+  interface
+    subroutine  psb_s_base_clean_zeros(a, info)
+      import :: psb_ipk_, psb_s_base_sparse_mat
+      class(psb_s_base_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)              :: info
+    end subroutine psb_s_base_clean_zeros
+  end interface
   
   !
   !> Function  transp:
@@ -1235,6 +1249,18 @@ module psb_s_base_mat_mod
       class(psb_s_coo_sparse_mat), intent(inout) :: a
     end subroutine psb_s_coo_trim
   end interface
+  !
+  !> 
+  !! \memberof  psb_s_coo_sparse_mat
+  !! \see psb_s_base_mat_mod::psb_s_base_clean_zeros
+  !
+  interface
+    subroutine  psb_s_coo_clean_zeros(a,info)
+      import :: psb_ipk_, psb_s_coo_sparse_mat
+      class(psb_s_coo_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)             :: info
+    end subroutine psb_s_coo_clean_zeros
+  end interface
   
   !
   !> 
diff --git a/base/modules/psb_s_comm_mod.f90 b/base/modules/psb_s_comm_mod.f90
index 15bb7bdc4..4a3bd0184 100644
--- a/base/modules/psb_s_comm_mod.f90
+++ b/base/modules/psb_s_comm_mod.f90
@@ -97,7 +97,7 @@ module psb_s_comm_mod
   interface psb_scatter
     subroutine  psb_sscatterm(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      real(psb_spk_), intent(out) :: locx(:,:)
+      real(psb_spk_), intent(out), allocatable :: locx(:,:)
       real(psb_spk_), intent(in)  :: globx(:,:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
@@ -105,7 +105,7 @@ module psb_s_comm_mod
     end subroutine psb_sscatterm
     subroutine  psb_sscatterv(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      real(psb_spk_), intent(out) :: locx(:)
+      real(psb_spk_), intent(out), allocatable :: locx(:)
       real(psb_spk_), intent(in)  :: globx(:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
diff --git a/base/modules/psb_s_mat_mod.f90 b/base/modules/psb_s_mat_mod.f90
index 1b4994a75..515c3ad32 100644
--- a/base/modules/psb_s_mat_mod.f90
+++ b/base/modules/psb_s_mat_mod.f90
@@ -44,7 +44,7 @@
 ! the functionalities to have the encapsulated class change its
 ! type dynamically, and to extract/input an inner object.
 !
-! A sparse matric has a state corresponding to its progression
+! A sparse matrix has a state corresponding to its progression
 ! through the application life.
 ! In particular, computational methods can only be invoked when
 ! the matrix is in the ASSEMBLED state, whereas the other states are
@@ -133,6 +133,7 @@ module psb_s_mat_mod
     procedure, pass(a) :: m_csclip    => psb_s_csclip
     procedure, pass(a) :: b_csclip    => psb_s_b_csclip
     generic, public    :: csclip      => b_csclip, m_csclip
+    procedure, pass(a) :: clean_zeros => psb_s_clean_zeros
     procedure, pass(a) :: reall       => psb_s_reallocate_nz
     procedure, pass(a) :: get_neigh   => psb_s_get_neigh
     procedure, pass(a) :: reinit      => psb_s_reinit
@@ -1177,5 +1178,15 @@ contains
 
   end function psb_s_get_nz_row
 
+  subroutine psb_s_clean_zeros(a,info)
+    implicit none 
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_sspmat_type), intent(inout) :: a
+
+    info = 0 
+    if (allocated(a%a)) call a%a%clean_zeros(info)
+
+  end subroutine psb_s_clean_zeros
+
 
 end module psb_s_mat_mod
diff --git a/base/modules/psb_z_base_mat_mod.f90 b/base/modules/psb_z_base_mat_mod.f90
index fcc375fef..4e4f22c7c 100644
--- a/base/modules/psb_z_base_mat_mod.f90
+++ b/base/modules/psb_z_base_mat_mod.f90
@@ -42,7 +42,7 @@ module psb_z_base_mat_mod
   !! The psb_z_base_sparse_mat type, extending psb_base_sparse_mat,
   !! defines a middle level  complex(psb_dpk_) sparse matrix object.
   !! This class object itself does not have any additional members
-  !! with respect to those of the base class. No methods can be fully
+  !! with respect to those of the base class. Most methods cannot be fully
   !! implemented at this level, but we can define the interface for the
   !! computational methods requiring the knowledge of the underlying
   !! field, such as the matrix-vector product; this interface is defined,
@@ -78,6 +78,7 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: mold          => psb_z_base_mold 
     procedure, pass(a) :: clone         => psb_z_base_clone
     procedure, pass(a) :: make_nonunit  => psb_z_base_make_nonunit
+    procedure, pass(a) :: clean_zeros   => psb_z_base_clean_zeros
     
     !
     ! Transpose methods: defined here but not implemented. 
@@ -161,6 +162,7 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: get_nz_row   => psb_z_coo_get_nz_row
     procedure, pass(a) :: fix          => psb_z_fix_coo
     procedure, pass(a) :: trim         => psb_z_coo_trim
+    procedure, pass(a) :: clean_zeros  => psb_z_coo_clean_zeros
     procedure, pass(a) :: print        => psb_z_coo_print
     procedure, pass(a) :: free         => z_coo_free
     procedure, pass(a) :: mold         => psb_z_coo_mold
@@ -697,6 +699,18 @@ module psb_z_base_mat_mod
       integer(psb_ipk_), intent(out)            :: info
     end subroutine psb_z_base_mv_from_fmt
   end interface
+  !
+  !> 
+  !! \memberof  psb_z_base_sparse_mat
+  !! \see psb_z_base_mat_mod::psb_z_base_clean_zeros
+  !
+  interface
+    subroutine  psb_z_base_clean_zeros(a, info)
+      import :: psb_ipk_, psb_z_base_sparse_mat
+      class(psb_z_base_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)              :: info
+    end subroutine psb_z_base_clean_zeros
+  end interface
   
   !
   !> Function  transp:
@@ -1235,6 +1249,18 @@ module psb_z_base_mat_mod
       class(psb_z_coo_sparse_mat), intent(inout) :: a
     end subroutine psb_z_coo_trim
   end interface
+  !
+  !> 
+  !! \memberof  psb_z_coo_sparse_mat
+  !! \see psb_z_base_mat_mod::psb_z_base_clean_zeros
+  !
+  interface
+    subroutine  psb_z_coo_clean_zeros(a,info)
+      import :: psb_ipk_, psb_z_coo_sparse_mat
+      class(psb_z_coo_sparse_mat), intent(inout) :: a
+      integer(psb_ipk_), intent(out)             :: info
+    end subroutine psb_z_coo_clean_zeros
+  end interface
   
   !
   !> 
diff --git a/base/modules/psb_z_comm_mod.f90 b/base/modules/psb_z_comm_mod.f90
index e0c8b7972..4a01c9087 100644
--- a/base/modules/psb_z_comm_mod.f90
+++ b/base/modules/psb_z_comm_mod.f90
@@ -97,7 +97,7 @@ module psb_z_comm_mod
   interface psb_scatter
     subroutine  psb_zscatterm(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      complex(psb_dpk_), intent(out) :: locx(:,:)
+      complex(psb_dpk_), intent(out), allocatable :: locx(:,:)
       complex(psb_dpk_), intent(in)  :: globx(:,:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
@@ -105,7 +105,7 @@ module psb_z_comm_mod
     end subroutine psb_zscatterm
     subroutine  psb_zscatterv(globx, locx, desc_a, info, root)
       use psb_desc_mod
-      complex(psb_dpk_), intent(out) :: locx(:)
+      complex(psb_dpk_), intent(out), allocatable :: locx(:)
       complex(psb_dpk_), intent(in)  :: globx(:)
       type(psb_desc_type), intent(in)  :: desc_a
       integer(psb_ipk_), intent(out)             :: info
diff --git a/base/modules/psb_z_mat_mod.f90 b/base/modules/psb_z_mat_mod.f90
index 4b32deb56..ebda6475b 100644
--- a/base/modules/psb_z_mat_mod.f90
+++ b/base/modules/psb_z_mat_mod.f90
@@ -44,7 +44,7 @@
 ! the functionalities to have the encapsulated class change its
 ! type dynamically, and to extract/input an inner object.
 !
-! A sparse matric has a state corresponding to its progression
+! A sparse matrix has a state corresponding to its progression
 ! through the application life.
 ! In particular, computational methods can only be invoked when
 ! the matrix is in the ASSEMBLED state, whereas the other states are
@@ -133,6 +133,7 @@ module psb_z_mat_mod
     procedure, pass(a) :: m_csclip    => psb_z_csclip
     procedure, pass(a) :: b_csclip    => psb_z_b_csclip
     generic, public    :: csclip      => b_csclip, m_csclip
+    procedure, pass(a) :: clean_zeros => psb_z_clean_zeros
     procedure, pass(a) :: reall       => psb_z_reallocate_nz
     procedure, pass(a) :: get_neigh   => psb_z_get_neigh
     procedure, pass(a) :: reinit      => psb_z_reinit
@@ -1177,5 +1178,15 @@ contains
 
   end function psb_z_get_nz_row
 
+  subroutine psb_z_clean_zeros(a,info)
+    implicit none 
+    integer(psb_ipk_), intent(out)        :: info
+    class(psb_zspmat_type), intent(inout) :: a
+
+    info = 0 
+    if (allocated(a%a)) call a%a%clean_zeros(info)
+
+  end subroutine psb_z_clean_zeros
+
 
 end module psb_z_mat_mod
diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90
index 2d3882900..d99ec5fc0 100644
--- a/base/serial/impl/psb_c_base_mat_impl.F90
+++ b/base/serial/impl/psb_c_base_mat_impl.F90
@@ -318,6 +318,23 @@ subroutine psb_c_base_mv_from_fmt(a,b,info)
 
 end subroutine psb_c_base_mv_from_fmt
 
+subroutine  psb_c_base_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_base_clean_zeros
+  implicit none 
+  class(psb_c_base_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_), intent(out)              :: info
+  !
+  type(psb_c_coo_sparse_mat) :: tmpcoo
+
+  call a%mv_to_coo(tmpcoo,info)
+  if (info == 0) call tmpcoo%clean_zeros(info)  
+  if (info == 0) call a%mv_from_coo(tmpcoo,info)
+  
+end subroutine psb_c_base_clean_zeros
+
+
+
 subroutine psb_c_base_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
   use psb_error_mod
   use psb_c_base_mat_mod, psb_protect_name => psb_c_base_csput_a
diff --git a/base/serial/impl/psb_c_coo_impl.f90 b/base/serial/impl/psb_c_coo_impl.f90
index 511faba7a..016f1f220 100644
--- a/base/serial/impl/psb_c_coo_impl.f90
+++ b/base/serial/impl/psb_c_coo_impl.f90
@@ -320,6 +320,31 @@ subroutine  psb_c_coo_trim(a)
 
 end subroutine psb_c_coo_trim
 
+subroutine  psb_c_coo_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_clean_zeros
+  implicit none 
+  class(psb_c_coo_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_) :: info
+  !
+  integer(psb_ipk_) :: i,j,k, nzin
+
+  info = 0 
+  nzin = a%get_nzeros()
+  j = 0
+  do i=1, nzin
+    if (a%val(i) /= czero) then
+      j = j + 1
+      a%val(j) = a%val(i)
+      a%ia(j)  = a%ia(i)
+      a%ja(j)  = a%ja(i)
+    end if
+  end do
+  call a%set_nzeros(j) 
+  call a%trim()
+end subroutine psb_c_coo_clean_zeros
+
+
 
 subroutine  psb_c_coo_allocate_mnnz(m,n,a,nz) 
   use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_allocate_mnnz
diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90
index 804372ccc..e485751c2 100644
--- a/base/serial/impl/psb_d_base_mat_impl.F90
+++ b/base/serial/impl/psb_d_base_mat_impl.F90
@@ -318,6 +318,23 @@ subroutine psb_d_base_mv_from_fmt(a,b,info)
 
 end subroutine psb_d_base_mv_from_fmt
 
+subroutine  psb_d_base_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_base_clean_zeros
+  implicit none 
+  class(psb_d_base_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_), intent(out)              :: info
+  !
+  type(psb_d_coo_sparse_mat) :: tmpcoo
+
+  call a%mv_to_coo(tmpcoo,info)
+  if (info == 0) call tmpcoo%clean_zeros(info)  
+  if (info == 0) call a%mv_from_coo(tmpcoo,info)
+  
+end subroutine psb_d_base_clean_zeros
+
+
+
 subroutine psb_d_base_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
   use psb_error_mod
   use psb_d_base_mat_mod, psb_protect_name => psb_d_base_csput_a
diff --git a/base/serial/impl/psb_d_coo_impl.f90 b/base/serial/impl/psb_d_coo_impl.f90
index e33be37f3..27a68ea77 100644
--- a/base/serial/impl/psb_d_coo_impl.f90
+++ b/base/serial/impl/psb_d_coo_impl.f90
@@ -320,6 +320,31 @@ subroutine  psb_d_coo_trim(a)
 
 end subroutine psb_d_coo_trim
 
+subroutine  psb_d_coo_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_clean_zeros
+  implicit none 
+  class(psb_d_coo_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_) :: info
+  !
+  integer(psb_ipk_) :: i,j,k, nzin
+
+  info = 0 
+  nzin = a%get_nzeros()
+  j = 0
+  do i=1, nzin
+    if (a%val(i) /= dzero) then
+      j = j + 1
+      a%val(j) = a%val(i)
+      a%ia(j)  = a%ia(i)
+      a%ja(j)  = a%ja(i)
+    end if
+  end do
+  call a%set_nzeros(j) 
+  call a%trim()
+end subroutine psb_d_coo_clean_zeros
+
+
 
 subroutine  psb_d_coo_allocate_mnnz(m,n,a,nz) 
   use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_allocate_mnnz
diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90
index 288ae128e..f7ed3384c 100644
--- a/base/serial/impl/psb_s_base_mat_impl.F90
+++ b/base/serial/impl/psb_s_base_mat_impl.F90
@@ -318,6 +318,23 @@ subroutine psb_s_base_mv_from_fmt(a,b,info)
 
 end subroutine psb_s_base_mv_from_fmt
 
+subroutine  psb_s_base_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_base_clean_zeros
+  implicit none 
+  class(psb_s_base_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_), intent(out)              :: info
+  !
+  type(psb_s_coo_sparse_mat) :: tmpcoo
+
+  call a%mv_to_coo(tmpcoo,info)
+  if (info == 0) call tmpcoo%clean_zeros(info)  
+  if (info == 0) call a%mv_from_coo(tmpcoo,info)
+  
+end subroutine psb_s_base_clean_zeros
+
+
+
 subroutine psb_s_base_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
   use psb_error_mod
   use psb_s_base_mat_mod, psb_protect_name => psb_s_base_csput_a
diff --git a/base/serial/impl/psb_s_coo_impl.f90 b/base/serial/impl/psb_s_coo_impl.f90
index 4753fe570..6e36d2945 100644
--- a/base/serial/impl/psb_s_coo_impl.f90
+++ b/base/serial/impl/psb_s_coo_impl.f90
@@ -320,6 +320,31 @@ subroutine  psb_s_coo_trim(a)
 
 end subroutine psb_s_coo_trim
 
+subroutine  psb_s_coo_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_clean_zeros
+  implicit none 
+  class(psb_s_coo_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_) :: info
+  !
+  integer(psb_ipk_) :: i,j,k, nzin
+
+  info = 0 
+  nzin = a%get_nzeros()
+  j = 0
+  do i=1, nzin
+    if (a%val(i) /= szero) then
+      j = j + 1
+      a%val(j) = a%val(i)
+      a%ia(j)  = a%ia(i)
+      a%ja(j)  = a%ja(i)
+    end if
+  end do
+  call a%set_nzeros(j) 
+  call a%trim()
+end subroutine psb_s_coo_clean_zeros
+
+
 
 subroutine  psb_s_coo_allocate_mnnz(m,n,a,nz) 
   use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_allocate_mnnz
diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90
index 42e72e5d0..f56bd024e 100644
--- a/base/serial/impl/psb_z_base_mat_impl.F90
+++ b/base/serial/impl/psb_z_base_mat_impl.F90
@@ -318,6 +318,23 @@ subroutine psb_z_base_mv_from_fmt(a,b,info)
 
 end subroutine psb_z_base_mv_from_fmt
 
+subroutine  psb_z_base_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_base_clean_zeros
+  implicit none 
+  class(psb_z_base_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_), intent(out)              :: info
+  !
+  type(psb_z_coo_sparse_mat) :: tmpcoo
+
+  call a%mv_to_coo(tmpcoo,info)
+  if (info == 0) call tmpcoo%clean_zeros(info)  
+  if (info == 0) call a%mv_from_coo(tmpcoo,info)
+  
+end subroutine psb_z_base_clean_zeros
+
+
+
 subroutine psb_z_base_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
   use psb_error_mod
   use psb_z_base_mat_mod, psb_protect_name => psb_z_base_csput_a
diff --git a/base/serial/impl/psb_z_coo_impl.f90 b/base/serial/impl/psb_z_coo_impl.f90
index d7136a129..79b22b09d 100644
--- a/base/serial/impl/psb_z_coo_impl.f90
+++ b/base/serial/impl/psb_z_coo_impl.f90
@@ -320,6 +320,31 @@ subroutine  psb_z_coo_trim(a)
 
 end subroutine psb_z_coo_trim
 
+subroutine  psb_z_coo_clean_zeros(a, info)
+  use psb_error_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_clean_zeros
+  implicit none 
+  class(psb_z_coo_sparse_mat), intent(inout) :: a
+  integer(psb_ipk_) :: info
+  !
+  integer(psb_ipk_) :: i,j,k, nzin
+
+  info = 0 
+  nzin = a%get_nzeros()
+  j = 0
+  do i=1, nzin
+    if (a%val(i) /= zzero) then
+      j = j + 1
+      a%val(j) = a%val(i)
+      a%ia(j)  = a%ia(i)
+      a%ja(j)  = a%ja(i)
+    end if
+  end do
+  call a%set_nzeros(j) 
+  call a%trim()
+end subroutine psb_z_coo_clean_zeros
+
+
 
 subroutine  psb_z_coo_allocate_mnnz(m,n,a,nz) 
   use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_allocate_mnnz