diff --git a/base/modules/psb_d_mat_mod.f03 b/base/modules/psb_d_mat_mod.f03
index ae57ad89..521f8a79 100644
--- a/base/modules/psb_d_mat_mod.f03
+++ b/base/modules/psb_d_mat_mod.f03
@@ -56,6 +56,9 @@ module psb_d_mat_mod
     procedure, pass(a) :: d_csclip
     procedure, pass(a) :: d_b_csclip
     generic, public    :: csclip => d_b_csclip, d_csclip
+    procedure, pass(a) :: d_clip_d_ip
+    procedure, pass(a) :: d_clip_d
+    generic, public    :: clip_diag => d_clip_d_ip, d_clip_d
     procedure, pass(a) :: reall => reallocate_nz
     procedure, pass(a) :: get_neigh
     procedure, pass(a) :: d_cscnv
@@ -99,7 +102,7 @@ module psb_d_mat_mod
   private :: get_nrows, get_ncols, get_nzeros, get_size, &
        & get_state, get_dupl, is_null, is_bld, is_upd, &
        & is_asb, is_sorted, is_upper, is_lower, is_triangle, &
-       & is_unit, get_neigh, csall, csput, d_csgetrow,&
+       & is_unit, get_neigh, csall, csput, d_csgetrow,  d_clip_d_ip, d_clip_d,&
        & d_csgetblk, d_csclip, d_b_csclip, d_cscnv, d_cscnv_ip, &
        & reallocate_nz, free, trim, &
        & sparse_print, reinit, &
@@ -1599,6 +1602,125 @@ contains
   end subroutine d_cscnv_base
 
 
+  subroutine d_clip_d(a,b,info)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    use psb_d_base_mat_mod
+    implicit none
+    
+    class(psb_d_sparse_mat), intent(in) :: a
+    class(psb_d_sparse_mat), intent(out) :: b
+    integer,intent(out)                  :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='clip_diag'
+    logical, parameter :: debug=.false.
+    type(psb_d_coo_sparse_mat), allocatable  :: acoo
+    integer :: i, j, nz
+    
+    info = 0
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    allocate(acoo,stat=info)    
+    if (info == 0) call a%a%cp_to_coo(acoo,info)
+    if (info /= 0) then 
+      info = 4000
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    nz = acoo%get_nzeros()
+    j = 0
+    do i=1, nz
+      if (acoo%ia(i) /= acoo%ja(i)) then 
+        j = j + 1 
+        acoo%ia(j)  = acoo%ia(i)
+        acoo%ja(j)  = acoo%ja(i)
+        acoo%val(j) = acoo%val(i)
+      end if
+    end do
+    call acoo%set_nzeros(j)
+    call acoo%trim()
+    call b%mv_from(acoo)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 continue
+    call psb_erractionrestore(err_act)
+    
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    
+  end subroutine d_clip_d
+
+
+  subroutine d_clip_d_ip(a,info)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    use psb_d_base_mat_mod
+    implicit none
+    
+    class(psb_d_sparse_mat), intent(inout) :: a
+    integer,intent(out)                  :: info
+
+    Integer :: err_act
+    character(len=20)  :: name='clip_diag'
+    logical, parameter :: debug=.false.
+    type(psb_d_coo_sparse_mat), allocatable  :: acoo
+    integer :: i, j, nz
+
+    info = 0
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    allocate(acoo,stat=info)    
+    if (info == 0) call a%a%mv_to_coo(acoo,info)
+    if (info /= 0) then
+      info = 4000
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    nz = acoo%get_nzeros()
+    j = 0
+    do i=1, nz
+      if (acoo%ia(i) /= acoo%ja(i)) then 
+        j = j + 1 
+        acoo%ia(j)  = acoo%ia(i)
+        acoo%ja(j)  = acoo%ja(i)
+        acoo%val(j) = acoo%val(i)
+      end if
+    end do
+    call acoo%set_nzeros(j)
+    call acoo%trim()
+    call a%mv_from(acoo)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine d_clip_d_ip
 
   subroutine d_mv_from(a,b)
     use psb_error_mod
diff --git a/base/modules/psb_s_mat_mod.f03 b/base/modules/psb_s_mat_mod.f03
index c2543086..7ea026d3 100644
--- a/base/modules/psb_s_mat_mod.f03
+++ b/base/modules/psb_s_mat_mod.f03
@@ -54,6 +54,9 @@ module psb_s_mat_mod
     procedure, pass(a) :: s_csgetblk
     generic, public    :: csget => s_csgetptn, s_csgetrow, s_csgetblk 
     procedure, pass(a) :: csclip
+    procedure, pass(a) :: s_clip_d_ip
+    procedure, pass(a) :: s_clip_d
+    generic, public    :: clip_diag => s_clip_d_ip, s_clip_d
     procedure, pass(a) :: reall => reallocate_nz
     procedure, pass(a) :: get_neigh
     procedure, pass(a) :: s_cscnv
@@ -63,8 +66,19 @@ module psb_s_mat_mod
     procedure, pass(a) :: print => sparse_print
     procedure, pass(a) :: s_mv_from
     generic, public    :: mv_from => s_mv_from
+    procedure, pass(a) :: s_mv_to
+    generic, public    :: mv_to => s_mv_to
     procedure, pass(a) :: s_cp_from
     generic, public    :: cp_from => s_cp_from
+    procedure, pass(a) :: s_cp_to
+    generic, public    :: cp_to => s_cp_to
+    procedure, pass(a) :: s_transp_1mat
+    procedure, pass(a) :: s_transp_2mat
+    generic, public    :: transp => s_transp_1mat, s_transp_2mat
+    procedure, pass(a) :: s_transc_1mat
+    procedure, pass(a) :: s_transc_2mat
+    generic, public    :: transc => s_transc_1mat, s_transc_2mat
+
     
     
     ! Computational routines 
@@ -85,7 +99,7 @@ module psb_s_mat_mod
   private :: get_nrows, get_ncols, get_nzeros, get_size, &
        & get_state, get_dupl, is_null, is_bld, is_upd, &
        & is_asb, is_sorted, is_upper, is_lower, is_triangle, &
-       & is_unit, get_neigh, csall, csput, s_csgetrow,&
+       & is_unit, get_neigh, csall, csput, s_csgetrow,  s_clip_d_ip, s_clip_d,&
        & s_csgetblk, csclip, s_cscnv, s_cscnv_ip, &
        & reallocate_nz, free, trim, &
        & sparse_print, reinit, &
@@ -94,7 +108,9 @@ module psb_s_mat_mod
        & set_upd, set_asb, set_sorted, &
        & set_upper, set_lower, set_triangle, &
        & set_unit, get_diag, get_nz_row, s_csgetptn, &
-       & s_mv_from, s_cp_from
+       & s_mv_from, s_mv_to, s_cp_from, s_cp_to,&
+       & s_transp_1mat, s_transp_2mat, &
+       & s_transc_1mat, s_transc_2mat
 
   interface psb_sizeof
     module procedure s_sizeof
@@ -1486,6 +1502,126 @@ contains
 
   end subroutine s_cscnv_ip
 
+  subroutine s_clip_d(a,b,info)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    use psb_s_base_mat_mod
+    implicit none
+    
+    class(psb_s_sparse_mat), intent(in) :: a
+    class(psb_s_sparse_mat), intent(out) :: b
+    integer,intent(out)                  :: info
+    
+    Integer :: err_act
+    character(len=20)  :: name='clip_diag'
+    logical, parameter :: debug=.false.
+    type(psb_s_coo_sparse_mat), allocatable  :: acoo
+    integer :: i, j, nz
+    
+    info = 0
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    allocate(acoo,stat=info)    
+    if (info == 0) call a%a%cp_to_coo(acoo,info)
+    if (info /= 0) then 
+      info = 4000
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    nz = acoo%get_nzeros()
+    j = 0
+    do i=1, nz
+      if (acoo%ia(i) /= acoo%ja(i)) then 
+        j = j + 1 
+        acoo%ia(j)  = acoo%ia(i)
+        acoo%ja(j)  = acoo%ja(i)
+        acoo%val(j) = acoo%val(i)
+      end if
+    end do
+    call acoo%set_nzeros(j)
+    call acoo%trim()
+    call b%mv_from(acoo)
+    
+    call psb_erractionrestore(err_act)
+    return
+    
+9999 continue
+    call psb_erractionrestore(err_act)
+    
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    
+  end subroutine s_clip_d
+
+
+  subroutine s_clip_d_ip(a,info)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    use psb_s_base_mat_mod
+    implicit none
+    
+    class(psb_s_sparse_mat), intent(inout) :: a
+    integer,intent(out)                  :: info
+
+    Integer :: err_act
+    character(len=20)  :: name='clip_diag'
+    logical, parameter :: debug=.false.
+    type(psb_s_coo_sparse_mat), allocatable  :: acoo
+    integer :: i, j, nz
+
+    info = 0
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    allocate(acoo,stat=info)    
+    if (info == 0) call a%a%mv_to_coo(acoo,info)
+    if (info /= 0) then
+      info = 4000
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+    
+    nz = acoo%get_nzeros()
+    j = 0
+    do i=1, nz
+      if (acoo%ia(i) /= acoo%ja(i)) then 
+        j = j + 1 
+        acoo%ia(j)  = acoo%ia(i)
+        acoo%ja(j)  = acoo%ja(i)
+        acoo%val(j) = acoo%val(i)
+      end if
+    end do
+    call acoo%set_nzeros(j)
+    call acoo%trim()
+    call a%mv_from(acoo)
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine s_clip_d_ip
+
   subroutine s_mv_from(a,b)
     use psb_error_mod
     use psb_string_mod
@@ -1530,6 +1666,32 @@ contains
     end if
   end subroutine s_cp_from
 
+  subroutine s_mv_to(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(inout) :: a
+    class(psb_s_base_sparse_mat), intent(out) :: b
+    integer :: info
+
+    call b%mv_from_fmt(a%a,info)
+    
+    return
+  end subroutine s_mv_to
+
+  subroutine s_cp_to(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(in) :: a
+    class(psb_s_base_sparse_mat), intent(out) :: b
+    integer :: info
+
+    call b%cp_from_fmt(a%a,info)
+    
+    return
+  end subroutine s_cp_to
+
   subroutine s_sparse_mat_move(a,b,info)
     use psb_error_mod
     use psb_string_mod
@@ -1582,6 +1744,153 @@ contains
   end subroutine s_sparse_mat_clone
 
 
+  subroutine s_transp_1mat(a)
+    use psb_error_mod 
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(inout) :: a
+
+    Integer :: err_act, info 
+    character(len=20) :: name='transp'
+    logical, parameter :: debug=.false.
+
+
+    call psb_erractionsave(err_act) 
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name) 
+      goto 9999 
+    endif
+
+    call a%a%transp()
+
+    call psb_erractionrestore(err_act) 
+    return
+
+9999 continue 
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error() 
+      return
+    end if
+    
+  end subroutine s_transp_1mat
+
+
+  subroutine s_transp_2mat(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(out) :: a
+    class(psb_s_sparse_mat), intent(in)  :: b
+
+    Integer :: err_act, info
+    character(len=20)  :: name='transp'
+    logical, parameter :: debug=.false.
+
+
+    call psb_erractionsave(err_act)
+    if (b%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    allocate(a%a,source=b%a,stat=info)
+    if (info /= 0) then 
+      info = 4000
+      goto 9999
+    end if
+    call a%a%transp(b%a)    
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine s_transp_2mat
+
+  subroutine s_transc_1mat(a)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(inout) :: a
+
+    Integer :: err_act, info
+    character(len=20)  :: name='transc'
+    logical, parameter :: debug=.false.
+
+
+    call psb_erractionsave(err_act)
+    if (a%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    call a%a%transc()
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine s_transc_1mat
+
+
+  subroutine s_transc_2mat(a,b)
+    use psb_error_mod
+    use psb_string_mod
+    implicit none 
+    class(psb_s_sparse_mat), intent(out) :: a
+    class(psb_s_sparse_mat), intent(in)  :: b
+
+    Integer :: err_act, info
+    character(len=20)  :: name='transc'
+    logical, parameter :: debug=.false.
+
+
+    call psb_erractionsave(err_act)
+    if (b%is_null()) then 
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    endif
+
+    allocate(a%a,source=b%a,stat=info)
+    if (info /= 0) then 
+      info = 4000
+      goto 9999
+    end if
+    call a%a%transc(b%a)    
+    
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+
+  end subroutine s_transc_2mat
+
+
   subroutine reinit(a,clear)
     use psb_error_mod
     implicit none 
diff --git a/test/fileread/runs/dfs.inp b/test/fileread/runs/dfs.inp
index cc084fb9..daf4debc 100644
--- a/test/fileread/runs/dfs.inp
+++ b/test/fileread/runs/dfs.inp
@@ -7,7 +7,7 @@ BJAC                Preconditioner NONE  DIAG  BJAC
 CSR                 Storage format CSR COO JAD 
 0                   IPART: Partition method 0: BLK   2: graph (with Metis) 
 2                   ISTOPC
-00010               ITMAX
+00200               ITMAX
 01                  ITRACE
 30                  IRST     (restart for RGMRES and BiCGSTABL) 
 1.d-6               EPS
diff --git a/test/serial/d_coo_matgen.f03 b/test/serial/d_coo_matgen.f03
index a7641bea..27cb7be5 100644
--- a/test/serial/d_coo_matgen.f03
+++ b/test/serial/d_coo_matgen.f03
@@ -17,7 +17,7 @@ program d_coo_matgen
   real(psb_dpk_) :: t1, t2, tprec 
 
   ! sparse matrix and preconditioner
-  type(psb_dspmat_type) :: a
+  type(psb_d_sparse_mat) :: a
   type(psb_dprec_type)  :: prec
   ! descriptor
   type(psb_desc_type)   :: desc_a
@@ -145,7 +145,7 @@ contains
     type(psb_desc_type)            :: desc_a
     integer                        :: ictxt, info
     character                      :: afmt*5
-    type(psb_dspmat_type)    :: a
+    type(psb_d_sparse_mat)   :: a
     real(psb_dpk_)           :: zt(nb),glob_x,glob_y,glob_z
     integer                  :: m,n,nnz,glob_row,nlr,i,ii,ib,k
     integer                  :: x,y,z,ia,indx_owner
@@ -161,10 +161,10 @@ contains
     real(psb_dpk_),parameter   :: rhs=0.d0,one=1.d0,zero=0.d0
     real(psb_dpk_)   :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcpy, tmov
     real(psb_dpk_)   :: a1, a2, a3, a4, b1, b2, b3 
-    external           :: a1, a2, a3, a4, b1, b2, b3
-    integer            :: err_act
+    external         :: a1, a2, a3, a4, b1, b2, b3
+    integer          :: err_act
 
-    character(len=20)  :: name, ch_err
+    character(len=20)  :: name, ch_err, asbfmt
 
     info = 0
     name = 'create_matrix'
@@ -404,8 +404,9 @@ contains
     tmov = psb_wtime()-t1
 !!$    call acsr%print(22)
     if(iam == psb_root_) then
+      asbfmt = a%get_fmt()
       write(*,'("The matrix has been generated and assembled in ",a3," format.")')&
-           &   a%fida(1:3)
+           &   asbfmt
       write(*,'("-allocation  time : ",es12.5)') talc
       write(*,'("-coeff. gen. time : ",es12.5)') tgen
       write(*,'("-assembly    time : ",es12.5)') tasb
diff --git a/test/serial/d_matgen.f03 b/test/serial/d_matgen.f03
index b92cd23c..af31936c 100644
--- a/test/serial/d_matgen.f03
+++ b/test/serial/d_matgen.f03
@@ -18,7 +18,7 @@ program d_matgen
   real(psb_dpk_) :: t1, t2, tprec 
 
   ! sparse matrix and preconditioner
-  type(psb_dspmat_type) :: a
+  type(psb_d_sparse_mat) :: a
   type(psb_dprec_type)  :: prec
   ! descriptor
   type(psb_desc_type)   :: desc_a
@@ -147,7 +147,7 @@ contains
     type(psb_desc_type)            :: desc_a
     integer                        :: ictxt, info
     character                      :: afmt*5
-    type(psb_dspmat_type)    :: a
+    type(psb_d_sparse_mat)    :: a
     real(psb_dpk_)           :: zt(nb),glob_x,glob_y,glob_z
     integer                  :: m,n,nnz,glob_row,nlr,i,ii,ib,k
     integer                  :: x,y,z,ia,indx_owner
@@ -165,8 +165,8 @@ contains
     real(psb_dpk_),parameter   :: rhs=0.d0,one=1.d0,zero=0.d0
     real(psb_dpk_)   :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcpy, tmov
     real(psb_dpk_)   :: a1, a2, a3, a4, b1, b2, b3 
-    external           :: a1, a2, a3, a4, b1, b2, b3
-    integer            :: err_act
+    external         :: a1, a2, a3, a4, b1, b2, b3
+    integer          :: err_act
 
     character(len=20)  :: name, ch_err
 
diff --git a/test/serial/psb_d_cxx_impl.f03 b/test/serial/psb_d_cxx_impl.f03
index b07990e8..cdca88a6 100644
--- a/test/serial/psb_d_cxx_impl.f03
+++ b/test/serial/psb_d_cxx_impl.f03
@@ -14,6 +14,7 @@
 
 subroutine d_cxx_csmv_impl(alpha,a,x,beta,y,info,trans) 
   use psb_error_mod
+  use psb_string_mod
   use psb_d_cxx_mat_mod, psb_protect_name => d_cxx_csmv_impl
   implicit none 
   class(psb_d_cxx_sparse_mat), intent(in) :: a
@@ -40,14 +41,13 @@ subroutine d_cxx_csmv_impl(alpha,a,x,beta,y,info,trans)
   end if
 
   if (.not.a%is_asb()) then 
-    write(0,*) 'Error: csmv called on an unassembled mat'
     info = 1121
     call psb_errpush(info,name)
     goto 9999
   endif
 
 
-  tra = ((trans_=='T').or.(trans_=='t'))
+  tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
 
   if (tra) then 
     m = a%get_ncols()
@@ -57,227 +57,249 @@ subroutine d_cxx_csmv_impl(alpha,a,x,beta,y,info,trans)
     m = a%get_nrows()
   end if
 
+  call d_cxx_csmv_inner(m,n,alpha,a%irp,a%ja,a%val,&
+       & a%is_triangle(),a%is_unit(),&
+       & x,beta,y,tra) 
 
-  if (alpha == dzero) then
-    if (beta == dzero) then
-      do i = 1, m
-        y(i) = dzero
-      enddo
-    else
-      do  i = 1, m
-        y(i) = beta*y(i)
-      end do
-    endif
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+
+  if (err_act == psb_act_abort_) then
+    call psb_error()
     return
   end if
+  return
 
-  if (.not.tra) then 
+contains
+  subroutine d_cxx_csmv_inner(m,n,alpha,irp,ja,val,is_triangle,is_unit,&
+       & x,beta,y,tra) 
+    integer, intent(in)             :: m,n,irp(*),ja(*)
+    real(psb_dpk_), intent(in)      :: alpha, beta, x(*),val(*)
+    real(psb_dpk_), intent(inout)   :: y(*)
+    logical, intent(in)             :: is_triangle,is_unit,tra
 
-    if (beta == dzero) then 
 
-      if (alpha == done) then 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = acc
+    integer   :: i,j,k, ir, jc
+    real(psb_dpk_) :: acc
+
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i) = beta*y(i)
         end do
+      endif
+      return
+    end if
 
-      else if (alpha == -done) then 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = -acc
-        end do
+    if (.not.tra) then 
 
-      else 
+      if (beta == dzero) then 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = alpha*acc
-        end do
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = acc
+          end do
 
-      end if
+        else if (alpha == -done) then 
 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc - val(j) * x(ja(j))          
+            enddo
+            y(i) = acc
+          end do
 
-    else if (beta == done) then 
+        else 
 
-      if (alpha == done) then 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = y(i) + acc
-        end do
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = alpha*acc
+          end do
 
-      else if (alpha == -done) then 
+        end if
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = y(i) -acc
-        end do
 
-      else 
+      else if (beta == done) then 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = y(i) + alpha*acc
-        end do
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = y(i) + acc
+          end do
 
-      end if
+        else if (alpha == -done) then 
 
-    else if (beta == -done) then 
+          do i=1,m 
+            acc  = y(i)
+            do j=irp(i), irp(i+1)-1
+              acc  = acc - val(j) * x(ja(j))          
+            enddo
+            y(i) = acc
+          end do
 
-      if (alpha == done) then 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = -y(i) + acc
-        end do
+        else 
 
-      else if (alpha == -done) then 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = y(i) + alpha*acc
+          end do
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = -y(i) -acc
-        end do
+        end if
 
-      else 
+      else if (beta == -done) then 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = -y(i) + alpha*acc
-        end do
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = -y(i)
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = acc
+          end do
 
-      end if
+        else if (alpha == -done) then 
 
-    else 
+          do i=1,m 
+            acc  = y(i)
+            do j=irp(i), irp(i+1)-1
+              acc  = acc - val(j) * x(ja(j))          
+            enddo
+            y(i) =  acc
+          end do
 
-      if (alpha == done) then 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = beta*y(i) + acc
-        end do
+        else 
 
-      else if (alpha == -done) then 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = -y(i) + alpha*acc
+          end do
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = beta*y(i) - acc
-        end do
+        end if
 
       else 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j))          
-          enddo
-          y(i) = beta*y(i) + alpha*acc
-        end do
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = beta*y(i) + acc
+          end do
 
-      end if
+        else if (alpha == -done) then 
 
-    end if
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = beta*y(i) - acc
+          end do
 
-  else if (tra) then 
+        else 
 
-    if (beta == dzero) then 
-      do i=1, m
-        y(i) = dzero
-      end do
-    else if (beta == done) then 
-      ! Do nothing
-    else if (beta == -done) then 
-      do i=1, m
-        y(i) = -y(i) 
-      end do
-    else
-      do i=1, m
-        y(i) = beta*y(i) 
-      end do
-    end if
+          do i=1,m 
+            acc  = dzero
+            do j=irp(i), irp(i+1)-1
+              acc  = acc + val(j) * x(ja(j))          
+            enddo
+            y(i) = beta*y(i) + alpha*acc
+          end do
 
-    if (alpha.eq.done) then
+        end if
 
-      do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir) = y(ir) +  a%val(j)*x(i)
-        end do
-      enddo
+      end if
 
-    else if (alpha.eq.-done) then
+    else if (tra) then 
 
-      do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir) = y(ir) -  a%val(j)*x(i)
+      if (beta == dzero) then 
+        do i=1, m
+          y(i) = dzero
         end do
-      enddo
+      else if (beta == done) then 
+        ! Do nothing
+      else if (beta == -done) then 
+        do i=1, m
+          y(i) = -y(i) 
+        end do
+      else
+        do i=1, m
+          y(i) = beta*y(i) 
+        end do
+      end if
 
-    else                    
+      if (alpha == done) then
 
-      do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir) = y(ir) + alpha*a%val(j)*x(i)
-        end do
-      enddo
+        do i=1,n
+          do j=irp(i), irp(i+1)-1
+            ir = ja(j)
+            y(ir) = y(ir) +  val(j)*x(i)
+          end do
+        enddo
 
-    end if
+      else if (alpha == -done) then
 
-  endif
+        do i=1,n
+          do j=irp(i), irp(i+1)-1
+            ir = ja(j)
+            y(ir) = y(ir) -  val(j)*x(i)
+          end do
+        enddo
 
-  if (a%is_unit()) then 
-    do i=1, min(m,n)
-      y(i) = y(i) + alpha*x(i)
-    end do
-  end if
+      else                    
 
-  call psb_erractionrestore(err_act)
-  return
-9999 continue
-  call psb_erractionrestore(err_act)
+        do i=1,n
+          do j=irp(i), irp(i+1)-1
+            ir = ja(j)
+            y(ir) = y(ir) + alpha*val(j)*x(i)
+          end do
+        enddo
+
+      end if
+
+    endif
+
+    if (is_triangle.and.is_unit) then 
+      do i=1, min(m,n)
+        y(i) = y(i) + alpha*x(i)
+      end do
+    end if
+
+
+  end subroutine d_cxx_csmv_inner
 
-  if (err_act == psb_act_abort_) then
-    call psb_error()
-    return
-  end if
-  return
 
 end subroutine d_cxx_csmv_impl
 
 subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans) 
   use psb_error_mod
+  use psb_string_mod
   use psb_d_cxx_mat_mod, psb_protect_name => d_cxx_csmm_impl
   implicit none 
   class(psb_d_cxx_sparse_mat), intent(in) :: a
@@ -303,12 +325,12 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
     trans_ = 'N'
   end if
 
-  tra = ((trans_=='T').or.(trans_=='t'))
   if (.not.a%is_asb()) then 
     info = 1121
     call psb_errpush(info,name)
     goto 9999
   endif
+  tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
 
   if (tra) then 
     m = a%get_ncols()
@@ -326,15 +348,43 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
     call psb_errpush(info,name,a_err='allocate')
     goto 9999
   end if
+  
+  call  d_cxx_csmm_inner(m,n,nc,alpha,a%irp,a%ja,a%val, &
+       & a%is_triangle(),a%is_unit(),x,size(x,1), &
+       & beta,y,size(y,1),tra,acc) 
+
+
+  call psb_erractionrestore(err_act)
+  return
+9999 continue
+  call psb_erractionrestore(err_act)
+
+  if (err_act == psb_act_abort_) then
+    call psb_error()
+    return
+  end if
+  return
 
+contains
+  subroutine d_cxx_csmm_inner(m,n,nc,alpha,irp,ja,val,&
+       & is_triangle,is_unit,x,ldx,beta,y,ldy,tra,acc) 
+    integer, intent(in)             :: m,n,ldx,ldy,nc,irp(*),ja(*)
+    real(psb_dpk_), intent(in)      :: alpha, beta, x(ldx,*),val(*)
+    real(psb_dpk_), intent(inout)   :: y(ldy,*)
+    logical, intent(in)             :: is_triangle,is_unit,tra
+
+    real(psb_dpk_), intent(inout)   :: acc(*)
+    integer   :: i,j,k, ir, jc
+    
+    
   if (alpha == dzero) then
     if (beta == dzero) then
       do i = 1, m
-        y(i,:) = dzero
+        y(i,1:nc) = dzero
       enddo
     else
       do  i = 1, m
-        y(i,:) = beta*y(i,:)
+        y(i,1:nc) = beta*y(i,1:nc)
       end do
     endif
     return
@@ -346,31 +396,31 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
       if (alpha == done) then 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = acc
+          y(i,1:nc) = acc(1:nc)
         end do
 
       else if (alpha == -done) then 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) - val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = -acc
+          y(i,1:nc) = acc(1:nc)
         end do
 
       else 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = alpha*acc
+          y(i,1:nc) = alpha*acc(1:nc)
         end do
 
       end if
@@ -380,31 +430,31 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
       if (alpha == done) then 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = y(i,1:nc)
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = y(i,:) + acc
+          y(i,1:nc) = acc(1:nc)
         end do
 
       else if (alpha == -done) then 
 
-        do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+        do i=1,m
+          acc(1:nc)  = y(i,1:nc)
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) - val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = y(i,:) -acc
+          y(i,1:nc) = acc(1:nc)
         end do
 
       else 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = y(i,:) + alpha*acc
+          y(i,1:nc) = y(i,1:nc) + alpha*acc(1:nc)
         end do
 
       end if
@@ -413,31 +463,31 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
       if (alpha == done) then 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = -y(i,:) + acc
+          y(i,1:nc) = -y(i,1:nc) + acc(1:nc)
         end do
 
       else if (alpha == -done) then 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = -y(i,:) -acc
+          y(i,1:nc) = -y(i,1:nc) -acc(1:nc)
         end do
 
       else 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = -y(i,:) + alpha*acc
+          y(i,1:nc) = -y(i,1:nc) + alpha*acc(1:nc)
         end do
 
       end if
@@ -446,31 +496,31 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
       if (alpha == done) then 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = beta*y(i,:) + acc
+          y(i,1:nc) = beta*y(i,1:nc) + acc(1:nc)
         end do
 
       else if (alpha == -done) then 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = beta*y(i,:) - acc
+          y(i,1:nc) = beta*y(i,1:nc) - acc(1:nc)
         end do
 
       else 
 
         do i=1,m 
-          acc  = dzero
-          do j=a%irp(i), a%irp(i+1)-1
-            acc  = acc + a%val(j) * x(a%ja(j),:)          
+          acc(1:nc)  = dzero
+          do j=irp(i), irp(i+1)-1
+            acc(1:nc)  = acc(1:nc) + val(j) * x(ja(j),1:nc)          
           enddo
-          y(i,:) = beta*y(i,:) + alpha*acc
+          y(i,1:nc) = beta*y(i,1:nc) + alpha*acc(1:nc)
         end do
 
       end if
@@ -481,44 +531,44 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
     if (beta == dzero) then 
       do i=1, m
-        y(i,:) = dzero
+        y(i,1:nc) = dzero
       end do
     else if (beta == done) then 
       ! Do nothing
     else if (beta == -done) then 
       do i=1, m
-        y(i,:) = -y(i,:) 
+        y(i,1:nc) = -y(i,1:nc) 
       end do
     else
       do i=1, m
-        y(i,:) = beta*y(i,:) 
+        y(i,1:nc) = beta*y(i,1:nc) 
       end do
     end if
 
-    if (alpha.eq.done) then
+    if (alpha == done) then
 
       do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir,:) = y(ir,:) +  a%val(j)*x(i,:)
+        do j=irp(i), irp(i+1)-1
+          ir = ja(j)
+          y(ir,1:nc) = y(ir,1:nc) +  val(j)*x(i,1:nc)
         end do
       enddo
 
-    else if (alpha.eq.-done) then
+    else if (alpha == -done) then
 
       do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir,:) = y(ir,:) -  a%val(j)*x(i,:)
+        do j=irp(i), irp(i+1)-1
+          ir = ja(j)
+          y(ir,1:nc) = y(ir,1:nc) -  val(j)*x(i,1:nc)
         end do
       enddo
 
     else                    
 
       do i=1,n
-        do j=a%irp(i), a%irp(i+1)-1
-          ir = a%ja(j)
-          y(ir,:) = y(ir,:) + alpha*a%val(j)*x(i,:)
+        do j=irp(i), irp(i+1)-1
+          ir = ja(j)
+          y(ir,1:nc) = y(ir,1:nc) + alpha*val(j)*x(i,1:nc)
         end do
       enddo
 
@@ -526,28 +576,20 @@ subroutine d_cxx_csmm_impl(alpha,a,x,beta,y,info,trans)
 
   endif
 
-  if (a%is_unit()) then 
+  if (is_triangle.and.is_unit) then 
     do i=1, min(m,n)
-      y(i,:) = y(i,:) + alpha*x(i,:)
+      y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc)
     end do
   end if
 
-  call psb_erractionrestore(err_act)
-  return
-9999 continue
-  call psb_erractionrestore(err_act)
-
-  if (err_act == psb_act_abort_) then
-    call psb_error()
-    return
-  end if
-  return
+end subroutine d_cxx_csmm_inner
 
 end subroutine d_cxx_csmm_impl
 
 
 subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans) 
   use psb_error_mod
+  use psb_string_mod
   use psb_d_cxx_mat_mod, psb_protect_name => d_cxx_cssv_impl
   implicit none 
   class(psb_d_cxx_sparse_mat), intent(in) :: a
@@ -567,7 +609,6 @@ subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans)
 
   info = 0
   call psb_erractionsave(err_act)
-
   if (present(trans)) then
     trans_ = trans
   else
@@ -579,7 +620,7 @@ subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans)
     goto 9999
   endif
 
-  tra = ((trans_=='T').or.(trans_=='t'))
+  tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
   m = a%get_nrows()
 
   if (.not. (a%is_triangle())) then 
@@ -588,7 +629,18 @@ subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans)
     goto 9999
   end if
 
+  if (size(x)<m) then 
+    info = 36
+    call psb_errpush(info,name,i_err=(/3,m,0,0,0/))
+    goto 9999
+  end if
 
+  if (size(y)<m) then 
+    info = 36
+    call psb_errpush(info,name,i_err=(/5,m,0,0,0/))
+    goto 9999
+  end if
+  
   if (alpha == dzero) then
     if (beta == dzero) then
       do i = 1, m
@@ -603,21 +655,31 @@ subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans)
   end if
 
   if (beta == dzero) then 
-    call inner_cxxsv(tra,a,x,y)
-    do  i = 1, m
-      y(i) = alpha*y(i)
-    end do
+
+    call inner_cxxsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
+         & a%irp,a%ja,a%val,x,y) 
+    if (alpha == done) then 
+      ! do nothing
+    else if (alpha == -done) then 
+      do  i = 1, m
+        y(i) = -y(i)
+      end do
+    else
+      do  i = 1, m
+        y(i) = alpha*y(i)
+      end do
+    end if
   else 
     allocate(tmp(m), stat=info) 
     if (info /= 0) then 
-      write(0,*) 'Memory allocation error in CXXSV '
       return
     end if
-    tmp(1:m) = x(1:m)
-    call inner_cxxsv(tra,a,tmp,y)
-    do  i = 1, m
-      y(i) = alpha*tmp(i) + beta*y(i)
-    end do
+
+    call inner_cxxsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),&
+         & a%irp,a%ja,a%val,x,tmp) 
+
+    call psb_geaxpby(m,alpha,tmp,beta,y,info)
+
   end if
 
   call psb_erractionrestore(err_act)
@@ -634,53 +696,54 @@ subroutine d_cxx_cssv_impl(alpha,a,x,beta,y,info,trans)
 
 contains 
 
-  subroutine inner_cxxsv(tra,a,x,y) 
+  subroutine inner_cxxsv(tra,lower,unit,n,irp,ja,val,x,y) 
     implicit none 
-    logical, intent(in)                 :: tra  
-    class(psb_d_cxx_sparse_mat), intent(in) :: a
-    real(psb_dpk_), intent(in)          :: x(:)
-    real(psb_dpk_), intent(out)         :: y(:)
+    logical, intent(in)                 :: tra,lower,unit  
+    integer, intent(in)                 :: irp(*), ja(*),n
+    real(psb_dpk_), intent(in)          :: val(*)
+    real(psb_dpk_), intent(in)          :: x(*)
+    real(psb_dpk_), intent(out)         :: y(*)
 
     integer :: i,j,k,m, ir, jc
     real(psb_dpk_) :: acc
 
     if (.not.tra) then 
 
-      if (a%is_lower()) then 
-        if (a%is_unit()) then 
-          do i=1, a%get_nrows()
+      if (lower) then 
+        if (unit) then 
+          do i=1, n
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j))
+            do j=irp(i), irp(i+1)-1
+              acc = acc + val(j)*y(ja(j))
             end do
             y(i) = x(i) - acc
           end do
-        else if (.not.a%is_unit()) then 
-          do i=1, a%get_nrows()
+        else if (.not.unit) then 
+          do i=1, n
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-2
-              acc = acc + a%val(j)*x(a%ja(j))
+            do j=irp(i), irp(i+1)-2
+              acc = acc + val(j)*y(ja(j))
             end do
-            y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1)
+            y(i) = (x(i) - acc)/val(irp(i+1)-1)
           end do
         end if
-      else if (a%is_upper()) then 
+      else if (.not.lower) then 
 
-        if (a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1 
+        if (unit) then 
+          do i=n, 1, -1 
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j))
+            do j=irp(i), irp(i+1)-1
+              acc = acc + val(j)*y(ja(j))
             end do
             y(i) = x(i) - acc
           end do
-        else if (.not.a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1 
+        else if (.not.unit) then 
+          do i=n, 1, -1 
             acc = dzero 
-            do j=a%irp(i)+1, a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j))
+            do j=irp(i)+1, irp(i+1)-1
+              acc = acc + val(j)*y(ja(j))
             end do
-            y(i) = (x(i) - acc)/a%val(a%irp(i))
+            y(i) = (x(i) - acc)/val(irp(i))
           end do
         end if
 
@@ -688,46 +751,46 @@ contains
 
     else if (tra) then 
 
-      do i=1, a%get_nrows()
+      do i=1, n
         y(i) = x(i)
       end do
 
-      if (a%is_lower()) then 
-        if (a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1
+      if (lower) then 
+        if (unit) then 
+          do i=n, 1, -1
             acc = y(i) 
-            do j=a%irp(i), a%irp(i+1)-1
-              jc    = a%ja(j)
-              y(jc) = y(jc) - a%val(j)*acc 
+            do j=irp(i), irp(i+1)-1
+              jc    = ja(j)
+              y(jc) = y(jc) - val(j)*acc 
             end do
           end do
-        else if (.not.a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1
-            y(i) = y(i)/a%val(a%irp(i+1)-1)
+        else if (.not.unit) then 
+          do i=n, 1, -1
+            y(i) = y(i)/val(irp(i+1)-1)
             acc  = y(i) 
-            do j=a%irp(i), a%irp(i+1)-2
-              jc    = a%ja(j)
-              y(jc) = y(jc) - a%val(j)*acc 
+            do j=irp(i), irp(i+1)-2
+              jc    = ja(j)
+              y(jc) = y(jc) - val(j)*acc 
             end do
           end do
         end if
-      else if (a%is_upper()) then 
+      else if (.not.lower) then 
 
-        if (a%is_unit()) then 
-          do i=1, a%get_nrows()
+        if (unit) then 
+          do i=1, n
             acc  = y(i) 
-            do j=a%irp(i), a%irp(i+1)-1
-              jc    = a%ja(j)
-              y(jc) = y(jc) - a%val(j)*acc 
+            do j=irp(i), irp(i+1)-1
+              jc    = ja(j)
+              y(jc) = y(jc) - val(j)*acc 
             end do
           end do
-        else if (.not.a%is_unit()) then 
-          do i=1, a%get_nrows()
-            y(i) = y(i)/a%val(a%irp(i))
+        else if (.not.unit) then 
+          do i=1, n
+            y(i) = y(i)/val(irp(i))
             acc  = y(i) 
-            do j=a%irp(i)+1, a%irp(i+1)-1
-              jc    = a%ja(j)
-              y(jc) = y(jc) - a%val(j)*acc 
+            do j=irp(i)+1, irp(i+1)-1
+              jc    = ja(j)
+              y(jc) = y(jc) - val(j)*acc 
             end do
           end do
         end if
@@ -742,6 +805,7 @@ end subroutine d_cxx_cssv_impl
 
 subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans) 
   use psb_error_mod
+  use psb_string_mod
   use psb_d_cxx_mat_mod, psb_protect_name => d_cxx_cssm_impl
   implicit none 
   class(psb_d_cxx_sparse_mat), intent(in) :: a
@@ -756,7 +820,7 @@ subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans)
   real(psb_dpk_), allocatable :: tmp(:,:)
   logical   :: tra
   Integer :: err_act
-  character(len=20)  :: name='d_base_csmm'
+  character(len=20)  :: name='d_cxx_cssm'
   logical, parameter :: debug=.false.
 
   info = 0
@@ -774,12 +838,12 @@ subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans)
   endif
 
 
-  tra = ((trans_=='T').or.(trans_=='t'))
+  tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C')
+
   m   = a%get_nrows()
   nc  = min(size(x,2) , size(y,2)) 
 
   if (.not. (a%is_triangle())) then 
-    write(0,*) 'Called SM on a non-triangular mat!'
     info = 1121
     call psb_errpush(info,name)
     goto 9999
@@ -800,9 +864,10 @@ subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans)
   end if
 
   if (beta == dzero) then 
-    call inner_cxxsm(tra,a,x,y,info)
+    call inner_cxxsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
+         & a%irp,a%ja,a%val,x,size(x,1),y,size(y,1),info) 
     do  i = 1, m
-      y(i,:) = alpha*y(i,:)
+      y(i,1:nc) = alpha*y(i,1:nc)
     end do
   else 
     allocate(tmp(m,nc), stat=info) 
@@ -812,10 +877,10 @@ subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans)
       goto 9999
     end if
 
-    tmp(1:m,:) = x(1:m,:)
-    call inner_cxxsm(tra,a,tmp,y,info)
+    call inner_cxxsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,&
+         & a%irp,a%ja,a%val,x,size(x,1),tmp,size(tmp,1),info) 
     do  i = 1, m
-      y(i,:) = alpha*tmp(i,:) + beta*y(i,:)
+      y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc)
     end do
   end if
 
@@ -841,18 +906,19 @@ subroutine d_cxx_cssm_impl(alpha,a,x,beta,y,info,trans)
 
 contains 
 
-  subroutine inner_cxxsm(tra,a,x,y,info) 
+  subroutine inner_cxxsm(tra,lower,unit,nr,nc,&
+       & irp,ja,val,x,ldx,y,ldy,info) 
     implicit none 
-    logical, intent(in)                 :: tra  
-    class(psb_d_cxx_sparse_mat), intent(in) :: a
-    real(psb_dpk_), intent(in)          :: x(:,:)
-    real(psb_dpk_), intent(out)         :: y(:,:)
+    logical, intent(in)                 :: tra,lower,unit
+    integer, intent(in)                 :: nr,nc,ldx,ldy,irp(*),ja(*)
+    real(psb_dpk_), intent(in)          :: val(*), x(ldx,*)
+    real(psb_dpk_), intent(out)         :: y(ldy,*)
     integer, intent(out)                :: info
     integer :: i,j,k,m, ir, jc
     real(psb_dpk_), allocatable  :: acc(:)
 
     info = 0
-    allocate(acc(size(x,2)), stat=info)
+    allocate(acc(nc), stat=info)
     if(info /= 0) then
       info=4010
       return
@@ -861,41 +927,41 @@ contains
 
     if (.not.tra) then 
 
-      if (a%is_lower()) then 
-        if (a%is_unit()) then 
-          do i=1, a%get_nrows()
+      if (lower) then 
+        if (unit) then 
+          do i=1, nr
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j),:)
+            do j=irp(i), irp(i+1)-1
+              acc = acc + val(j)*y(ja(j),1:nc)
             end do
-            y(i,:) = x(i,:) - acc
+            y(i,1:nc) = x(i,1:nc) - acc
           end do
-        else if (.not.a%is_unit()) then 
-          do i=1, a%get_nrows()
+        else if (.not.unit) then 
+          do i=1, nr
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-2
-              acc = acc + a%val(j)*x(a%ja(j),:)
+            do j=irp(i), irp(i+1)-2
+              acc = acc + val(j)*y(ja(j),1:nc)
             end do
-            y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1)
+            y(i,1:nc) = (x(i,1:nc) - acc)/val(irp(i+1)-1)
           end do
         end if
-      else if (a%is_upper()) then 
+      else if (.not.lower) then 
 
-        if (a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1 
+        if (unit) then 
+          do i=nr, 1, -1 
             acc = dzero 
-            do j=a%irp(i), a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j),:)
+            do j=irp(i), irp(i+1)-1
+              acc = acc + val(j)*y(ja(j),1:nc)
             end do
-            y(i,:) = x(i,:) - acc
+            y(i,1:nc) = x(i,1:nc) - acc
           end do
-        else if (.not.a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1 
+        else if (.not.unit) then 
+          do i=nr, 1, -1 
             acc = dzero 
-            do j=a%irp(i)+1, a%irp(i+1)-1
-              acc = acc + a%val(j)*x(a%ja(j),:)
+            do j=irp(i)+1, irp(i+1)-1
+              acc = acc + val(j)*y(ja(j),1:nc)
             end do
-            y(i,:) = (x(i,:) - acc)/a%val(a%irp(i))
+            y(i,1:nc) = (x(i,1:nc) - acc)/val(irp(i))
           end do
         end if
 
@@ -903,46 +969,46 @@ contains
 
     else if (tra) then 
 
-      do i=1, a%get_nrows()
-        y(i,:) = x(i,:)
+      do i=1, nr
+        y(i,1:nc) = x(i,1:nc)
       end do
 
-      if (a%is_lower()) then 
-        if (a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1
-            acc = y(i,:) 
-            do j=a%irp(i), a%irp(i+1)-1
-              jc    = a%ja(j)
-              y(jc,:) = y(jc,:) - a%val(j)*acc 
+      if (lower) then 
+        if (unit) then 
+          do i=nr, 1, -1
+            acc = y(i,1:nc) 
+            do j=irp(i), irp(i+1)-1
+              jc    = ja(j)
+              y(jc,1:nc) = y(jc,1:nc) - val(j)*acc 
             end do
           end do
-        else if (.not.a%is_unit()) then 
-          do i=a%get_nrows(), 1, -1
-            y(i,:) = y(i,:)/a%val(a%irp(i+1)-1)
-            acc  = y(i,:) 
-            do j=a%irp(i), a%irp(i+1)-2
-              jc    = a%ja(j)
-              y(jc,:) = y(jc,:) - a%val(j)*acc 
+        else if (.not.unit) then 
+          do i=nr, 1, -1
+            y(i,1:nc) = y(i,1:nc)/val(irp(i+1)-1)
+            acc  = y(i,1:nc) 
+            do j=irp(i), irp(i+1)-2
+              jc    = ja(j)
+              y(jc,1:nc) = y(jc,1:nc) - val(j)*acc 
             end do
           end do
         end if
-      else if (a%is_upper()) then 
-
-        if (a%is_unit()) then 
-          do i=1, a%get_nrows()
-            acc  = y(i,:) 
-            do j=a%irp(i), a%irp(i+1)-1
-              jc    = a%ja(j)
-              y(jc,:) = y(jc,:) - a%val(j)*acc 
+      else if (.not.lower) then 
+
+        if (unit) then 
+          do i=1, nr
+            acc  = y(i,1:nc) 
+            do j=irp(i), irp(i+1)-1
+              jc    = ja(j)
+              y(jc,1:nc) = y(jc,1:nc) - val(j)*acc 
             end do
           end do
-        else if (.not.a%is_unit()) then 
-          do i=1, a%get_nrows()
-            y(i,:) = y(i,:)/a%val(a%irp(i))
-            acc    = y(i,:) 
-            do j=a%irp(i)+1, a%irp(i+1)-1
-              jc      = a%ja(j)
-              y(jc,:) = y(jc,:) - a%val(j)*acc 
+        else if (.not.unit) then 
+          do i=1, nr
+            y(i,1:nc) = y(i,1:nc)/val(irp(i))
+            acc    = y(i,1:nc) 
+            do j=irp(i)+1, irp(i+1)-1
+              jc      = ja(j)
+              y(jc,1:nc) = y(jc,1:nc) - val(j)*acc 
             end do
           end do
         end if
@@ -993,6 +1059,181 @@ end function d_cxx_csnmi_impl
 !=====================================   
 
 
+subroutine d_cxx_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
+     & jmin,jmax,iren,append,nzin,rscale,cscale)
+  ! Output is always in  COO format 
+  use psb_error_mod
+  use psb_const_mod
+  use psb_error_mod
+  use psb_d_base_mat_mod
+  use psb_d_cxx_mat_mod, psb_protect_name => d_cxx_csgetptn_impl
+  implicit none
+
+  class(psb_d_cxx_sparse_mat), intent(in) :: a
+  integer, intent(in)                  :: imin,imax
+  integer, intent(out)                 :: nz
+  integer, allocatable, intent(inout)  :: ia(:), ja(:)
+  integer,intent(out)                  :: info
+  logical, intent(in), optional        :: append
+  integer, intent(in), optional        :: iren(:)
+  integer, intent(in), optional        :: jmin,jmax, nzin
+  logical, intent(in), optional        :: rscale,cscale
+
+  logical :: append_, rscale_, cscale_ 
+  integer :: nzin_, jmin_, jmax_, err_act, i
+  character(len=20)  :: name='csget'
+  logical, parameter :: debug=.false.
+
+  call psb_erractionsave(err_act)
+  info = 0
+
+  if (present(jmin)) then
+    jmin_ = jmin
+  else
+    jmin_ = 1
+  endif
+  if (present(jmax)) then
+    jmax_ = jmax
+  else
+    jmax_ = a%get_ncols()
+  endif
+
+  if ((imax<imin).or.(jmax_<jmin_)) then 
+    nz = 0
+    return
+  end if
+
+  if (present(append)) then
+    append_=append
+  else
+    append_=.false.
+  endif
+  if ((append_).and.(present(nzin))) then 
+    nzin_ = nzin
+  else
+    nzin_ = 0
+  endif
+  if (present(rscale)) then 
+    rscale_ = rscale
+  else
+    rscale_ = .false.
+  endif
+  if (present(cscale)) then 
+    cscale_ = cscale
+  else
+    cscale_ = .false.
+  endif
+  if ((rscale_.or.cscale_).and.(present(iren))) then 
+    info = 583
+    call psb_errpush(info,name,a_err='iren (rscale.or.cscale)')
+    goto 9999
+  end if
+
+  call cxx_getptn(imin,imax,jmin_,jmax_,a,nz,ia,ja,nzin_,append_,info,iren)
+  
+  if (rscale_) then 
+    do i=nzin_+1, nzin_+nz
+      ia(i) = ia(i) - imin + 1
+    end do
+  end if
+  if (cscale_) then 
+    do i=nzin_+1, nzin_+nz
+      ja(i) = ja(i) - jmin_ + 1
+    end do
+  end if
+
+  if (info /= 0) goto 9999
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_erractionrestore(err_act)
+
+  if (err_act == psb_act_abort_) then
+    call psb_error()
+    return
+  end if
+  return
+
+contains
+
+  subroutine cxx_getptn(imin,imax,jmin,jmax,a,nz,ia,ja,nzin,append,info,&
+       & iren)
+
+    use psb_const_mod
+    use psb_error_mod
+    use psb_realloc_mod
+    use psb_sort_mod
+    implicit none
+
+    class(psb_d_cxx_sparse_mat), intent(in)    :: a
+    integer                              :: imin,imax,jmin,jmax
+    integer, intent(out)                 :: nz
+    integer, allocatable, intent(inout)  :: ia(:), ja(:)
+    integer, intent(in)                  :: nzin
+    logical, intent(in)                  :: append
+    integer                              :: info
+    integer, optional                    :: iren(:)
+    integer  :: nzin_, nza, idx,i,j,k, nzt, irw, lrw
+    integer  :: debug_level, debug_unit
+    character(len=20) :: name='cxx_getptn'
+
+    debug_unit  = psb_get_debug_unit()
+    debug_level = psb_get_debug_level()
+
+    nza = a%get_nzeros()
+    irw = imin
+    lrw = min(imax,a%get_nrows())
+    if (irw<0) then 
+      info = 2
+      return
+    end if
+
+    if (append) then 
+      nzin_ = nzin
+    else
+      nzin_ = 0
+    endif
+
+    nzt = a%irp(lrw+1)-a%irp(irw)
+    nz = 0 
+
+
+    call psb_ensure_size(nzin_+nzt,ia,info)
+    if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
+
+    if (info /= 0) return
+    
+    if (present(iren)) then 
+      do i=irw, lrw
+        do j=a%irp(i), a%irp(i+1) - 1
+          if ((jmin <= a%ja(j)).and.(a%ja(j)<=jmax)) then 
+            nzin_ = nzin_ + 1
+            nz    = nz + 1
+            ia(nzin_)  = iren(i)
+            ja(nzin_)  = iren(a%ja(j))
+          end if
+        enddo
+      end do
+    else
+      do i=irw, lrw
+        do j=a%irp(i), a%irp(i+1) - 1
+          if ((jmin <= a%ja(j)).and.(a%ja(j)<=jmax)) then 
+            nzin_ = nzin_ + 1
+            nz    = nz + 1
+            ia(nzin_)  = (i)
+            ja(nzin_)  = (a%ja(j))
+          end if
+        enddo
+      end do
+    end if
+
+  end subroutine cxx_getptn
+  
+end subroutine d_cxx_csgetptn_impl
+
+
 subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
      & jmin,jmax,iren,append,nzin,rscale,cscale)
   ! Output is always in  COO format 
@@ -1021,7 +1262,7 @@ subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
 
   call psb_erractionsave(err_act)
   info = 0
-
+  
   if (present(jmin)) then
     jmin_ = jmin
   else
@@ -1033,7 +1274,10 @@ subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
     jmax_ = a%get_ncols()
   endif
 
-  if ((imax<imin).or.(jmax_<jmin_)) return
+  if ((imax<imin).or.(jmax_<jmin_)) then 
+    nz = 0
+    return
+  end if
 
   if (present(append)) then
     append_=append
@@ -1120,7 +1364,6 @@ contains
     irw = imin
     lrw = min(imax,a%get_nrows())
     if (irw<0) then 
-      write(debug_unit,*) ' spgtrow Error : idx no good ',irw
       info = 2
       return
     end if
@@ -1131,13 +1374,14 @@ contains
       nzin_ = 0
     endif
 
-    nzt = a%irp(lrw)-a%irp(irw)
+    nzt = a%irp(lrw+1)-a%irp(irw)
     nz = 0 
 
-    
+
     call psb_ensure_size(nzin_+nzt,ia,info)
     if (info==0) call psb_ensure_size(nzin_+nzt,ja,info)
     if (info==0) call psb_ensure_size(nzin_+nzt,val,info)
+
     if (info /= 0) return
     
     if (present(iren)) then 
@@ -1246,7 +1490,7 @@ contains
     integer, intent(out) :: info
     integer, intent(in), optional  :: gtl(:)
     integer  :: i,ir,ic, ilr, ilc, ip, &
-         & i1,i2,nc,nnz,dupl,ng
+         & i1,i2,nr,nc,nnz,dupl,ng
     integer              :: debug_level, debug_unit
     character(len=20)    :: name='d_cxx_srch_upd'
 
@@ -1264,6 +1508,8 @@ contains
     ilr = -1 
     ilc = -1 
     nnz = a%get_nzeros()
+    nr  = a%get_nrows()
+    nc  = a%get_ncols()
 
     if (present(gtl)) then 
       ng = size(gtl)
@@ -1281,7 +1527,7 @@ contains
           if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
             ir = gtl(ir)
             ic = gtl(ic)
-            if ((ir > 0).and.(ir <= a%m)) then 
+            if ((ir > 0).and.(ir <= nr)) then 
               i1 = a%irp(ir)
               i2 = a%irp(ir+1)
               nc=i2-i1
@@ -1317,7 +1563,7 @@ contains
           if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then 
             ir = gtl(ir)
             ic = gtl(ic)
-            if ((ir > 0).and.(ir <= a%m)) then 
+            if ((ir > 0).and.(ir <= nr)) then 
               i1 = a%irp(ir)
               i2 = a%irp(ir+1)
               nc = i2-i1
@@ -1361,7 +1607,7 @@ contains
           ir = ia(i)
           ic = ja(i) 
 
-          if ((ir > 0).and.(ir <= a%m)) then 
+          if ((ir > 0).and.(ir <= nr)) then 
 
             i1 = a%irp(ir)
             i2 = a%irp(ir+1)
@@ -1394,7 +1640,7 @@ contains
         do i=1, nz
           ir = ia(i)
           ic = ja(i) 
-          if ((ir > 0).and.(ir <= a%m)) then 
+          if ((ir > 0).and.(ir <= nr)) then 
             i1 = a%irp(ir)
             i2 = a%irp(ir+1)
             nc = i2-i1
@@ -1480,7 +1726,8 @@ subroutine d_cp_cxx_to_coo_impl(a,b,info)
   nc  = a%get_ncols()
   nza = a%get_nzeros()
 
-  call b%allocate(nr,nr,nza)
+  call b%allocate(nr,nc,nza)
+  call b%psb_d_base_sparse_mat%cp_from(a%psb_d_base_sparse_mat)
 
   do i=1, nr
     do j=a%irp(i),a%irp(i+1)-1
@@ -1489,14 +1736,7 @@ subroutine d_cp_cxx_to_coo_impl(a,b,info)
       b%val(j) = a%val(j)
     end do
   end do
-  
   call b%set_nzeros(a%get_nzeros())
-  call b%set_nrows(a%get_nrows())
-  call b%set_ncols(a%get_ncols())
-  call b%set_dupl(a%get_dupl())
-  call b%set_state(a%get_state())
-  call b%set_triangle(a%is_triangle())
-  call b%set_upper(a%is_upper())
   call b%fix(info)
 
 
@@ -1528,14 +1768,8 @@ subroutine d_mv_cxx_to_coo_impl(a,b,info)
   nc  = a%get_ncols()
   nza = a%get_nzeros()
 
+  call b%psb_d_base_sparse_mat%mv_from(a%psb_d_base_sparse_mat)
   call b%set_nzeros(a%get_nzeros())
-  call b%set_nrows(a%get_nrows())
-  call b%set_ncols(a%get_ncols())
-  call b%set_dupl(a%get_dupl())
-  call b%set_state(a%get_state())
-  call b%set_triangle(a%is_triangle())
-  call b%set_upper(a%is_upper())
-
   call move_alloc(a%ja,b%ja)
   call move_alloc(a%val,b%val)
   call psb_realloc(nza,b%ia,info)
@@ -1580,13 +1814,9 @@ subroutine d_mv_cxx_from_coo_impl(a,b,info)
   nr  = b%get_nrows()
   nc  = b%get_ncols()
   nza = b%get_nzeros()
+  
+  call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat)
 
-  call a%set_nrows(b%get_nrows())
-  call a%set_ncols(b%get_ncols())
-  call a%set_dupl(b%get_dupl())
-  call a%set_state(b%get_state())
-  call a%set_triangle(b%is_triangle())
-  call a%set_upper(b%is_upper())
   ! Dirty trick: call move_alloc to have the new data allocated just once.
   call move_alloc(b%ia,itemp)
   call move_alloc(b%ja,a%ja)
@@ -1669,8 +1899,16 @@ subroutine d_mv_cxx_to_fmt_impl(a,b,info)
   info = 0
 
   select type (b)
-  class is (psb_d_coo_sparse_mat) 
+  type is (psb_d_coo_sparse_mat) 
     call a%mv_to_coo(b,info)
+    ! Need to fix trivial copies! 
+  type is (psb_d_cxx_sparse_mat) 
+    call b%psb_d_base_sparse_mat%mv_from(a%psb_d_base_sparse_mat)
+    call move_alloc(a%irp, b%irp)
+    call move_alloc(a%ja,  b%ja)
+    call move_alloc(a%val, b%val)
+    call a%free()
+    
   class default
     call tmp%mv_from_fmt(a,info)
     if (info == 0) call b%mv_from_coo(tmp,info)
@@ -1702,8 +1940,15 @@ subroutine d_cp_cxx_to_fmt_impl(a,b,info)
 
 
   select type (b)
-  class is (psb_d_coo_sparse_mat) 
+  type is (psb_d_coo_sparse_mat) 
     call a%cp_to_coo(b,info)
+
+  type is (psb_d_cxx_sparse_mat) 
+    call b%psb_d_base_sparse_mat%cp_from(a%psb_d_base_sparse_mat)
+    b%irp = a%irp
+    b%ja  = a%ja
+    b%val = a%val
+
   class default
     call tmp%cp_from_fmt(a,info)
     if (info == 0) call b%mv_from_coo(tmp,info)
@@ -1734,8 +1979,16 @@ subroutine d_mv_cxx_from_fmt_impl(a,b,info)
   info = 0
 
   select type (b)
-  class is (psb_d_coo_sparse_mat) 
+  type is (psb_d_coo_sparse_mat) 
     call a%mv_from_coo(b,info)
+
+  type is (psb_d_cxx_sparse_mat) 
+    call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat)
+    call move_alloc(b%irp, a%irp)
+    call move_alloc(b%ja,  a%ja)
+    call move_alloc(b%val, a%val)
+    call b%free()
+
   class default
     call tmp%mv_from_fmt(b,info)
     if (info == 0) call a%mv_from_coo(tmp,info)
@@ -1759,7 +2012,7 @@ subroutine d_cp_cxx_from_fmt_impl(a,b,info)
   !locals
   type(psb_d_coo_sparse_mat) :: tmp
   logical             :: rwshr_
-  Integer             :: nza, nr, i,j,irw, idl,err_act, nc
+  Integer             :: nz, nr, i,j,irw, idl,err_act, nc
   Integer, Parameter  :: maxtry=8
   integer              :: debug_level, debug_unit
   character(len=20)   :: name
@@ -1767,8 +2020,15 @@ subroutine d_cp_cxx_from_fmt_impl(a,b,info)
   info = 0
 
   select type (b)
-  class is (psb_d_coo_sparse_mat) 
+  type is (psb_d_coo_sparse_mat) 
     call a%cp_from_coo(b,info)
+
+  type is (psb_d_cxx_sparse_mat) 
+    call a%psb_d_base_sparse_mat%cp_from(b%psb_d_base_sparse_mat)
+    a%irp = b%irp
+    a%ja  = b%ja
+    a%val = b%val
+
   class default
     call tmp%cp_from_fmt(b,info)
     if (info == 0) call a%mv_from_coo(tmp,info)
diff --git a/test/serial/psb_d_cxx_mat_mod.f03 b/test/serial/psb_d_cxx_mat_mod.f03
index a978207b..0b9c129b 100644
--- a/test/serial/psb_d_cxx_mat_mod.f03
+++ b/test/serial/psb_d_cxx_mat_mod.f03
@@ -29,12 +29,22 @@ module psb_d_cxx_mat_mod
     procedure, pass(a) :: mv_from_coo => d_mv_cxx_from_coo
     procedure, pass(a) :: mv_to_fmt => d_mv_cxx_to_fmt
     procedure, pass(a) :: mv_from_fmt => d_mv_cxx_from_fmt
+    procedure, pass(a) :: csgetptn => d_cxx_csgetptn
     procedure, pass(a) :: d_csgetrow => d_cxx_csgetrow
+    procedure, pass(a) :: get_nz_row => d_cxx_get_nz_row
     procedure, pass(a) :: get_size => d_cxx_get_size
     procedure, pass(a) :: free => d_cxx_free
     procedure, pass(a) :: trim => d_cxx_trim
     procedure, pass(a) :: print => d_cxx_print
+    procedure, pass(a) :: sizeof => d_cxx_sizeof
+    procedure, pass(a) :: reinit => d_cxx_reinit
+    procedure, pass(a) :: d_cxx_cp_from
+    generic, public    :: cp_from => d_cxx_cp_from
+    procedure, pass(a) :: d_cxx_mv_from
+    generic, public    :: mv_from => d_cxx_mv_from
+
   end type psb_d_cxx_sparse_mat
+
   private :: d_cxx_get_nzeros, d_cxx_csmm, d_cxx_csmv, d_cxx_cssm, d_cxx_cssv, &
        & d_cxx_csput, d_cxx_reallocate_nz, d_cxx_allocate_mnnz, &
        & d_cxx_free,  d_cxx_print, d_cxx_get_fmt, d_cxx_csnmi, get_diag, &
@@ -42,7 +52,8 @@ module psb_d_cxx_mat_mod
        & d_mv_cxx_to_coo, d_mv_cxx_from_coo, &
        & d_cp_cxx_to_fmt, d_cp_cxx_from_fmt, &
        & d_mv_cxx_to_fmt, d_mv_cxx_from_fmt, &
-       & d_cxx_scals, d_cxx_scal, d_cxx_trim, d_cxx_csgetrow, d_cxx_get_size
+       & d_cxx_scals, d_cxx_scal, d_cxx_trim, d_cxx_csgetrow, d_cxx_get_size, &
+       & d_cxx_sizeof, d_cxx_csgetptn, d_cxx_get_nz_row, d_cxx_reinit
 
 
   interface 
@@ -147,6 +158,25 @@ module psb_d_cxx_mat_mod
     end subroutine d_cxx_csput_impl
   end interface
 
+  interface 
+    subroutine d_cxx_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      use psb_const_mod
+      import psb_d_cxx_sparse_mat
+      implicit none
+      
+      class(psb_d_cxx_sparse_mat), intent(in) :: a
+      integer, intent(in)                  :: imin,imax
+      integer, intent(out)                 :: nz
+      integer, allocatable, intent(inout)  :: ia(:), ja(:)
+      integer,intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer, intent(in), optional        :: iren(:)
+      integer, intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine d_cxx_csgetptn_impl
+  end interface
+
   interface 
     subroutine d_cxx_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,&
          & jmin,jmax,iren,append,nzin,rscale,cscale)
@@ -234,6 +264,18 @@ contains
   !
   !=====================================
 
+  
+  function d_cxx_sizeof(a) result(res)
+    implicit none 
+    class(psb_d_cxx_sparse_mat), intent(in) :: a
+    integer(psb_long_int_k_) :: res
+    res = 8 
+    res = res + psb_sizeof_dp  * size(a%val)
+    res = res + psb_sizeof_int * size(a%irp)
+    res = res + psb_sizeof_int * size(a%ja)
+      
+  end function d_cxx_sizeof
+
   function d_cxx_get_fmt(a) result(res)
     implicit none 
     class(psb_d_cxx_sparse_mat), intent(in) :: a
@@ -245,7 +287,7 @@ contains
     implicit none 
     class(psb_d_cxx_sparse_mat), intent(in) :: a
     integer :: res
-    res = a%irp(a%m+1)-1
+    res = a%irp(a%get_nrows()+1)-1
   end function d_cxx_get_nzeros
 
   function d_cxx_get_size(a) result(res)
@@ -272,6 +314,26 @@ contains
 
   end function d_cxx_get_size
 
+
+
+  function  d_cxx_get_nz_row(idx,a) result(res)
+    use psb_const_mod
+    implicit none
+    
+    class(psb_d_cxx_sparse_mat), intent(in) :: a
+    integer, intent(in)                  :: idx
+    integer                              :: res
+    
+    res = 0 
+ 
+    if ((1<=idx).and.(idx<=a%get_nrows())) then 
+      res = a%irp(idx+1)-a%irp(idx)
+    end if
+    
+  end function d_cxx_get_nz_row
+
+
+
   !=====================================
   !
   !
@@ -299,7 +361,8 @@ contains
 
     call psb_realloc(nz,a%ja,info)
     if (info == 0) call psb_realloc(nz,a%val,info)
-    if (info == 0) call psb_realloc(max(nz,a%m+1,a%n+1),a%irp,info)
+    if (info == 0) call psb_realloc(&
+         & max(nz,a%get_nrows()+1,a%get_ncols()+1),a%irp,info)
     if (info /= 0) then 
       call psb_errpush(4000,name)
       goto 9999
@@ -382,6 +445,49 @@ contains
     return
   end subroutine d_cxx_csput
 
+  subroutine d_cxx_csgetptn(imin,imax,a,nz,ia,ja,info,&
+       & jmin,jmax,iren,append,nzin,rscale,cscale)
+    ! Output is always in  COO format 
+    use psb_error_mod
+    use psb_const_mod
+    implicit none
+    
+    class(psb_d_cxx_sparse_mat), intent(in) :: a
+    integer, intent(in)                  :: imin,imax
+    integer, intent(out)                 :: nz
+    integer, allocatable, intent(inout)  :: ia(:), ja(:)
+    integer,intent(out)                  :: info
+    logical, intent(in), optional        :: append
+    integer, intent(in), optional        :: iren(:)
+    integer, intent(in), optional        :: jmin,jmax, nzin
+    logical, intent(in), optional        :: rscale,cscale
+    Integer :: err_act
+    character(len=20)  :: name='csget'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 0
+
+    call d_cxx_csgetptn_impl(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+
+    if (info /= 0) goto 9999
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+
+  end subroutine d_cxx_csgetptn
+
+
   subroutine d_cxx_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
        & jmin,jmax,iren,append,nzin,rscale,cscale)
     ! Output is always in  COO format 
@@ -591,6 +697,55 @@ contains
 
   end subroutine d_cxx_free
 
+  subroutine d_cxx_reinit(a,clear)
+    use psb_error_mod
+    implicit none 
+
+    class(psb_d_cxx_sparse_mat), intent(inout) :: a   
+    logical, intent(in), optional :: clear
+
+    Integer :: err_act, info
+    character(len=20)  :: name='reinit'
+    logical  :: clear_
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 0
+
+
+    if (present(clear)) then 
+      clear_ = clear
+    else
+      clear_ = .true.
+    end if
+
+    if (a%is_bld() .or. a%is_upd()) then 
+      ! do nothing
+      return
+    else if (a%is_asb()) then 
+      if (clear_) a%val(:) = dzero
+      call a%set_upd()
+    else
+      info = 1121
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+
+  end subroutine d_cxx_reinit
+
+
   subroutine  d_cxx_trim(a)
     use psb_realloc_mod
     use psb_error_mod
@@ -927,6 +1082,7 @@ contains
       call a%set_ncols(n)
       call a%set_bld()
       call a%set_triangle(.false.)
+      call a%set_unit(.false.)
     end if
 
     call psb_erractionrestore(err_act)
@@ -1025,6 +1181,81 @@ contains
   end subroutine d_cxx_print
 
 
+  subroutine d_cxx_cp_from(a,b)
+    use psb_error_mod
+    implicit none 
+
+    class(psb_d_cxx_sparse_mat), intent(out) :: a
+    type(psb_d_cxx_sparse_mat), intent(in)   :: b
+
+
+    Integer :: err_act, info
+    character(len=20)  :: name='cp_from'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+
+    info = 0
+
+    call a%allocate(b%get_nrows(),b%get_ncols(),b%get_nzeros())
+    call a%psb_d_base_sparse_mat%cp_from(b%psb_d_base_sparse_mat)
+    a%irp = b%irp 
+    a%ja  = b%ja
+    a%val = b%val 
+
+    if (info /= 0) goto 9999
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    call psb_errpush(info,name)
+
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine d_cxx_cp_from
+
+  subroutine d_cxx_mv_from(a,b)
+    use psb_error_mod
+    implicit none 
+
+    class(psb_d_cxx_sparse_mat), intent(out)  :: a
+    type(psb_d_cxx_sparse_mat), intent(inout) :: b
+    
+
+    Integer :: err_act, info
+    character(len=20)  :: name='mv_from'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 0
+    call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat)
+    call move_alloc(b%irp, a%irp)
+    call move_alloc(b%ja,  a%ja)
+    call move_alloc(b%val, a%val)
+    call b%free()
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+
+    call psb_errpush(info,name)
+
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine d_cxx_mv_from
+
+
+
   !=====================================
   !
   !
@@ -1151,7 +1382,6 @@ contains
 
     
     if (.not. (a%is_triangle())) then 
-      write(0,*) 'Called SM on a non-triangular mat!'
       info = 1121
       call psb_errpush(info,name)
       goto 9999
@@ -1205,7 +1435,6 @@ contains
 
     
     if (.not. (a%is_triangle())) then 
-      write(0,*) 'Called SM on a non-triangular mat!'
       info = 1121
       call psb_errpush(info,name)
       goto 9999