diff --git a/base/modules/psblas/psb_c_psblas_mod.F90 b/base/modules/psblas/psb_c_psblas_mod.F90
index b5f9d97a..b4f1fe28 100644
--- a/base/modules/psblas/psb_c_psblas_mod.F90
+++ b/base/modules/psblas/psb_c_psblas_mod.F90
@@ -574,6 +574,26 @@ module psb_c_psblas_mod
       type(psb_desc_type), intent (in)        :: desc_a
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_ccmp_vect
+    subroutine psb_ccmp_spmatval(a,val,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_cspmat_type, psb_spk_
+      type(psb_cspmat_type), intent(inout)  :: a
+      complex(psb_spk_), intent(in)             :: val
+      real(psb_spk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_ccmp_spmatval
+    subroutine psb_ccmp_spmat(a,b,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_cspmat_type, psb_spk_
+      type(psb_cspmat_type), intent(inout)  :: a
+      type(psb_cspmat_type), intent(inout)  :: b
+      real(psb_spk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_ccmp_spmat
   end interface
   interface psb_geaddconst
     subroutine psb_caddconst_vect(x,b,z,desc_a,info)
diff --git a/base/modules/psblas/psb_d_psblas_mod.F90 b/base/modules/psblas/psb_d_psblas_mod.F90
index 05e4b343..a7748116 100644
--- a/base/modules/psblas/psb_d_psblas_mod.F90
+++ b/base/modules/psblas/psb_d_psblas_mod.F90
@@ -585,6 +585,26 @@ module psb_d_psblas_mod
       type(psb_desc_type), intent (in)        :: desc_a
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_dcmp_vect
+    subroutine psb_dcmp_spmatval(a,val,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_dspmat_type, psb_dpk_
+      type(psb_dspmat_type), intent(inout)  :: a
+      real(psb_dpk_), intent(in)             :: val
+      real(psb_dpk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_dcmp_spmatval
+    subroutine psb_dcmp_spmat(a,b,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_dspmat_type, psb_dpk_
+      type(psb_dspmat_type), intent(inout)  :: a
+      type(psb_dspmat_type), intent(inout)  :: b
+      real(psb_dpk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_dcmp_spmat
   end interface
   interface psb_geaddconst
     subroutine psb_daddconst_vect(x,b,z,desc_a,info)
diff --git a/base/modules/psblas/psb_s_psblas_mod.F90 b/base/modules/psblas/psb_s_psblas_mod.F90
index 651a0c6e..803c7328 100644
--- a/base/modules/psblas/psb_s_psblas_mod.F90
+++ b/base/modules/psblas/psb_s_psblas_mod.F90
@@ -585,6 +585,26 @@ module psb_s_psblas_mod
       type(psb_desc_type), intent (in)        :: desc_a
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_scmp_vect
+    subroutine psb_scmp_spmatval(a,val,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_sspmat_type, psb_spk_
+      type(psb_sspmat_type), intent(inout)  :: a
+      real(psb_spk_), intent(in)             :: val
+      real(psb_spk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_scmp_spmatval
+    subroutine psb_scmp_spmat(a,b,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_sspmat_type, psb_spk_
+      type(psb_sspmat_type), intent(inout)  :: a
+      type(psb_sspmat_type), intent(inout)  :: b
+      real(psb_spk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_scmp_spmat
   end interface
   interface psb_geaddconst
     subroutine psb_saddconst_vect(x,b,z,desc_a,info)
diff --git a/base/modules/psblas/psb_z_psblas_mod.F90 b/base/modules/psblas/psb_z_psblas_mod.F90
index 63fed5ae..60a373e4 100644
--- a/base/modules/psblas/psb_z_psblas_mod.F90
+++ b/base/modules/psblas/psb_z_psblas_mod.F90
@@ -574,6 +574,26 @@ module psb_z_psblas_mod
       type(psb_desc_type), intent (in)        :: desc_a
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_zcmp_vect
+    subroutine psb_zcmp_spmatval(a,val,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_zspmat_type, psb_dpk_
+      type(psb_zspmat_type), intent(inout)  :: a
+      complex(psb_dpk_), intent(in)             :: val
+      real(psb_dpk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_zcmp_spmatval
+    subroutine psb_zcmp_spmat(a,b,tol,desc_a,res,info)
+      import :: psb_desc_type, psb_ipk_, &
+           & psb_lpk_, psb_zspmat_type, psb_dpk_
+      type(psb_zspmat_type), intent(inout)  :: a
+      type(psb_zspmat_type), intent(inout)  :: b
+      real(psb_dpk_), intent(in)            :: tol
+      type(psb_desc_type), intent (in)        :: desc_a
+      integer(psb_ipk_), intent(out)          :: info
+      logical, intent(out)                    :: res
+    end subroutine psb_zcmp_spmat
   end interface
   interface psb_geaddconst
     subroutine psb_zaddconst_vect(x,b,z,desc_a,info)
diff --git a/base/modules/serial/psb_c_base_mat_mod.F90 b/base/modules/serial/psb_c_base_mat_mod.F90
index c7ba1a56..a52f7a7c 100644
--- a/base/modules/serial/psb_c_base_mat_mod.F90
+++ b/base/modules/serial/psb_c_base_mat_mod.F90
@@ -127,6 +127,9 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_c_base_aclsum
     procedure, pass(a) :: scalpid     => psb_c_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_c_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_c_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_c_base_cmpmat
+    generic, public    :: spcmp       => cmpval, cmpmat
   end type psb_c_base_sparse_mat
 
   private :: c_base_mat_sync, c_base_mat_is_host, c_base_mat_is_dev, &
@@ -230,6 +233,8 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_c_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_c_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_c_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_c_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_c_coo_cmpmat
   end type psb_c_coo_sparse_mat
 
   private :: c_coo_get_nzeros, c_coo_set_nzeros, &
@@ -295,6 +300,9 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_lc_base_aclsum
     procedure, pass(a) :: scalpid     => psb_lc_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_lc_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_lc_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_lc_base_cmpmat
+    generic, public    :: spcmp         => cmpval, cmpmat
     !
     ! Convert internal indices
     !
@@ -397,6 +405,8 @@ module psb_c_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_lc_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_lc_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_lc_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_lc_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_lc_coo_cmpmat
     !
     ! This is COO specific
     !
@@ -1499,6 +1509,52 @@ module psb_c_base_mat_mod
     end subroutine psb_c_base_spaxpby
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_c_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_c_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_c_base_sparse_mat), intent(inout) :: a
+          complex(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_c_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_c_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_c_base_sparse_mat), intent(inout) :: a
+          class(psb_c_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_c_base_sparse_mat
@@ -2169,6 +2225,34 @@ module psb_c_base_mat_mod
     end subroutine psb_c_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_c_coo_sparse_mat
+  !! \see psb_c_base_mat_mod::psb_c_base_cmpval
+  interface
+      function psb_c_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_c_coo_sparse_mat), intent(inout) :: a
+          complex(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_c_coo_sparse_mat
+  !! \see psb_c_base_mat_mod::psb_c_base_cmpmat
+  interface
+      function psb_c_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_c_coo_sparse_mat), intent(inout) :: a
+          class(psb_c_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_coo_cmpmat
+  end interface
+
   ! == =================
   !
   ! BASE interfaces
@@ -2966,6 +3050,52 @@ module psb_c_base_mat_mod
     end subroutine psb_lc_base_scal
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_lc_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_lc_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_lc_base_sparse_mat), intent(inout) :: a
+          complex(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lc_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_lc_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_lc_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_lc_base_sparse_mat), intent(inout) :: a
+          class(psb_lc_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lc_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_lc_base_sparse_mat
@@ -3569,6 +3699,34 @@ module psb_c_base_mat_mod
     end subroutine psb_lc_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_lc_coo_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cmpval
+  interface
+      function psb_lc_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_lc_coo_sparse_mat), intent(inout) :: a
+          complex(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lc_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_lc_coo_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cmpmat
+  interface
+      function psb_lc_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_lc_coo_sparse_mat), intent(inout)  :: a
+          class(psb_lc_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lc_coo_cmpmat
+  end interface
+
 contains
 
 
diff --git a/base/modules/serial/psb_c_mat_mod.F90 b/base/modules/serial/psb_c_mat_mod.F90
index 3beb188a..82d565ec 100644
--- a/base/modules/serial/psb_c_mat_mod.F90
+++ b/base/modules/serial/psb_c_mat_mod.F90
@@ -235,6 +235,9 @@ module psb_c_mat_mod
     generic, public    :: spsm     => cssm, cssv, cssv_v
     procedure, pass(a) :: scalpid  => psb_c_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_c_spaxpby
+    procedure, pass(a) :: cmpval   => psb_c_cmpval
+    procedure, pass(a) :: cmpmat   => psb_c_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_cspmat_type
 
@@ -421,6 +424,9 @@ module psb_c_mat_mod
     generic, public    :: scal     => scals, scalv
     procedure, pass(a) :: scalpid  => psb_lc_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_lc_spaxpby
+    procedure, pass(a) :: cmpval   => psb_lc_cmpval
+    procedure, pass(a) :: cmpmat   => psb_lc_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_lcspmat_type
 
@@ -1181,6 +1187,28 @@ module psb_c_mat_mod
       end subroutine psb_c_spaxpby
   end interface
 
+  interface
+      function psb_c_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_
+          class(psb_cspmat_type), intent(inout) :: a
+          complex(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_cmpval
+  end interface
+
+  interface
+      function psb_c_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_
+          class(psb_cspmat_type), intent(inout) :: a
+          class(psb_cspmat_type), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_c_cmpmat
+  end interface
+
   ! == ===================================
   !
   !
@@ -1856,6 +1884,25 @@ module psb_c_mat_mod
     end function psb_lc_aclsum
   end interface
 
+  interface  psb_cmpmat
+      function psb_lc_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_
+          class(psb_lcspmat_type), intent(inout) :: a
+          complex(psb_spk_), intent(in)              :: val
+          real(psb_spk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_lc_cmpval
+      function psb_lc_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_
+          class(psb_lcspmat_type), intent(inout) :: a
+          class(psb_lcspmat_type), intent(inout) :: b
+          real(psb_spk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_lc_cmpmat
+  end interface
+
 contains
 
   subroutine  psb_c_set_mat_default(a)
diff --git a/base/modules/serial/psb_d_base_mat_mod.F90 b/base/modules/serial/psb_d_base_mat_mod.F90
index 6e3c7a46..a4b78e82 100644
--- a/base/modules/serial/psb_d_base_mat_mod.F90
+++ b/base/modules/serial/psb_d_base_mat_mod.F90
@@ -127,6 +127,9 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_d_base_aclsum
     procedure, pass(a) :: scalpid     => psb_d_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_d_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_d_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_d_base_cmpmat
+    generic, public    :: spcmp       => cmpval, cmpmat
   end type psb_d_base_sparse_mat
 
   private :: d_base_mat_sync, d_base_mat_is_host, d_base_mat_is_dev, &
@@ -230,6 +233,8 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_d_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_d_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_d_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_d_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_d_coo_cmpmat
   end type psb_d_coo_sparse_mat
 
   private :: d_coo_get_nzeros, d_coo_set_nzeros, &
@@ -295,6 +300,9 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_ld_base_aclsum
     procedure, pass(a) :: scalpid     => psb_ld_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_ld_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_ld_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_ld_base_cmpmat
+    generic, public    :: spcmp         => cmpval, cmpmat
     !
     ! Convert internal indices
     !
@@ -397,6 +405,8 @@ module psb_d_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_ld_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_ld_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_ld_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_ld_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_ld_coo_cmpmat
     !
     ! This is COO specific
     !
@@ -1499,6 +1509,52 @@ module psb_d_base_mat_mod
     end subroutine psb_d_base_spaxpby
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_d_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_d_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_d_base_sparse_mat), intent(inout) :: a
+          real(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_d_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_d_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_d_base_sparse_mat), intent(inout) :: a
+          class(psb_d_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_d_base_sparse_mat
@@ -2169,6 +2225,34 @@ module psb_d_base_mat_mod
     end subroutine psb_d_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_d_coo_sparse_mat
+  !! \see psb_d_base_mat_mod::psb_d_base_cmpval
+  interface
+      function psb_d_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_d_coo_sparse_mat), intent(inout) :: a
+          real(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_d_coo_sparse_mat
+  !! \see psb_d_base_mat_mod::psb_d_base_cmpmat
+  interface
+      function psb_d_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_d_coo_sparse_mat), intent(inout) :: a
+          class(psb_d_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_coo_cmpmat
+  end interface
+
   ! == =================
   !
   ! BASE interfaces
@@ -2966,6 +3050,52 @@ module psb_d_base_mat_mod
     end subroutine psb_ld_base_scal
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_ld_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_ld_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_ld_base_sparse_mat), intent(inout) :: a
+          real(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ld_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_ld_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_ld_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_ld_base_sparse_mat), intent(inout) :: a
+          class(psb_ld_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ld_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_ld_base_sparse_mat
@@ -3569,6 +3699,34 @@ module psb_d_base_mat_mod
     end subroutine psb_ld_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_ld_coo_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cmpval
+  interface
+      function psb_ld_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_ld_coo_sparse_mat), intent(inout) :: a
+          real(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ld_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_ld_coo_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cmpmat
+  interface
+      function psb_ld_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_ld_coo_sparse_mat), intent(inout)  :: a
+          class(psb_ld_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ld_coo_cmpmat
+  end interface
+
 contains
 
 
diff --git a/base/modules/serial/psb_d_mat_mod.F90 b/base/modules/serial/psb_d_mat_mod.F90
index c0f1577b..26319239 100644
--- a/base/modules/serial/psb_d_mat_mod.F90
+++ b/base/modules/serial/psb_d_mat_mod.F90
@@ -235,6 +235,9 @@ module psb_d_mat_mod
     generic, public    :: spsm     => cssm, cssv, cssv_v
     procedure, pass(a) :: scalpid  => psb_d_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_d_spaxpby
+    procedure, pass(a) :: cmpval   => psb_d_cmpval
+    procedure, pass(a) :: cmpmat   => psb_d_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_dspmat_type
 
@@ -421,6 +424,9 @@ module psb_d_mat_mod
     generic, public    :: scal     => scals, scalv
     procedure, pass(a) :: scalpid  => psb_ld_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_ld_spaxpby
+    procedure, pass(a) :: cmpval   => psb_ld_cmpval
+    procedure, pass(a) :: cmpmat   => psb_ld_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_ldspmat_type
 
@@ -1181,6 +1187,28 @@ module psb_d_mat_mod
       end subroutine psb_d_spaxpby
   end interface
 
+  interface
+      function psb_d_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_
+          class(psb_dspmat_type), intent(inout) :: a
+          real(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_cmpval
+  end interface
+
+  interface
+      function psb_d_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_
+          class(psb_dspmat_type), intent(inout) :: a
+          class(psb_dspmat_type), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_d_cmpmat
+  end interface
+
   ! == ===================================
   !
   !
@@ -1856,6 +1884,25 @@ module psb_d_mat_mod
     end function psb_ld_aclsum
   end interface
 
+  interface  psb_cmpmat
+      function psb_ld_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_
+          class(psb_ldspmat_type), intent(inout) :: a
+          real(psb_dpk_), intent(in)              :: val
+          real(psb_dpk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_ld_cmpval
+      function psb_ld_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_
+          class(psb_ldspmat_type), intent(inout) :: a
+          class(psb_ldspmat_type), intent(inout) :: b
+          real(psb_dpk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_ld_cmpmat
+  end interface
+
 contains
 
   subroutine  psb_d_set_mat_default(a)
diff --git a/base/modules/serial/psb_s_base_mat_mod.F90 b/base/modules/serial/psb_s_base_mat_mod.F90
index 27f053a4..9069d48e 100644
--- a/base/modules/serial/psb_s_base_mat_mod.F90
+++ b/base/modules/serial/psb_s_base_mat_mod.F90
@@ -127,6 +127,9 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_s_base_aclsum
     procedure, pass(a) :: scalpid     => psb_s_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_s_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_s_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_s_base_cmpmat
+    generic, public    :: spcmp       => cmpval, cmpmat
   end type psb_s_base_sparse_mat
 
   private :: s_base_mat_sync, s_base_mat_is_host, s_base_mat_is_dev, &
@@ -230,6 +233,8 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_s_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_s_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_s_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_s_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_s_coo_cmpmat
   end type psb_s_coo_sparse_mat
 
   private :: s_coo_get_nzeros, s_coo_set_nzeros, &
@@ -295,6 +300,9 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_ls_base_aclsum
     procedure, pass(a) :: scalpid     => psb_ls_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_ls_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_ls_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_ls_base_cmpmat
+    generic, public    :: spcmp         => cmpval, cmpmat
     !
     ! Convert internal indices
     !
@@ -397,6 +405,8 @@ module psb_s_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_ls_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_ls_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_ls_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_ls_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_ls_coo_cmpmat
     !
     ! This is COO specific
     !
@@ -1499,6 +1509,52 @@ module psb_s_base_mat_mod
     end subroutine psb_s_base_spaxpby
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_s_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_s_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_s_base_sparse_mat), intent(inout) :: a
+          real(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_s_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_s_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_s_base_sparse_mat), intent(inout) :: a
+          class(psb_s_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_s_base_sparse_mat
@@ -2169,6 +2225,34 @@ module psb_s_base_mat_mod
     end subroutine psb_s_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_s_coo_sparse_mat
+  !! \see psb_s_base_mat_mod::psb_s_base_cmpval
+  interface
+      function psb_s_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_s_coo_sparse_mat), intent(inout) :: a
+          real(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_s_coo_sparse_mat
+  !! \see psb_s_base_mat_mod::psb_s_base_cmpmat
+  interface
+      function psb_s_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_s_coo_sparse_mat), intent(inout) :: a
+          class(psb_s_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_coo_cmpmat
+  end interface
+
   ! == =================
   !
   ! BASE interfaces
@@ -2966,6 +3050,52 @@ module psb_s_base_mat_mod
     end subroutine psb_ls_base_scal
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_ls_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_ls_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_ls_base_sparse_mat), intent(inout) :: a
+          real(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ls_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_ls_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_ls_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_ls_base_sparse_mat), intent(inout) :: a
+          class(psb_ls_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ls_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_ls_base_sparse_mat
@@ -3569,6 +3699,34 @@ module psb_s_base_mat_mod
     end subroutine psb_ls_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_ls_coo_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cmpval
+  interface
+      function psb_ls_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_ls_coo_sparse_mat), intent(inout) :: a
+          real(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ls_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_ls_coo_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cmpmat
+  interface
+      function psb_ls_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_ls_coo_sparse_mat), intent(inout)  :: a
+          class(psb_ls_base_sparse_mat), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_ls_coo_cmpmat
+  end interface
+
 contains
 
 
diff --git a/base/modules/serial/psb_s_mat_mod.F90 b/base/modules/serial/psb_s_mat_mod.F90
index c3b19516..32a2f378 100644
--- a/base/modules/serial/psb_s_mat_mod.F90
+++ b/base/modules/serial/psb_s_mat_mod.F90
@@ -235,6 +235,9 @@ module psb_s_mat_mod
     generic, public    :: spsm     => cssm, cssv, cssv_v
     procedure, pass(a) :: scalpid  => psb_s_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_s_spaxpby
+    procedure, pass(a) :: cmpval   => psb_s_cmpval
+    procedure, pass(a) :: cmpmat   => psb_s_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_sspmat_type
 
@@ -421,6 +424,9 @@ module psb_s_mat_mod
     generic, public    :: scal     => scals, scalv
     procedure, pass(a) :: scalpid  => psb_ls_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_ls_spaxpby
+    procedure, pass(a) :: cmpval   => psb_ls_cmpval
+    procedure, pass(a) :: cmpmat   => psb_ls_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_lsspmat_type
 
@@ -1181,6 +1187,28 @@ module psb_s_mat_mod
       end subroutine psb_s_spaxpby
   end interface
 
+  interface
+      function psb_s_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_
+          class(psb_sspmat_type), intent(inout) :: a
+          real(psb_spk_), intent(in)             :: val
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_cmpval
+  end interface
+
+  interface
+      function psb_s_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_
+          class(psb_sspmat_type), intent(inout) :: a
+          class(psb_sspmat_type), intent(inout) :: b
+          real(psb_spk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_s_cmpmat
+  end interface
+
   ! == ===================================
   !
   !
@@ -1856,6 +1884,25 @@ module psb_s_mat_mod
     end function psb_ls_aclsum
   end interface
 
+  interface  psb_cmpmat
+      function psb_ls_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_
+          class(psb_lsspmat_type), intent(inout) :: a
+          real(psb_spk_), intent(in)              :: val
+          real(psb_spk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_ls_cmpval
+      function psb_ls_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_
+          class(psb_lsspmat_type), intent(inout) :: a
+          class(psb_lsspmat_type), intent(inout) :: b
+          real(psb_spk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_ls_cmpmat
+  end interface
+
 contains
 
   subroutine  psb_s_set_mat_default(a)
diff --git a/base/modules/serial/psb_z_base_mat_mod.F90 b/base/modules/serial/psb_z_base_mat_mod.F90
index c43373db..cc7b2e22 100644
--- a/base/modules/serial/psb_z_base_mat_mod.F90
+++ b/base/modules/serial/psb_z_base_mat_mod.F90
@@ -127,6 +127,9 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_z_base_aclsum
     procedure, pass(a) :: scalpid     => psb_z_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_z_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_z_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_z_base_cmpmat
+    generic, public    :: spcmp       => cmpval, cmpmat
   end type psb_z_base_sparse_mat
 
   private :: z_base_mat_sync, z_base_mat_is_host, z_base_mat_is_dev, &
@@ -230,6 +233,8 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_z_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_z_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_z_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_z_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_z_coo_cmpmat
   end type psb_z_coo_sparse_mat
 
   private :: z_coo_get_nzeros, z_coo_set_nzeros, &
@@ -295,6 +300,9 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: aclsum      => psb_lz_base_aclsum
     procedure, pass(a) :: scalpid     => psb_lz_base_scalplusidentity
     procedure, pass(a) :: spaxpby     => psb_lz_base_spaxpby
+    procedure, pass(a) :: cmpval      => psb_lz_base_cmpval
+    procedure, pass(a) :: cmpmat      => psb_lz_base_cmpmat
+    generic, public    :: spcmp         => cmpval, cmpmat
     !
     ! Convert internal indices
     !
@@ -397,6 +405,8 @@ module psb_z_base_mat_mod
     procedure, pass(a) :: aclsum     => psb_lz_coo_aclsum
     procedure, pass(a) :: scalpid    => psb_lz_coo_scalplusidentity
     procedure, pass(a) :: spaxpby    => psb_lz_coo_spaxpby
+    procedure, pass(a) :: cmpval     => psb_lz_coo_cmpval
+    procedure, pass(a) :: cmpmat     => psb_lz_coo_cmpmat
     !
     ! This is COO specific
     !
@@ -1499,6 +1509,52 @@ module psb_z_base_mat_mod
     end subroutine psb_z_base_spaxpby
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_z_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_z_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_z_base_sparse_mat), intent(inout) :: a
+          complex(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_z_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_z_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_z_base_sparse_mat), intent(inout) :: a
+          class(psb_z_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_z_base_sparse_mat
@@ -2169,6 +2225,34 @@ module psb_z_base_mat_mod
     end subroutine psb_z_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_z_coo_sparse_mat
+  !! \see psb_z_base_mat_mod::psb_z_base_cmpval
+  interface
+      function psb_z_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_z_coo_sparse_mat), intent(inout) :: a
+          complex(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_z_coo_sparse_mat
+  !! \see psb_z_base_mat_mod::psb_z_base_cmpmat
+  interface
+      function psb_z_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_z_coo_sparse_mat), intent(inout) :: a
+          class(psb_z_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_coo_cmpmat
+  end interface
+
   ! == =================
   !
   ! BASE interfaces
@@ -2966,6 +3050,52 @@ module psb_z_base_mat_mod
     end subroutine psb_lz_base_scal
   end interface
 
+  !
+  !> Function  base_cmpval:
+  !! \memberof  psb_lz_base_sparse_mat
+  !! \brief Compare the element of A with the value val |A(i,j) -val| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param val    comparing element for the entries of A
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_lz_base_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_lz_base_sparse_mat), intent(inout) :: a
+          complex(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lz_base_cmpval
+  end interface
+
+  !
+  !> Function  base_cmpmat:
+  !! \memberof  psb_lz_base_sparse_mat
+  !! \brief Compare the element of A with the ones of B |A(i,j) - B(i,j)| < tol
+  !!
+  !! \param alpha  scaling for A
+  !! \param A      sparse matrix A (intent inout)
+  !! \param A      sparse matrix B (intent inout)
+  !! \param tol    tolerance to which the comparison is done
+  !! \param res    return logical
+  !! \param info   return code
+  !
+  interface
+      function psb_lz_base_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_lz_base_sparse_mat), intent(inout) :: a
+          class(psb_lz_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lz_base_cmpmat
+  end interface
+
   !
   !> Function  base_maxval:
   !! \memberof  psb_lz_base_sparse_mat
@@ -3569,6 +3699,34 @@ module psb_z_base_mat_mod
     end subroutine psb_lz_coo_spaxpby
   end interface
 
+  !
+  !! \memberof  psb_lz_coo_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cmpval
+  interface
+      function psb_lz_coo_cmpval(a,val,tol,info) result(res)
+          import
+          class(psb_lz_coo_sparse_mat), intent(inout) :: a
+          complex(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lz_coo_cmpval
+  end interface
+
+  !
+  !! \memberof  psb_lz_coo_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cmpmat
+  interface
+      function psb_lz_coo_cmpmat(a,b,tol,info) result(res)
+          import
+          class(psb_lz_coo_sparse_mat), intent(inout)  :: a
+          class(psb_lz_base_sparse_mat), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_lz_coo_cmpmat
+  end interface
+
 contains
 
 
diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90
index da7776a3..7a32cde8 100644
--- a/base/modules/serial/psb_z_mat_mod.F90
+++ b/base/modules/serial/psb_z_mat_mod.F90
@@ -235,6 +235,9 @@ module psb_z_mat_mod
     generic, public    :: spsm     => cssm, cssv, cssv_v
     procedure, pass(a) :: scalpid  => psb_z_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_z_spaxpby
+    procedure, pass(a) :: cmpval   => psb_z_cmpval
+    procedure, pass(a) :: cmpmat   => psb_z_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_zspmat_type
 
@@ -421,6 +424,9 @@ module psb_z_mat_mod
     generic, public    :: scal     => scals, scalv
     procedure, pass(a) :: scalpid  => psb_lz_scalplusidentity
     procedure, pass(a) :: spaxpby  => psb_lz_spaxpby
+    procedure, pass(a) :: cmpval   => psb_lz_cmpval
+    procedure, pass(a) :: cmpmat   => psb_lz_cmpmat
+    generic, public    :: spcmp    => cmpval, cmpmat
 
   end type psb_lzspmat_type
 
@@ -1181,6 +1187,28 @@ module psb_z_mat_mod
       end subroutine psb_z_spaxpby
   end interface
 
+  interface
+      function psb_z_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_
+          class(psb_zspmat_type), intent(inout) :: a
+          complex(psb_dpk_), intent(in)             :: val
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_cmpval
+  end interface
+
+  interface
+      function psb_z_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_
+          class(psb_zspmat_type), intent(inout) :: a
+          class(psb_zspmat_type), intent(inout) :: b
+          real(psb_dpk_), intent(in)            :: tol
+          logical                                 :: res
+          integer(psb_ipk_), intent(out)          :: info
+      end function psb_z_cmpmat
+  end interface
+
   ! == ===================================
   !
   !
@@ -1856,6 +1884,25 @@ module psb_z_mat_mod
     end function psb_lz_aclsum
   end interface
 
+  interface  psb_cmpmat
+      function psb_lz_cmpval(a,val,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_
+          class(psb_lzspmat_type), intent(inout) :: a
+          complex(psb_dpk_), intent(in)              :: val
+          real(psb_dpk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_lz_cmpval
+      function psb_lz_cmpmat(a,b,tol,info) result(res)
+          import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_
+          class(psb_lzspmat_type), intent(inout) :: a
+          class(psb_lzspmat_type), intent(inout) :: b
+          real(psb_dpk_), intent(in)             :: tol
+          logical                                  :: res
+          integer(psb_ipk_), intent(out)           :: info
+      end function psb_lz_cmpmat
+  end interface
+
 contains
 
   subroutine  psb_z_set_mat_default(a)
diff --git a/base/psblas/psb_ccmp_vect.f90 b/base/psblas/psb_ccmp_vect.f90
index a0429852..64298a74 100644
--- a/base/psblas/psb_ccmp_vect.f90
+++ b/base/psblas/psb_ccmp_vect.f90
@@ -104,3 +104,140 @@ subroutine psb_ccmp_vect(x,c,z,desc_a,info)
 
 end subroutine psb_ccmp_vect
 
+subroutine psb_ccmp_spmatval(a,val,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_ccmp_spmatval
+  implicit none
+  type(psb_cspmat_type), intent(inout)  :: a
+  complex(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja
+  integer(psb_ipk_) :: iia,jja
+
+  name='psb_ccmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(val,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+end subroutine psb_ccmp_spmatval
+
+subroutine psb_ccmp_spmat(a,b,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_ccmp_spmat
+  implicit none
+  type(psb_cspmat_type), intent(inout)  :: a
+  type(psb_cspmat_type), intent(inout)  :: b
+  real(psb_spk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja,ib,jb
+  integer(psb_ipk_) :: iia,jja,iib,jjb
+
+  name='psb_ccmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  ib = 1
+  jb = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_1'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+  call psb_chkmat(m,n,ib,jb,desc_a,info,iib,jjb)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_2'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(b,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+
+end subroutine psb_ccmp_spmat
diff --git a/base/psblas/psb_dcmp_vect.f90 b/base/psblas/psb_dcmp_vect.f90
index 77147497..98d4251e 100644
--- a/base/psblas/psb_dcmp_vect.f90
+++ b/base/psblas/psb_dcmp_vect.f90
@@ -226,3 +226,140 @@ subroutine psb_dmask_vect(c,x,m,t,desc_a,info)
 end subroutine psb_dmask_vect
 
 
+subroutine psb_dcmp_spmatval(a,val,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_dcmp_spmatval
+  implicit none
+  type(psb_dspmat_type), intent(inout)  :: a
+  real(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja
+  integer(psb_ipk_) :: iia,jja
+
+  name='psb_dcmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(val,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+end subroutine psb_dcmp_spmatval
+
+subroutine psb_dcmp_spmat(a,b,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_dcmp_spmat
+  implicit none
+  type(psb_dspmat_type), intent(inout)  :: a
+  type(psb_dspmat_type), intent(inout)  :: b
+  real(psb_dpk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja,ib,jb
+  integer(psb_ipk_) :: iia,jja,iib,jjb
+
+  name='psb_dcmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  ib = 1
+  jb = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_1'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+  call psb_chkmat(m,n,ib,jb,desc_a,info,iib,jjb)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_2'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(b,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+
+end subroutine psb_dcmp_spmat
diff --git a/base/psblas/psb_scmp_vect.f90 b/base/psblas/psb_scmp_vect.f90
index a7b65b1d..beb2067a 100644
--- a/base/psblas/psb_scmp_vect.f90
+++ b/base/psblas/psb_scmp_vect.f90
@@ -226,3 +226,140 @@ subroutine psb_smask_vect(c,x,m,t,desc_a,info)
 end subroutine psb_smask_vect
 
 
+subroutine psb_scmp_spmatval(a,val,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_scmp_spmatval
+  implicit none
+  type(psb_sspmat_type), intent(inout)  :: a
+  real(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja
+  integer(psb_ipk_) :: iia,jja
+
+  name='psb_scmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(val,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+end subroutine psb_scmp_spmatval
+
+subroutine psb_scmp_spmat(a,b,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_scmp_spmat
+  implicit none
+  type(psb_sspmat_type), intent(inout)  :: a
+  type(psb_sspmat_type), intent(inout)  :: b
+  real(psb_spk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja,ib,jb
+  integer(psb_ipk_) :: iia,jja,iib,jjb
+
+  name='psb_scmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  ib = 1
+  jb = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_1'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+  call psb_chkmat(m,n,ib,jb,desc_a,info,iib,jjb)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_2'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(b,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+
+end subroutine psb_scmp_spmat
diff --git a/base/psblas/psb_zcmp_vect.f90 b/base/psblas/psb_zcmp_vect.f90
index 50d9bce5..b607e0b8 100644
--- a/base/psblas/psb_zcmp_vect.f90
+++ b/base/psblas/psb_zcmp_vect.f90
@@ -104,3 +104,140 @@ subroutine psb_zcmp_vect(x,c,z,desc_a,info)
 
 end subroutine psb_zcmp_vect
 
+subroutine psb_zcmp_spmatval(a,val,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_zcmp_spmatval
+  implicit none
+  type(psb_zspmat_type), intent(inout)  :: a
+  complex(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja
+  integer(psb_ipk_) :: iia,jja
+
+  name='psb_zcmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(val,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+end subroutine psb_zcmp_spmatval
+
+subroutine psb_zcmp_spmat(a,b,tol,desc_a,res,info)
+  use psb_base_mod, psb_protect_name => psb_zcmp_spmat
+  implicit none
+  type(psb_zspmat_type), intent(inout)  :: a
+  type(psb_zspmat_type), intent(inout)  :: b
+  real(psb_dpk_), intent(in)            :: tol
+  type(psb_desc_type), intent (in)        :: desc_a
+  integer(psb_ipk_), intent(out)          :: info
+  logical, intent(out)                    :: res
+
+  ! Local
+  integer(psb_ipk_) :: ictxt, np, me
+  integer(psb_ipk_) :: err_act
+  character(len=20) :: name, ch_err
+  integer(psb_ipk_) :: debug_level, debug_unit
+  integer(psb_lpk_) :: m,n,ia,ja,ib,jb
+  integer(psb_ipk_) :: iia,jja,iib,jjb
+
+  name='psb_zcmp_spmatval'
+  info=psb_success_
+  call psb_erractionsave(err_act)
+  if  (psb_errstatus_fatal()) then
+    info = psb_err_internal_error_ ;    goto 9999
+  end if
+  debug_unit  = psb_get_debug_unit()
+  debug_level = psb_get_debug_level()
+
+  ictxt=desc_a%get_context()
+  call psb_info(ictxt, me, np)
+  if (np == -1) then
+    info = psb_err_context_error_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  ia = 1
+  ja = 1
+  ib = 1
+  jb = 1
+  m    = desc_a%get_global_rows()
+  n    = desc_a%get_global_cols()
+  ! checking for matrix correctness
+  call psb_chkmat(m,n,ia,ja,desc_a,info,iia,jja)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_1'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+  call psb_chkmat(m,n,ib,jb,desc_a,info,iib,jjb)
+  if(info /= psb_success_) then
+    info=psb_err_from_subroutine_
+    ch_err='psb_chkmat_2'
+    call psb_errpush(info,name,a_err=ch_err)
+    goto 9999
+  end if
+
+  res = a%spcmp(b,tol,info)
+
+  call psb_lallreduceand(ictxt,res)
+
+  call psb_erractionrestore(err_act)
+  if (debug_level >= psb_debug_comp_) then
+    call psb_barrier(ictxt)
+    write(debug_unit,*) me,' ',trim(name),' Returning '
+  endif
+  return
+
+9999 call psb_error_handler(ictxt,err_act)
+
+  return
+
+end subroutine psb_zcmp_spmat
diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90
index b9d37636..5ac480e6 100644
--- a/base/serial/impl/psb_c_base_mat_impl.F90
+++ b/base/serial/impl/psb_c_base_mat_impl.F90
@@ -1865,6 +1865,92 @@ subroutine psb_c_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_c_base_spaxpby
 
+function psb_c_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_base_cmpval
+
+  class(psb_c_base_sparse_mat), intent(inout) :: a
+  complex(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_c_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_c_base_cmpval
+
+function psb_c_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_base_cmpmat
+
+  class(psb_c_base_sparse_mat), intent(inout) :: a
+  class(psb_c_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_c_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_c_base_cmpmat
+
 ! == ==================================
 !
 !
@@ -3819,6 +3905,92 @@ subroutine psb_lc_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_lc_base_spaxpby
 
+function psb_lc_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_lc_base_cmpval
+
+  class(psb_lc_base_sparse_mat), intent(inout) :: a
+  complex(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_lc_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lc_base_cmpval
+
+function psb_lc_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_lc_base_cmpmat
+
+  class(psb_lc_base_sparse_mat), intent(inout) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_lc_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lc_base_cmpmat
+
 subroutine psb_lc_base_get_diag(a,d,info)
   use psb_error_mod
   use psb_const_mod
diff --git a/base/serial/impl/psb_c_coo_impl.f90 b/base/serial/impl/psb_c_coo_impl.f90
index 733f576e..644c3638 100644
--- a/base/serial/impl/psb_c_coo_impl.f90
+++ b/base/serial/impl/psb_c_coo_impl.f90
@@ -277,6 +277,106 @@ subroutine psb_c_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_c_coo_spaxpby
 
+function psb_c_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_cmpval
+
+  class(psb_c_coo_sparse_mat),  intent(inout) :: a
+  complex(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_ipk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_c_coo_cmpval
+
+function psb_c_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_cmpmat
+
+  class(psb_c_coo_sparse_mat),  intent(inout) :: a
+  class(psb_c_base_sparse_mat),  intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_ipk_)              :: nza, nzb, M, N
+  type(psb_c_coo_sparse_mat) :: tcoo, bcoo
+  real(psb_spk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_spk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = maxval(abs(tcoo%val));
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_c_coo_cmpmat
+
 subroutine  psb_c_coo_reallocate_nz(nz,a)
   use psb_c_base_mat_mod, psb_protect_name => psb_c_coo_reallocate_nz
   use psb_error_mod
@@ -4681,6 +4781,106 @@ subroutine psb_lc_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_lc_coo_spaxpby
 
+function psb_lc_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_lc_coo_cmpval
+
+  class(psb_lc_coo_sparse_mat),  intent(inout) :: a
+  complex(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_lpk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lc_coo_cmpval
+
+function psb_lc_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_base_mat_mod, psb_protect_name => psb_lc_coo_cmpmat
+
+  class(psb_lc_coo_sparse_mat),  intent(inout) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_lpk_)              :: nza, nzb, M, N
+  type(psb_lc_coo_sparse_mat)  :: tcoo, bcoo
+  real(psb_spk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_spk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = tcoo%spnmi()
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lc_coo_cmpmat
+
 subroutine  psb_lc_coo_reallocate_nz(nz,a)
   use psb_c_base_mat_mod, psb_protect_name => psb_lc_coo_reallocate_nz
   use psb_error_mod
diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90
index 22adbd3a..5663bc29 100644
--- a/base/serial/impl/psb_c_mat_impl.F90
+++ b/base/serial/impl/psb_c_mat_impl.F90
@@ -2521,6 +2521,78 @@ subroutine psb_c_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_c_spaxpby
 
+function psb_c_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cmpval
+  implicit none
+  class(psb_cspmat_type), intent(inout) :: a
+  complex(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_c_cmpval
+
+function psb_c_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cmpmat
+  implicit none
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_cspmat_type), intent(inout) :: b
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_c_cmpmat
+
 subroutine psb_c_mv_from_lb(a,b)
   use psb_error_mod
   use psb_const_mod
@@ -4720,6 +4792,78 @@ subroutine psb_lc_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_lc_spaxpby
 
+function psb_lc_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cmpval
+  implicit none
+  class(psb_lcspmat_type), intent(inout) :: a
+  complex(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_lc_cmpval
+
+function psb_lc_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cmpmat
+  implicit none
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_lcspmat_type), intent(inout) :: b
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_lc_cmpmat
+
 function psb_lc_maxval(a) result(res)
   use psb_c_mat_mod, psb_protect_name => psb_lc_maxval
   use psb_error_mod
diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90
index c432808a..65bc1bff 100644
--- a/base/serial/impl/psb_d_base_mat_impl.F90
+++ b/base/serial/impl/psb_d_base_mat_impl.F90
@@ -1865,6 +1865,92 @@ subroutine psb_d_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_d_base_spaxpby
 
+function psb_d_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_base_cmpval
+
+  class(psb_d_base_sparse_mat), intent(inout) :: a
+  real(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_d_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_d_base_cmpval
+
+function psb_d_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_base_cmpmat
+
+  class(psb_d_base_sparse_mat), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_d_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_d_base_cmpmat
+
 ! == ==================================
 !
 !
@@ -3819,6 +3905,92 @@ subroutine psb_ld_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_ld_base_spaxpby
 
+function psb_ld_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_ld_base_cmpval
+
+  class(psb_ld_base_sparse_mat), intent(inout) :: a
+  real(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_ld_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ld_base_cmpval
+
+function psb_ld_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_ld_base_cmpmat
+
+  class(psb_ld_base_sparse_mat), intent(inout) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_ld_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ld_base_cmpmat
+
 subroutine psb_ld_base_get_diag(a,d,info)
   use psb_error_mod
   use psb_const_mod
diff --git a/base/serial/impl/psb_d_coo_impl.f90 b/base/serial/impl/psb_d_coo_impl.f90
index 29ff8d89..7255cb7a 100644
--- a/base/serial/impl/psb_d_coo_impl.f90
+++ b/base/serial/impl/psb_d_coo_impl.f90
@@ -277,6 +277,106 @@ subroutine psb_d_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_d_coo_spaxpby
 
+function psb_d_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_cmpval
+
+  class(psb_d_coo_sparse_mat),  intent(inout) :: a
+  real(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_ipk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_d_coo_cmpval
+
+function psb_d_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_cmpmat
+
+  class(psb_d_coo_sparse_mat),  intent(inout) :: a
+  class(psb_d_base_sparse_mat),  intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_ipk_)              :: nza, nzb, M, N
+  type(psb_d_coo_sparse_mat) :: tcoo, bcoo
+  real(psb_dpk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_dpk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = maxval(abs(tcoo%val));
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_d_coo_cmpmat
+
 subroutine  psb_d_coo_reallocate_nz(nz,a)
   use psb_d_base_mat_mod, psb_protect_name => psb_d_coo_reallocate_nz
   use psb_error_mod
@@ -4681,6 +4781,106 @@ subroutine psb_ld_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_ld_coo_spaxpby
 
+function psb_ld_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_cmpval
+
+  class(psb_ld_coo_sparse_mat),  intent(inout) :: a
+  real(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_lpk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ld_coo_cmpval
+
+function psb_ld_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_cmpmat
+
+  class(psb_ld_coo_sparse_mat),  intent(inout) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_lpk_)              :: nza, nzb, M, N
+  type(psb_ld_coo_sparse_mat)  :: tcoo, bcoo
+  real(psb_dpk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_dpk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = tcoo%spnmi()
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ld_coo_cmpmat
+
 subroutine  psb_ld_coo_reallocate_nz(nz,a)
   use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_reallocate_nz
   use psb_error_mod
diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90
index 1ef8305d..c8288a71 100644
--- a/base/serial/impl/psb_d_mat_impl.F90
+++ b/base/serial/impl/psb_d_mat_impl.F90
@@ -2521,6 +2521,78 @@ subroutine psb_d_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_d_spaxpby
 
+function psb_d_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cmpval
+  implicit none
+  class(psb_dspmat_type), intent(inout) :: a
+  real(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_d_cmpval
+
+function psb_d_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cmpmat
+  implicit none
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_dspmat_type), intent(inout) :: b
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_d_cmpmat
+
 subroutine psb_d_mv_from_lb(a,b)
   use psb_error_mod
   use psb_const_mod
@@ -4720,6 +4792,78 @@ subroutine psb_ld_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_ld_spaxpby
 
+function psb_ld_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cmpval
+  implicit none
+  class(psb_ldspmat_type), intent(inout) :: a
+  real(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_ld_cmpval
+
+function psb_ld_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cmpmat
+  implicit none
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_ldspmat_type), intent(inout) :: b
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_ld_cmpmat
+
 function psb_ld_maxval(a) result(res)
   use psb_d_mat_mod, psb_protect_name => psb_ld_maxval
   use psb_error_mod
diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90
index cffc4d8a..1e791b36 100644
--- a/base/serial/impl/psb_s_base_mat_impl.F90
+++ b/base/serial/impl/psb_s_base_mat_impl.F90
@@ -1865,6 +1865,92 @@ subroutine psb_s_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_s_base_spaxpby
 
+function psb_s_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_base_cmpval
+
+  class(psb_s_base_sparse_mat), intent(inout) :: a
+  real(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_s_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_s_base_cmpval
+
+function psb_s_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_base_cmpmat
+
+  class(psb_s_base_sparse_mat), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_s_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_s_base_cmpmat
+
 ! == ==================================
 !
 !
@@ -3819,6 +3905,92 @@ subroutine psb_ls_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_ls_base_spaxpby
 
+function psb_ls_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_ls_base_cmpval
+
+  class(psb_ls_base_sparse_mat), intent(inout) :: a
+  real(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_ls_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ls_base_cmpval
+
+function psb_ls_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_ls_base_cmpmat
+
+  class(psb_ls_base_sparse_mat), intent(inout) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_ls_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ls_base_cmpmat
+
 subroutine psb_ls_base_get_diag(a,d,info)
   use psb_error_mod
   use psb_const_mod
diff --git a/base/serial/impl/psb_s_coo_impl.f90 b/base/serial/impl/psb_s_coo_impl.f90
index 05cafaa8..db8d908c 100644
--- a/base/serial/impl/psb_s_coo_impl.f90
+++ b/base/serial/impl/psb_s_coo_impl.f90
@@ -277,6 +277,106 @@ subroutine psb_s_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_s_coo_spaxpby
 
+function psb_s_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_cmpval
+
+  class(psb_s_coo_sparse_mat),  intent(inout) :: a
+  real(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_ipk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_s_coo_cmpval
+
+function psb_s_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_cmpmat
+
+  class(psb_s_coo_sparse_mat),  intent(inout) :: a
+  class(psb_s_base_sparse_mat),  intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_ipk_)              :: nza, nzb, M, N
+  type(psb_s_coo_sparse_mat) :: tcoo, bcoo
+  real(psb_spk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_spk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = maxval(abs(tcoo%val));
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_s_coo_cmpmat
+
 subroutine  psb_s_coo_reallocate_nz(nz,a)
   use psb_s_base_mat_mod, psb_protect_name => psb_s_coo_reallocate_nz
   use psb_error_mod
@@ -4681,6 +4781,106 @@ subroutine psb_ls_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_ls_coo_spaxpby
 
+function psb_ls_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_ls_coo_cmpval
+
+  class(psb_ls_coo_sparse_mat),  intent(inout) :: a
+  real(psb_spk_), intent(in)                   :: val
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_lpk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ls_coo_cmpval
+
+function psb_ls_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_base_mat_mod, psb_protect_name => psb_ls_coo_cmpmat
+
+  class(psb_ls_coo_sparse_mat),  intent(inout) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  real(psb_spk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_lpk_)              :: nza, nzb, M, N
+  type(psb_ls_coo_sparse_mat)  :: tcoo, bcoo
+  real(psb_spk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_spk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = tcoo%spnmi()
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_ls_coo_cmpmat
+
 subroutine  psb_ls_coo_reallocate_nz(nz,a)
   use psb_s_base_mat_mod, psb_protect_name => psb_ls_coo_reallocate_nz
   use psb_error_mod
diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90
index 9b3d1343..65deb5a7 100644
--- a/base/serial/impl/psb_s_mat_impl.F90
+++ b/base/serial/impl/psb_s_mat_impl.F90
@@ -2521,6 +2521,78 @@ subroutine psb_s_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_s_spaxpby
 
+function psb_s_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cmpval
+  implicit none
+  class(psb_sspmat_type), intent(inout) :: a
+  real(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_s_cmpval
+
+function psb_s_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cmpmat
+  implicit none
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_sspmat_type), intent(inout) :: b
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_s_cmpmat
+
 subroutine psb_s_mv_from_lb(a,b)
   use psb_error_mod
   use psb_const_mod
@@ -4720,6 +4792,78 @@ subroutine psb_ls_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_ls_spaxpby
 
+function psb_ls_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cmpval
+  implicit none
+  class(psb_lsspmat_type), intent(inout) :: a
+  real(psb_spk_), intent(in)             :: val
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_ls_cmpval
+
+function psb_ls_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cmpmat
+  implicit none
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_lsspmat_type), intent(inout) :: b
+  real(psb_spk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_ls_cmpmat
+
 function psb_ls_maxval(a) result(res)
   use psb_s_mat_mod, psb_protect_name => psb_ls_maxval
   use psb_error_mod
diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90
index 6490657d..d54000ff 100644
--- a/base/serial/impl/psb_z_base_mat_impl.F90
+++ b/base/serial/impl/psb_z_base_mat_impl.F90
@@ -1865,6 +1865,92 @@ subroutine psb_z_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_z_base_spaxpby
 
+function psb_z_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_base_cmpval
+
+  class(psb_z_base_sparse_mat), intent(inout) :: a
+  complex(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_z_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_z_base_cmpval
+
+function psb_z_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_base_cmpmat
+
+  class(psb_z_base_sparse_mat), intent(inout) :: a
+  class(psb_z_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_z_coo_sparse_mat) :: acoo
+
+  call a%cp_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_z_base_cmpmat
+
 ! == ==================================
 !
 !
@@ -3819,6 +3905,92 @@ subroutine psb_lz_base_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_lz_base_spaxpby
 
+function psb_lz_base_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_lz_base_cmpval
+
+  class(psb_lz_base_sparse_mat), intent(inout) :: a
+  complex(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  type(psb_lz_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(val,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpval')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lz_base_cmpval
+
+function psb_lz_base_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_lz_base_cmpmat
+
+  class(psb_lz_base_sparse_mat), intent(inout) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+  type(psb_lz_coo_sparse_mat) :: acoo
+
+  call a%mv_to_coo(acoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='mv_to_coo')
+    goto 9999
+  end if
+
+  ! Fix the indexes
+  call acoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  res = acoo%spcmp(b,tol,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cmpmat')
+    goto 9999
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lz_base_cmpmat
+
 subroutine psb_lz_base_get_diag(a,d,info)
   use psb_error_mod
   use psb_const_mod
diff --git a/base/serial/impl/psb_z_coo_impl.f90 b/base/serial/impl/psb_z_coo_impl.f90
index 5384d507..73667152 100644
--- a/base/serial/impl/psb_z_coo_impl.f90
+++ b/base/serial/impl/psb_z_coo_impl.f90
@@ -277,6 +277,106 @@ subroutine psb_z_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_z_coo_spaxpby
 
+function psb_z_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_cmpval
+
+  class(psb_z_coo_sparse_mat),  intent(inout) :: a
+  complex(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_ipk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_z_coo_cmpval
+
+function psb_z_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_cmpmat
+
+  class(psb_z_coo_sparse_mat),  intent(inout) :: a
+  class(psb_z_base_sparse_mat),  intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_ipk_)              :: nza, nzb, M, N
+  type(psb_z_coo_sparse_mat) :: tcoo, bcoo
+  real(psb_dpk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_dpk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = maxval(abs(tcoo%val));
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_z_coo_cmpmat
+
 subroutine  psb_z_coo_reallocate_nz(nz,a)
   use psb_z_base_mat_mod, psb_protect_name => psb_z_coo_reallocate_nz
   use psb_error_mod
@@ -4681,6 +4781,106 @@ subroutine psb_lz_coo_spaxpby(alpha,a,beta,b,info)
   return
 end subroutine psb_lz_coo_spaxpby
 
+function psb_lz_coo_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_lz_coo_cmpval
+
+  class(psb_lz_coo_sparse_mat),  intent(inout) :: a
+  complex(psb_dpk_), intent(in)                   :: val
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpval'
+  logical, parameter           :: debug=.false.
+  integer(psb_lpk_)            :: nza
+
+  nza = a%get_nzeros()
+
+  if (any(abs(a%val(1:nza)-val) > tol)) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  info = psb_success_
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lz_coo_cmpval
+
+function psb_lz_coo_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_base_mat_mod, psb_protect_name => psb_lz_coo_cmpmat
+
+  class(psb_lz_coo_sparse_mat),  intent(inout) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  real(psb_dpk_), intent(in)                  :: tol
+  integer(psb_ipk_), intent(out)                :: info
+  logical                                       :: res
+
+  ! Auxiliary
+  integer(psb_ipk_)            :: err_act
+  character(len=20)            :: name='cmpmat'
+  logical, parameter           :: debug=.false.
+
+  integer(psb_lpk_)              :: nza, nzb, M, N
+  type(psb_lz_coo_sparse_mat)  :: tcoo, bcoo
+  real(psb_dpk_)               :: normval
+
+  ! Copy (whatever) b format to coo
+  call b%cp_to_coo(bcoo,info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='cp_to_coo')
+    goto 9999
+  end if
+  ! Get information on the matrix
+  M = a%get_nrows()
+  N = a%get_ncols()
+  nza = a%get_nzeros()
+  nzb = b%get_nzeros()
+  ! Allocate (temporary) space for the solution
+  call tcoo%allocate(M,N,(nza+nzb))
+  ! Compute the sum
+  tcoo%ia(1:nza) = a%ia(1:nza)
+  tcoo%ja(1:nza) = a%ja(1:nza)
+  tcoo%val(1:nza) = a%val(1:nza)
+  tcoo%ia(nza+1:nza+nzb) = bcoo%ia(1:nzb)
+  tcoo%ja(nza+1:nza+nzb) = bcoo%ja(1:nzb)
+  tcoo%val(nza+1:nza+nzb) = (-1_psb_dpk_)*bcoo%val(1:nzb)
+  ! Fix the indexes
+  call tcoo%fix(info)
+  if (info /= psb_success_) then
+    info = psb_err_from_subroutine_
+    call psb_errpush(info,name, a_err='fix')
+    goto 9999
+  end if
+
+  normval = tcoo%spnmi()
+
+  if ( normval > tol) then
+    res = .false.
+  else
+    res = .true.
+  end if
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+end function psb_lz_coo_cmpmat
+
 subroutine  psb_lz_coo_reallocate_nz(nz,a)
   use psb_z_base_mat_mod, psb_protect_name => psb_lz_coo_reallocate_nz
   use psb_error_mod
diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90
index 424f3700..ed16fdec 100644
--- a/base/serial/impl/psb_z_mat_impl.F90
+++ b/base/serial/impl/psb_z_mat_impl.F90
@@ -2521,6 +2521,78 @@ subroutine psb_z_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_z_spaxpby
 
+function psb_z_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cmpval
+  implicit none
+  class(psb_zspmat_type), intent(inout) :: a
+  complex(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_z_cmpval
+
+function psb_z_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cmpmat
+  implicit none
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_zspmat_type), intent(inout) :: b
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_z_cmpmat
+
 subroutine psb_z_mv_from_lb(a,b)
   use psb_error_mod
   use psb_const_mod
@@ -4720,6 +4792,78 @@ subroutine psb_lz_spaxpby(alpha,a,beta,b,info)
 
 end subroutine psb_lz_spaxpby
 
+function psb_lz_cmpval(a,val,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cmpval
+  implicit none
+  class(psb_lzspmat_type), intent(inout) :: a
+  complex(psb_dpk_), intent(in)             :: val
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpval'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(val,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_lz_cmpval
+
+function psb_lz_cmpmat(a,b,tol,info) result(res)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cmpmat
+  implicit none
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_lzspmat_type), intent(inout) :: b
+  real(psb_dpk_), intent(in)            :: tol
+  logical                                 :: res
+  integer(psb_ipk_), intent(out)          :: info
+
+  integer(psb_ipk_) :: err_act
+  character(len=20)  :: name='cmpmat'
+  logical, parameter :: debug=.false.
+
+  res = .false.
+  info = psb_success_
+  call psb_erractionsave(err_act)
+  if (.not.allocated(a%a)) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  res = a%a%spcmp(b%a,tol,info)
+  if (info /= psb_success_) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end function psb_lz_cmpmat
+
 function psb_lz_maxval(a) result(res)
   use psb_z_mat_mod, psb_protect_name => psb_lz_maxval
   use psb_error_mod
diff --git a/cbind/base/psb_c_psblas_cbind_mod.f90 b/cbind/base/psb_c_psblas_cbind_mod.f90
index dbf8a88a..703f86a2 100644
--- a/cbind/base/psb_c_psblas_cbind_mod.f90
+++ b/cbind/base/psb_c_psblas_cbind_mod.f90
@@ -463,6 +463,66 @@ contains
 
   end function psb_c_cgecmp
 
+  function psb_c_cgecmpmat(ah,bh,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_cspmat)  :: ah,bh
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_cspmat_type), pointer :: ap,bp
+    integer(psb_c_ipk_)          :: info
+    real(c_float_complex), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+    if (c_associated(bh%item)) then
+      call c_f_pointer(bh%item,bp)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,bp,tol,descp,res,info)
+
+  end function psb_c_cgecmpmat
+
+  function psb_c_cgecmpmat_val(ah,val,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_cspmat)  :: ah
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_cspmat_type), pointer :: ap
+    integer(psb_c_ipk_)          :: info
+    complex(c_float_complex), value :: val
+    real(c_float_complex), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,val,tol,descp,res,info)
+
+  end function psb_c_cgecmpmat_val
+
   function psb_c_cgeaddconst(xh,bh,zh,cdh) bind(c) result(res)
     implicit none
     integer(psb_c_ipk_)    :: res
diff --git a/cbind/base/psb_d_psblas_cbind_mod.f90 b/cbind/base/psb_d_psblas_cbind_mod.f90
index aef9636f..198c6f79 100644
--- a/cbind/base/psb_d_psblas_cbind_mod.f90
+++ b/cbind/base/psb_d_psblas_cbind_mod.f90
@@ -463,6 +463,66 @@ contains
 
   end function psb_c_dgecmp
 
+  function psb_c_dgecmpmat(ah,bh,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_dspmat)  :: ah,bh
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_dspmat_type), pointer :: ap,bp
+    integer(psb_c_ipk_)          :: info
+    real(c_double), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+    if (c_associated(bh%item)) then
+      call c_f_pointer(bh%item,bp)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,bp,tol,descp,res,info)
+
+  end function psb_c_dgecmpmat
+
+  function psb_c_dgecmpmat_val(ah,val,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_dspmat)  :: ah
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_dspmat_type), pointer :: ap
+    integer(psb_c_ipk_)          :: info
+    real(c_double), value :: val
+    real(c_double), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,val,tol,descp,res,info)
+
+  end function psb_c_dgecmpmat_val
+
   function psb_c_dgeaddconst(xh,bh,zh,cdh) bind(c) result(res)
     implicit none
     integer(psb_c_ipk_)    :: res
diff --git a/cbind/base/psb_s_psblas_cbind_mod.f90 b/cbind/base/psb_s_psblas_cbind_mod.f90
index ab773819..96a0f304 100644
--- a/cbind/base/psb_s_psblas_cbind_mod.f90
+++ b/cbind/base/psb_s_psblas_cbind_mod.f90
@@ -463,6 +463,66 @@ contains
 
   end function psb_c_sgecmp
 
+  function psb_c_sgecmpmat(ah,bh,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_sspmat)  :: ah,bh
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_sspmat_type), pointer :: ap,bp
+    integer(psb_c_ipk_)          :: info
+    real(c_float), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+    if (c_associated(bh%item)) then
+      call c_f_pointer(bh%item,bp)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,bp,tol,descp,res,info)
+
+  end function psb_c_sgecmpmat
+
+  function psb_c_sgecmpmat_val(ah,val,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_sspmat)  :: ah
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_sspmat_type), pointer :: ap
+    integer(psb_c_ipk_)          :: info
+    real(c_float), value :: val
+    real(c_float), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,val,tol,descp,res,info)
+
+  end function psb_c_sgecmpmat_val
+
   function psb_c_sgeaddconst(xh,bh,zh,cdh) bind(c) result(res)
     implicit none
     integer(psb_c_ipk_)    :: res
diff --git a/cbind/base/psb_z_psblas_cbind_mod.f90 b/cbind/base/psb_z_psblas_cbind_mod.f90
index b7bf0038..fbd5797a 100644
--- a/cbind/base/psb_z_psblas_cbind_mod.f90
+++ b/cbind/base/psb_z_psblas_cbind_mod.f90
@@ -463,6 +463,66 @@ contains
 
   end function psb_c_zgecmp
 
+  function psb_c_zgecmpmat(ah,bh,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_zspmat)  :: ah,bh
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_zspmat_type), pointer :: ap,bp
+    integer(psb_c_ipk_)          :: info
+    real(c_double_complex), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+    if (c_associated(bh%item)) then
+      call c_f_pointer(bh%item,bp)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,bp,tol,descp,res,info)
+
+  end function psb_c_zgecmpmat
+
+  function psb_c_zgecmpmat_val(ah,val,tol,cdh) bind(c) result(res)
+    implicit none
+    logical               :: res
+
+    type(psb_c_zspmat)  :: ah
+    type(psb_c_descriptor) :: cdh
+
+    type(psb_desc_type), pointer :: descp
+    type(psb_zspmat_type), pointer :: ap
+    integer(psb_c_ipk_)          :: info
+    complex(c_double_complex), value :: val
+    real(c_double_complex), value :: tol
+
+    if (c_associated(cdh%item)) then
+      call c_f_pointer(cdh%item,descp)
+    else
+      return
+    end if
+    if (c_associated(ah%item)) then
+      call c_f_pointer(ah%item,ap)
+    else
+      return
+    end if
+
+    call psb_gecmp(ap,val,tol,descp,res,info)
+
+  end function psb_c_zgecmpmat_val
+
   function psb_c_zgeaddconst(xh,bh,zh,cdh) bind(c) result(res)
     implicit none
     integer(psb_c_ipk_)    :: res