Implementations for @LX@ spmat.

base/comm/psb_cspgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_c_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_csp_allgather
+
+
+subroutine  psb_lcsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_cspmat_type), intent(inout) :: loca
+  type(psb_lcspmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_lc_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_c_spk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_c_spk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_lcsp_allgather

base/comm/psb_dspgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_d_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_dsp_allgather
+
+
+subroutine  psb_ldsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_dspmat_type), intent(inout) :: loca
+  type(psb_ldspmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_ld_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_r_dpk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_r_dpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_ldsp_allgather

base/comm/psb_ispgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_i_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_isp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_isp_allgather
+
+
+subroutine  psb_@LX@sp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_ispmat_type), intent(inout) :: loca
+  type(psb_@LX@spmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_@LX@_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_ipk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_ipk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_@LX@sp_allgather

base/comm/psb_lspgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_l_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_lsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_lsp_allgather
+
+
+subroutine  psb_@LX@sp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_lspmat_type), intent(inout) :: loca
+  type(psb_@LX@spmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_@LX@_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_lpk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_@LX@sp_allgather

base/comm/psb_sspgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_s_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_ssp_allgather
+
+
+subroutine  psb_lssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_sspmat_type), intent(inout) :: loca
+  type(psb_lsspmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_ls_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_r_spk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_r_spk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_lssp_allgather

base/comm/psb_zspgather.F90

@@ -51,13 +51,14 @@ subroutine  psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   logical, intent(in), optional   :: keepnum,keeploc
 
   type(psb_z_coo_sparse_mat)      :: loc_coo, glob_coo
-  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
+  integer(psb_ipk_) :: nrg, ncg, nzg, nzl
-  integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
-  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k
   logical :: keepnum_, keeploc_
   integer(psb_mpk_) :: ictxt,np,me
   integer(psb_mpk_) :: icomm, minfo, ndx
   integer(psb_mpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_lpk_), allocatable :: locia(:), locja(:), glbia(:), glbja(:)
   integer(psb_ipk_) :: ierr(5)
   character(len=20) :: name
   integer(psb_ipk_) :: debug_level, debug_unit
@@ -169,3 +170,136 @@ subroutine  psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
   return
  
 end subroutine psb_zsp_allgather
+
+
+subroutine  psb_lzsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+  use psb_desc_mod
+  use psb_error_mod
+  use psb_penv_mod
+  use psb_mat_mod
+  use psb_tools_mod
+#ifdef MPI_MOD
+  use mpi
+#endif
+  implicit none
+#ifdef MPI_H
+  include 'mpif.h'
+#endif
+  type(psb_zspmat_type), intent(inout) :: loca
+  type(psb_lzspmat_type), intent(inout) :: globa
+  type(psb_desc_type), intent(in) :: desc_a
+  integer(psb_ipk_), intent(out)            :: info
+  integer(psb_ipk_), intent(in), optional   :: root, dupl
+  logical, intent(in), optional   :: keepnum,keeploc
+
+  type(psb_lz_coo_sparse_mat)     :: loc_coo, glob_coo
+  integer(psb_lpk_) :: nrg, ncg, nzg
+  integer(psb_ipk_) :: err_act, dupl_
+  integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
+  logical :: keepnum_, keeploc_
+  integer(psb_mpk_) :: ictxt,np,me
+  integer(psb_mpk_) :: icomm, minfo, ndx
+  integer(psb_lpk_), allocatable :: nzbr(:), idisp(:)
+  integer(psb_ipk_) :: ierr(5)
+  character(len=20) :: name
+  integer(psb_ipk_) :: debug_level, debug_unit
+
+  name='psb_gather'
+  if (psb_get_errstatus().ne.0) return 
+  info=psb_success_
+
+  call psb_erractionsave(err_act)
+  ictxt = desc_a%get_context()
+  icomm = desc_a%get_mpic()
+  call psb_info(ictxt, me, np)
+
+  if (present(keepnum)) then 
+    keepnum_ = keepnum
+  else
+    keepnum_ = .true.
+  end if
+  if (present(keeploc)) then 
+    keeploc_ = keeploc
+  else
+    keeploc_ = .true.
+  end if
+  call globa%free()
+
+  if (keepnum_) then 
+    nrg = desc_a%get_global_rows()
+    ncg = desc_a%get_global_rows()
+
+    allocate(nzbr(np), idisp(np),stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_request_
+      ierr(1) = 2*np
+      call psb_errpush(info,name,i_err=ierr,a_err='integer')
+      goto 9999      
+    end if
+
+    
+    if (keeploc_) then
+      call loca%cp_to(loc_coo)
+    else
+      call loca%mv_to(loc_coo)
+    end if
+    nzl = loc_coo%get_nzeros()
+    call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
+    call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I')
+    nzbr(:) = 0
+    nzbr(me+1) = nzl
+    call psb_sum(ictxt,nzbr(1:np))
+    nzg = sum(nzbr)
+    if (info == psb_success_) call glob_coo%allocate(nrg,ncg,nzg)
+    if (info /= psb_success_) goto 9999
+
+    do ip=1,np
+      idisp(ip) = sum(nzbr(1:ip-1))
+    enddo
+    ndx = nzbr(me+1) 
+    call mpi_allgatherv(loc_coo%val,ndx,psb_mpi_c_dpk_,&
+         & glob_coo%val,nzbr,idisp,&
+         & psb_mpi_c_dpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ia,ndx,psb_mpi_lpk_,&
+         & glob_coo%ia,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    if (minfo == psb_success_) call &
+         & mpi_allgatherv(loc_coo%ja,ndx,psb_mpi_lpk_,&
+         & glob_coo%ja,nzbr,idisp,&
+         & psb_mpi_lpk_,icomm,minfo)
+    
+    if (minfo /= psb_success_) then 
+      info  = minfo
+      call psb_errpush(psb_err_internal_error_,name,a_err=' from mpi_allgatherv')
+      goto 9999
+    end if    
+    call loc_coo%free()
+    !
+    ! Is the code below safe? For very large cases
+    ! the indices in glob_coo will overflow. But then,
+    ! for very large cases it does not make sense to
+    ! gather the matrix on a single procecss anyway...
+    !
+    call glob_coo%set_nzeros(nzg)
+    if (present(dupl)) call glob_coo%set_dupl(dupl)
+    call globa%mv_from(glob_coo)
+
+  else
+    write(psb_err_unit,*) 'SP_ALLGATHER: Not implemented yet with keepnum ',keepnum_
+    info = -1
+    goto 9999
+  end if
+
+
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 continue
+  call psb_errpush(info,name)
+  call psb_error_handler(ione*ictxt,err_act)
+
+  return
+ 
+end subroutine psb_lzsp_allgather

base/modules/comm/psb_c_comm_mod.f90

@@ -31,7 +31,7 @@
 !    
 module psb_c_comm_mod
   use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_spk_
-  use psb_mat_mod, only  : psb_cspmat_type
+  use psb_mat_mod, only  : psb_cspmat_type, psb_lcspmat_type
   
   use psb_c_vect_mod, only : psb_c_vect_type, psb_c_base_vect_type
   use psb_c_multivect_mod, only : psb_c_multivect_type, psb_c_base_multivect_type
@@ -105,6 +105,16 @@ module psb_c_comm_mod
       integer(psb_ipk_), intent(in), optional   :: root,dupl
       logical, intent(in), optional   :: keepnum,keeploc
     end subroutine psb_csp_allgather
+    subroutine psb_lcsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+      import
+      implicit none
+      type(psb_cspmat_type), intent(inout) :: loca
+      type(psb_lcspmat_type), intent(out)   :: globa
+      type(psb_desc_type), intent(in) :: desc_a
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_ipk_), intent(in), optional   :: root,dupl
+      logical, intent(in), optional   :: keepnum,keeploc
+    end subroutine psb_lcsp_allgather
     subroutine psb_cgather_vect(globx, locx, desc_a, info, root)
       import
       implicit none

base/modules/comm/psb_d_comm_mod.f90

@@ -31,7 +31,7 @@
 !    
 module psb_d_comm_mod
   use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_dpk_
-  use psb_mat_mod, only  : psb_dspmat_type
+  use psb_mat_mod, only  : psb_dspmat_type, psb_ldspmat_type
   
   use psb_d_vect_mod, only : psb_d_vect_type, psb_d_base_vect_type
   use psb_d_multivect_mod, only : psb_d_multivect_type, psb_d_base_multivect_type
@@ -105,6 +105,16 @@ module psb_d_comm_mod
       integer(psb_ipk_), intent(in), optional   :: root,dupl
       logical, intent(in), optional   :: keepnum,keeploc
     end subroutine psb_dsp_allgather
+    subroutine psb_ldsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+      import
+      implicit none
+      type(psb_dspmat_type), intent(inout) :: loca
+      type(psb_ldspmat_type), intent(out)   :: globa
+      type(psb_desc_type), intent(in) :: desc_a
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_ipk_), intent(in), optional   :: root,dupl
+      logical, intent(in), optional   :: keepnum,keeploc
+    end subroutine psb_ldsp_allgather
     subroutine psb_dgather_vect(globx, locx, desc_a, info, root)
       import
       implicit none

base/modules/comm/psb_s_comm_mod.f90

@@ -31,7 +31,7 @@
 !    
 module psb_s_comm_mod
   use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_spk_
-  use psb_mat_mod, only  : psb_sspmat_type
+  use psb_mat_mod, only  : psb_sspmat_type, psb_lsspmat_type
   
   use psb_s_vect_mod, only : psb_s_vect_type, psb_s_base_vect_type
   use psb_s_multivect_mod, only : psb_s_multivect_type, psb_s_base_multivect_type
@@ -105,6 +105,16 @@ module psb_s_comm_mod
       integer(psb_ipk_), intent(in), optional   :: root,dupl
       logical, intent(in), optional   :: keepnum,keeploc
     end subroutine psb_ssp_allgather
+    subroutine psb_lssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+      import
+      implicit none
+      type(psb_sspmat_type), intent(inout) :: loca
+      type(psb_lsspmat_type), intent(out)   :: globa
+      type(psb_desc_type), intent(in) :: desc_a
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_ipk_), intent(in), optional   :: root,dupl
+      logical, intent(in), optional   :: keepnum,keeploc
+    end subroutine psb_lssp_allgather
     subroutine psb_sgather_vect(globx, locx, desc_a, info, root)
       import
       implicit none

base/modules/comm/psb_z_comm_mod.f90

@@ -31,7 +31,7 @@
 !    
 module psb_z_comm_mod
   use psb_desc_mod, only : psb_desc_type, psb_ipk_, psb_dpk_
-  use psb_mat_mod, only  : psb_zspmat_type
+  use psb_mat_mod, only  : psb_zspmat_type, psb_lzspmat_type
   
   use psb_z_vect_mod, only : psb_z_vect_type, psb_z_base_vect_type
   use psb_z_multivect_mod, only : psb_z_multivect_type, psb_z_base_multivect_type
@@ -105,6 +105,16 @@ module psb_z_comm_mod
       integer(psb_ipk_), intent(in), optional   :: root,dupl
       logical, intent(in), optional   :: keepnum,keeploc
     end subroutine psb_zsp_allgather
+    subroutine psb_lzsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc)
+      import
+      implicit none
+      type(psb_zspmat_type), intent(inout) :: loca
+      type(psb_lzspmat_type), intent(out)   :: globa
+      type(psb_desc_type), intent(in) :: desc_a
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_ipk_), intent(in), optional   :: root,dupl
+      logical, intent(in), optional   :: keepnum,keeploc
+    end subroutine psb_lzsp_allgather
     subroutine psb_zgather_vect(globx, locx, desc_a, info, root)
       import
       implicit none

base/modules/serial/psb_c_mat_mod.F90

@@ -194,7 +194,22 @@ module psb_c_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_c_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_cspmat_clone
-
+    !
+    ! To/from lc
+    !
+    procedure, pass(a) :: mv_from_lb  => psb_c_mv_from_lb
+    procedure, pass(a) :: mv_to_lb    => psb_c_mv_to_lb
+    procedure, pass(a) :: cp_from_lb  => psb_c_cp_from_lb
+    procedure, pass(a) :: cp_to_lb    => psb_c_cp_to_lb
+    procedure, pass(a) :: mv_from_l   => psb_c_mv_from_l 
+    procedure, pass(a) :: mv_to_l     => psb_c_mv_to_l 
+    procedure, pass(a) :: cp_from_l   => psb_c_cp_from_l 
+    procedure, pass(a) :: cp_to_l     => psb_c_cp_to_l 
+    generic, public    :: mv_from     => mv_from_lb, mv_from_l
+    generic, public    :: mv_to       => mv_to_lb, mv_to_l
+    generic, public    :: cp_from     => cp_from_lb, cp_from_l
+    generic, public    :: cp_to       => cp_to_lb, cp_to_l
+    
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_c_get_diag
     procedure, pass(a) :: maxval   => psb_c_maxval
@@ -355,6 +370,21 @@ module psb_c_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_lc_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_lcspmat_clone
+    !
+    ! To/from c
+    !
+    procedure, pass(a) :: mv_from_ib  => psb_lc_mv_from_ib
+    procedure, pass(a) :: mv_to_ib    => psb_lc_mv_to_ib
+    procedure, pass(a) :: cp_from_ib  => psb_lc_cp_from_ib
+    procedure, pass(a) :: cp_to_ib    => psb_lc_cp_to_ib
+    procedure, pass(a) :: mv_from_i   => psb_lc_mv_from_i 
+    procedure, pass(a) :: mv_to_i     => psb_lc_mv_to_i 
+    procedure, pass(a) :: cp_from_i   => psb_lc_cp_from_i 
+    procedure, pass(a) :: cp_to_i     => psb_lc_cp_to_i 
+    generic, public    :: mv_from     => mv_from_ib, mv_from_i
+    generic, public    :: mv_to       => mv_to_ib, mv_to_i
+    generic, public    :: cp_from     => cp_from_ib, cp_from_i
+    generic, public    :: cp_to       => cp_to_ib, cp_to_i
 
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_lc_get_diag
@@ -845,7 +875,73 @@ module psb_c_mat_mod
       class(psb_c_base_sparse_mat), intent(inout) :: b
     end subroutine psb_c_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_c_mv_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat
+      class(psb_cspmat_type), intent(inout) :: a
+      class(psb_lc_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_c_mv_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_c_cp_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat
+      class(psb_cspmat_type), intent(out) :: a
+      class(psb_lc_base_sparse_mat), intent(in) :: b
+    end subroutine psb_c_cp_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_c_mv_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat
+      class(psb_cspmat_type), intent(inout) :: a
+      class(psb_lc_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_c_mv_to_lb
+  end interface
   
+  interface 
+    subroutine psb_c_cp_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lc_base_sparse_mat    
+      class(psb_cspmat_type), intent(in) :: a
+      class(psb_lc_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_c_cp_to_lb
+  end interface
+
+    interface 
+    subroutine psb_c_mv_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type
+      class(psb_cspmat_type), intent(inout) :: a
+      class(psb_lcspmat_type), intent(inout) :: b
+    end subroutine psb_c_mv_from_l
+  end interface
+  
+  interface 
+    subroutine psb_c_cp_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type
+      class(psb_cspmat_type), intent(out) :: a
+      class(psb_lcspmat_type), intent(in) :: b
+    end subroutine psb_c_cp_from_l
+  end interface
+  
+  interface 
+    subroutine psb_c_mv_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type
+      class(psb_cspmat_type), intent(inout) :: a
+      class(psb_lcspmat_type), intent(inout) :: b
+    end subroutine psb_c_mv_to_l
+  end interface
+  
+  interface 
+    subroutine psb_c_cp_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_spk_, psb_lcspmat_type
+      class(psb_cspmat_type), intent(in) :: a
+      class(psb_lcspmat_type), intent(inout) :: b
+    end subroutine psb_c_cp_to_l
+  end interface
+
   !
   ! Transfer the internal allocation to the target.
   !  
@@ -1467,6 +1563,73 @@ module psb_c_mat_mod
       class(psb_lc_base_sparse_mat), intent(inout) :: b
     end subroutine psb_lc_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_lc_mv_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat
+      class(psb_lcspmat_type), intent(inout) :: a
+      class(psb_c_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lc_mv_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_lc_cp_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat
+      class(psb_lcspmat_type), intent(out) :: a
+      class(psb_c_base_sparse_mat), intent(in) :: b
+    end subroutine psb_lc_cp_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_lc_mv_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat
+      class(psb_lcspmat_type), intent(inout) :: a
+      class(psb_c_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lc_mv_to_ib
+  end interface
+  
+  interface 
+    subroutine psb_lc_cp_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_c_base_sparse_mat    
+      class(psb_lcspmat_type), intent(in) :: a
+      class(psb_c_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lc_cp_to_ib
+  end interface
+
+    interface 
+    subroutine psb_lc_mv_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type
+      class(psb_lcspmat_type), intent(inout) :: a
+      class(psb_cspmat_type), intent(inout) :: b
+    end subroutine psb_lc_mv_from_i
+  end interface
+  
+  interface 
+    subroutine psb_lc_cp_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type
+      class(psb_lcspmat_type), intent(out) :: a
+      class(psb_cspmat_type), intent(in) :: b
+    end subroutine psb_lc_cp_from_i
+  end interface
+  
+  interface 
+    subroutine psb_lc_mv_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type
+      class(psb_lcspmat_type), intent(inout) :: a
+      class(psb_cspmat_type), intent(inout) :: b
+    end subroutine psb_lc_mv_to_i
+  end interface
+  
+  interface 
+    subroutine psb_lc_cp_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_spk_, psb_cspmat_type
+      class(psb_lcspmat_type), intent(in) :: a
+      class(psb_cspmat_type), intent(inout) :: b
+    end subroutine psb_lc_cp_to_i
+  end interface
+
   
   !
   ! Transfer the internal allocation to the target.

base/modules/serial/psb_c_serial_mod.f90

@@ -204,14 +204,14 @@ module psb_c_serial_mod
   end interface psb_aspxpby
 
   interface psb_spspmm
-!!$    subroutine psb_cspspmm(a,b,c,info)
+    subroutine psb_lcspspmm(a,b,c,info)
-!!$      use psb_c_mat_mod, only : psb_cspmat_type
+      use psb_c_mat_mod, only : psb_lcspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_cspmat_type), intent(in)  :: a,b
+      type(psb_lcspmat_type), intent(in)  :: a,b
-!!$      type(psb_cspmat_type), intent(out) :: c
+      type(psb_lcspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_cspspmm
+    end subroutine psb_lcspspmm
     subroutine psb_lccsrspspmm(a,b,c,info)
       use psb_c_mat_mod, only : psb_lc_csr_sparse_mat
       import :: psb_ipk_
@@ -231,14 +231,14 @@ module psb_c_serial_mod
   end interface psb_spspmm
 
   interface psb_symbmm
-!!$    subroutine psb_csymbmm(a,b,c,info)
+    subroutine psb_lcsymbmm(a,b,c,info)
-!!$      use psb_c_mat_mod, only : psb_cspmat_type
+      use psb_c_mat_mod, only : psb_lcspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_cspmat_type), intent(in)  :: a,b
+      type(psb_lcspmat_type), intent(in)  :: a,b
-!!$      type(psb_cspmat_type), intent(out) :: c
+      type(psb_lcspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_csymbmm
+    end subroutine psb_lcsymbmm
     subroutine psb_lcbase_symbmm(a,b,c,info)
       use psb_c_mat_mod, only : psb_lc_base_sparse_mat, psb_lc_csr_sparse_mat
       import :: psb_ipk_
@@ -250,13 +250,13 @@ module psb_c_serial_mod
   end interface psb_symbmm
 
   interface psb_numbmm
-!!$    subroutine psb_cnumbmm(a,b,c)
+    subroutine psb_lcnumbmm(a,b,c)
-!!$      use psb_c_mat_mod, only : psb_cspmat_type
+      use psb_c_mat_mod, only : psb_lcspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_cspmat_type), intent(in) :: a,b
+      type(psb_lcspmat_type), intent(in) :: a,b
-!!$      type(psb_cspmat_type), intent(inout)  :: c
+      type(psb_lcspmat_type), intent(inout)  :: c
-!!$    end subroutine psb_cnumbmm
+    end subroutine psb_lcnumbmm
     subroutine psb_lcbase_numbmm(a,b,c)
       use psb_c_mat_mod, only : psb_lc_base_sparse_mat, psb_lc_csr_sparse_mat
       import :: psb_ipk_

base/modules/serial/psb_d_mat_mod.F90

@@ -194,7 +194,22 @@ module psb_d_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_d_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_dspmat_clone
-
+    !
+    ! To/from ld
+    !
+    procedure, pass(a) :: mv_from_lb  => psb_d_mv_from_lb
+    procedure, pass(a) :: mv_to_lb    => psb_d_mv_to_lb
+    procedure, pass(a) :: cp_from_lb  => psb_d_cp_from_lb
+    procedure, pass(a) :: cp_to_lb    => psb_d_cp_to_lb
+    procedure, pass(a) :: mv_from_l   => psb_d_mv_from_l 
+    procedure, pass(a) :: mv_to_l     => psb_d_mv_to_l 
+    procedure, pass(a) :: cp_from_l   => psb_d_cp_from_l 
+    procedure, pass(a) :: cp_to_l     => psb_d_cp_to_l 
+    generic, public    :: mv_from     => mv_from_lb, mv_from_l
+    generic, public    :: mv_to       => mv_to_lb, mv_to_l
+    generic, public    :: cp_from     => cp_from_lb, cp_from_l
+    generic, public    :: cp_to       => cp_to_lb, cp_to_l
+    
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_d_get_diag
     procedure, pass(a) :: maxval   => psb_d_maxval
@@ -355,6 +370,21 @@ module psb_d_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_ld_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_ldspmat_clone
+    !
+    ! To/from d
+    !
+    procedure, pass(a) :: mv_from_ib  => psb_ld_mv_from_ib
+    procedure, pass(a) :: mv_to_ib    => psb_ld_mv_to_ib
+    procedure, pass(a) :: cp_from_ib  => psb_ld_cp_from_ib
+    procedure, pass(a) :: cp_to_ib    => psb_ld_cp_to_ib
+    procedure, pass(a) :: mv_from_i   => psb_ld_mv_from_i 
+    procedure, pass(a) :: mv_to_i     => psb_ld_mv_to_i 
+    procedure, pass(a) :: cp_from_i   => psb_ld_cp_from_i 
+    procedure, pass(a) :: cp_to_i     => psb_ld_cp_to_i 
+    generic, public    :: mv_from     => mv_from_ib, mv_from_i
+    generic, public    :: mv_to       => mv_to_ib, mv_to_i
+    generic, public    :: cp_from     => cp_from_ib, cp_from_i
+    generic, public    :: cp_to       => cp_to_ib, cp_to_i
 
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_ld_get_diag
@@ -845,7 +875,73 @@ module psb_d_mat_mod
       class(psb_d_base_sparse_mat), intent(inout) :: b
     end subroutine psb_d_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_d_mv_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat
+      class(psb_dspmat_type), intent(inout) :: a
+      class(psb_ld_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_d_mv_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_d_cp_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat
+      class(psb_dspmat_type), intent(out) :: a
+      class(psb_ld_base_sparse_mat), intent(in) :: b
+    end subroutine psb_d_cp_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_d_mv_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat
+      class(psb_dspmat_type), intent(inout) :: a
+      class(psb_ld_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_d_mv_to_lb
+  end interface
   
+  interface 
+    subroutine psb_d_cp_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ld_base_sparse_mat    
+      class(psb_dspmat_type), intent(in) :: a
+      class(psb_ld_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_d_cp_to_lb
+  end interface
+
+    interface 
+    subroutine psb_d_mv_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type
+      class(psb_dspmat_type), intent(inout) :: a
+      class(psb_ldspmat_type), intent(inout) :: b
+    end subroutine psb_d_mv_from_l
+  end interface
+  
+  interface 
+    subroutine psb_d_cp_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type
+      class(psb_dspmat_type), intent(out) :: a
+      class(psb_ldspmat_type), intent(in) :: b
+    end subroutine psb_d_cp_from_l
+  end interface
+  
+  interface 
+    subroutine psb_d_mv_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type
+      class(psb_dspmat_type), intent(inout) :: a
+      class(psb_ldspmat_type), intent(inout) :: b
+    end subroutine psb_d_mv_to_l
+  end interface
+  
+  interface 
+    subroutine psb_d_cp_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_dpk_, psb_ldspmat_type
+      class(psb_dspmat_type), intent(in) :: a
+      class(psb_ldspmat_type), intent(inout) :: b
+    end subroutine psb_d_cp_to_l
+  end interface
+
   !
   ! Transfer the internal allocation to the target.
   !  
@@ -1467,6 +1563,73 @@ module psb_d_mat_mod
       class(psb_ld_base_sparse_mat), intent(inout) :: b
     end subroutine psb_ld_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_ld_mv_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat
+      class(psb_ldspmat_type), intent(inout) :: a
+      class(psb_d_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ld_mv_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_ld_cp_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat
+      class(psb_ldspmat_type), intent(out) :: a
+      class(psb_d_base_sparse_mat), intent(in) :: b
+    end subroutine psb_ld_cp_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_ld_mv_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat
+      class(psb_ldspmat_type), intent(inout) :: a
+      class(psb_d_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ld_mv_to_ib
+  end interface
+  
+  interface 
+    subroutine psb_ld_cp_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_d_base_sparse_mat    
+      class(psb_ldspmat_type), intent(in) :: a
+      class(psb_d_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ld_cp_to_ib
+  end interface
+
+    interface 
+    subroutine psb_ld_mv_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type
+      class(psb_ldspmat_type), intent(inout) :: a
+      class(psb_dspmat_type), intent(inout) :: b
+    end subroutine psb_ld_mv_from_i
+  end interface
+  
+  interface 
+    subroutine psb_ld_cp_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type
+      class(psb_ldspmat_type), intent(out) :: a
+      class(psb_dspmat_type), intent(in) :: b
+    end subroutine psb_ld_cp_from_i
+  end interface
+  
+  interface 
+    subroutine psb_ld_mv_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type
+      class(psb_ldspmat_type), intent(inout) :: a
+      class(psb_dspmat_type), intent(inout) :: b
+    end subroutine psb_ld_mv_to_i
+  end interface
+  
+  interface 
+    subroutine psb_ld_cp_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_dpk_, psb_dspmat_type
+      class(psb_ldspmat_type), intent(in) :: a
+      class(psb_dspmat_type), intent(inout) :: b
+    end subroutine psb_ld_cp_to_i
+  end interface
+
   
   !
   ! Transfer the internal allocation to the target.

base/modules/serial/psb_d_serial_mod.f90

@@ -204,14 +204,14 @@ module psb_d_serial_mod
   end interface psb_aspxpby
 
   interface psb_spspmm
-!!$    subroutine psb_dspspmm(a,b,c,info)
+    subroutine psb_ldspspmm(a,b,c,info)
-!!$      use psb_d_mat_mod, only : psb_dspmat_type
+      use psb_d_mat_mod, only : psb_ldspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_dspmat_type), intent(in)  :: a,b
+      type(psb_ldspmat_type), intent(in)  :: a,b
-!!$      type(psb_dspmat_type), intent(out) :: c
+      type(psb_ldspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_dspspmm
+    end subroutine psb_ldspspmm
     subroutine psb_ldcsrspspmm(a,b,c,info)
       use psb_d_mat_mod, only : psb_ld_csr_sparse_mat
       import :: psb_ipk_
@@ -231,14 +231,14 @@ module psb_d_serial_mod
   end interface psb_spspmm
 
   interface psb_symbmm
-!!$    subroutine psb_dsymbmm(a,b,c,info)
+    subroutine psb_ldsymbmm(a,b,c,info)
-!!$      use psb_d_mat_mod, only : psb_dspmat_type
+      use psb_d_mat_mod, only : psb_ldspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_dspmat_type), intent(in)  :: a,b
+      type(psb_ldspmat_type), intent(in)  :: a,b
-!!$      type(psb_dspmat_type), intent(out) :: c
+      type(psb_ldspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_dsymbmm
+    end subroutine psb_ldsymbmm
     subroutine psb_ldbase_symbmm(a,b,c,info)
       use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat
       import :: psb_ipk_
@@ -250,13 +250,13 @@ module psb_d_serial_mod
   end interface psb_symbmm
 
   interface psb_numbmm
-!!$    subroutine psb_dnumbmm(a,b,c)
+    subroutine psb_ldnumbmm(a,b,c)
-!!$      use psb_d_mat_mod, only : psb_dspmat_type
+      use psb_d_mat_mod, only : psb_ldspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_dspmat_type), intent(in) :: a,b
+      type(psb_ldspmat_type), intent(in) :: a,b
-!!$      type(psb_dspmat_type), intent(inout)  :: c
+      type(psb_ldspmat_type), intent(inout)  :: c
-!!$    end subroutine psb_dnumbmm
+    end subroutine psb_ldnumbmm
     subroutine psb_ldbase_numbmm(a,b,c)
       use psb_d_mat_mod, only : psb_ld_base_sparse_mat, psb_ld_csr_sparse_mat
       import :: psb_ipk_

base/modules/serial/psb_s_mat_mod.F90

@@ -194,7 +194,22 @@ module psb_s_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_s_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_sspmat_clone
-
+    !
+    ! To/from ls
+    !
+    procedure, pass(a) :: mv_from_lb  => psb_s_mv_from_lb
+    procedure, pass(a) :: mv_to_lb    => psb_s_mv_to_lb
+    procedure, pass(a) :: cp_from_lb  => psb_s_cp_from_lb
+    procedure, pass(a) :: cp_to_lb    => psb_s_cp_to_lb
+    procedure, pass(a) :: mv_from_l   => psb_s_mv_from_l 
+    procedure, pass(a) :: mv_to_l     => psb_s_mv_to_l 
+    procedure, pass(a) :: cp_from_l   => psb_s_cp_from_l 
+    procedure, pass(a) :: cp_to_l     => psb_s_cp_to_l 
+    generic, public    :: mv_from     => mv_from_lb, mv_from_l
+    generic, public    :: mv_to       => mv_to_lb, mv_to_l
+    generic, public    :: cp_from     => cp_from_lb, cp_from_l
+    generic, public    :: cp_to       => cp_to_lb, cp_to_l
+    
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_s_get_diag
     procedure, pass(a) :: maxval   => psb_s_maxval
@@ -355,6 +370,21 @@ module psb_s_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_ls_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_lsspmat_clone
+    !
+    ! To/from s
+    !
+    procedure, pass(a) :: mv_from_ib  => psb_ls_mv_from_ib
+    procedure, pass(a) :: mv_to_ib    => psb_ls_mv_to_ib
+    procedure, pass(a) :: cp_from_ib  => psb_ls_cp_from_ib
+    procedure, pass(a) :: cp_to_ib    => psb_ls_cp_to_ib
+    procedure, pass(a) :: mv_from_i   => psb_ls_mv_from_i 
+    procedure, pass(a) :: mv_to_i     => psb_ls_mv_to_i 
+    procedure, pass(a) :: cp_from_i   => psb_ls_cp_from_i 
+    procedure, pass(a) :: cp_to_i     => psb_ls_cp_to_i 
+    generic, public    :: mv_from     => mv_from_ib, mv_from_i
+    generic, public    :: mv_to       => mv_to_ib, mv_to_i
+    generic, public    :: cp_from     => cp_from_ib, cp_from_i
+    generic, public    :: cp_to       => cp_to_ib, cp_to_i
 
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_ls_get_diag
@@ -845,7 +875,73 @@ module psb_s_mat_mod
       class(psb_s_base_sparse_mat), intent(inout) :: b
     end subroutine psb_s_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_s_mv_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat
+      class(psb_sspmat_type), intent(inout) :: a
+      class(psb_ls_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_s_mv_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_s_cp_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat
+      class(psb_sspmat_type), intent(out) :: a
+      class(psb_ls_base_sparse_mat), intent(in) :: b
+    end subroutine psb_s_cp_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_s_mv_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat
+      class(psb_sspmat_type), intent(inout) :: a
+      class(psb_ls_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_s_mv_to_lb
+  end interface
   
+  interface 
+    subroutine psb_s_cp_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_ls_base_sparse_mat    
+      class(psb_sspmat_type), intent(in) :: a
+      class(psb_ls_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_s_cp_to_lb
+  end interface
+
+    interface 
+    subroutine psb_s_mv_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type
+      class(psb_sspmat_type), intent(inout) :: a
+      class(psb_lsspmat_type), intent(inout) :: b
+    end subroutine psb_s_mv_from_l
+  end interface
+  
+  interface 
+    subroutine psb_s_cp_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type
+      class(psb_sspmat_type), intent(out) :: a
+      class(psb_lsspmat_type), intent(in) :: b
+    end subroutine psb_s_cp_from_l
+  end interface
+  
+  interface 
+    subroutine psb_s_mv_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type
+      class(psb_sspmat_type), intent(inout) :: a
+      class(psb_lsspmat_type), intent(inout) :: b
+    end subroutine psb_s_mv_to_l
+  end interface
+  
+  interface 
+    subroutine psb_s_cp_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_spk_, psb_lsspmat_type
+      class(psb_sspmat_type), intent(in) :: a
+      class(psb_lsspmat_type), intent(inout) :: b
+    end subroutine psb_s_cp_to_l
+  end interface
+
   !
   ! Transfer the internal allocation to the target.
   !  
@@ -1467,6 +1563,73 @@ module psb_s_mat_mod
       class(psb_ls_base_sparse_mat), intent(inout) :: b
     end subroutine psb_ls_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_ls_mv_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat
+      class(psb_lsspmat_type), intent(inout) :: a
+      class(psb_s_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ls_mv_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_ls_cp_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat
+      class(psb_lsspmat_type), intent(out) :: a
+      class(psb_s_base_sparse_mat), intent(in) :: b
+    end subroutine psb_ls_cp_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_ls_mv_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat
+      class(psb_lsspmat_type), intent(inout) :: a
+      class(psb_s_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ls_mv_to_ib
+  end interface
+  
+  interface 
+    subroutine psb_ls_cp_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_s_base_sparse_mat    
+      class(psb_lsspmat_type), intent(in) :: a
+      class(psb_s_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_ls_cp_to_ib
+  end interface
+
+    interface 
+    subroutine psb_ls_mv_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type
+      class(psb_lsspmat_type), intent(inout) :: a
+      class(psb_sspmat_type), intent(inout) :: b
+    end subroutine psb_ls_mv_from_i
+  end interface
+  
+  interface 
+    subroutine psb_ls_cp_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type
+      class(psb_lsspmat_type), intent(out) :: a
+      class(psb_sspmat_type), intent(in) :: b
+    end subroutine psb_ls_cp_from_i
+  end interface
+  
+  interface 
+    subroutine psb_ls_mv_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type
+      class(psb_lsspmat_type), intent(inout) :: a
+      class(psb_sspmat_type), intent(inout) :: b
+    end subroutine psb_ls_mv_to_i
+  end interface
+  
+  interface 
+    subroutine psb_ls_cp_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_spk_, psb_sspmat_type
+      class(psb_lsspmat_type), intent(in) :: a
+      class(psb_sspmat_type), intent(inout) :: b
+    end subroutine psb_ls_cp_to_i
+  end interface
+
   
   !
   ! Transfer the internal allocation to the target.

base/modules/serial/psb_s_serial_mod.f90

@@ -204,14 +204,14 @@ module psb_s_serial_mod
   end interface psb_aspxpby
 
   interface psb_spspmm
-!!$    subroutine psb_sspspmm(a,b,c,info)
+    subroutine psb_lsspspmm(a,b,c,info)
-!!$      use psb_s_mat_mod, only : psb_sspmat_type
+      use psb_s_mat_mod, only : psb_lsspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_sspmat_type), intent(in)  :: a,b
+      type(psb_lsspmat_type), intent(in)  :: a,b
-!!$      type(psb_sspmat_type), intent(out) :: c
+      type(psb_lsspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_sspspmm
+    end subroutine psb_lsspspmm
     subroutine psb_lscsrspspmm(a,b,c,info)
       use psb_s_mat_mod, only : psb_ls_csr_sparse_mat
       import :: psb_ipk_
@@ -231,14 +231,14 @@ module psb_s_serial_mod
   end interface psb_spspmm
 
   interface psb_symbmm
-!!$    subroutine psb_ssymbmm(a,b,c,info)
+    subroutine psb_lssymbmm(a,b,c,info)
-!!$      use psb_s_mat_mod, only : psb_sspmat_type
+      use psb_s_mat_mod, only : psb_lsspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_sspmat_type), intent(in)  :: a,b
+      type(psb_lsspmat_type), intent(in)  :: a,b
-!!$      type(psb_sspmat_type), intent(out) :: c
+      type(psb_lsspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_ssymbmm
+    end subroutine psb_lssymbmm
     subroutine psb_lsbase_symbmm(a,b,c,info)
       use psb_s_mat_mod, only : psb_ls_base_sparse_mat, psb_ls_csr_sparse_mat
       import :: psb_ipk_
@@ -250,13 +250,13 @@ module psb_s_serial_mod
   end interface psb_symbmm
 
   interface psb_numbmm
-!!$    subroutine psb_snumbmm(a,b,c)
+    subroutine psb_lsnumbmm(a,b,c)
-!!$      use psb_s_mat_mod, only : psb_sspmat_type
+      use psb_s_mat_mod, only : psb_lsspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_sspmat_type), intent(in) :: a,b
+      type(psb_lsspmat_type), intent(in) :: a,b
-!!$      type(psb_sspmat_type), intent(inout)  :: c
+      type(psb_lsspmat_type), intent(inout)  :: c
-!!$    end subroutine psb_snumbmm
+    end subroutine psb_lsnumbmm
     subroutine psb_lsbase_numbmm(a,b,c)
       use psb_s_mat_mod, only : psb_ls_base_sparse_mat, psb_ls_csr_sparse_mat
       import :: psb_ipk_

base/modules/serial/psb_serial_mod.f90

@@ -68,7 +68,7 @@ module psb_serial_mod
     end subroutine psb_d_nspaxpby
   end interface psb_nspaxpby
 
-  interface symbmm
+  interface 
     subroutine symbmm (n, m, l, ia, ja, diaga, &
          & ib, jb, diagb, ic, jc, diagc, index)
       import :: psb_ipk_
@@ -76,6 +76,13 @@ module psb_serial_mod
            & diagc,  index(*)
       integer(psb_ipk_), allocatable :: ic(:),jc(:)
     end subroutine symbmm
+    subroutine lsymbmm (n, m, l, ia, ja, diaga, &
+         & ib, jb, diagb, ic, jc, diagc, index)
+      import :: psb_ipk_, psb_lpk_
+      integer(psb_lpk_) :: n,m,l,  ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
+           & diagc,  index(*)
+      integer(psb_lpk_), allocatable :: ic(:),jc(:)
+    end subroutine lsymbmm
   end interface
   
 

base/modules/serial/psb_z_mat_mod.F90

@@ -194,7 +194,22 @@ module psb_z_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_z_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_zspmat_clone
-
+    !
+    ! To/from lz
+    !
+    procedure, pass(a) :: mv_from_lb  => psb_z_mv_from_lb
+    procedure, pass(a) :: mv_to_lb    => psb_z_mv_to_lb
+    procedure, pass(a) :: cp_from_lb  => psb_z_cp_from_lb
+    procedure, pass(a) :: cp_to_lb    => psb_z_cp_to_lb
+    procedure, pass(a) :: mv_from_l   => psb_z_mv_from_l 
+    procedure, pass(a) :: mv_to_l     => psb_z_mv_to_l 
+    procedure, pass(a) :: cp_from_l   => psb_z_cp_from_l 
+    procedure, pass(a) :: cp_to_l     => psb_z_cp_to_l 
+    generic, public    :: mv_from     => mv_from_lb, mv_from_l
+    generic, public    :: mv_to       => mv_to_lb, mv_to_l
+    generic, public    :: cp_from     => cp_from_lb, cp_from_l
+    generic, public    :: cp_to       => cp_to_lb, cp_to_l
+    
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_z_get_diag
     procedure, pass(a) :: maxval   => psb_z_maxval
@@ -355,6 +370,21 @@ module psb_z_mat_mod
     procedure, pass(a) :: cscnv_base  => psb_lz_cscnv_base
     generic, public    :: cscnv       => cscnv_np, cscnv_ip, cscnv_base
     procedure, pass(a) :: clone       => psb_lzspmat_clone
+    !
+    ! To/from z
+    !
+    procedure, pass(a) :: mv_from_ib  => psb_lz_mv_from_ib
+    procedure, pass(a) :: mv_to_ib    => psb_lz_mv_to_ib
+    procedure, pass(a) :: cp_from_ib  => psb_lz_cp_from_ib
+    procedure, pass(a) :: cp_to_ib    => psb_lz_cp_to_ib
+    procedure, pass(a) :: mv_from_i   => psb_lz_mv_from_i 
+    procedure, pass(a) :: mv_to_i     => psb_lz_mv_to_i 
+    procedure, pass(a) :: cp_from_i   => psb_lz_cp_from_i 
+    procedure, pass(a) :: cp_to_i     => psb_lz_cp_to_i 
+    generic, public    :: mv_from     => mv_from_ib, mv_from_i
+    generic, public    :: mv_to       => mv_to_ib, mv_to_i
+    generic, public    :: cp_from     => cp_from_ib, cp_from_i
+    generic, public    :: cp_to       => cp_to_ib, cp_to_i
 
     ! Computational routines 
     procedure, pass(a) :: get_diag => psb_lz_get_diag
@@ -845,7 +875,73 @@ module psb_z_mat_mod
       class(psb_z_base_sparse_mat), intent(inout) :: b
     end subroutine psb_z_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_z_mv_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat
+      class(psb_zspmat_type), intent(inout) :: a
+      class(psb_lz_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_z_mv_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_z_cp_from_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat
+      class(psb_zspmat_type), intent(out) :: a
+      class(psb_lz_base_sparse_mat), intent(in) :: b
+    end subroutine psb_z_cp_from_lb
+  end interface
+  
+  interface 
+    subroutine psb_z_mv_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat
+      class(psb_zspmat_type), intent(inout) :: a
+      class(psb_lz_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_z_mv_to_lb
+  end interface
   
+  interface 
+    subroutine psb_z_cp_to_lb(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lz_base_sparse_mat    
+      class(psb_zspmat_type), intent(in) :: a
+      class(psb_lz_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_z_cp_to_lb
+  end interface
+
+    interface 
+    subroutine psb_z_mv_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type
+      class(psb_zspmat_type), intent(inout) :: a
+      class(psb_lzspmat_type), intent(inout) :: b
+    end subroutine psb_z_mv_from_l
+  end interface
+  
+  interface 
+    subroutine psb_z_cp_from_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type
+      class(psb_zspmat_type), intent(out) :: a
+      class(psb_lzspmat_type), intent(in) :: b
+    end subroutine psb_z_cp_from_l
+  end interface
+  
+  interface 
+    subroutine psb_z_mv_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type
+      class(psb_zspmat_type), intent(inout) :: a
+      class(psb_lzspmat_type), intent(inout) :: b
+    end subroutine psb_z_mv_to_l
+  end interface
+  
+  interface 
+    subroutine psb_z_cp_to_l(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_dpk_, psb_lzspmat_type
+      class(psb_zspmat_type), intent(in) :: a
+      class(psb_lzspmat_type), intent(inout) :: b
+    end subroutine psb_z_cp_to_l
+  end interface
+
   !
   ! Transfer the internal allocation to the target.
   !  
@@ -1467,6 +1563,73 @@ module psb_z_mat_mod
       class(psb_lz_base_sparse_mat), intent(inout) :: b
     end subroutine psb_lz_cp_to
   end interface
+  !
+  ! Mixed type conversions
+  !
+  interface 
+    subroutine psb_lz_mv_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat
+      class(psb_lzspmat_type), intent(inout) :: a
+      class(psb_z_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lz_mv_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_lz_cp_from_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat
+      class(psb_lzspmat_type), intent(out) :: a
+      class(psb_z_base_sparse_mat), intent(in) :: b
+    end subroutine psb_lz_cp_from_ib
+  end interface
+  
+  interface 
+    subroutine psb_lz_mv_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat
+      class(psb_lzspmat_type), intent(inout) :: a
+      class(psb_z_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lz_mv_to_ib
+  end interface
+  
+  interface 
+    subroutine psb_lz_cp_to_ib(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_z_base_sparse_mat    
+      class(psb_lzspmat_type), intent(in) :: a
+      class(psb_z_base_sparse_mat), intent(inout) :: b
+    end subroutine psb_lz_cp_to_ib
+  end interface
+
+    interface 
+    subroutine psb_lz_mv_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type
+      class(psb_lzspmat_type), intent(inout) :: a
+      class(psb_zspmat_type), intent(inout) :: b
+    end subroutine psb_lz_mv_from_i
+  end interface
+  
+  interface 
+    subroutine psb_lz_cp_from_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type
+      class(psb_lzspmat_type), intent(out) :: a
+      class(psb_zspmat_type), intent(in) :: b
+    end subroutine psb_lz_cp_from_i
+  end interface
+  
+  interface 
+    subroutine psb_lz_mv_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type
+      class(psb_lzspmat_type), intent(inout) :: a
+      class(psb_zspmat_type), intent(inout) :: b
+    end subroutine psb_lz_mv_to_i
+  end interface
+  
+  interface 
+    subroutine psb_lz_cp_to_i(a,b)
+      import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_dpk_, psb_zspmat_type
+      class(psb_lzspmat_type), intent(in) :: a
+      class(psb_zspmat_type), intent(inout) :: b
+    end subroutine psb_lz_cp_to_i
+  end interface
+
   
   !
   ! Transfer the internal allocation to the target.

base/modules/serial/psb_z_serial_mod.f90

@@ -204,14 +204,14 @@ module psb_z_serial_mod
   end interface psb_aspxpby
 
   interface psb_spspmm
-!!$    subroutine psb_zspspmm(a,b,c,info)
+    subroutine psb_lzspspmm(a,b,c,info)
-!!$      use psb_z_mat_mod, only : psb_zspmat_type
+      use psb_z_mat_mod, only : psb_lzspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_zspmat_type), intent(in)  :: a,b
+      type(psb_lzspmat_type), intent(in)  :: a,b
-!!$      type(psb_zspmat_type), intent(out) :: c
+      type(psb_lzspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_zspspmm
+    end subroutine psb_lzspspmm
     subroutine psb_lzcsrspspmm(a,b,c,info)
       use psb_z_mat_mod, only : psb_lz_csr_sparse_mat
       import :: psb_ipk_
@@ -231,14 +231,14 @@ module psb_z_serial_mod
   end interface psb_spspmm
 
   interface psb_symbmm
-!!$    subroutine psb_zsymbmm(a,b,c,info)
+    subroutine psb_lzsymbmm(a,b,c,info)
-!!$      use psb_z_mat_mod, only : psb_zspmat_type
+      use psb_z_mat_mod, only : psb_lzspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_zspmat_type), intent(in)  :: a,b
+      type(psb_lzspmat_type), intent(in)  :: a,b
-!!$      type(psb_zspmat_type), intent(out) :: c
+      type(psb_lzspmat_type), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)                :: info
+      integer(psb_ipk_), intent(out)                :: info
-!!$    end subroutine psb_zsymbmm
+    end subroutine psb_lzsymbmm
     subroutine psb_lzbase_symbmm(a,b,c,info)
       use psb_z_mat_mod, only : psb_lz_base_sparse_mat, psb_lz_csr_sparse_mat
       import :: psb_ipk_
@@ -250,13 +250,13 @@ module psb_z_serial_mod
   end interface psb_symbmm
 
   interface psb_numbmm
-!!$    subroutine psb_znumbmm(a,b,c)
+    subroutine psb_lznumbmm(a,b,c)
-!!$      use psb_z_mat_mod, only : psb_zspmat_type
+      use psb_z_mat_mod, only : psb_lzspmat_type
-!!$      import :: psb_ipk_
+      import :: psb_ipk_
-!!$      implicit none 
+      implicit none 
-!!$      type(psb_zspmat_type), intent(in) :: a,b
+      type(psb_lzspmat_type), intent(in) :: a,b
-!!$      type(psb_zspmat_type), intent(inout)  :: c
+      type(psb_lzspmat_type), intent(inout)  :: c
-!!$    end subroutine psb_znumbmm
+    end subroutine psb_lznumbmm
     subroutine psb_lzbase_numbmm(a,b,c)
       use psb_z_mat_mod, only : psb_lz_base_sparse_mat, psb_lz_csr_sparse_mat
       import :: psb_ipk_

base/serial/Makefile

@@ -8,7 +8,7 @@ FOBJS = psb_lsame.o psi_m_serial_impl.o psi_e_serial_impl.o \
 	psb_sspspmm.o psb_dspspmm.o psb_cspspmm.o psb_zspspmm.o \
 	psb_ssymbmm.o psb_dsymbmm.o psb_csymbmm.o psb_zsymbmm.o \
 	psb_snumbmm.o psb_dnumbmm.o psb_cnumbmm.o psb_znumbmm.o \
-	smmp.o \
+	smmp.o lsmmp.o \
 	psb_sgeprt.o psb_dgeprt.o psb_cgeprt.o psb_zgeprt.o\
 	psb_spdot_srtd.o psb_aspxpby.o psb_spge_dot.o\
 	psb_sgelp.o psb_dgelp.o psb_cgelp.o psb_zgelp.o \

base/serial/impl/psb_c_mat_impl.F90

@@ -2432,6 +2432,155 @@ subroutine psb_c_scals(d,a,info)
 
 end subroutine psb_c_scals
 
+subroutine psb_c_mv_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_mv_from_lb
+  implicit none 
+  
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_lfmt(b,info)
+  
+end subroutine psb_c_mv_from_lb
+
+  
+subroutine psb_c_cp_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cp_from_lb
+  implicit none 
+  
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_lfmt(b,info)
+
+end subroutine psb_c_cp_from_lb
+
+subroutine psb_c_mv_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_mv_to_lb
+  implicit none 
+  
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_lfmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_c_mv_to_lb
+
+subroutine psb_c_cp_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cp_to_lb
+  implicit none 
+  class(psb_cspmat_type), intent(in) :: a
+  class(psb_lc_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_lfmt(b,info)
+  end if
+
+end subroutine psb_c_cp_to_lb
+
+subroutine psb_c_mv_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_mv_from_l
+  implicit none 
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_lcspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_c_mv_from_l
+
+  
+subroutine psb_c_cp_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cp_from_l
+  implicit none 
+
+  class(psb_cspmat_type), intent(out) :: a
+  class(psb_lcspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_c_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_c_cp_from_l
+
+subroutine psb_c_mv_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_mv_to_l
+  implicit none 
+
+  class(psb_cspmat_type), intent(inout) :: a
+  class(psb_lcspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_lc_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_c_mv_to_l
+
+subroutine psb_c_cp_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_c_cp_to_l
+  implicit none 
+  
+  class(psb_cspmat_type), intent(in) :: a
+  class(psb_lcspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_lc_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_c_cp_to_l
+
+
+!
+!
+! lc versions
+!
+
 
 subroutine  psb_lc_set_nrows(m,a) 
   use psb_c_mat_mod, psb_protect_name => psb_lc_set_nrows
@@ -4331,5 +4480,148 @@ subroutine psb_lc_scals(d,a,info)
 
 end subroutine psb_lc_scals
 
+subroutine psb_lc_mv_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_mv_from_ib
+  implicit none 
+  
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_c_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_ifmt(b,info)
+  
+end subroutine psb_lc_mv_from_ib
+  
+subroutine psb_lc_cp_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cp_from_ib
+  implicit none 
+  
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_c_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_ifmt(b,info)
+
+end subroutine psb_lc_cp_from_ib
+
+subroutine psb_lc_mv_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_mv_to_ib
+  implicit none 
+  
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_c_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_ifmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_lc_mv_to_ib
+
+subroutine psb_lc_cp_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cp_to_ib
+  implicit none 
+  class(psb_lcspmat_type), intent(in) :: a
+  class(psb_c_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_ifmt(b,info)
+  end if
+
+end subroutine psb_lc_cp_to_ib
+
+subroutine psb_lc_mv_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_mv_from_i
+  implicit none 
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_cspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_lc_mv_from_i
+
+  
+subroutine psb_lc_cp_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cp_from_i
+  implicit none 
+
+  class(psb_lcspmat_type), intent(out) :: a
+  class(psb_cspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_lc_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_lc_cp_from_i
+
+subroutine psb_lc_mv_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_mv_to_i
+  implicit none 
+
+  class(psb_lcspmat_type), intent(inout) :: a
+  class(psb_cspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_c_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_lc_mv_to_i
+
+subroutine psb_lc_cp_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_c_mat_mod, psb_protect_name => psb_lc_cp_to_i
+  implicit none 
+  
+  class(psb_lcspmat_type), intent(in) :: a
+  class(psb_cspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_c_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_lc_cp_to_i
 
 
+
+  

base/serial/impl/psb_d_mat_impl.F90

@@ -2432,6 +2432,155 @@ subroutine psb_d_scals(d,a,info)
 
 end subroutine psb_d_scals
 
+subroutine psb_d_mv_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_mv_from_lb
+  implicit none 
+  
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_lfmt(b,info)
+  
+end subroutine psb_d_mv_from_lb
+
+  
+subroutine psb_d_cp_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cp_from_lb
+  implicit none 
+  
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_lfmt(b,info)
+
+end subroutine psb_d_cp_from_lb
+
+subroutine psb_d_mv_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_mv_to_lb
+  implicit none 
+  
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_lfmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_d_mv_to_lb
+
+subroutine psb_d_cp_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cp_to_lb
+  implicit none 
+  class(psb_dspmat_type), intent(in) :: a
+  class(psb_ld_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_lfmt(b,info)
+  end if
+
+end subroutine psb_d_cp_to_lb
+
+subroutine psb_d_mv_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_mv_from_l
+  implicit none 
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_ldspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_d_mv_from_l
+
+  
+subroutine psb_d_cp_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cp_from_l
+  implicit none 
+
+  class(psb_dspmat_type), intent(out) :: a
+  class(psb_ldspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_d_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_d_cp_from_l
+
+subroutine psb_d_mv_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_mv_to_l
+  implicit none 
+
+  class(psb_dspmat_type), intent(inout) :: a
+  class(psb_ldspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_ld_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_d_mv_to_l
+
+subroutine psb_d_cp_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_d_cp_to_l
+  implicit none 
+  
+  class(psb_dspmat_type), intent(in) :: a
+  class(psb_ldspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_ld_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_d_cp_to_l
+
+
+!
+!
+! ld versions
+!
+
 
 subroutine  psb_ld_set_nrows(m,a) 
   use psb_d_mat_mod, psb_protect_name => psb_ld_set_nrows
@@ -4331,5 +4480,148 @@ subroutine psb_ld_scals(d,a,info)
 
 end subroutine psb_ld_scals
 
+subroutine psb_ld_mv_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_mv_from_ib
+  implicit none 
+  
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_ifmt(b,info)
+  
+end subroutine psb_ld_mv_from_ib
+  
+subroutine psb_ld_cp_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cp_from_ib
+  implicit none 
+  
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_ifmt(b,info)
+
+end subroutine psb_ld_cp_from_ib
+
+subroutine psb_ld_mv_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_mv_to_ib
+  implicit none 
+  
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_ifmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_ld_mv_to_ib
+
+subroutine psb_ld_cp_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cp_to_ib
+  implicit none 
+  class(psb_ldspmat_type), intent(in) :: a
+  class(psb_d_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_ifmt(b,info)
+  end if
+
+end subroutine psb_ld_cp_to_ib
+
+subroutine psb_ld_mv_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_mv_from_i
+  implicit none 
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_dspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_ld_mv_from_i
+
+  
+subroutine psb_ld_cp_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cp_from_i
+  implicit none 
+
+  class(psb_ldspmat_type), intent(out) :: a
+  class(psb_dspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_ld_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_ld_cp_from_i
+
+subroutine psb_ld_mv_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_mv_to_i
+  implicit none 
+
+  class(psb_ldspmat_type), intent(inout) :: a
+  class(psb_dspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_d_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_ld_mv_to_i
+
+subroutine psb_ld_cp_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_d_mat_mod, psb_protect_name => psb_ld_cp_to_i
+  implicit none 
+  
+  class(psb_ldspmat_type), intent(in) :: a
+  class(psb_dspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_d_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_ld_cp_to_i
 
 
+
+  

base/serial/impl/psb_s_mat_impl.F90

@@ -2432,6 +2432,155 @@ subroutine psb_s_scals(d,a,info)
 
 end subroutine psb_s_scals
 
+subroutine psb_s_mv_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_mv_from_lb
+  implicit none 
+  
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_lfmt(b,info)
+  
+end subroutine psb_s_mv_from_lb
+
+  
+subroutine psb_s_cp_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cp_from_lb
+  implicit none 
+  
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_lfmt(b,info)
+
+end subroutine psb_s_cp_from_lb
+
+subroutine psb_s_mv_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_mv_to_lb
+  implicit none 
+  
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_lfmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_s_mv_to_lb
+
+subroutine psb_s_cp_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cp_to_lb
+  implicit none 
+  class(psb_sspmat_type), intent(in) :: a
+  class(psb_ls_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_lfmt(b,info)
+  end if
+
+end subroutine psb_s_cp_to_lb
+
+subroutine psb_s_mv_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_mv_from_l
+  implicit none 
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_lsspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_s_mv_from_l
+
+  
+subroutine psb_s_cp_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cp_from_l
+  implicit none 
+
+  class(psb_sspmat_type), intent(out) :: a
+  class(psb_lsspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_s_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_s_cp_from_l
+
+subroutine psb_s_mv_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_mv_to_l
+  implicit none 
+
+  class(psb_sspmat_type), intent(inout) :: a
+  class(psb_lsspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_ls_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_s_mv_to_l
+
+subroutine psb_s_cp_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_s_cp_to_l
+  implicit none 
+  
+  class(psb_sspmat_type), intent(in) :: a
+  class(psb_lsspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_ls_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_s_cp_to_l
+
+
+!
+!
+! ls versions
+!
+
 
 subroutine  psb_ls_set_nrows(m,a) 
   use psb_s_mat_mod, psb_protect_name => psb_ls_set_nrows
@@ -4331,5 +4480,148 @@ subroutine psb_ls_scals(d,a,info)
 
 end subroutine psb_ls_scals
 
+subroutine psb_ls_mv_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_mv_from_ib
+  implicit none 
+  
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_ifmt(b,info)
+  
+end subroutine psb_ls_mv_from_ib
+  
+subroutine psb_ls_cp_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cp_from_ib
+  implicit none 
+  
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_ifmt(b,info)
+
+end subroutine psb_ls_cp_from_ib
+
+subroutine psb_ls_mv_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_mv_to_ib
+  implicit none 
+  
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_ifmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_ls_mv_to_ib
+
+subroutine psb_ls_cp_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cp_to_ib
+  implicit none 
+  class(psb_lsspmat_type), intent(in) :: a
+  class(psb_s_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_ifmt(b,info)
+  end if
+
+end subroutine psb_ls_cp_to_ib
+
+subroutine psb_ls_mv_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_mv_from_i
+  implicit none 
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_sspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_ls_mv_from_i
+
+  
+subroutine psb_ls_cp_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cp_from_i
+  implicit none 
+
+  class(psb_lsspmat_type), intent(out) :: a
+  class(psb_sspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_ls_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_ls_cp_from_i
+
+subroutine psb_ls_mv_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_mv_to_i
+  implicit none 
+
+  class(psb_lsspmat_type), intent(inout) :: a
+  class(psb_sspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_s_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_ls_mv_to_i
+
+subroutine psb_ls_cp_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_s_mat_mod, psb_protect_name => psb_ls_cp_to_i
+  implicit none 
+  
+  class(psb_lsspmat_type), intent(in) :: a
+  class(psb_sspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_s_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_ls_cp_to_i
 
 
+
+  

base/serial/impl/psb_z_mat_impl.F90

@@ -2432,6 +2432,155 @@ subroutine psb_z_scals(d,a,info)
 
 end subroutine psb_z_scals
 
+subroutine psb_z_mv_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_mv_from_lb
+  implicit none 
+  
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_lfmt(b,info)
+  
+end subroutine psb_z_mv_from_lb
+
+  
+subroutine psb_z_cp_from_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cp_from_lb
+  implicit none 
+  
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_lfmt(b,info)
+
+end subroutine psb_z_cp_from_lb
+
+subroutine psb_z_mv_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_mv_to_lb
+  implicit none 
+  
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_lfmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_z_mv_to_lb
+
+subroutine psb_z_cp_to_lb(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cp_to_lb
+  implicit none 
+  class(psb_zspmat_type), intent(in) :: a
+  class(psb_lz_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_lfmt(b,info)
+  end if
+
+end subroutine psb_z_cp_to_lb
+
+subroutine psb_z_mv_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_mv_from_l
+  implicit none 
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_lzspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_z_mv_from_l
+
+  
+subroutine psb_z_cp_from_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cp_from_l
+  implicit none 
+
+  class(psb_zspmat_type), intent(out) :: a
+  class(psb_lzspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_z_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_lfmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_z_cp_from_l
+
+subroutine psb_z_mv_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_mv_to_l
+  implicit none 
+
+  class(psb_zspmat_type), intent(inout) :: a
+  class(psb_lzspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_lz_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_z_mv_to_l
+
+subroutine psb_z_cp_to_l(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_z_cp_to_l
+  implicit none 
+  
+  class(psb_zspmat_type), intent(in) :: a
+  class(psb_lzspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_lz_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_lfmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_z_cp_to_l
+
+
+!
+!
+! lz versions
+!
+
 
 subroutine  psb_lz_set_nrows(m,a) 
   use psb_z_mat_mod, psb_protect_name => psb_lz_set_nrows
@@ -4331,5 +4480,148 @@ subroutine psb_lz_scals(d,a,info)
 
 end subroutine psb_lz_scals
 
+subroutine psb_lz_mv_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_mv_from_ib
+  implicit none 
+  
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_z_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%mv_from_ifmt(b,info)
+  
+end subroutine psb_lz_mv_from_ib
+  
+subroutine psb_lz_cp_from_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cp_from_ib
+  implicit none 
+  
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_z_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info)
+  if (info == psb_success_) call a%a%cp_from_ifmt(b,info)
+
+end subroutine psb_lz_cp_from_ib
+
+subroutine psb_lz_mv_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_mv_to_ib
+  implicit none 
+  
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_z_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%mv_to_ifmt(b,info)
+    call a%free()
+  end if
+  
+end subroutine psb_lz_mv_to_ib
+
+subroutine psb_lz_cp_to_ib(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cp_to_ib
+  implicit none 
+  class(psb_lzspmat_type), intent(in) :: a
+  class(psb_z_base_sparse_mat), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (.not.allocated(a%a)) then
+    call b%free()
+  else
+    call a%a%cp_to_ifmt(b,info)
+  end if
+
+end subroutine psb_lz_cp_to_ib
+
+subroutine psb_lz_mv_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_mv_from_i
+  implicit none 
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_zspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info)
+    call a%a%mv_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+  call b%free()
+  
+end subroutine psb_lz_mv_from_i
+
+  
+subroutine psb_lz_cp_from_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cp_from_i
+  implicit none 
+
+  class(psb_lzspmat_type), intent(out) :: a
+  class(psb_zspmat_type), intent(in) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(b%a)) then
+    if (.not.allocated(a%a)) allocate(psb_lz_csr_sparse_mat :: a%a, stat=info)
+    call a%a%cp_from_ifmt(b%a,info)
+  else
+    call a%free()
+  end if
+end subroutine psb_lz_cp_from_i
+
+subroutine psb_lz_mv_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_mv_to_i
+  implicit none 
+
+  class(psb_lzspmat_type), intent(inout) :: a
+  class(psb_zspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_z_csr_sparse_mat :: b%a, stat=info)
+    call a%a%mv_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  call a%free()
+  
+end subroutine psb_lz_mv_to_i
+
+subroutine psb_lz_cp_to_i(a,b)
+  use psb_error_mod
+  use psb_const_mod
+  use psb_z_mat_mod, psb_protect_name => psb_lz_cp_to_i
+  implicit none 
+  
+  class(psb_lzspmat_type), intent(in) :: a
+  class(psb_zspmat_type), intent(inout) :: b
+  integer(psb_ipk_) :: info
+  
+  if (allocated(a%a)) then
+    if (.not.allocated(b%a)) allocate(psb_z_csr_sparse_mat :: b%a, stat=info)
+    call a%a%cp_to_ifmt(b%a,info)
+  else
+    call b%free()
+  end if
+  
+end subroutine psb_lz_cp_to_i
 
 
+
+  

base/serial/psb_cnumbmm.f90

@@ -39,7 +39,6 @@
 !       rewritten in Fortran 95/2003 making use of our sparse matrix facilities.
 !
 !
-
 subroutine psb_cnumbmm(a,b,c)
   use psb_base_mod, psb_protect_name => psb_cnumbmm
   implicit none 
@@ -234,3 +233,200 @@ contains
   end subroutine gen_numbmm
 
 end subroutine psb_cbase_numbmm
+
+
+
+subroutine psb_lcnumbmm(a,b,c)
+  use psb_base_mod, psb_protect_name => psb_lcnumbmm
+  implicit none 
+
+  type(psb_lcspmat_type), intent(in) :: a,b
+  type(psb_lcspmat_type), intent(inout)  :: c
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_numbmm'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  select type(aa=>c%a)
+  type is (psb_lc_csr_sparse_mat)
+    call psb_numbmm(a%a,b%a,aa)
+  class default
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end select
+
+  call c%set_asb()
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lcnumbmm
+
+subroutine psb_lcbase_numbmm(a,b,c)
+  use psb_mat_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_lcbase_numbmm
+  implicit none 
+
+  class(psb_lc_base_sparse_mat), intent(in) :: a,b
+  type(psb_lc_csr_sparse_mat), intent(inout)  :: c
+  integer(psb_ipk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb  
+  character(len=20)     :: name
+  complex(psb_spk_), allocatable :: temp(:)
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  name='psb_numbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+  endif
+  allocate(temp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  call psb_ensure_size(ione*size(c%ja),c%val,info)
+  select type(a)
+  type is (psb_lc_csr_sparse_mat) 
+    select type(b)
+    type is (psb_lc_csr_sparse_mat) 
+      call csr_numbmm(a,b,c,temp,info)
+    class default
+      call gen_numbmm(a,b,c,temp,info)
+    end select
+  class default
+    call gen_numbmm(a,b,c,temp,info)
+  end select
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call c%set_asb()
+  deallocate(temp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains 
+
+  subroutine csr_numbmm(a,b,c,temp,info)
+    type(psb_lc_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_lc_csr_sparse_mat), intent(inout) :: c
+    complex(psb_spk_)                          :: temp(:)
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    call lcnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,&
+         & b%irp,b%ja,lzero,b%val,&
+         & c%irp,c%ja,lzero,c%val,temp)
+
+    
+  end subroutine csr_numbmm
+
+  subroutine gen_numbmm(a,b,c,temp,info)
+    class(psb_lc_base_sparse_mat), intent(in)  :: a,b
+    type(psb_lc_csr_sparse_mat), intent(inout) :: c
+    integer(psb_ipk_) :: info
+    complex(psb_spk_)      :: temp(:)
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    complex(psb_spk_), allocatable :: aval(:),bval(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
+    complex(psb_spk_)      :: ajj
+
+    n = a%get_nrows()
+    m = a%get_ncols() 
+    l = b%get_ncols()
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = czero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m
+            info = 1
+            return
+          
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+            info = psb_err_pivot_too_small_
+            return
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%irp(i),c%irp(i+1)-1
+        if((c%ja(j)<1).or. (c%ja(j) > maxlmn))  then 
+          write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
+            info = psb_err_invalid_ovr_num_
+            return
+        else
+          c%val(j) = temp(c%ja(j))
+          temp(c%ja(j)) = czero
+        endif
+      end do
+    end do
+
+    
+  end subroutine gen_numbmm
+
+end subroutine psb_lcbase_numbmm

base/serial/psb_csymbmm.f90

@@ -255,3 +255,223 @@ contains
   end subroutine gen_symbmm
 
 end subroutine psb_cbase_symbmm
+
+
+
+subroutine psb_lcsymbmm(a,b,c,info)
+  use psb_base_mod, psb_protect_name => psb_lcsymbmm
+  implicit none 
+
+  type(psb_lcspmat_type), intent(in)    :: a,b
+  type(psb_lcspmat_type), intent(out)   :: c
+  integer(psb_ipk_), intent(out)                  :: info
+  type(psb_lc_csr_sparse_mat), allocatable :: ccsr
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  allocate(ccsr,stat=info)    
+
+  if (info == psb_success_) then 
+      call psb_symbmm(a%a,b%a,ccsr,info)
+  else
+    info = psb_err_alloc_dealloc_
+  end if
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call move_alloc(ccsr,c%a)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lcsymbmm
+
+subroutine psb_lcbase_symbmm(a,b,c,info)
+  use psb_mat_mod
+  use psb_serial_mod, psb_protect_name => psb_lcbase_symbmm
+  implicit none 
+
+  class(psb_lc_base_sparse_mat), intent(in) :: a,b
+  type(psb_lc_csr_sparse_mat), intent(out)  :: c
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_lpk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb
+  character(len=20)     :: name
+  integer(psb_ipk_) :: err_act
+  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+    info = psb_err_invalid_matrix_sizes_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  allocate(itemp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
+  select type(a)
+  type is (psb_lc_csr_sparse_mat) 
+    select type(b)
+    type is (psb_lc_csr_sparse_mat) 
+      call csr_symbmm(a,b,c,itemp,info)
+    class default
+      call gen_symbmm(a,b,c,itemp,info)
+    end select
+  class default
+    call gen_symbmm(a,b,c,itemp,info)
+  end select
+
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+
+  call psb_realloc(size(c%ja),c%val,info)
+  deallocate(itemp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains
+
+  subroutine csr_symbmm(a,b,c,itemp,info)
+    type(psb_lc_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_lc_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: itemp(:) 
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+    call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,&
+         & b%irp,b%ja,lzero,&
+         & c%irp,c%ja,lzero,itemp)
+    
+  end subroutine csr_symbmm
+  subroutine gen_symbmm(a,b,c,index,info)
+    class(psb_lc_base_sparse_mat), intent(in)  :: a,b
+    type(psb_lc_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: index(:)
+    integer(psb_ipk_) :: info
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+
+    n = ma
+    m = na 
+    l = nb 
+    maxlmn = max(l,m,n)
+
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & stat=info)
+    if (info /= psb_success_) then 
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i=1,maxlmn
+      index(i)=0
+    end do
+
+    c%irp(1)=1
+    minlm = min(l,m)
+    minmn = min(m,n)
+
+    main: do  i=1,n
+      istart=-1
+      length=0
+      call a%csget(i,i,nazr,iarw,iacl,info)
+      do jj=1, nazr
+
+        j=iacl(jj)
+
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m
+          info = 1
+          return
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
+            info=psb_err_pivot_too_small_
+            return
+          else
+            if(index(ibcl(k)) == 0) then
+              index(ibcl(k))=istart
+              istart=ibcl(k)
+              length=length+1
+            endif
+          endif
+        end do
+      end do
+
+      c%irp(i+1)=c%irp(i)+length
+
+      if (c%irp(i+1) > size(c%ja)) then 
+        if (n > (2*i)) then 
+          nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
+        else
+          nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i)))   )
+        endif
+        call psb_realloc(nze,c%ja,info)
+      end if
+      do j= c%irp(i),c%irp(i+1)-1
+        c%ja(j)=istart
+        istart=index(istart)
+        index(c%ja(j))=0
+      end do
+      call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
+      index(i) = 0
+    end do main
+
+  end subroutine gen_symbmm
+
+end subroutine psb_lcbase_symbmm

base/serial/psb_dnumbmm.f90

@@ -39,7 +39,6 @@
 !       rewritten in Fortran 95/2003 making use of our sparse matrix facilities.
 !
 !
-
 subroutine psb_dnumbmm(a,b,c)
   use psb_base_mod, psb_protect_name => psb_dnumbmm
   implicit none 
@@ -234,3 +233,200 @@ contains
   end subroutine gen_numbmm
 
 end subroutine psb_dbase_numbmm
+
+
+
+subroutine psb_ldnumbmm(a,b,c)
+  use psb_base_mod, psb_protect_name => psb_ldnumbmm
+  implicit none 
+
+  type(psb_ldspmat_type), intent(in) :: a,b
+  type(psb_ldspmat_type), intent(inout)  :: c
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_numbmm'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  select type(aa=>c%a)
+  type is (psb_ld_csr_sparse_mat)
+    call psb_numbmm(a%a,b%a,aa)
+  class default
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end select
+
+  call c%set_asb()
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_ldnumbmm
+
+subroutine psb_ldbase_numbmm(a,b,c)
+  use psb_mat_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_ldbase_numbmm
+  implicit none 
+
+  class(psb_ld_base_sparse_mat), intent(in) :: a,b
+  type(psb_ld_csr_sparse_mat), intent(inout)  :: c
+  integer(psb_ipk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb  
+  character(len=20)     :: name
+  real(psb_dpk_), allocatable :: temp(:)
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  name='psb_numbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+  endif
+  allocate(temp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  call psb_ensure_size(ione*size(c%ja),c%val,info)
+  select type(a)
+  type is (psb_ld_csr_sparse_mat) 
+    select type(b)
+    type is (psb_ld_csr_sparse_mat) 
+      call csr_numbmm(a,b,c,temp,info)
+    class default
+      call gen_numbmm(a,b,c,temp,info)
+    end select
+  class default
+    call gen_numbmm(a,b,c,temp,info)
+  end select
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call c%set_asb()
+  deallocate(temp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains 
+
+  subroutine csr_numbmm(a,b,c,temp,info)
+    type(psb_ld_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_ld_csr_sparse_mat), intent(inout) :: c
+    real(psb_dpk_)                          :: temp(:)
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    call ldnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,&
+         & b%irp,b%ja,lzero,b%val,&
+         & c%irp,c%ja,lzero,c%val,temp)
+
+    
+  end subroutine csr_numbmm
+
+  subroutine gen_numbmm(a,b,c,temp,info)
+    class(psb_ld_base_sparse_mat), intent(in)  :: a,b
+    type(psb_ld_csr_sparse_mat), intent(inout) :: c
+    integer(psb_ipk_) :: info
+    real(psb_dpk_)      :: temp(:)
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    real(psb_dpk_), allocatable :: aval(:),bval(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
+    real(psb_dpk_)      :: ajj
+
+    n = a%get_nrows()
+    m = a%get_ncols() 
+    l = b%get_ncols()
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = dzero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m
+            info = 1
+            return
+          
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+            info = psb_err_pivot_too_small_
+            return
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%irp(i),c%irp(i+1)-1
+        if((c%ja(j)<1).or. (c%ja(j) > maxlmn))  then 
+          write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
+            info = psb_err_invalid_ovr_num_
+            return
+        else
+          c%val(j) = temp(c%ja(j))
+          temp(c%ja(j)) = dzero
+        endif
+      end do
+    end do
+
+    
+  end subroutine gen_numbmm
+
+end subroutine psb_ldbase_numbmm

base/serial/psb_dsymbmm.f90

@@ -255,3 +255,223 @@ contains
   end subroutine gen_symbmm
 
 end subroutine psb_dbase_symbmm
+
+
+
+subroutine psb_ldsymbmm(a,b,c,info)
+  use psb_base_mod, psb_protect_name => psb_ldsymbmm
+  implicit none 
+
+  type(psb_ldspmat_type), intent(in)    :: a,b
+  type(psb_ldspmat_type), intent(out)   :: c
+  integer(psb_ipk_), intent(out)                  :: info
+  type(psb_ld_csr_sparse_mat), allocatable :: ccsr
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  allocate(ccsr,stat=info)    
+
+  if (info == psb_success_) then 
+      call psb_symbmm(a%a,b%a,ccsr,info)
+  else
+    info = psb_err_alloc_dealloc_
+  end if
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call move_alloc(ccsr,c%a)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_ldsymbmm
+
+subroutine psb_ldbase_symbmm(a,b,c,info)
+  use psb_mat_mod
+  use psb_serial_mod, psb_protect_name => psb_ldbase_symbmm
+  implicit none 
+
+  class(psb_ld_base_sparse_mat), intent(in) :: a,b
+  type(psb_ld_csr_sparse_mat), intent(out)  :: c
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_lpk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb
+  character(len=20)     :: name
+  integer(psb_ipk_) :: err_act
+  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+    info = psb_err_invalid_matrix_sizes_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  allocate(itemp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
+  select type(a)
+  type is (psb_ld_csr_sparse_mat) 
+    select type(b)
+    type is (psb_ld_csr_sparse_mat) 
+      call csr_symbmm(a,b,c,itemp,info)
+    class default
+      call gen_symbmm(a,b,c,itemp,info)
+    end select
+  class default
+    call gen_symbmm(a,b,c,itemp,info)
+  end select
+
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+
+  call psb_realloc(size(c%ja),c%val,info)
+  deallocate(itemp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains
+
+  subroutine csr_symbmm(a,b,c,itemp,info)
+    type(psb_ld_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_ld_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: itemp(:) 
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+    call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,&
+         & b%irp,b%ja,lzero,&
+         & c%irp,c%ja,lzero,itemp)
+    
+  end subroutine csr_symbmm
+  subroutine gen_symbmm(a,b,c,index,info)
+    class(psb_ld_base_sparse_mat), intent(in)  :: a,b
+    type(psb_ld_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: index(:)
+    integer(psb_ipk_) :: info
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+
+    n = ma
+    m = na 
+    l = nb 
+    maxlmn = max(l,m,n)
+
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & stat=info)
+    if (info /= psb_success_) then 
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i=1,maxlmn
+      index(i)=0
+    end do
+
+    c%irp(1)=1
+    minlm = min(l,m)
+    minmn = min(m,n)
+
+    main: do  i=1,n
+      istart=-1
+      length=0
+      call a%csget(i,i,nazr,iarw,iacl,info)
+      do jj=1, nazr
+
+        j=iacl(jj)
+
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m
+          info = 1
+          return
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
+            info=psb_err_pivot_too_small_
+            return
+          else
+            if(index(ibcl(k)) == 0) then
+              index(ibcl(k))=istart
+              istart=ibcl(k)
+              length=length+1
+            endif
+          endif
+        end do
+      end do
+
+      c%irp(i+1)=c%irp(i)+length
+
+      if (c%irp(i+1) > size(c%ja)) then 
+        if (n > (2*i)) then 
+          nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
+        else
+          nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i)))   )
+        endif
+        call psb_realloc(nze,c%ja,info)
+      end if
+      do j= c%irp(i),c%irp(i+1)-1
+        c%ja(j)=istart
+        istart=index(istart)
+        index(c%ja(j))=0
+      end do
+      call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
+      index(i) = 0
+    end do main
+
+  end subroutine gen_symbmm
+
+end subroutine psb_ldbase_symbmm

base/serial/psb_snumbmm.f90

@@ -39,7 +39,6 @@
 !       rewritten in Fortran 95/2003 making use of our sparse matrix facilities.
 !
 !
-
 subroutine psb_snumbmm(a,b,c)
   use psb_base_mod, psb_protect_name => psb_snumbmm
   implicit none 
@@ -234,3 +233,200 @@ contains
   end subroutine gen_numbmm
 
 end subroutine psb_sbase_numbmm
+
+
+
+subroutine psb_lsnumbmm(a,b,c)
+  use psb_base_mod, psb_protect_name => psb_lsnumbmm
+  implicit none 
+
+  type(psb_lsspmat_type), intent(in) :: a,b
+  type(psb_lsspmat_type), intent(inout)  :: c
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_numbmm'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  select type(aa=>c%a)
+  type is (psb_ls_csr_sparse_mat)
+    call psb_numbmm(a%a,b%a,aa)
+  class default
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end select
+
+  call c%set_asb()
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lsnumbmm
+
+subroutine psb_lsbase_numbmm(a,b,c)
+  use psb_mat_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_lsbase_numbmm
+  implicit none 
+
+  class(psb_ls_base_sparse_mat), intent(in) :: a,b
+  type(psb_ls_csr_sparse_mat), intent(inout)  :: c
+  integer(psb_ipk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb  
+  character(len=20)     :: name
+  real(psb_spk_), allocatable :: temp(:)
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  name='psb_numbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+  endif
+  allocate(temp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  call psb_ensure_size(ione*size(c%ja),c%val,info)
+  select type(a)
+  type is (psb_ls_csr_sparse_mat) 
+    select type(b)
+    type is (psb_ls_csr_sparse_mat) 
+      call csr_numbmm(a,b,c,temp,info)
+    class default
+      call gen_numbmm(a,b,c,temp,info)
+    end select
+  class default
+    call gen_numbmm(a,b,c,temp,info)
+  end select
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call c%set_asb()
+  deallocate(temp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains 
+
+  subroutine csr_numbmm(a,b,c,temp,info)
+    type(psb_ls_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_ls_csr_sparse_mat), intent(inout) :: c
+    real(psb_spk_)                          :: temp(:)
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    call lsnumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,&
+         & b%irp,b%ja,lzero,b%val,&
+         & c%irp,c%ja,lzero,c%val,temp)
+
+    
+  end subroutine csr_numbmm
+
+  subroutine gen_numbmm(a,b,c,temp,info)
+    class(psb_ls_base_sparse_mat), intent(in)  :: a,b
+    type(psb_ls_csr_sparse_mat), intent(inout) :: c
+    integer(psb_ipk_) :: info
+    real(psb_spk_)      :: temp(:)
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    real(psb_spk_), allocatable :: aval(:),bval(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
+    real(psb_spk_)      :: ajj
+
+    n = a%get_nrows()
+    m = a%get_ncols() 
+    l = b%get_ncols()
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = szero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m
+            info = 1
+            return
+          
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+            info = psb_err_pivot_too_small_
+            return
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%irp(i),c%irp(i+1)-1
+        if((c%ja(j)<1).or. (c%ja(j) > maxlmn))  then 
+          write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
+            info = psb_err_invalid_ovr_num_
+            return
+        else
+          c%val(j) = temp(c%ja(j))
+          temp(c%ja(j)) = szero
+        endif
+      end do
+    end do
+
+    
+  end subroutine gen_numbmm
+
+end subroutine psb_lsbase_numbmm

base/serial/psb_ssymbmm.f90

@@ -255,3 +255,223 @@ contains
   end subroutine gen_symbmm
 
 end subroutine psb_sbase_symbmm
+
+
+
+subroutine psb_lssymbmm(a,b,c,info)
+  use psb_base_mod, psb_protect_name => psb_lssymbmm
+  implicit none 
+
+  type(psb_lsspmat_type), intent(in)    :: a,b
+  type(psb_lsspmat_type), intent(out)   :: c
+  integer(psb_ipk_), intent(out)                  :: info
+  type(psb_ls_csr_sparse_mat), allocatable :: ccsr
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  allocate(ccsr,stat=info)    
+
+  if (info == psb_success_) then 
+      call psb_symbmm(a%a,b%a,ccsr,info)
+  else
+    info = psb_err_alloc_dealloc_
+  end if
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call move_alloc(ccsr,c%a)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lssymbmm
+
+subroutine psb_lsbase_symbmm(a,b,c,info)
+  use psb_mat_mod
+  use psb_serial_mod, psb_protect_name => psb_lsbase_symbmm
+  implicit none 
+
+  class(psb_ls_base_sparse_mat), intent(in) :: a,b
+  type(psb_ls_csr_sparse_mat), intent(out)  :: c
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_lpk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb
+  character(len=20)     :: name
+  integer(psb_ipk_) :: err_act
+  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+    info = psb_err_invalid_matrix_sizes_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  allocate(itemp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
+  select type(a)
+  type is (psb_ls_csr_sparse_mat) 
+    select type(b)
+    type is (psb_ls_csr_sparse_mat) 
+      call csr_symbmm(a,b,c,itemp,info)
+    class default
+      call gen_symbmm(a,b,c,itemp,info)
+    end select
+  class default
+    call gen_symbmm(a,b,c,itemp,info)
+  end select
+
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+
+  call psb_realloc(size(c%ja),c%val,info)
+  deallocate(itemp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains
+
+  subroutine csr_symbmm(a,b,c,itemp,info)
+    type(psb_ls_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_ls_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: itemp(:) 
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+    call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,&
+         & b%irp,b%ja,lzero,&
+         & c%irp,c%ja,lzero,itemp)
+    
+  end subroutine csr_symbmm
+  subroutine gen_symbmm(a,b,c,index,info)
+    class(psb_ls_base_sparse_mat), intent(in)  :: a,b
+    type(psb_ls_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: index(:)
+    integer(psb_ipk_) :: info
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+
+    n = ma
+    m = na 
+    l = nb 
+    maxlmn = max(l,m,n)
+
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & stat=info)
+    if (info /= psb_success_) then 
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i=1,maxlmn
+      index(i)=0
+    end do
+
+    c%irp(1)=1
+    minlm = min(l,m)
+    minmn = min(m,n)
+
+    main: do  i=1,n
+      istart=-1
+      length=0
+      call a%csget(i,i,nazr,iarw,iacl,info)
+      do jj=1, nazr
+
+        j=iacl(jj)
+
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m
+          info = 1
+          return
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
+            info=psb_err_pivot_too_small_
+            return
+          else
+            if(index(ibcl(k)) == 0) then
+              index(ibcl(k))=istart
+              istart=ibcl(k)
+              length=length+1
+            endif
+          endif
+        end do
+      end do
+
+      c%irp(i+1)=c%irp(i)+length
+
+      if (c%irp(i+1) > size(c%ja)) then 
+        if (n > (2*i)) then 
+          nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
+        else
+          nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i)))   )
+        endif
+        call psb_realloc(nze,c%ja,info)
+      end if
+      do j= c%irp(i),c%irp(i+1)-1
+        c%ja(j)=istart
+        istart=index(istart)
+        index(c%ja(j))=0
+      end do
+      call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
+      index(i) = 0
+    end do main
+
+  end subroutine gen_symbmm
+
+end subroutine psb_lsbase_symbmm

base/serial/psb_znumbmm.f90

@@ -39,7 +39,6 @@
 !       rewritten in Fortran 95/2003 making use of our sparse matrix facilities.
 !
 !
-
 subroutine psb_znumbmm(a,b,c)
   use psb_base_mod, psb_protect_name => psb_znumbmm
   implicit none 
@@ -234,3 +233,200 @@ contains
   end subroutine gen_numbmm
 
 end subroutine psb_zbase_numbmm
+
+
+
+subroutine psb_lznumbmm(a,b,c)
+  use psb_base_mod, psb_protect_name => psb_lznumbmm
+  implicit none 
+
+  type(psb_lzspmat_type), intent(in) :: a,b
+  type(psb_lzspmat_type), intent(inout)  :: c
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_numbmm'
+
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null()).or.(c%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  select type(aa=>c%a)
+  type is (psb_lz_csr_sparse_mat)
+    call psb_numbmm(a%a,b%a,aa)
+  class default
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  end select
+
+  call c%set_asb()
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lznumbmm
+
+subroutine psb_lzbase_numbmm(a,b,c)
+  use psb_mat_mod
+  use psb_string_mod
+  use psb_serial_mod, psb_protect_name => psb_lzbase_numbmm
+  implicit none 
+
+  class(psb_lz_base_sparse_mat), intent(in) :: a,b
+  type(psb_lz_csr_sparse_mat), intent(inout)  :: c
+  integer(psb_ipk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb  
+  character(len=20)     :: name
+  complex(psb_dpk_), allocatable :: temp(:)
+  integer(psb_ipk_) :: info
+  integer(psb_ipk_) :: err_act
+  name='psb_numbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+  endif
+  allocate(temp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+
+  !
+  ! Note: we still have to test about possible performance hits. 
+  !
+  !
+  call psb_ensure_size(ione*size(c%ja),c%val,info)
+  select type(a)
+  type is (psb_lz_csr_sparse_mat) 
+    select type(b)
+    type is (psb_lz_csr_sparse_mat) 
+      call csr_numbmm(a,b,c,temp,info)
+    class default
+      call gen_numbmm(a,b,c,temp,info)
+    end select
+  class default
+    call gen_numbmm(a,b,c,temp,info)
+  end select
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name)
+    goto 9999
+  end if
+
+  call c%set_asb()
+  deallocate(temp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains 
+
+  subroutine csr_numbmm(a,b,c,temp,info)
+    type(psb_lz_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_lz_csr_sparse_mat), intent(inout) :: c
+    complex(psb_dpk_)                          :: temp(:)
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    call lznumbmm(ma,na,nb,a%irp,a%ja,lzero,a%val,&
+         & b%irp,b%ja,lzero,b%val,&
+         & c%irp,c%ja,lzero,c%val,temp)
+
+    
+  end subroutine csr_numbmm
+
+  subroutine gen_numbmm(a,b,c,temp,info)
+    class(psb_lz_base_sparse_mat), intent(in)  :: a,b
+    type(psb_lz_csr_sparse_mat), intent(inout) :: c
+    integer(psb_ipk_) :: info
+    complex(psb_dpk_)      :: temp(:)
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    complex(psb_dpk_), allocatable :: aval(:),bval(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,nazr,nbzr,jj,minlm,minmn,minln
+    complex(psb_dpk_)      :: ajj
+
+    n = a%get_nrows()
+    m = a%get_ncols() 
+    l = b%get_ncols()
+    maxlmn = max(l,m,n)
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & aval(maxlmn),bval(maxlmn), stat=info)
+    if (info /= psb_success_) then
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i = 1,maxlmn
+      temp(i) = zzero
+    end do
+    minlm = min(l,m)
+    minln = min(l,n)
+    minmn = min(m,n)
+    do  i = 1,n
+
+      call a%csget(i,i,nazr,iarw,iacl,aval,info)
+      do jj=1, nazr
+        j=iacl(jj)
+        ajj = aval(jj)
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' NUMBMM: Problem with A ',i,jj,j,m
+            info = 1
+            return
+          
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,bval,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in NUMBM 1:',j,k,ibcl(k),maxlmn
+            info = psb_err_pivot_too_small_
+            return
+          else
+            temp(ibcl(k)) = temp(ibcl(k)) + ajj * bval(k)
+          endif
+        enddo
+      end do
+      do  j = c%irp(i),c%irp(i+1)-1
+        if((c%ja(j)<1).or. (c%ja(j) > maxlmn))  then 
+          write(psb_err_unit,*) ' NUMBMM: output problem',i,j,c%ja(j),maxlmn
+            info = psb_err_invalid_ovr_num_
+            return
+        else
+          c%val(j) = temp(c%ja(j))
+          temp(c%ja(j)) = zzero
+        endif
+      end do
+    end do
+
+    
+  end subroutine gen_numbmm
+
+end subroutine psb_lzbase_numbmm

base/serial/psb_zsymbmm.f90

@@ -255,3 +255,223 @@ contains
   end subroutine gen_symbmm
 
 end subroutine psb_zbase_symbmm
+
+
+
+subroutine psb_lzsymbmm(a,b,c,info)
+  use psb_base_mod, psb_protect_name => psb_lzsymbmm
+  implicit none 
+
+  type(psb_lzspmat_type), intent(in)    :: a,b
+  type(psb_lzspmat_type), intent(out)   :: c
+  integer(psb_ipk_), intent(out)                  :: info
+  type(psb_lz_csr_sparse_mat), allocatable :: ccsr
+  integer(psb_ipk_) :: err_act
+  character(len=*), parameter ::  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  if ((a%is_null()) .or.(b%is_null())) then
+    info = psb_err_invalid_mat_state_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+
+  allocate(ccsr,stat=info)    
+
+  if (info == psb_success_) then 
+      call psb_symbmm(a%a,b%a,ccsr,info)
+  else
+    info = psb_err_alloc_dealloc_
+  end if
+  
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+  call move_alloc(ccsr,c%a)
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+end subroutine psb_lzsymbmm
+
+subroutine psb_lzbase_symbmm(a,b,c,info)
+  use psb_mat_mod
+  use psb_serial_mod, psb_protect_name => psb_lzbase_symbmm
+  implicit none 
+
+  class(psb_lz_base_sparse_mat), intent(in) :: a,b
+  type(psb_lz_csr_sparse_mat), intent(out)  :: c
+  integer(psb_ipk_), intent(out)                     :: info
+  integer(psb_lpk_), allocatable  :: itemp(:)
+  integer(psb_lpk_) :: nze, ma,na,mb,nb
+  character(len=20)     :: name
+  integer(psb_ipk_) :: err_act
+  name='psb_symbmm'
+  call psb_erractionsave(err_act)
+  info = psb_success_
+
+  ma = a%get_nrows()
+  na = a%get_ncols()
+  mb = b%get_nrows()
+  nb = b%get_ncols()
+
+
+  if ( mb /= na ) then 
+    write(psb_err_unit,*) 'Mismatch in SYMBMM: ',ma,na,mb,nb
+    info = psb_err_invalid_matrix_sizes_
+    call psb_errpush(info,name)
+    goto 9999
+  endif
+  allocate(itemp(max(ma,na,mb,nb)),stat=info)    
+  if (info /= psb_success_) then 
+    info = psb_err_alloc_dealloc_ 
+    call psb_Errpush(info,name)
+    goto 9999
+  endif
+  !
+  ! Note: we need to test whether there is a performance impact 
+  !       in not using the original Douglas & Bank code. 
+  !
+  select type(a)
+  type is (psb_lz_csr_sparse_mat) 
+    select type(b)
+    type is (psb_lz_csr_sparse_mat) 
+      call csr_symbmm(a,b,c,itemp,info)
+    class default
+      call gen_symbmm(a,b,c,itemp,info)
+    end select
+  class default
+    call gen_symbmm(a,b,c,itemp,info)
+  end select
+
+  if (info /= psb_success_) then 
+    call psb_errpush(info,name) 
+    goto 9999
+  end if
+
+  call psb_realloc(size(c%ja),c%val,info)
+  deallocate(itemp) 
+
+  call psb_erractionrestore(err_act)
+  return
+
+9999 call psb_error_handler(err_act)
+
+  return
+
+contains
+
+  subroutine csr_symbmm(a,b,c,itemp,info)
+    type(psb_lz_csr_sparse_mat), intent(in)  :: a,b
+    type(psb_lz_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: itemp(:) 
+    integer(psb_ipk_), intent(out)                    :: info
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    info = psb_success_
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+    
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+    call lsymbmm(ma,na,nb,a%irp,a%ja,lzero,&
+         & b%irp,b%ja,lzero,&
+         & c%irp,c%ja,lzero,itemp)
+    
+  end subroutine csr_symbmm
+  subroutine gen_symbmm(a,b,c,index,info)
+    class(psb_lz_base_sparse_mat), intent(in)  :: a,b
+    type(psb_lz_csr_sparse_mat), intent(out) :: c
+    integer(psb_lpk_) :: index(:)
+    integer(psb_ipk_) :: info
+    integer(psb_lpk_), allocatable  :: iarw(:), iacl(:),ibrw(:),ibcl(:)
+    integer(psb_lpk_) :: maxlmn,i,j,m,n,k,l,istart,length,nazr,nbzr,jj,minlm,minmn
+    integer(psb_lpk_) :: nze, ma,na,mb,nb
+
+    ma = a%get_nrows()
+    na = a%get_ncols()
+    mb = b%get_nrows()
+    nb = b%get_ncols()
+
+    nze = max(ma+1,2*ma)
+    call c%allocate(ma,nb,nze)
+
+    n = ma
+    m = na 
+    l = nb 
+    maxlmn = max(l,m,n)
+
+    allocate(iarw(maxlmn),iacl(maxlmn),ibrw(maxlmn),ibcl(maxlmn),&
+         & stat=info)
+    if (info /= psb_success_) then 
+      info = psb_err_alloc_dealloc_
+      return
+    endif
+
+    do i=1,maxlmn
+      index(i)=0
+    end do
+
+    c%irp(1)=1
+    minlm = min(l,m)
+    minmn = min(m,n)
+
+    main: do  i=1,n
+      istart=-1
+      length=0
+      call a%csget(i,i,nazr,iarw,iacl,info)
+      do jj=1, nazr
+
+        j=iacl(jj)
+
+        if ((j<1).or.(j>m)) then 
+          write(psb_err_unit,*) ' SymbMM: Problem with A ',i,jj,j,m
+          info = 1
+          return
+        endif
+        call b%csget(j,j,nbzr,ibrw,ibcl,info)
+        do k=1,nbzr
+          if ((ibcl(k)<1).or.(ibcl(k)>maxlmn)) then 
+            write(psb_err_unit,*) 'Problem in SYMBMM 1:',j,k,ibcl(k),maxlmn
+            info=psb_err_pivot_too_small_
+            return
+          else
+            if(index(ibcl(k)) == 0) then
+              index(ibcl(k))=istart
+              istart=ibcl(k)
+              length=length+1
+            endif
+          endif
+        end do
+      end do
+
+      c%irp(i+1)=c%irp(i)+length
+
+      if (c%irp(i+1) > size(c%ja)) then 
+        if (n > (2*i)) then 
+          nze = max(c%irp(i+1), c%irp(i)*((n+i-1)/i))
+        else
+          nze = max(c%irp(i+1), nint((dble(c%irp(i))*(dble(n)/i)))   )
+        endif
+        call psb_realloc(nze,c%ja,info)
+      end if
+      do j= c%irp(i),c%irp(i+1)-1
+        c%ja(j)=istart
+        istart=index(istart)
+        index(c%ja(j))=0
+      end do
+      call psb_msort(c%ja(c%irp(i):c%irp(i)+length-1))
+      index(i) = 0
+    end do main
+
+  end subroutine gen_symbmm
+
+end subroutine psb_lzbase_symbmm