New LCSC stuff.

New LSMMP
New HBIO with LSPX.

base/modules/serial/psb_c_csc_mat_mod.f90

@@ -100,9 +100,58 @@ module psb_c_csc_mat_mod
 
   end type psb_c_csc_sparse_mat
 
- private :: c_csc_get_nzeros, c_csc_free,  c_csc_get_fmt, &
+  private :: c_csc_get_nzeros, c_csc_free,  c_csc_get_fmt, &
        & c_csc_get_size, c_csc_sizeof, c_csc_get_nz_col
 
+
+  !> \namespace  psb_base_mod  \class  psb_c_csc_sparse_mat
+  !! \extends psb_c_base_mat_mod::psb_c_base_sparse_mat
+  !! 
+  !! psb_c_csc_sparse_mat type and the related methods.
+  !! 
+  type, extends(psb_lc_base_sparse_mat) :: psb_lc_csc_sparse_mat
+
+    !> Pointers to beginning of cols in IA and VAL. 
+    integer(psb_lpk_), allocatable :: icp(:)
+    !> Row indices.
+    integer(psb_lpk_), allocatable :: ia(:)
+    !> Coefficient values. 
+    complex(psb_spk_), allocatable :: val(:)
+
+  contains
+    procedure, pass(a) :: is_by_cols  => lc_csc_is_by_cols
+    procedure, pass(a) :: get_size    => lc_csc_get_size
+    procedure, pass(a) :: get_nzeros  => lc_csc_get_nzeros
+    procedure, nopass  :: get_fmt     => lc_csc_get_fmt
+    procedure, pass(a) :: sizeof      => lc_csc_sizeof
+    procedure, pass(a) :: scals       => psb_lc_csc_scals
+    procedure, pass(a) :: scalv       => psb_lc_csc_scal
+    procedure, pass(a) :: reallocate_nz => psb_lc_csc_reallocate_nz
+    procedure, pass(a) :: allocate_mnnz => psb_lc_csc_allocate_mnnz
+    procedure, pass(a) :: cp_to_coo   => psb_lc_cp_csc_to_coo
+    procedure, pass(a) :: cp_from_coo => psb_lc_cp_csc_from_coo
+    procedure, pass(a) :: cp_to_fmt   => psb_lc_cp_csc_to_fmt
+    procedure, pass(a) :: cp_from_fmt => psb_lc_cp_csc_from_fmt
+    procedure, pass(a) :: mv_to_coo   => psb_lc_mv_csc_to_coo
+    procedure, pass(a) :: mv_from_coo => psb_lc_mv_csc_from_coo
+    procedure, pass(a) :: mv_to_fmt   => psb_lc_mv_csc_to_fmt
+    procedure, pass(a) :: mv_from_fmt => psb_lc_mv_csc_from_fmt
+    procedure, pass(a) :: csput_a      => psb_lc_csc_csput_a
+    procedure, pass(a) :: get_diag    => psb_lc_csc_get_diag
+    procedure, pass(a) :: csgetptn    => psb_lc_csc_csgetptn
+    procedure, pass(a) :: csgetrow   => psb_lc_csc_csgetrow
+    procedure, pass(a) :: get_nz_col  => lc_csc_get_nz_col
+    procedure, pass(a) :: reinit      => psb_lc_csc_reinit
+    procedure, pass(a) :: trim        => psb_lc_csc_trim
+    procedure, pass(a) :: print       => psb_lc_csc_print
+    procedure, pass(a) :: free        => lc_csc_free
+    procedure, pass(a) :: mold        => psb_lc_csc_mold
+
+  end type psb_lc_csc_sparse_mat
+
+  private :: lc_csc_get_nzeros, lc_csc_free,  lc_csc_get_fmt, &
+       & lc_csc_get_size, lc_csc_sizeof, lc_csc_get_nz_col
+
   !> \memberof psb_c_csc_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
   interface
@@ -492,6 +541,287 @@ module psb_c_csc_mat_mod
   end interface
   
 
+  !
+  ! lc
+  ! 
+  !> \memberof psb_lc_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reallocate_nz
+  interface
+    subroutine  psb_lc_csc_reallocate_nz(nz,a) 
+      import
+      integer(psb_lpk_), intent(in) :: nz
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_lc_csc_reallocate_nz
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reinit
+  interface 
+    subroutine psb_lc_csc_reinit(a,clear)
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a   
+      logical, intent(in), optional :: clear
+    end subroutine psb_lc_csc_reinit
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_trim
+  interface
+    subroutine  psb_lc_csc_trim(a)
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_lc_csc_trim
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_mold
+  interface 
+    subroutine psb_lc_csc_mold(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(in)                  :: a
+      class(psb_lc_base_sparse_mat), intent(inout), allocatable :: b
+      integer(psb_ipk_), intent(out)                           :: info
+    end subroutine psb_lc_csc_mold
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_allocate_mnnz
+  interface
+    subroutine  psb_lc_csc_allocate_mnnz(m,n,a,nz) 
+      import
+      integer(psb_lpk_), intent(in) :: m,n
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      integer(psb_lpk_), intent(in), optional :: nz
+    end subroutine psb_lc_csc_allocate_mnnz
+  end interface
+
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_print
+  interface
+    subroutine psb_lc_csc_print(iout,a,iv,head,ivr,ivc)
+      import
+      integer(psb_ipk_), intent(in)               :: iout
+      class(psb_lc_csc_sparse_mat), intent(in) :: a   
+      integer(psb_lpk_), intent(in), optional     :: iv(:)
+      character(len=*), optional        :: head
+      integer(psb_lpk_), intent(in), optional     :: ivr(:), ivc(:)
+    end subroutine psb_lc_csc_print
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cp_to_coo
+  interface 
+    subroutine psb_lc_cp_csc_to_coo(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(in) :: a
+      class(psb_lc_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lc_cp_csc_to_coo
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cp_from_coo
+  interface 
+    subroutine psb_lc_cp_csc_from_coo(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      class(psb_lc_coo_sparse_mat), intent(in)    :: b
+      integer(psb_ipk_), intent(out)               :: info
+    end subroutine psb_lc_cp_csc_from_coo
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cp_to_fmt
+  interface 
+    subroutine psb_lc_cp_csc_to_fmt(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(in)   :: a
+      class(psb_lc_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                       :: info
+    end subroutine psb_lc_cp_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cp_from_fmt
+  interface 
+    subroutine psb_lc_cp_csc_from_fmt(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      class(psb_lc_base_sparse_mat), intent(in)   :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lc_cp_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_mv_to_coo
+  interface 
+    subroutine psb_lc_mv_csc_to_coo(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      class(psb_lc_coo_sparse_mat), intent(inout)   :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lc_mv_csc_to_coo
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_mv_from_coo
+  interface 
+    subroutine psb_lc_mv_csc_from_coo(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      class(psb_lc_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lc_mv_csc_from_coo
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_mv_to_fmt
+  interface 
+    subroutine psb_lc_mv_csc_to_fmt(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      class(psb_lc_base_sparse_mat), intent(inout)  :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lc_mv_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_mv_from_fmt
+  interface 
+    subroutine psb_lc_mv_csc_from_fmt(a,b,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout)  :: a
+      class(psb_lc_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                         :: info
+    end subroutine psb_lc_mv_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_cp_from
+  interface 
+    subroutine psb_lc_csc_cp_from(a,b)
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      type(psb_lc_csc_sparse_mat), intent(in)   :: b
+    end subroutine psb_lc_csc_cp_from
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_mv_from
+  interface 
+    subroutine psb_lc_csc_mv_from(a,b)
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout)  :: a
+      type(psb_lc_csc_sparse_mat), intent(inout) :: b
+    end subroutine psb_lc_csc_mv_from
+  end interface
+  
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_csput_a
+  interface 
+    subroutine psb_lc_csc_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      complex(psb_spk_), intent(in)      :: val(:)
+      integer(psb_lpk_), intent(in)             :: nz,ia(:), ja(:),&
+           &  imin,imax,jmin,jmax
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_lpk_), intent(in), optional   :: gtl(:)
+    end subroutine psb_lc_csc_csput_a
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_base_mat_mod::psb_base_csgetptn
+  interface 
+    subroutine psb_lc_csc_csgetptn(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_lc_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lc_csc_csgetptn
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_csgetrow
+  interface 
+    subroutine psb_lc_csc_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_lc_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      complex(psb_spk_), allocatable,  intent(inout)    :: val(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lc_csc_csgetrow
+  end interface
+
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_csgetblk
+  interface 
+    subroutine psb_lc_csc_csgetblk(imin,imax,a,b,info,&
+       & jmin,jmax,iren,append,rscale,cscale)
+      import
+      class(psb_lc_csc_sparse_mat), intent(in) :: a
+      class(psb_lc_coo_sparse_mat), intent(inout) :: b
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lc_csc_csgetblk
+  end interface
+        
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_get_diag
+  interface 
+    subroutine psb_lc_csc_get_diag(a,d,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(in) :: a
+      complex(psb_spk_), intent(out)     :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lc_csc_get_diag
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_scal
+  interface 
+    subroutine psb_lc_csc_scal(d,a,info,side) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      complex(psb_spk_), intent(in)      :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+      character, intent(in), optional :: side
+    end subroutine psb_lc_csc_scal
+  end interface
+  
+  !> \memberof psb_lc_csc_sparse_mat
+  !! \see psb_lc_base_mat_mod::psb_lc_base_scals
+  interface
+    subroutine psb_lc_csc_scals(d,a,info) 
+      import
+      class(psb_lc_csc_sparse_mat), intent(inout) :: a
+      complex(psb_spk_), intent(in)      :: d
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lc_csc_scals
+  end interface
+  
+
+  
 contains 
 
   ! == ===================================
@@ -520,7 +850,7 @@ contains
     implicit none 
     class(psb_c_csc_sparse_mat), intent(in) :: a
     integer(psb_epk_) :: res
-    res = 8 
+    res = 2 * psb_sizeof_ip
     res = res + (2*psb_sizeof_sp)  * psb_size(a%val)
     res = res + psb_sizeof_ip * psb_size(a%icp)
     res = res + psb_sizeof_ip * psb_size(a%ia)
@@ -609,4 +939,126 @@ contains
 
   end subroutine c_csc_free
 
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Getters 
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================
+
+  
+  function lc_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_lc_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function lc_csc_is_by_cols
+
+  !
+  ! lc
+  !
+  
+  function lc_csc_sizeof(a) result(res)
+    implicit none 
+    class(psb_lc_csc_sparse_mat), intent(in) :: a
+    integer(psb_epk_) :: res
+    res = 2*psb_sizeof_lp
+    res = res + (2*psb_sizeof_sp)  * psb_size(a%val)
+    res = res + psb_sizeof_lp * psb_size(a%icp)
+    res = res + psb_sizeof_lp * psb_size(a%ia)
+      
+  end function lc_csc_sizeof
+
+  function lc_csc_get_fmt() result(res)
+    implicit none 
+    character(len=5) :: res
+    res = 'CSC'
+  end function lc_csc_get_fmt
+  
+  function lc_csc_get_nzeros(a) result(res)
+    implicit none 
+    class(psb_lc_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+    res = a%icp(a%get_ncols()+1)-1
+  end function lc_csc_get_nzeros
+
+  function lc_csc_get_size(a) result(res)
+    implicit none 
+    class(psb_lc_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+
+    res = -1
+    
+    if (allocated(a%ia)) then 
+      res = size(a%ia)
+    end if
+    if (allocated(a%val)) then 
+      if (res >= 0) then 
+        res = min(res,size(a%val))
+      else 
+        res = size(a%val)
+      end if
+    end if
+
+  end function lc_csc_get_size
+
+
+
+  function  lc_csc_get_nz_col(idx,a) result(res)
+    use psb_const_mod
+    implicit none
+    
+    class(psb_lc_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_), intent(in)                  :: idx
+    integer(psb_lpk_) :: res
+    
+    res = 0 
+ 
+    if ((1<=idx).and.(idx<=a%get_ncols())) then 
+      res = a%icp(idx+1)-a%icp(idx)
+    end if
+    
+  end function lc_csc_get_nz_col
+
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Data management
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================  
+
+
+  subroutine  lc_csc_free(a) 
+    implicit none 
+
+    class(psb_lc_csc_sparse_mat), intent(inout) :: a
+
+    if (allocated(a%icp)) deallocate(a%icp)
+    if (allocated(a%ia)) deallocate(a%ia)
+    if (allocated(a%val)) deallocate(a%val)
+    call a%set_null()
+    call a%set_nrows(0_psb_lpk_)
+    call a%set_ncols(0_psb_lpk_)
+    
+    return
+
+  end subroutine lc_csc_free
+
+
+
 end module psb_c_csc_mat_mod

base/modules/serial/psb_c_mat_mod.F90

@@ -80,7 +80,7 @@ module psb_c_mat_mod
 
   use psb_c_base_mat_mod
   use psb_c_csr_mat_mod,  only : psb_c_csr_sparse_mat, psb_lc_csr_sparse_mat
-  use psb_c_csc_mat_mod,  only : psb_c_csc_sparse_mat
+  use psb_c_csc_mat_mod,  only : psb_c_csc_sparse_mat, psb_lc_csc_sparse_mat
 
   type :: psb_cspmat_type
 

base/modules/serial/psb_c_serial_mod.f90

@@ -220,14 +220,14 @@ module psb_c_serial_mod
       type(psb_lc_csr_sparse_mat), intent(out) :: c
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_lccsrspspmm
-!!$    subroutine psb_ccscspspmm(a,b,c,info)
-!!$      use psb_c_mat_mod, only : psb_c_csc_sparse_mat
-!!$      import :: psb_ipk_
-!!$      implicit none 
-!!$      class(psb_c_csc_sparse_mat), intent(in) :: a,b
-!!$      type(psb_c_csc_sparse_mat), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)          :: info
-!!$    end subroutine psb_ccscspspmm
+    subroutine psb_lccscspspmm(a,b,c,info)
+      use psb_c_mat_mod, only : psb_lc_csc_sparse_mat
+      import :: psb_ipk_
+      implicit none 
+      class(psb_lc_csc_sparse_mat), intent(in) :: a,b
+      type(psb_lc_csc_sparse_mat), intent(out) :: c
+      integer(psb_ipk_), intent(out)          :: info
+    end subroutine psb_lccscspspmm
   end interface psb_spspmm
 
   interface psb_symbmm

base/modules/serial/psb_d_csc_mat_mod.f90

@@ -100,9 +100,58 @@ module psb_d_csc_mat_mod
 
   end type psb_d_csc_sparse_mat
 
- private :: d_csc_get_nzeros, d_csc_free,  d_csc_get_fmt, &
+  private :: d_csc_get_nzeros, d_csc_free,  d_csc_get_fmt, &
        & d_csc_get_size, d_csc_sizeof, d_csc_get_nz_col
 
+
+  !> \namespace  psb_base_mod  \class  psb_d_csc_sparse_mat
+  !! \extends psb_d_base_mat_mod::psb_d_base_sparse_mat
+  !! 
+  !! psb_d_csc_sparse_mat type and the related methods.
+  !! 
+  type, extends(psb_ld_base_sparse_mat) :: psb_ld_csc_sparse_mat
+
+    !> Pointers to beginning of cols in IA and VAL. 
+    integer(psb_lpk_), allocatable :: icp(:)
+    !> Row indices.
+    integer(psb_lpk_), allocatable :: ia(:)
+    !> Coefficient values. 
+    real(psb_dpk_), allocatable :: val(:)
+
+  contains
+    procedure, pass(a) :: is_by_cols  => ld_csc_is_by_cols
+    procedure, pass(a) :: get_size    => ld_csc_get_size
+    procedure, pass(a) :: get_nzeros  => ld_csc_get_nzeros
+    procedure, nopass  :: get_fmt     => ld_csc_get_fmt
+    procedure, pass(a) :: sizeof      => ld_csc_sizeof
+    procedure, pass(a) :: scals       => psb_ld_csc_scals
+    procedure, pass(a) :: scalv       => psb_ld_csc_scal
+    procedure, pass(a) :: reallocate_nz => psb_ld_csc_reallocate_nz
+    procedure, pass(a) :: allocate_mnnz => psb_ld_csc_allocate_mnnz
+    procedure, pass(a) :: cp_to_coo   => psb_ld_cp_csc_to_coo
+    procedure, pass(a) :: cp_from_coo => psb_ld_cp_csc_from_coo
+    procedure, pass(a) :: cp_to_fmt   => psb_ld_cp_csc_to_fmt
+    procedure, pass(a) :: cp_from_fmt => psb_ld_cp_csc_from_fmt
+    procedure, pass(a) :: mv_to_coo   => psb_ld_mv_csc_to_coo
+    procedure, pass(a) :: mv_from_coo => psb_ld_mv_csc_from_coo
+    procedure, pass(a) :: mv_to_fmt   => psb_ld_mv_csc_to_fmt
+    procedure, pass(a) :: mv_from_fmt => psb_ld_mv_csc_from_fmt
+    procedure, pass(a) :: csput_a      => psb_ld_csc_csput_a
+    procedure, pass(a) :: get_diag    => psb_ld_csc_get_diag
+    procedure, pass(a) :: csgetptn    => psb_ld_csc_csgetptn
+    procedure, pass(a) :: csgetrow   => psb_ld_csc_csgetrow
+    procedure, pass(a) :: get_nz_col  => ld_csc_get_nz_col
+    procedure, pass(a) :: reinit      => psb_ld_csc_reinit
+    procedure, pass(a) :: trim        => psb_ld_csc_trim
+    procedure, pass(a) :: print       => psb_ld_csc_print
+    procedure, pass(a) :: free        => ld_csc_free
+    procedure, pass(a) :: mold        => psb_ld_csc_mold
+
+  end type psb_ld_csc_sparse_mat
+
+  private :: ld_csc_get_nzeros, ld_csc_free,  ld_csc_get_fmt, &
+       & ld_csc_get_size, ld_csc_sizeof, ld_csc_get_nz_col
+
   !> \memberof psb_d_csc_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
   interface
@@ -492,6 +541,287 @@ module psb_d_csc_mat_mod
   end interface
   
 
+  !
+  ! ld
+  ! 
+  !> \memberof psb_ld_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reallocate_nz
+  interface
+    subroutine  psb_ld_csc_reallocate_nz(nz,a) 
+      import
+      integer(psb_lpk_), intent(in) :: nz
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_ld_csc_reallocate_nz
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reinit
+  interface 
+    subroutine psb_ld_csc_reinit(a,clear)
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a   
+      logical, intent(in), optional :: clear
+    end subroutine psb_ld_csc_reinit
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_trim
+  interface
+    subroutine  psb_ld_csc_trim(a)
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_ld_csc_trim
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_mold
+  interface 
+    subroutine psb_ld_csc_mold(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(in)                  :: a
+      class(psb_ld_base_sparse_mat), intent(inout), allocatable :: b
+      integer(psb_ipk_), intent(out)                           :: info
+    end subroutine psb_ld_csc_mold
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_allocate_mnnz
+  interface
+    subroutine  psb_ld_csc_allocate_mnnz(m,n,a,nz) 
+      import
+      integer(psb_lpk_), intent(in) :: m,n
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      integer(psb_lpk_), intent(in), optional :: nz
+    end subroutine psb_ld_csc_allocate_mnnz
+  end interface
+
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_print
+  interface
+    subroutine psb_ld_csc_print(iout,a,iv,head,ivr,ivc)
+      import
+      integer(psb_ipk_), intent(in)               :: iout
+      class(psb_ld_csc_sparse_mat), intent(in) :: a   
+      integer(psb_lpk_), intent(in), optional     :: iv(:)
+      character(len=*), optional        :: head
+      integer(psb_lpk_), intent(in), optional     :: ivr(:), ivc(:)
+    end subroutine psb_ld_csc_print
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cp_to_coo
+  interface 
+    subroutine psb_ld_cp_csc_to_coo(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(in) :: a
+      class(psb_ld_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ld_cp_csc_to_coo
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cp_from_coo
+  interface 
+    subroutine psb_ld_cp_csc_from_coo(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      class(psb_ld_coo_sparse_mat), intent(in)    :: b
+      integer(psb_ipk_), intent(out)               :: info
+    end subroutine psb_ld_cp_csc_from_coo
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cp_to_fmt
+  interface 
+    subroutine psb_ld_cp_csc_to_fmt(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(in)   :: a
+      class(psb_ld_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                       :: info
+    end subroutine psb_ld_cp_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cp_from_fmt
+  interface 
+    subroutine psb_ld_cp_csc_from_fmt(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      class(psb_ld_base_sparse_mat), intent(in)   :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ld_cp_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_mv_to_coo
+  interface 
+    subroutine psb_ld_mv_csc_to_coo(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      class(psb_ld_coo_sparse_mat), intent(inout)   :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ld_mv_csc_to_coo
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_mv_from_coo
+  interface 
+    subroutine psb_ld_mv_csc_from_coo(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      class(psb_ld_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ld_mv_csc_from_coo
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_mv_to_fmt
+  interface 
+    subroutine psb_ld_mv_csc_to_fmt(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      class(psb_ld_base_sparse_mat), intent(inout)  :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ld_mv_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_mv_from_fmt
+  interface 
+    subroutine psb_ld_mv_csc_from_fmt(a,b,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout)  :: a
+      class(psb_ld_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                         :: info
+    end subroutine psb_ld_mv_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_cp_from
+  interface 
+    subroutine psb_ld_csc_cp_from(a,b)
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      type(psb_ld_csc_sparse_mat), intent(in)   :: b
+    end subroutine psb_ld_csc_cp_from
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_mv_from
+  interface 
+    subroutine psb_ld_csc_mv_from(a,b)
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout)  :: a
+      type(psb_ld_csc_sparse_mat), intent(inout) :: b
+    end subroutine psb_ld_csc_mv_from
+  end interface
+  
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_csput_a
+  interface 
+    subroutine psb_ld_csc_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      real(psb_dpk_), intent(in)      :: val(:)
+      integer(psb_lpk_), intent(in)             :: nz,ia(:), ja(:),&
+           &  imin,imax,jmin,jmax
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_lpk_), intent(in), optional   :: gtl(:)
+    end subroutine psb_ld_csc_csput_a
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_base_mat_mod::psb_base_csgetptn
+  interface 
+    subroutine psb_ld_csc_csgetptn(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_ld_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ld_csc_csgetptn
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_csgetrow
+  interface 
+    subroutine psb_ld_csc_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_ld_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      real(psb_dpk_), allocatable,  intent(inout)    :: val(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ld_csc_csgetrow
+  end interface
+
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_csgetblk
+  interface 
+    subroutine psb_ld_csc_csgetblk(imin,imax,a,b,info,&
+       & jmin,jmax,iren,append,rscale,cscale)
+      import
+      class(psb_ld_csc_sparse_mat), intent(in) :: a
+      class(psb_ld_coo_sparse_mat), intent(inout) :: b
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ld_csc_csgetblk
+  end interface
+        
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_get_diag
+  interface 
+    subroutine psb_ld_csc_get_diag(a,d,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(in) :: a
+      real(psb_dpk_), intent(out)     :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ld_csc_get_diag
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_scal
+  interface 
+    subroutine psb_ld_csc_scal(d,a,info,side) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      real(psb_dpk_), intent(in)      :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+      character, intent(in), optional :: side
+    end subroutine psb_ld_csc_scal
+  end interface
+  
+  !> \memberof psb_ld_csc_sparse_mat
+  !! \see psb_ld_base_mat_mod::psb_ld_base_scals
+  interface
+    subroutine psb_ld_csc_scals(d,a,info) 
+      import
+      class(psb_ld_csc_sparse_mat), intent(inout) :: a
+      real(psb_dpk_), intent(in)      :: d
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ld_csc_scals
+  end interface
+  
+
+  
 contains 
 
   ! == ===================================
@@ -520,7 +850,7 @@ contains
     implicit none 
     class(psb_d_csc_sparse_mat), intent(in) :: a
     integer(psb_epk_) :: res
-    res = 8 
+    res = 2 * psb_sizeof_ip
     res = res + psb_sizeof_dp  * psb_size(a%val)
     res = res + psb_sizeof_ip * psb_size(a%icp)
     res = res + psb_sizeof_ip * psb_size(a%ia)
@@ -609,4 +939,126 @@ contains
 
   end subroutine d_csc_free
 
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Getters 
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================
+
+  
+  function ld_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_ld_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function ld_csc_is_by_cols
+
+  !
+  ! ld
+  !
+  
+  function ld_csc_sizeof(a) result(res)
+    implicit none 
+    class(psb_ld_csc_sparse_mat), intent(in) :: a
+    integer(psb_epk_) :: res
+    res = 2*psb_sizeof_lp
+    res = res + psb_sizeof_dp  * psb_size(a%val)
+    res = res + psb_sizeof_lp * psb_size(a%icp)
+    res = res + psb_sizeof_lp * psb_size(a%ia)
+      
+  end function ld_csc_sizeof
+
+  function ld_csc_get_fmt() result(res)
+    implicit none 
+    character(len=5) :: res
+    res = 'CSC'
+  end function ld_csc_get_fmt
+  
+  function ld_csc_get_nzeros(a) result(res)
+    implicit none 
+    class(psb_ld_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+    res = a%icp(a%get_ncols()+1)-1
+  end function ld_csc_get_nzeros
+
+  function ld_csc_get_size(a) result(res)
+    implicit none 
+    class(psb_ld_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+
+    res = -1
+    
+    if (allocated(a%ia)) then 
+      res = size(a%ia)
+    end if
+    if (allocated(a%val)) then 
+      if (res >= 0) then 
+        res = min(res,size(a%val))
+      else 
+        res = size(a%val)
+      end if
+    end if
+
+  end function ld_csc_get_size
+
+
+
+  function  ld_csc_get_nz_col(idx,a) result(res)
+    use psb_const_mod
+    implicit none
+    
+    class(psb_ld_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_), intent(in)                  :: idx
+    integer(psb_lpk_) :: res
+    
+    res = 0 
+ 
+    if ((1<=idx).and.(idx<=a%get_ncols())) then 
+      res = a%icp(idx+1)-a%icp(idx)
+    end if
+    
+  end function ld_csc_get_nz_col
+
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Data management
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================  
+
+
+  subroutine  ld_csc_free(a) 
+    implicit none 
+
+    class(psb_ld_csc_sparse_mat), intent(inout) :: a
+
+    if (allocated(a%icp)) deallocate(a%icp)
+    if (allocated(a%ia)) deallocate(a%ia)
+    if (allocated(a%val)) deallocate(a%val)
+    call a%set_null()
+    call a%set_nrows(0_psb_lpk_)
+    call a%set_ncols(0_psb_lpk_)
+    
+    return
+
+  end subroutine ld_csc_free
+
+
+
 end module psb_d_csc_mat_mod

base/modules/serial/psb_d_mat_mod.F90

@@ -80,7 +80,7 @@ module psb_d_mat_mod
 
   use psb_d_base_mat_mod
   use psb_d_csr_mat_mod,  only : psb_d_csr_sparse_mat, psb_ld_csr_sparse_mat
-  use psb_d_csc_mat_mod,  only : psb_d_csc_sparse_mat
+  use psb_d_csc_mat_mod,  only : psb_d_csc_sparse_mat, psb_ld_csc_sparse_mat
 
   type :: psb_dspmat_type
 

base/modules/serial/psb_d_serial_mod.f90

@@ -220,14 +220,14 @@ module psb_d_serial_mod
       type(psb_ld_csr_sparse_mat), intent(out) :: c
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_ldcsrspspmm
-!!$    subroutine psb_dcscspspmm(a,b,c,info)
-!!$      use psb_d_mat_mod, only : psb_d_csc_sparse_mat
-!!$      import :: psb_ipk_
-!!$      implicit none 
-!!$      class(psb_d_csc_sparse_mat), intent(in) :: a,b
-!!$      type(psb_d_csc_sparse_mat), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)          :: info
-!!$    end subroutine psb_dcscspspmm
+    subroutine psb_ldcscspspmm(a,b,c,info)
+      use psb_d_mat_mod, only : psb_ld_csc_sparse_mat
+      import :: psb_ipk_
+      implicit none 
+      class(psb_ld_csc_sparse_mat), intent(in) :: a,b
+      type(psb_ld_csc_sparse_mat), intent(out) :: c
+      integer(psb_ipk_), intent(out)          :: info
+    end subroutine psb_ldcscspspmm
   end interface psb_spspmm
 
   interface psb_symbmm

base/modules/serial/psb_s_csc_mat_mod.f90

@@ -100,9 +100,58 @@ module psb_s_csc_mat_mod
 
   end type psb_s_csc_sparse_mat
 
- private :: s_csc_get_nzeros, s_csc_free,  s_csc_get_fmt, &
+  private :: s_csc_get_nzeros, s_csc_free,  s_csc_get_fmt, &
        & s_csc_get_size, s_csc_sizeof, s_csc_get_nz_col
 
+
+  !> \namespace  psb_base_mod  \class  psb_s_csc_sparse_mat
+  !! \extends psb_s_base_mat_mod::psb_s_base_sparse_mat
+  !! 
+  !! psb_s_csc_sparse_mat type and the related methods.
+  !! 
+  type, extends(psb_ls_base_sparse_mat) :: psb_ls_csc_sparse_mat
+
+    !> Pointers to beginning of cols in IA and VAL. 
+    integer(psb_lpk_), allocatable :: icp(:)
+    !> Row indices.
+    integer(psb_lpk_), allocatable :: ia(:)
+    !> Coefficient values. 
+    real(psb_spk_), allocatable :: val(:)
+
+  contains
+    procedure, pass(a) :: is_by_cols  => ls_csc_is_by_cols
+    procedure, pass(a) :: get_size    => ls_csc_get_size
+    procedure, pass(a) :: get_nzeros  => ls_csc_get_nzeros
+    procedure, nopass  :: get_fmt     => ls_csc_get_fmt
+    procedure, pass(a) :: sizeof      => ls_csc_sizeof
+    procedure, pass(a) :: scals       => psb_ls_csc_scals
+    procedure, pass(a) :: scalv       => psb_ls_csc_scal
+    procedure, pass(a) :: reallocate_nz => psb_ls_csc_reallocate_nz
+    procedure, pass(a) :: allocate_mnnz => psb_ls_csc_allocate_mnnz
+    procedure, pass(a) :: cp_to_coo   => psb_ls_cp_csc_to_coo
+    procedure, pass(a) :: cp_from_coo => psb_ls_cp_csc_from_coo
+    procedure, pass(a) :: cp_to_fmt   => psb_ls_cp_csc_to_fmt
+    procedure, pass(a) :: cp_from_fmt => psb_ls_cp_csc_from_fmt
+    procedure, pass(a) :: mv_to_coo   => psb_ls_mv_csc_to_coo
+    procedure, pass(a) :: mv_from_coo => psb_ls_mv_csc_from_coo
+    procedure, pass(a) :: mv_to_fmt   => psb_ls_mv_csc_to_fmt
+    procedure, pass(a) :: mv_from_fmt => psb_ls_mv_csc_from_fmt
+    procedure, pass(a) :: csput_a      => psb_ls_csc_csput_a
+    procedure, pass(a) :: get_diag    => psb_ls_csc_get_diag
+    procedure, pass(a) :: csgetptn    => psb_ls_csc_csgetptn
+    procedure, pass(a) :: csgetrow   => psb_ls_csc_csgetrow
+    procedure, pass(a) :: get_nz_col  => ls_csc_get_nz_col
+    procedure, pass(a) :: reinit      => psb_ls_csc_reinit
+    procedure, pass(a) :: trim        => psb_ls_csc_trim
+    procedure, pass(a) :: print       => psb_ls_csc_print
+    procedure, pass(a) :: free        => ls_csc_free
+    procedure, pass(a) :: mold        => psb_ls_csc_mold
+
+  end type psb_ls_csc_sparse_mat
+
+  private :: ls_csc_get_nzeros, ls_csc_free,  ls_csc_get_fmt, &
+       & ls_csc_get_size, ls_csc_sizeof, ls_csc_get_nz_col
+
   !> \memberof psb_s_csc_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
   interface
@@ -492,6 +541,287 @@ module psb_s_csc_mat_mod
   end interface
   
 
+  !
+  ! ls
+  ! 
+  !> \memberof psb_ls_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reallocate_nz
+  interface
+    subroutine  psb_ls_csc_reallocate_nz(nz,a) 
+      import
+      integer(psb_lpk_), intent(in) :: nz
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_ls_csc_reallocate_nz
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reinit
+  interface 
+    subroutine psb_ls_csc_reinit(a,clear)
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a   
+      logical, intent(in), optional :: clear
+    end subroutine psb_ls_csc_reinit
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_trim
+  interface
+    subroutine  psb_ls_csc_trim(a)
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_ls_csc_trim
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_mold
+  interface 
+    subroutine psb_ls_csc_mold(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(in)                  :: a
+      class(psb_ls_base_sparse_mat), intent(inout), allocatable :: b
+      integer(psb_ipk_), intent(out)                           :: info
+    end subroutine psb_ls_csc_mold
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_allocate_mnnz
+  interface
+    subroutine  psb_ls_csc_allocate_mnnz(m,n,a,nz) 
+      import
+      integer(psb_lpk_), intent(in) :: m,n
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      integer(psb_lpk_), intent(in), optional :: nz
+    end subroutine psb_ls_csc_allocate_mnnz
+  end interface
+
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_print
+  interface
+    subroutine psb_ls_csc_print(iout,a,iv,head,ivr,ivc)
+      import
+      integer(psb_ipk_), intent(in)               :: iout
+      class(psb_ls_csc_sparse_mat), intent(in) :: a   
+      integer(psb_lpk_), intent(in), optional     :: iv(:)
+      character(len=*), optional        :: head
+      integer(psb_lpk_), intent(in), optional     :: ivr(:), ivc(:)
+    end subroutine psb_ls_csc_print
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cp_to_coo
+  interface 
+    subroutine psb_ls_cp_csc_to_coo(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(in) :: a
+      class(psb_ls_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ls_cp_csc_to_coo
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cp_from_coo
+  interface 
+    subroutine psb_ls_cp_csc_from_coo(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      class(psb_ls_coo_sparse_mat), intent(in)    :: b
+      integer(psb_ipk_), intent(out)               :: info
+    end subroutine psb_ls_cp_csc_from_coo
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cp_to_fmt
+  interface 
+    subroutine psb_ls_cp_csc_to_fmt(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(in)   :: a
+      class(psb_ls_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                       :: info
+    end subroutine psb_ls_cp_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cp_from_fmt
+  interface 
+    subroutine psb_ls_cp_csc_from_fmt(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      class(psb_ls_base_sparse_mat), intent(in)   :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ls_cp_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_mv_to_coo
+  interface 
+    subroutine psb_ls_mv_csc_to_coo(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      class(psb_ls_coo_sparse_mat), intent(inout)   :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ls_mv_csc_to_coo
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_mv_from_coo
+  interface 
+    subroutine psb_ls_mv_csc_from_coo(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      class(psb_ls_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ls_mv_csc_from_coo
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_mv_to_fmt
+  interface 
+    subroutine psb_ls_mv_csc_to_fmt(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      class(psb_ls_base_sparse_mat), intent(inout)  :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_ls_mv_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_mv_from_fmt
+  interface 
+    subroutine psb_ls_mv_csc_from_fmt(a,b,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout)  :: a
+      class(psb_ls_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                         :: info
+    end subroutine psb_ls_mv_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_cp_from
+  interface 
+    subroutine psb_ls_csc_cp_from(a,b)
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      type(psb_ls_csc_sparse_mat), intent(in)   :: b
+    end subroutine psb_ls_csc_cp_from
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_mv_from
+  interface 
+    subroutine psb_ls_csc_mv_from(a,b)
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout)  :: a
+      type(psb_ls_csc_sparse_mat), intent(inout) :: b
+    end subroutine psb_ls_csc_mv_from
+  end interface
+  
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_csput_a
+  interface 
+    subroutine psb_ls_csc_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      real(psb_spk_), intent(in)      :: val(:)
+      integer(psb_lpk_), intent(in)             :: nz,ia(:), ja(:),&
+           &  imin,imax,jmin,jmax
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_lpk_), intent(in), optional   :: gtl(:)
+    end subroutine psb_ls_csc_csput_a
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_base_mat_mod::psb_base_csgetptn
+  interface 
+    subroutine psb_ls_csc_csgetptn(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_ls_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ls_csc_csgetptn
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_csgetrow
+  interface 
+    subroutine psb_ls_csc_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_ls_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      real(psb_spk_), allocatable,  intent(inout)    :: val(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ls_csc_csgetrow
+  end interface
+
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_csgetblk
+  interface 
+    subroutine psb_ls_csc_csgetblk(imin,imax,a,b,info,&
+       & jmin,jmax,iren,append,rscale,cscale)
+      import
+      class(psb_ls_csc_sparse_mat), intent(in) :: a
+      class(psb_ls_coo_sparse_mat), intent(inout) :: b
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_ls_csc_csgetblk
+  end interface
+        
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_get_diag
+  interface 
+    subroutine psb_ls_csc_get_diag(a,d,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(in) :: a
+      real(psb_spk_), intent(out)     :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ls_csc_get_diag
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_scal
+  interface 
+    subroutine psb_ls_csc_scal(d,a,info,side) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      real(psb_spk_), intent(in)      :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+      character, intent(in), optional :: side
+    end subroutine psb_ls_csc_scal
+  end interface
+  
+  !> \memberof psb_ls_csc_sparse_mat
+  !! \see psb_ls_base_mat_mod::psb_ls_base_scals
+  interface
+    subroutine psb_ls_csc_scals(d,a,info) 
+      import
+      class(psb_ls_csc_sparse_mat), intent(inout) :: a
+      real(psb_spk_), intent(in)      :: d
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_ls_csc_scals
+  end interface
+  
+
+  
 contains 
 
   ! == ===================================
@@ -520,7 +850,7 @@ contains
     implicit none 
     class(psb_s_csc_sparse_mat), intent(in) :: a
     integer(psb_epk_) :: res
-    res = 8 
+    res = 2 * psb_sizeof_ip
     res = res + psb_sizeof_sp  * psb_size(a%val)
     res = res + psb_sizeof_ip * psb_size(a%icp)
     res = res + psb_sizeof_ip * psb_size(a%ia)
@@ -609,4 +939,126 @@ contains
 
   end subroutine s_csc_free
 
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Getters 
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================
+
+  
+  function ls_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_ls_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function ls_csc_is_by_cols
+
+  !
+  ! ls
+  !
+  
+  function ls_csc_sizeof(a) result(res)
+    implicit none 
+    class(psb_ls_csc_sparse_mat), intent(in) :: a
+    integer(psb_epk_) :: res
+    res = 2*psb_sizeof_lp
+    res = res + psb_sizeof_sp  * psb_size(a%val)
+    res = res + psb_sizeof_lp * psb_size(a%icp)
+    res = res + psb_sizeof_lp * psb_size(a%ia)
+      
+  end function ls_csc_sizeof
+
+  function ls_csc_get_fmt() result(res)
+    implicit none 
+    character(len=5) :: res
+    res = 'CSC'
+  end function ls_csc_get_fmt
+  
+  function ls_csc_get_nzeros(a) result(res)
+    implicit none 
+    class(psb_ls_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+    res = a%icp(a%get_ncols()+1)-1
+  end function ls_csc_get_nzeros
+
+  function ls_csc_get_size(a) result(res)
+    implicit none 
+    class(psb_ls_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+
+    res = -1
+    
+    if (allocated(a%ia)) then 
+      res = size(a%ia)
+    end if
+    if (allocated(a%val)) then 
+      if (res >= 0) then 
+        res = min(res,size(a%val))
+      else 
+        res = size(a%val)
+      end if
+    end if
+
+  end function ls_csc_get_size
+
+
+
+  function  ls_csc_get_nz_col(idx,a) result(res)
+    use psb_const_mod
+    implicit none
+    
+    class(psb_ls_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_), intent(in)                  :: idx
+    integer(psb_lpk_) :: res
+    
+    res = 0 
+ 
+    if ((1<=idx).and.(idx<=a%get_ncols())) then 
+      res = a%icp(idx+1)-a%icp(idx)
+    end if
+    
+  end function ls_csc_get_nz_col
+
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Data management
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================  
+
+
+  subroutine  ls_csc_free(a) 
+    implicit none 
+
+    class(psb_ls_csc_sparse_mat), intent(inout) :: a
+
+    if (allocated(a%icp)) deallocate(a%icp)
+    if (allocated(a%ia)) deallocate(a%ia)
+    if (allocated(a%val)) deallocate(a%val)
+    call a%set_null()
+    call a%set_nrows(0_psb_lpk_)
+    call a%set_ncols(0_psb_lpk_)
+    
+    return
+
+  end subroutine ls_csc_free
+
+
+
 end module psb_s_csc_mat_mod

base/modules/serial/psb_s_mat_mod.F90

@@ -80,7 +80,7 @@ module psb_s_mat_mod
 
   use psb_s_base_mat_mod
   use psb_s_csr_mat_mod,  only : psb_s_csr_sparse_mat, psb_ls_csr_sparse_mat
-  use psb_s_csc_mat_mod,  only : psb_s_csc_sparse_mat
+  use psb_s_csc_mat_mod,  only : psb_s_csc_sparse_mat, psb_ls_csc_sparse_mat
 
   type :: psb_sspmat_type
 

base/modules/serial/psb_s_serial_mod.f90

@@ -220,14 +220,14 @@ module psb_s_serial_mod
       type(psb_ls_csr_sparse_mat), intent(out) :: c
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_lscsrspspmm
-!!$    subroutine psb_scscspspmm(a,b,c,info)
-!!$      use psb_s_mat_mod, only : psb_s_csc_sparse_mat
-!!$      import :: psb_ipk_
-!!$      implicit none 
-!!$      class(psb_s_csc_sparse_mat), intent(in) :: a,b
-!!$      type(psb_s_csc_sparse_mat), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)          :: info
-!!$    end subroutine psb_scscspspmm
+    subroutine psb_lscscspspmm(a,b,c,info)
+      use psb_s_mat_mod, only : psb_ls_csc_sparse_mat
+      import :: psb_ipk_
+      implicit none 
+      class(psb_ls_csc_sparse_mat), intent(in) :: a,b
+      type(psb_ls_csc_sparse_mat), intent(out) :: c
+      integer(psb_ipk_), intent(out)          :: info
+    end subroutine psb_lscscspspmm
   end interface psb_spspmm
 
   interface psb_symbmm

base/modules/serial/psb_z_csc_mat_mod.f90

@@ -100,9 +100,58 @@ module psb_z_csc_mat_mod
 
   end type psb_z_csc_sparse_mat
 
- private :: z_csc_get_nzeros, z_csc_free,  z_csc_get_fmt, &
+  private :: z_csc_get_nzeros, z_csc_free,  z_csc_get_fmt, &
        & z_csc_get_size, z_csc_sizeof, z_csc_get_nz_col
 
+
+  !> \namespace  psb_base_mod  \class  psb_z_csc_sparse_mat
+  !! \extends psb_z_base_mat_mod::psb_z_base_sparse_mat
+  !! 
+  !! psb_z_csc_sparse_mat type and the related methods.
+  !! 
+  type, extends(psb_lz_base_sparse_mat) :: psb_lz_csc_sparse_mat
+
+    !> Pointers to beginning of cols in IA and VAL. 
+    integer(psb_lpk_), allocatable :: icp(:)
+    !> Row indices.
+    integer(psb_lpk_), allocatable :: ia(:)
+    !> Coefficient values. 
+    complex(psb_dpk_), allocatable :: val(:)
+
+  contains
+    procedure, pass(a) :: is_by_cols  => lz_csc_is_by_cols
+    procedure, pass(a) :: get_size    => lz_csc_get_size
+    procedure, pass(a) :: get_nzeros  => lz_csc_get_nzeros
+    procedure, nopass  :: get_fmt     => lz_csc_get_fmt
+    procedure, pass(a) :: sizeof      => lz_csc_sizeof
+    procedure, pass(a) :: scals       => psb_lz_csc_scals
+    procedure, pass(a) :: scalv       => psb_lz_csc_scal
+    procedure, pass(a) :: reallocate_nz => psb_lz_csc_reallocate_nz
+    procedure, pass(a) :: allocate_mnnz => psb_lz_csc_allocate_mnnz
+    procedure, pass(a) :: cp_to_coo   => psb_lz_cp_csc_to_coo
+    procedure, pass(a) :: cp_from_coo => psb_lz_cp_csc_from_coo
+    procedure, pass(a) :: cp_to_fmt   => psb_lz_cp_csc_to_fmt
+    procedure, pass(a) :: cp_from_fmt => psb_lz_cp_csc_from_fmt
+    procedure, pass(a) :: mv_to_coo   => psb_lz_mv_csc_to_coo
+    procedure, pass(a) :: mv_from_coo => psb_lz_mv_csc_from_coo
+    procedure, pass(a) :: mv_to_fmt   => psb_lz_mv_csc_to_fmt
+    procedure, pass(a) :: mv_from_fmt => psb_lz_mv_csc_from_fmt
+    procedure, pass(a) :: csput_a      => psb_lz_csc_csput_a
+    procedure, pass(a) :: get_diag    => psb_lz_csc_get_diag
+    procedure, pass(a) :: csgetptn    => psb_lz_csc_csgetptn
+    procedure, pass(a) :: csgetrow   => psb_lz_csc_csgetrow
+    procedure, pass(a) :: get_nz_col  => lz_csc_get_nz_col
+    procedure, pass(a) :: reinit      => psb_lz_csc_reinit
+    procedure, pass(a) :: trim        => psb_lz_csc_trim
+    procedure, pass(a) :: print       => psb_lz_csc_print
+    procedure, pass(a) :: free        => lz_csc_free
+    procedure, pass(a) :: mold        => psb_lz_csc_mold
+
+  end type psb_lz_csc_sparse_mat
+
+  private :: lz_csc_get_nzeros, lz_csc_free,  lz_csc_get_fmt, &
+       & lz_csc_get_size, lz_csc_sizeof, lz_csc_get_nz_col
+
   !> \memberof psb_z_csc_sparse_mat
   !| \see psb_base_mat_mod::psb_base_reallocate_nz
   interface
@@ -492,6 +541,287 @@ module psb_z_csc_mat_mod
   end interface
   
 
+  !
+  ! lz
+  ! 
+  !> \memberof psb_lz_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reallocate_nz
+  interface
+    subroutine  psb_lz_csc_reallocate_nz(nz,a) 
+      import
+      integer(psb_lpk_), intent(in) :: nz
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_lz_csc_reallocate_nz
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_reinit
+  interface 
+    subroutine psb_lz_csc_reinit(a,clear)
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a   
+      logical, intent(in), optional :: clear
+    end subroutine psb_lz_csc_reinit
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_trim
+  interface
+    subroutine  psb_lz_csc_trim(a)
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+    end subroutine psb_lz_csc_trim
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_mold
+  interface 
+    subroutine psb_lz_csc_mold(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(in)                  :: a
+      class(psb_lz_base_sparse_mat), intent(inout), allocatable :: b
+      integer(psb_ipk_), intent(out)                           :: info
+    end subroutine psb_lz_csc_mold
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !| \see psb_base_mat_mod::psb_base_allocate_mnnz
+  interface
+    subroutine  psb_lz_csc_allocate_mnnz(m,n,a,nz) 
+      import
+      integer(psb_lpk_), intent(in) :: m,n
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      integer(psb_lpk_), intent(in), optional :: nz
+    end subroutine psb_lz_csc_allocate_mnnz
+  end interface
+
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_print
+  interface
+    subroutine psb_lz_csc_print(iout,a,iv,head,ivr,ivc)
+      import
+      integer(psb_ipk_), intent(in)               :: iout
+      class(psb_lz_csc_sparse_mat), intent(in) :: a   
+      integer(psb_lpk_), intent(in), optional     :: iv(:)
+      character(len=*), optional        :: head
+      integer(psb_lpk_), intent(in), optional     :: ivr(:), ivc(:)
+    end subroutine psb_lz_csc_print
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cp_to_coo
+  interface 
+    subroutine psb_lz_cp_csc_to_coo(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(in) :: a
+      class(psb_lz_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lz_cp_csc_to_coo
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cp_from_coo
+  interface 
+    subroutine psb_lz_cp_csc_from_coo(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      class(psb_lz_coo_sparse_mat), intent(in)    :: b
+      integer(psb_ipk_), intent(out)               :: info
+    end subroutine psb_lz_cp_csc_from_coo
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cp_to_fmt
+  interface 
+    subroutine psb_lz_cp_csc_to_fmt(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(in)   :: a
+      class(psb_lz_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                       :: info
+    end subroutine psb_lz_cp_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cp_from_fmt
+  interface 
+    subroutine psb_lz_cp_csc_from_fmt(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      class(psb_lz_base_sparse_mat), intent(in)   :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lz_cp_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_mv_to_coo
+  interface 
+    subroutine psb_lz_mv_csc_to_coo(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      class(psb_lz_coo_sparse_mat), intent(inout)   :: b
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lz_mv_csc_to_coo
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_mv_from_coo
+  interface 
+    subroutine psb_lz_mv_csc_from_coo(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      class(psb_lz_coo_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lz_mv_csc_from_coo
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_mv_to_fmt
+  interface 
+    subroutine psb_lz_mv_csc_to_fmt(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      class(psb_lz_base_sparse_mat), intent(inout)  :: b
+      integer(psb_ipk_), intent(out)                        :: info
+    end subroutine psb_lz_mv_csc_to_fmt
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_mv_from_fmt
+  interface 
+    subroutine psb_lz_mv_csc_from_fmt(a,b,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout)  :: a
+      class(psb_lz_base_sparse_mat), intent(inout) :: b
+      integer(psb_ipk_), intent(out)                         :: info
+    end subroutine psb_lz_mv_csc_from_fmt
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_cp_from
+  interface 
+    subroutine psb_lz_csc_cp_from(a,b)
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      type(psb_lz_csc_sparse_mat), intent(in)   :: b
+    end subroutine psb_lz_csc_cp_from
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_mv_from
+  interface 
+    subroutine psb_lz_csc_mv_from(a,b)
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout)  :: a
+      type(psb_lz_csc_sparse_mat), intent(inout) :: b
+    end subroutine psb_lz_csc_mv_from
+  end interface
+  
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_csput_a
+  interface 
+    subroutine psb_lz_csc_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      complex(psb_dpk_), intent(in)      :: val(:)
+      integer(psb_lpk_), intent(in)             :: nz,ia(:), ja(:),&
+           &  imin,imax,jmin,jmax
+      integer(psb_ipk_), intent(out)            :: info
+      integer(psb_lpk_), intent(in), optional   :: gtl(:)
+    end subroutine psb_lz_csc_csput_a
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_base_mat_mod::psb_base_csgetptn
+  interface 
+    subroutine psb_lz_csc_csgetptn(imin,imax,a,nz,ia,ja,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_lz_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lz_csc_csgetptn
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_csgetrow
+  interface 
+    subroutine psb_lz_csc_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
+         & jmin,jmax,iren,append,nzin,rscale,cscale)
+      import
+      class(psb_lz_csc_sparse_mat), intent(in) :: a
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_lpk_), intent(out)                 :: nz
+      integer(psb_lpk_), allocatable, intent(inout)  :: ia(:), ja(:)
+      complex(psb_dpk_), allocatable,  intent(inout)    :: val(:)
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax, nzin
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lz_csc_csgetrow
+  end interface
+
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_csgetblk
+  interface 
+    subroutine psb_lz_csc_csgetblk(imin,imax,a,b,info,&
+       & jmin,jmax,iren,append,rscale,cscale)
+      import
+      class(psb_lz_csc_sparse_mat), intent(in) :: a
+      class(psb_lz_coo_sparse_mat), intent(inout) :: b
+      integer(psb_lpk_), intent(in)                  :: imin,imax
+      integer(psb_ipk_),intent(out)                  :: info
+      logical, intent(in), optional        :: append
+      integer(psb_lpk_), intent(in), optional        :: iren(:)
+      integer(psb_lpk_), intent(in), optional        :: jmin,jmax
+      logical, intent(in), optional        :: rscale,cscale
+    end subroutine psb_lz_csc_csgetblk
+  end interface
+        
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_get_diag
+  interface 
+    subroutine psb_lz_csc_get_diag(a,d,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(in) :: a
+      complex(psb_dpk_), intent(out)     :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lz_csc_get_diag
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_scal
+  interface 
+    subroutine psb_lz_csc_scal(d,a,info,side) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      complex(psb_dpk_), intent(in)      :: d(:)
+      integer(psb_ipk_), intent(out)            :: info
+      character, intent(in), optional :: side
+    end subroutine psb_lz_csc_scal
+  end interface
+  
+  !> \memberof psb_lz_csc_sparse_mat
+  !! \see psb_lz_base_mat_mod::psb_lz_base_scals
+  interface
+    subroutine psb_lz_csc_scals(d,a,info) 
+      import
+      class(psb_lz_csc_sparse_mat), intent(inout) :: a
+      complex(psb_dpk_), intent(in)      :: d
+      integer(psb_ipk_), intent(out)            :: info
+    end subroutine psb_lz_csc_scals
+  end interface
+  
+
+  
 contains 
 
   ! == ===================================
@@ -520,7 +850,7 @@ contains
     implicit none 
     class(psb_z_csc_sparse_mat), intent(in) :: a
     integer(psb_epk_) :: res
-    res = 8 
+    res = 2 * psb_sizeof_ip
     res = res + (2*psb_sizeof_dp)  * psb_size(a%val)
     res = res + psb_sizeof_ip * psb_size(a%icp)
     res = res + psb_sizeof_ip * psb_size(a%ia)
@@ -609,4 +939,126 @@ contains
 
   end subroutine z_csc_free
 
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Getters 
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================
+
+  
+  function lz_csc_is_by_cols(a) result(res)
+    implicit none 
+    class(psb_lz_csc_sparse_mat), intent(in) :: a
+    logical  :: res
+    res = .true.
+     
+  end function lz_csc_is_by_cols
+
+  !
+  ! lz
+  !
+  
+  function lz_csc_sizeof(a) result(res)
+    implicit none 
+    class(psb_lz_csc_sparse_mat), intent(in) :: a
+    integer(psb_epk_) :: res
+    res = 2*psb_sizeof_lp
+    res = res + (2*psb_sizeof_dp)  * psb_size(a%val)
+    res = res + psb_sizeof_lp * psb_size(a%icp)
+    res = res + psb_sizeof_lp * psb_size(a%ia)
+      
+  end function lz_csc_sizeof
+
+  function lz_csc_get_fmt() result(res)
+    implicit none 
+    character(len=5) :: res
+    res = 'CSC'
+  end function lz_csc_get_fmt
+  
+  function lz_csc_get_nzeros(a) result(res)
+    implicit none 
+    class(psb_lz_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+    res = a%icp(a%get_ncols()+1)-1
+  end function lz_csc_get_nzeros
+
+  function lz_csc_get_size(a) result(res)
+    implicit none 
+    class(psb_lz_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_) :: res
+
+    res = -1
+    
+    if (allocated(a%ia)) then 
+      res = size(a%ia)
+    end if
+    if (allocated(a%val)) then 
+      if (res >= 0) then 
+        res = min(res,size(a%val))
+      else 
+        res = size(a%val)
+      end if
+    end if
+
+  end function lz_csc_get_size
+
+
+
+  function  lz_csc_get_nz_col(idx,a) result(res)
+    use psb_const_mod
+    implicit none
+    
+    class(psb_lz_csc_sparse_mat), intent(in) :: a
+    integer(psb_lpk_), intent(in)                  :: idx
+    integer(psb_lpk_) :: res
+    
+    res = 0 
+ 
+    if ((1<=idx).and.(idx<=a%get_ncols())) then 
+      res = a%icp(idx+1)-a%icp(idx)
+    end if
+    
+  end function lz_csc_get_nz_col
+
+
+
+  ! == ===================================
+  !
+  !
+  !
+  ! Data management
+  !
+  !
+  !
+  !
+  !
+  ! == ===================================  
+
+
+  subroutine  lz_csc_free(a) 
+    implicit none 
+
+    class(psb_lz_csc_sparse_mat), intent(inout) :: a
+
+    if (allocated(a%icp)) deallocate(a%icp)
+    if (allocated(a%ia)) deallocate(a%ia)
+    if (allocated(a%val)) deallocate(a%val)
+    call a%set_null()
+    call a%set_nrows(0_psb_lpk_)
+    call a%set_ncols(0_psb_lpk_)
+    
+    return
+
+  end subroutine lz_csc_free
+
+
+
 end module psb_z_csc_mat_mod

base/modules/serial/psb_z_mat_mod.F90

@@ -80,7 +80,7 @@ module psb_z_mat_mod
 
   use psb_z_base_mat_mod
   use psb_z_csr_mat_mod,  only : psb_z_csr_sparse_mat, psb_lz_csr_sparse_mat
-  use psb_z_csc_mat_mod,  only : psb_z_csc_sparse_mat
+  use psb_z_csc_mat_mod,  only : psb_z_csc_sparse_mat, psb_lz_csc_sparse_mat
 
   type :: psb_zspmat_type
 

base/modules/serial/psb_z_serial_mod.f90

@@ -220,14 +220,14 @@ module psb_z_serial_mod
       type(psb_lz_csr_sparse_mat), intent(out) :: c
       integer(psb_ipk_), intent(out)          :: info
     end subroutine psb_lzcsrspspmm
-!!$    subroutine psb_zcscspspmm(a,b,c,info)
-!!$      use psb_z_mat_mod, only : psb_z_csc_sparse_mat
-!!$      import :: psb_ipk_
-!!$      implicit none 
-!!$      class(psb_z_csc_sparse_mat), intent(in) :: a,b
-!!$      type(psb_z_csc_sparse_mat), intent(out) :: c
-!!$      integer(psb_ipk_), intent(out)          :: info
-!!$    end subroutine psb_zcscspspmm
+    subroutine psb_lzcscspspmm(a,b,c,info)
+      use psb_z_mat_mod, only : psb_lz_csc_sparse_mat
+      import :: psb_ipk_
+      implicit none 
+      class(psb_lz_csc_sparse_mat), intent(in) :: a,b
+      type(psb_lz_csc_sparse_mat), intent(out) :: c
+      integer(psb_ipk_), intent(out)          :: info
+    end subroutine psb_lzcscspspmm
   end interface psb_spspmm
 
   interface psb_symbmm

base/serial/impl/psb_c_mat_impl.F90

@@ -3690,8 +3690,8 @@ subroutine psb_lc_cscnv(a,b,info,type,mold,upd,dupl)
       allocate(psb_lc_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_lc_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_lc_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_lc_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)
@@ -3791,8 +3791,8 @@ subroutine psb_lc_cscnv_ip(a,info,type,mold,dupl)
       allocate(psb_lc_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_lc_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_lc_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_lc_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)

base/serial/impl/psb_d_mat_impl.F90

@@ -3690,8 +3690,8 @@ subroutine psb_ld_cscnv(a,b,info,type,mold,upd,dupl)
       allocate(psb_ld_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_ld_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_ld_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_ld_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)
@@ -3791,8 +3791,8 @@ subroutine psb_ld_cscnv_ip(a,info,type,mold,dupl)
       allocate(psb_ld_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_ld_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_ld_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_ld_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)

base/serial/impl/psb_s_mat_impl.F90

@@ -3690,8 +3690,8 @@ subroutine psb_ls_cscnv(a,b,info,type,mold,upd,dupl)
       allocate(psb_ls_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_ls_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_ls_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_ls_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)
@@ -3791,8 +3791,8 @@ subroutine psb_ls_cscnv_ip(a,info,type,mold,dupl)
       allocate(psb_ls_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_ls_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_ls_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_ls_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)

base/serial/impl/psb_z_mat_impl.F90

@@ -3690,8 +3690,8 @@ subroutine psb_lz_cscnv(a,b,info,type,mold,upd,dupl)
       allocate(psb_lz_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_lz_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_lz_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_lz_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)
@@ -3791,8 +3791,8 @@ subroutine psb_lz_cscnv_ip(a,info,type,mold,dupl)
       allocate(psb_lz_csr_sparse_mat :: altmp, stat=info) 
     case ('COO')
       allocate(psb_lz_coo_sparse_mat :: altmp, stat=info) 
-!!$    case ('CSC')
-!!$      allocate(psb_lz_csc_sparse_mat :: altmp, stat=info) 
+    case ('CSC')
+      allocate(psb_lz_csc_sparse_mat :: altmp, stat=info) 
     case default
       info = psb_err_format_unknown_ 
       call psb_errpush(info,name,a_err=type)

base/serial/lsmmp.f90

@@ -0,0 +1,478 @@
+!   
+!                Parallel Sparse BLAS  version 3.5
+!      (C) Copyright 2006-2018
+!        Salvatore Filippone    
+!        Alfredo Buttari      
+!   
+!    Redistribution and use in source and binary forms, with or without
+!    modification, are permitted provided that the following conditions
+!    are met:
+!      1. Redistributions of source code must retain the above copyright
+!         notice, this list of conditions and the following disclaimer.
+!      2. Redistributions in binary form must reproduce the above copyright
+!         notice, this list of conditions, and the following disclaimer in the
+!         documentation and/or other materials provided with the distribution.
+!      3. The name of the PSBLAS group or the names of its contributors may
+!         not be used to endorse or promote products derived from this
+!         software without specific written permission.
+!   
+!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
+!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+!    POSSIBILITY OF SUCH DAMAGE.
+!   
+!
+! Original code adapted from:
+! == =====================================================================
+! Sparse Matrix Multiplication Package
+!
+! Randolph E. Bank and Craig C. Douglas
+!
+! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
+!
+! Compile this with the following command (or a similar one):
+!
+!     f77 -c -O smmp.f
+!
+! == =====================================================================
+subroutine lsymbmm(n, m, l, ia, ja, diaga, ib, jb, diagb,&
+     &  ic, jc, diagc,  index)
+  use psb_const_mod
+  use psb_realloc_mod
+  use psb_sort_mod, only: psb_msort
+  !
+  integer(psb_lpk_) ::       ia(*), ja(*), diaga, &
+       &  ib(*), jb(*), diagb,  diagc,  index(*)
+  integer(psb_lpk_), allocatable :: ic(:),jc(:)
+  integer(psb_lpk_) :: nze
+  integer(psb_ipk_) :: info
+
+  !
+  !       symbolic matrix multiply c=a*b
+  !
+  if (size(ic) < n+1) then 
+    write(psb_err_unit,*)&
+         &  'Called realloc in SYMBMM '
+    call psb_realloc(n+1,ic,info)
+    if (info /= psb_success_) then 
+      write(psb_err_unit,*)&
+           &  'realloc failed in SYMBMM ',info
+    end if
+  endif
+  maxlmn = max(l,m,n)
+  do i=1,maxlmn
+    index(i)=0
+  end do
+  if (diagc.eq.0) then
+    ic(1)=1
+  else
+    ic(1)=n+2
+  endif
+  minlm = min(l,m)
+  minmn = min(m,n)
+  !
+  !    main loop
+  !
+  do i=1,n
+    istart=-1
+    length=0
+    !
+    !    merge row lists
+    !
+    rowi: do jj=ia(i),ia(i+1)
+      !    a = d + ...
+      if (jj.eq.ia(i+1)) then
+        if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
+        j = i
+      else
+        j=ja(jj)
+      endif
+      !    b = d + ...
+      if (index(j).eq.0 .and. diagb.eq.1 .and. j.le.minlm)then
+        index(j)=istart
+        istart=j
+        length=length+1
+      endif
+      if ((j<1).or.(j>m)) then 
+        write(psb_err_unit,*)&
+             &  ' SymbMM: Problem with A ',i,jj,j,m
+      endif
+      do k=ib(j),ib(j+1)-1 
+        if ((jb(k)<1).or.(jb(k)>maxlmn)) then 
+          write(psb_err_unit,*)&
+               & 'Problem in SYMBMM 1:',j,k,jb(k),maxlmn
+        else
+          if(index(jb(k)).eq.0) then
+            index(jb(k))=istart
+            istart=jb(k)
+            length=length+1
+          endif
+        endif
+      end do
+    end do rowi
+
+    !
+    !   row i of jc
+    !
+    if (diagc.eq.1 .and. index(i).ne.0) length = length - 1
+    ic(i+1)=ic(i)+length
+
+    if (ic(i+1) > size(jc)) then 
+      if (n > (2*i)) then 
+        nze = max(ic(i+1), ic(i)*((n+i-1)/i))
+      else
+        nze = max(ic(i+1), nint((dble(ic(i))*(dble(n)/i)))   )
+      endif
+      call psb_realloc(nze,jc,info)
+    end if
+
+    do j= ic(i),ic(i+1)-1
+      if (diagc.eq.1 .and. istart.eq.i) then
+        istart = index(istart)
+        index(i) = 0
+      endif
+      jc(j)=istart
+      istart=index(istart)
+      index(jc(j))=0
+    end do
+    call psb_msort(jc(ic(i):ic(i)+length -1))
+    index(i) = 0
+  end do
+  return
+end subroutine lsymbmm
+! == =====================================================================
+! Sparse Matrix Multiplication Package
+!
+! Randolph E. Bank and Craig C. Douglas
+!
+! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
+!
+! Compile this with the following command (or a similar one):
+!
+!     f77 -c -O smmp.f
+!
+! == =====================================================================
+subroutine lcnumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
+     & ic, jc, diagc, c, temp)
+  !
+  use psb_const_mod
+  integer(psb_lpk_) :: ia(*), ja(*), diaga,&
+       &  ib(*), jb(*), diagb,  ic(*), jc(*), diagc 
+  !
+  complex(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
+  !
+  !       numeric matrix multiply c=a*b
+  !
+  maxlmn = max(l,m,n)
+  do i = 1,maxlmn
+    temp(i) = 0.
+  end do
+  minlm = min(l,m)
+  minln = min(l,n)
+  minmn = min(m,n)
+  !
+  !   c = a*b
+  !
+  do i = 1,n
+    rowi: do  jj = ia(i),ia(i+1)
+      !    a = d + ...
+      if (jj.eq.ia(i+1)) then
+        if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
+        j = i
+        ajj = a(i)
+      else
+        j=ja(jj)
+        ajj = a(jj)
+      endif
+      !    b = d + ...
+      if (diagb.eq.1 .and. j.le.minlm) &
+           &    temp(j) = temp(j) + ajj * b(j)
+      if ((j<1).or.(j>m)) then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: Problem with A ',i,jj,j,m
+      endif
+
+      do  k = ib(j),ib(j+1)-1
+        if((jb(k)<1).or. (jb(k) > maxlmn))  then 
+          write(psb_err_unit,*)&
+               & ' NUMBMM: jb  problem',j,k,jb(k),maxlmn
+        else
+          temp(jb(k)) = temp(jb(k)) + ajj * b(k)
+        endif
+      end do
+    end do rowi
+
+    !    c = d + ...
+    if (diagc.eq.1 .and. i.le.minln) then
+      c(i) = temp(i)
+      temp(i) = 0.
+    endif
+    !$$$        if (mod(i,100) == 1)
+    !$$$     +    write(psb_err_unit,*)
+    !$$$     ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
+    do j = ic(i),ic(i+1)-1
+      if((jc(j)<1).or. (jc(j) > maxlmn))  then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: output problem',i,j,jc(j),maxlmn
+      else
+        c(j) = temp(jc(j))
+        temp(jc(j)) = 0.
+      endif
+    end do
+  end do
+
+  return
+end subroutine lcnumbmm
+! == =====================================================================
+! Sparse Matrix Multiplication Package
+!
+! Randolph E. Bank and Craig C. Douglas
+!
+! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
+!
+! Compile this with the following command (or a similar one):
+!
+!     f77 -c -O smmp.f
+!
+! == =====================================================================
+subroutine ldnumbmm(n, m, l,  ia, ja, diaga, a, ib, jb, diagb, b,&
+     &  ic, jc, diagc, c, temp)
+  use psb_const_mod
+  !
+  integer(psb_lpk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
+       &  ic(*), jc(*), diagc 
+  !
+  real(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
+  !
+  !       numeric matrix multiply c=a*b
+  !
+  maxlmn = max(l,m,n)
+  do i = 1,maxlmn
+    temp(i) = 0.
+  end do
+  minlm = min(l,m)
+  minln = min(l,n)
+  minmn = min(m,n)
+  !
+  !   c = a*b
+  !
+  do i = 1,n
+    rowi: do jj = ia(i),ia(i+1)
+      !    a = d + ...
+      if (jj.eq.ia(i+1)) then
+        if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
+        j = i
+        ajj = a(i)
+      else
+        j=ja(jj)
+        ajj = a(jj)
+      endif
+      !    b = d + ...
+      if (diagb.eq.1 .and. j.le.minlm) &
+           &  temp(j) = temp(j) + ajj * b(j)
+      if ((j<1).or.(j>m)) then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: Problem with A ',i,jj,j,m
+      endif
+
+      do k = ib(j),ib(j+1)-1
+        if((jb(k)<1).or. (jb(k) > maxlmn))  then 
+          write(psb_err_unit,*)&
+               & ' NUMBMM: jb  problem',j,k,jb(k),maxlmn
+        else
+          temp(jb(k)) = temp(jb(k)) + ajj * b(k)
+        endif
+      end do
+    end do rowi
+
+    !    c = d + ...
+    if (diagc.eq.1 .and. i.le.minln) then
+      c(i) = temp(i)
+      temp(i) = 0.
+    endif
+    !$$$        if (mod(i,100) == 1)
+    !$$$     +    write(psb_err_unit,*)
+    !$$$     ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
+    do j = ic(i),ic(i+1)-1
+      if((jc(j)<1).or. (jc(j) > maxlmn))  then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: output problem',i,j,jc(j),maxlmn
+      else
+        c(j) = temp(jc(j))
+        temp(jc(j)) = 0.
+      endif
+    end do
+  end do
+
+  return
+end subroutine ldnumbmm
+! == =====================================================================
+! Sparse Matrix Multiplication Package
+!
+! Randolph E. Bank and Craig C. Douglas
+!
+! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
+!
+! Compile this with the following command (or a similar one):
+!
+!     f77 -c -O smmp.f
+!
+! == =====================================================================
+subroutine lsnumbmm(n, m, l,  ia, ja, diaga, a, ib, jb, diagb, b,&
+     &  ic, jc, diagc, c, temp)
+  use psb_const_mod
+  !
+  integer(psb_lpk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
+       &  ic(*), jc(*), diagc 
+  !
+  real(psb_spk_) :: a(*), b(*), c(*), temp(*),ajj
+  !
+  !       numeric matrix multiply c=a*b
+  !
+  maxlmn = max(l,m,n)
+  do i = 1,maxlmn
+    temp(i) = 0.
+  end do
+  minlm = min(l,m)
+  minln = min(l,n)
+  minmn = min(m,n)
+  !
+  !   c = a*b
+  !
+  do  i = 1,n
+    rowi: do jj = ia(i),ia(i+1)
+      !    a = d + ...
+      if (jj.eq.ia(i+1)) then
+        if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
+        j = i
+        ajj = a(i)
+      else
+        j=ja(jj)
+        ajj = a(jj)
+      endif
+      !    b = d + ...
+      if (diagb.eq.1 .and. j.le.minlm) &
+           &  temp(j) = temp(j) + ajj * b(j)
+      if ((j<1).or.(j>m)) then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: Problem with A ',i,jj,j,m
+      endif
+
+      do  k = ib(j),ib(j+1)-1
+        if((jb(k)<1).or. (jb(k) > maxlmn))  then 
+          write(psb_err_unit,*)&
+               & ' NUMBMM: jb  problem',j,k,jb(k),maxlmn
+        else
+          temp(jb(k)) = temp(jb(k)) + ajj * b(k)
+        endif
+      end do
+    end do rowi
+
+    !    c = d + ...
+    if (diagc.eq.1 .and. i.le.minln) then
+      c(i) = temp(i)
+      temp(i) = 0.
+    endif
+    !$$$        if (mod(i,100) == 1)
+    !$$$     +    write(psb_err_unit,*)
+    !$$$     ' NUMBMM: Fixing row ',i,ic(i),ic(i+1)-1
+    do j = ic(i),ic(i+1)-1
+      if((jc(j)<1).or. (jc(j) > maxlmn))  then 
+        write(psb_err_unit,*)&
+             &  ' NUMBMM: output problem',i,j,jc(j),maxlmn
+      else
+        c(j) = temp(jc(j))
+        temp(jc(j)) = 0.
+      endif
+    end do
+  end do
+
+  return
+end subroutine lsnumbmm
+! == =====================================================================
+! Sparse Matrix Multiplication Package
+!
+! Randolph E. Bank and Craig C. Douglas
+!
+! na.bank@na-net.ornl.gov and na.cdouglas@na-net.ornl.gov
+!
+! Compile this with the following command (or a similar one):
+!
+!     f77 -c -O smmp.f
+!
+! == =====================================================================
+subroutine lznumbmm(n, m, l, ia, ja, diaga, a, ib, jb, diagb, b,&
+     &  ic, jc, diagc, c, temp)
+  !
+  use psb_const_mod
+  integer(psb_lpk_) :: ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
+       &  ic(*), jc(*), diagc 
+  !
+  complex(psb_dpk_) :: a(*), b(*), c(*), temp(*),ajj
+  !
+  !       numeric matrix multiply c=a*b
+  !
+  maxlmn = max(l,m,n)
+  do i = 1,maxlmn
+    temp(i) = 0.
+  end do
+  minlm = min(l,m)
+  minln = min(l,n)
+  minmn = min(m,n)
+  !
+  !   c = a*b
+  !
+  do i = 1,n
+    rowi: do jj = ia(i),ia(i+1)
+      !    a = d + ...
+      if (jj.eq.ia(i+1)) then
+        if (diaga.eq.0 .or. i.gt.minmn) cycle rowi
+        j = i
+        ajj = a(i)
+      else
+        j=ja(jj)
+        ajj = a(jj)
+      endif
+      !    b = d + ...
+      if (diagb.eq.1 .and. j.le.minlm) &
+           & temp(j) = temp(j) + ajj * b(j)
+      if ((j<1).or.(j>m)) then 
+        write(psb_err_unit,*)&
+             &  ' NUMBMM: Problem with A ',i,jj,j,m
+      endif
+
+      do k = ib(j),ib(j+1)-1
+        if((jb(k)<1).or. (jb(k) > maxlmn))  then 
+          write(psb_err_unit,*)&
+               & ' NUMBMM: jb  problem',j,k,jb(k),maxlmn
+        else
+          temp(jb(k)) = temp(jb(k)) + ajj * b(k)
+        endif
+      end do
+    end do rowi
+
+    !    c = d + ...
+    if (diagc.eq.1 .and. i.le.minln) then
+      c(i) = temp(i)
+      temp(i) = 0.
+    endif
+    do j = ic(i),ic(i+1)-1
+      if((jc(j)<1).or. (jc(j) > maxlmn))  then 
+        write(psb_err_unit,*)&
+             & ' NUMBMM: output problem',i,j,jc(j),maxlmn
+      else
+        c(j) = temp(jc(j))
+        temp(jc(j)) = 0.
+      endif
+    end do
+  end do
+
+  return
+end subroutine lznumbmm

krylov/psb_cfcg.F90

@@ -121,9 +121,8 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err,cond
 ! =   Local data
-  type(psb_c_vect_type)  :: v, w
-  type(psb_c_vect_type), dimension(0:1) ::  d
-  complex(psb_spk_) :: alpha, tau, tau1, beta, delta
+  type(psb_c_vect_type)  :: v, w, d , q, r
+  complex(psb_spk_) :: alpha, beta, delta, gamma, theta
   real(psb_dpk_) :: derr
   integer(psb_ipk_) ::  i, idx, nc2l, it, itx, istop_, itmax_, itrace_
   integer(psb_ipk_) :: n_col, naux, err_act
@@ -131,6 +130,7 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_) :: debug_level, debug_unit
   integer(psb_ipk_) :: np, me, ictxt
   complex(psb_spk_), allocatable, target   :: aux(:)
+  complex(psb_spk_)   :: vres(3)
   character(len=20)           :: name
   type(psb_itconv_type)       :: stopdat
   character(len=*), parameter :: methdname='FCG'
@@ -192,112 +192,113 @@ subroutine psb_cfcg_vect(a,prec,b,x,eps,desc_a,info,&
   end if
 
 
-  !Assemble w, v
-
-  call psb_geasb(w,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-  call psb_geasb(v,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-
-  !Assemble d(0) and d(1)
-  call psb_geasb(d(0),&
-       & desc_a,info,&
+  !Assemble w, v, d, q, r, u 
+  call psb_geasb(w, desc_a,info,&
        & scratch=.true.,mold=x%v)
-  call psb_geasb(d(1),&
-       & desc_a,info,&
+  call psb_geasb(v, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(d, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(q, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(r, desc_a,info,&
        & scratch=.true.,mold=x%v)
 
-
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
-  itx=0
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,&
+       & a,x,b,eps,desc_a,stopdat,info)
+  itx = 0 
 
   restart: do 
     if (itx>= itmax_) exit restart 
 
-    ! w=b
-    call psb_geaxpby(cone,b,czero,w,&
-         &   desc_a,info)
+    ! r=b -Ax
+    call psb_geaxpby(cone,b,czero,r, desc_a,info)
+    if (info == psb_success_) call psb_spmm(-cone,a,x,cone,r,desc_a,info)
+    if (info /= psb_success_) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='Error during residue')
+      goto 9999
+    end if
+    
 
-      if (psb_errstatus_fatal()) then 
-        nc2l = desc_a%get_local_cols()
-        info=psb_err_alloc_request_
-        call psb_errpush(info,name,i_err=(/2*nc2l/),&
-             & a_err='complex(psb_spk_)')
-        goto 9999      
-      end if
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-   !Compute v = Ax  
 
-    call psb_spmm(cone,a,x,czero,v,desc_a,info)
+    ! Apply the preconditioner v=Pr
+    ! Compute w = Av  
+    call prec%apply(r,v,desc_a,info,work=aux)  
+    if (info == psb_success_) call psb_spmm(cone,a,v,czero,w,desc_a,info)
     if (info /= psb_success_) then
       call psb_errpush(psb_err_internal_error_,name,&
-        & a_err='Error during residue')
+           & a_err='Error during residue')
       goto 9999
     end if
 
-    !Compute w = -Ax + b
 
-    call psb_geaxpby(-cone, v, cone, w, desc_a, info)   
+    vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+    vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
 
-    !Apply the preconditioner
 
-    idx=0
+    call psb_sum(ictxt, vres(1:2))
 
-    call prec%apply(w,d(idx),desc_a,info,work=aux)  
+    alpha = vres(1)
+    beta  = vres(2)
 
-    delta = psb_gedot(d(idx), w, desc_a, info)
+    ! d = v
+    call psb_geaxpby(cone, v, czero, d, desc_a, info)   
+    ! q = w
+    call psb_geaxpby(cone, w, czero, q, desc_a, info)   
 
+    ! compute delta=beta
+    ! then 
+    ! x = x + (alpha/delta)*d
+    ! r = r - (alpha/delta)*q
 
-    !Loop
+    delta = beta
+    theta = alpha/delta
+    call psb_geaxpby(theta, d, cone, x, desc_a, info)   
+    call psb_geaxpby(-theta, q, cone, r, desc_a, info)   
 
-    if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
+    iteration: do 
 
-    if (info /= psb_success_) Then 
-      call psb_errpush(psb_err_from_subroutine_non_,name)
-      goto 9999
-    End If
+      itx = itx + 1
 
-    iteration: do 
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-      call psb_spmm(cone,a,d(idx),czero,v,desc_a,info)
+      ! Apply the preconditioner v = Pr
+      ! Compute w = Av  
+      call prec%apply(r,v,desc_a,info,work=aux)  
+      if (info == psb_success_) call psb_spmm(cone,a,v,czero,w,desc_a,info)
       if (info /= psb_success_) then
-       call psb_errpush(psb_err_internal_error_,name,&
-          & a_err='Error during residue')
-       goto 9999
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='Error during residual'); goto 9999
       end if
-      tau = psb_gedot(d(idx), v, desc_a, info) 
-
 
-      alpha = delta/tau
-      !Update solution x
-      call psb_geaxpby(alpha, d(idx), cone, x, desc_a, info)   
-      !Update residual w
-      call psb_geaxpby(-alpha, v, cone, w, desc_a, info) 
+      vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+      vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
+      vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) 
 
-      itx = itx + 1
-      idx=mod(itx ,2)
+      call psb_sum(ictxt, vres(1:3))
 
-      call d(idx)%set(czero)  
-      call prec%apply(w,d(idx),desc_a,info,work=aux)    
+      alpha = vres(1)
+      beta  = vres(2)
+      gamma = vres(3)
 
-      tau1= psb_gedot(d(idx), v, desc_a, info)
-      beta=tau1/tau
+      ! Compute d = v-(gamma/delta)*d
+      !         q = w-(gamma/delta)*q
+      theta= gamma/delta
+      call psb_geaxpby(cone, v, -theta, d, desc_a, info)   
+      call psb_geaxpby(cone, w, -theta, q , desc_a, info)   
 
-      if (idx == 1) then
-        call psb_geaxpby(-beta, d(idx - 1), cone, d(idx), desc_a, info)   
-      else
-        call psb_geaxpby(-beta, d(idx + 1), cone, d(idx), desc_a, info)         
-      endif
-    
-      delta = psb_gedot(w, d(idx), desc_a, info)
+      ! update delta
+      delta = beta - (gamma*gamma)/delta
 
-      if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
-      if (info /= psb_success_) Then 
-        call psb_errpush(psb_err_from_subroutine_non_,name)
-        goto 9999
-      End If
+      ! update u and r      
+      ! u = u + (alpha/delta)*d
+      ! r = r - (alpha/delta)*q
+      theta= alpha/delta
+      call psb_geaxpby(theta, d, cone, x, desc_a, info)   
+      call psb_geaxpby(-theta, q, cone, r, desc_a, info)   
 
     end do iteration
   end do restart

krylov/psb_dfcg.F90

@@ -121,9 +121,8 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err,cond
 ! =   Local data
-  type(psb_d_vect_type)  :: v, w
-  type(psb_d_vect_type), dimension(0:1) ::  d
-  real(psb_dpk_) :: alpha, tau, tau1, beta, delta
+  type(psb_d_vect_type)  :: v, w, d , q, r
+  real(psb_dpk_) :: alpha, beta, delta, gamma, theta
   real(psb_dpk_) :: derr
   integer(psb_ipk_) ::  i, idx, nc2l, it, itx, istop_, itmax_, itrace_
   integer(psb_ipk_) :: n_col, naux, err_act
@@ -131,6 +130,7 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_) :: debug_level, debug_unit
   integer(psb_ipk_) :: np, me, ictxt
   real(psb_dpk_), allocatable, target   :: aux(:)
+  real(psb_dpk_)   :: vres(3)
   character(len=20)           :: name
   type(psb_itconv_type)       :: stopdat
   character(len=*), parameter :: methdname='FCG'
@@ -192,112 +192,113 @@ subroutine psb_dfcg_vect(a,prec,b,x,eps,desc_a,info,&
   end if
 
 
-  !Assemble w, v
-
-  call psb_geasb(w,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-  call psb_geasb(v,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-
-  !Assemble d(0) and d(1)
-  call psb_geasb(d(0),&
-       & desc_a,info,&
+  !Assemble w, v, d, q, r, u 
+  call psb_geasb(w, desc_a,info,&
        & scratch=.true.,mold=x%v)
-  call psb_geasb(d(1),&
-       & desc_a,info,&
+  call psb_geasb(v, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(d, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(q, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(r, desc_a,info,&
        & scratch=.true.,mold=x%v)
 
-
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
-  itx=0
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,&
+       & a,x,b,eps,desc_a,stopdat,info)
+  itx = 0 
 
   restart: do 
     if (itx>= itmax_) exit restart 
 
-    ! w=b
-    call psb_geaxpby(done,b,dzero,w,&
-         &   desc_a,info)
+    ! r=b -Ax
+    call psb_geaxpby(done,b,dzero,r, desc_a,info)
+    if (info == psb_success_) call psb_spmm(-done,a,x,done,r,desc_a,info)
+    if (info /= psb_success_) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='Error during residue')
+      goto 9999
+    end if
+    
 
-      if (psb_errstatus_fatal()) then 
-        nc2l = desc_a%get_local_cols()
-        info=psb_err_alloc_request_
-        call psb_errpush(info,name,i_err=(/2*nc2l/),&
-             & a_err='real(psb_dpk_)')
-        goto 9999      
-      end if
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-   !Compute v = Ax  
 
-    call psb_spmm(done,a,x,dzero,v,desc_a,info)
+    ! Apply the preconditioner v=Pr
+    ! Compute w = Av  
+    call prec%apply(r,v,desc_a,info,work=aux)  
+    if (info == psb_success_) call psb_spmm(done,a,v,dzero,w,desc_a,info)
     if (info /= psb_success_) then
       call psb_errpush(psb_err_internal_error_,name,&
-        & a_err='Error during residue')
+           & a_err='Error during residue')
       goto 9999
     end if
 
-    !Compute w = -Ax + b
 
-    call psb_geaxpby(-done, v, done, w, desc_a, info)   
+    vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+    vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
 
-    !Apply the preconditioner
 
-    idx=0
+    call psb_sum(ictxt, vres(1:2))
 
-    call prec%apply(w,d(idx),desc_a,info,work=aux)  
+    alpha = vres(1)
+    beta  = vres(2)
 
-    delta = psb_gedot(d(idx), w, desc_a, info)
+    ! d = v
+    call psb_geaxpby(done, v, dzero, d, desc_a, info)   
+    ! q = w
+    call psb_geaxpby(done, w, dzero, q, desc_a, info)   
 
+    ! compute delta=beta
+    ! then 
+    ! x = x + (alpha/delta)*d
+    ! r = r - (alpha/delta)*q
 
-    !Loop
+    delta = beta
+    theta = alpha/delta
+    call psb_geaxpby(theta, d, done, x, desc_a, info)   
+    call psb_geaxpby(-theta, q, done, r, desc_a, info)   
 
-    if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
+    iteration: do 
 
-    if (info /= psb_success_) Then 
-      call psb_errpush(psb_err_from_subroutine_non_,name)
-      goto 9999
-    End If
+      itx = itx + 1
 
-    iteration: do 
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-      call psb_spmm(done,a,d(idx),dzero,v,desc_a,info)
+      ! Apply the preconditioner v = Pr
+      ! Compute w = Av  
+      call prec%apply(r,v,desc_a,info,work=aux)  
+      if (info == psb_success_) call psb_spmm(done,a,v,dzero,w,desc_a,info)
       if (info /= psb_success_) then
-       call psb_errpush(psb_err_internal_error_,name,&
-          & a_err='Error during residue')
-       goto 9999
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='Error during residual'); goto 9999
       end if
-      tau = psb_gedot(d(idx), v, desc_a, info) 
-
 
-      alpha = delta/tau
-      !Update solution x
-      call psb_geaxpby(alpha, d(idx), done, x, desc_a, info)   
-      !Update residual w
-      call psb_geaxpby(-alpha, v, done, w, desc_a, info) 
+      vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+      vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
+      vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) 
 
-      itx = itx + 1
-      idx=mod(itx ,2)
+      call psb_sum(ictxt, vres(1:3))
 
-      call d(idx)%set(dzero)  
-      call prec%apply(w,d(idx),desc_a,info,work=aux)    
+      alpha = vres(1)
+      beta  = vres(2)
+      gamma = vres(3)
 
-      tau1= psb_gedot(d(idx), v, desc_a, info)
-      beta=tau1/tau
+      ! Compute d = v-(gamma/delta)*d
+      !         q = w-(gamma/delta)*q
+      theta= gamma/delta
+      call psb_geaxpby(done, v, -theta, d, desc_a, info)   
+      call psb_geaxpby(done, w, -theta, q , desc_a, info)   
 
-      if (idx == 1) then
-        call psb_geaxpby(-beta, d(idx - 1), done, d(idx), desc_a, info)   
-      else
-        call psb_geaxpby(-beta, d(idx + 1), done, d(idx), desc_a, info)         
-      endif
-    
-      delta = psb_gedot(w, d(idx), desc_a, info)
+      ! update delta
+      delta = beta - (gamma*gamma)/delta
 
-      if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
-      if (info /= psb_success_) Then 
-        call psb_errpush(psb_err_from_subroutine_non_,name)
-        goto 9999
-      End If
+      ! update u and r      
+      ! u = u + (alpha/delta)*d
+      ! r = r - (alpha/delta)*q
+      theta= alpha/delta
+      call psb_geaxpby(theta, d, done, x, desc_a, info)   
+      call psb_geaxpby(-theta, q, done, r, desc_a, info)   
 
     end do iteration
   end do restart

krylov/psb_sfcg.F90

@@ -121,9 +121,8 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_spk_), Optional, Intent(out) :: err,cond
 ! =   Local data
-  type(psb_s_vect_type)  :: v, w
-  type(psb_s_vect_type), dimension(0:1) ::  d
-  real(psb_spk_) :: alpha, tau, tau1, beta, delta
+  type(psb_s_vect_type)  :: v, w, d , q, r
+  real(psb_spk_) :: alpha, beta, delta, gamma, theta
   real(psb_dpk_) :: derr
   integer(psb_ipk_) ::  i, idx, nc2l, it, itx, istop_, itmax_, itrace_
   integer(psb_ipk_) :: n_col, naux, err_act
@@ -131,6 +130,7 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_) :: debug_level, debug_unit
   integer(psb_ipk_) :: np, me, ictxt
   real(psb_spk_), allocatable, target   :: aux(:)
+  real(psb_spk_)   :: vres(3)
   character(len=20)           :: name
   type(psb_itconv_type)       :: stopdat
   character(len=*), parameter :: methdname='FCG'
@@ -192,112 +192,113 @@ subroutine psb_sfcg_vect(a,prec,b,x,eps,desc_a,info,&
   end if
 
 
-  !Assemble w, v
-
-  call psb_geasb(w,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-  call psb_geasb(v,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-
-  !Assemble d(0) and d(1)
-  call psb_geasb(d(0),&
-       & desc_a,info,&
+  !Assemble w, v, d, q, r, u 
+  call psb_geasb(w, desc_a,info,&
        & scratch=.true.,mold=x%v)
-  call psb_geasb(d(1),&
-       & desc_a,info,&
+  call psb_geasb(v, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(d, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(q, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(r, desc_a,info,&
        & scratch=.true.,mold=x%v)
 
-
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
-  itx=0
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,&
+       & a,x,b,eps,desc_a,stopdat,info)
+  itx = 0 
 
   restart: do 
     if (itx>= itmax_) exit restart 
 
-    ! w=b
-    call psb_geaxpby(sone,b,szero,w,&
-         &   desc_a,info)
+    ! r=b -Ax
+    call psb_geaxpby(sone,b,szero,r, desc_a,info)
+    if (info == psb_success_) call psb_spmm(-sone,a,x,sone,r,desc_a,info)
+    if (info /= psb_success_) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='Error during residue')
+      goto 9999
+    end if
+    
 
-      if (psb_errstatus_fatal()) then 
-        nc2l = desc_a%get_local_cols()
-        info=psb_err_alloc_request_
-        call psb_errpush(info,name,i_err=(/2*nc2l/),&
-             & a_err='real(psb_spk_)')
-        goto 9999      
-      end if
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-   !Compute v = Ax  
 
-    call psb_spmm(sone,a,x,szero,v,desc_a,info)
+    ! Apply the preconditioner v=Pr
+    ! Compute w = Av  
+    call prec%apply(r,v,desc_a,info,work=aux)  
+    if (info == psb_success_) call psb_spmm(sone,a,v,szero,w,desc_a,info)
     if (info /= psb_success_) then
       call psb_errpush(psb_err_internal_error_,name,&
-        & a_err='Error during residue')
+           & a_err='Error during residue')
       goto 9999
     end if
 
-    !Compute w = -Ax + b
 
-    call psb_geaxpby(-sone, v, sone, w, desc_a, info)   
+    vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+    vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
 
-    !Apply the preconditioner
 
-    idx=0
+    call psb_sum(ictxt, vres(1:2))
 
-    call prec%apply(w,d(idx),desc_a,info,work=aux)  
+    alpha = vres(1)
+    beta  = vres(2)
 
-    delta = psb_gedot(d(idx), w, desc_a, info)
+    ! d = v
+    call psb_geaxpby(sone, v, szero, d, desc_a, info)   
+    ! q = w
+    call psb_geaxpby(sone, w, szero, q, desc_a, info)   
 
+    ! compute delta=beta
+    ! then 
+    ! x = x + (alpha/delta)*d
+    ! r = r - (alpha/delta)*q
 
-    !Loop
+    delta = beta
+    theta = alpha/delta
+    call psb_geaxpby(theta, d, sone, x, desc_a, info)   
+    call psb_geaxpby(-theta, q, sone, r, desc_a, info)   
 
-    if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
+    iteration: do 
 
-    if (info /= psb_success_) Then 
-      call psb_errpush(psb_err_from_subroutine_non_,name)
-      goto 9999
-    End If
+      itx = itx + 1
 
-    iteration: do 
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-      call psb_spmm(sone,a,d(idx),szero,v,desc_a,info)
+      ! Apply the preconditioner v = Pr
+      ! Compute w = Av  
+      call prec%apply(r,v,desc_a,info,work=aux)  
+      if (info == psb_success_) call psb_spmm(sone,a,v,szero,w,desc_a,info)
       if (info /= psb_success_) then
-       call psb_errpush(psb_err_internal_error_,name,&
-          & a_err='Error during residue')
-       goto 9999
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='Error during residual'); goto 9999
       end if
-      tau = psb_gedot(d(idx), v, desc_a, info) 
-
 
-      alpha = delta/tau
-      !Update solution x
-      call psb_geaxpby(alpha, d(idx), sone, x, desc_a, info)   
-      !Update residual w
-      call psb_geaxpby(-alpha, v, sone, w, desc_a, info) 
+      vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+      vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
+      vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) 
 
-      itx = itx + 1
-      idx=mod(itx ,2)
+      call psb_sum(ictxt, vres(1:3))
 
-      call d(idx)%set(szero)  
-      call prec%apply(w,d(idx),desc_a,info,work=aux)    
+      alpha = vres(1)
+      beta  = vres(2)
+      gamma = vres(3)
 
-      tau1= psb_gedot(d(idx), v, desc_a, info)
-      beta=tau1/tau
+      ! Compute d = v-(gamma/delta)*d
+      !         q = w-(gamma/delta)*q
+      theta= gamma/delta
+      call psb_geaxpby(sone, v, -theta, d, desc_a, info)   
+      call psb_geaxpby(sone, w, -theta, q , desc_a, info)   
 
-      if (idx == 1) then
-        call psb_geaxpby(-beta, d(idx - 1), sone, d(idx), desc_a, info)   
-      else
-        call psb_geaxpby(-beta, d(idx + 1), sone, d(idx), desc_a, info)         
-      endif
-    
-      delta = psb_gedot(w, d(idx), desc_a, info)
+      ! update delta
+      delta = beta - (gamma*gamma)/delta
 
-      if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
-      if (info /= psb_success_) Then 
-        call psb_errpush(psb_err_from_subroutine_non_,name)
-        goto 9999
-      End If
+      ! update u and r      
+      ! u = u + (alpha/delta)*d
+      ! r = r - (alpha/delta)*q
+      theta= alpha/delta
+      call psb_geaxpby(theta, d, sone, x, desc_a, info)   
+      call psb_geaxpby(-theta, q, sone, r, desc_a, info)   
 
     end do iteration
   end do restart

krylov/psb_zfcg.F90

@@ -121,9 +121,8 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_), Optional, Intent(out)       :: iter
   real(psb_dpk_), Optional, Intent(out) :: err,cond
 ! =   Local data
-  type(psb_z_vect_type)  :: v, w
-  type(psb_z_vect_type), dimension(0:1) ::  d
-  complex(psb_dpk_) :: alpha, tau, tau1, beta, delta
+  type(psb_z_vect_type)  :: v, w, d , q, r
+  complex(psb_dpk_) :: alpha, beta, delta, gamma, theta
   real(psb_dpk_) :: derr
   integer(psb_ipk_) ::  i, idx, nc2l, it, itx, istop_, itmax_, itrace_
   integer(psb_ipk_) :: n_col, naux, err_act
@@ -131,6 +130,7 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
   integer(psb_ipk_) :: debug_level, debug_unit
   integer(psb_ipk_) :: np, me, ictxt
   complex(psb_dpk_), allocatable, target   :: aux(:)
+  complex(psb_dpk_)   :: vres(3)
   character(len=20)           :: name
   type(psb_itconv_type)       :: stopdat
   character(len=*), parameter :: methdname='FCG'
@@ -192,112 +192,113 @@ subroutine psb_zfcg_vect(a,prec,b,x,eps,desc_a,info,&
   end if
 
 
-  !Assemble w, v
-
-  call psb_geasb(w,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-  call psb_geasb(v,&
-       & desc_a,info,&
-       & scratch=.true.,mold=b%v)
-
-  !Assemble d(0) and d(1)
-  call psb_geasb(d(0),&
-       & desc_a,info,&
+  !Assemble w, v, d, q, r, u 
+  call psb_geasb(w, desc_a,info,&
        & scratch=.true.,mold=x%v)
-  call psb_geasb(d(1),&
-       & desc_a,info,&
+  call psb_geasb(v, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(d, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(q, desc_a,info,&
+       & scratch=.true.,mold=x%v)
+  call psb_geasb(r, desc_a,info,&
        & scratch=.true.,mold=x%v)
 
-
-    call psb_init_conv(methdname,istop_,itrace_,itmax_,a,x,b,eps,desc_a,stopdat,info)
-  itx=0
+  call psb_init_conv(methdname,istop_,itrace_,itmax_,&
+       & a,x,b,eps,desc_a,stopdat,info)
+  itx = 0 
 
   restart: do 
     if (itx>= itmax_) exit restart 
 
-    ! w=b
-    call psb_geaxpby(zone,b,zzero,w,&
-         &   desc_a,info)
+    ! r=b -Ax
+    call psb_geaxpby(zone,b,zzero,r, desc_a,info)
+    if (info == psb_success_) call psb_spmm(-zone,a,x,zone,r,desc_a,info)
+    if (info /= psb_success_) then
+      call psb_errpush(psb_err_internal_error_,name,&
+           & a_err='Error during residue')
+      goto 9999
+    end if
+    
 
-      if (psb_errstatus_fatal()) then 
-        nc2l = desc_a%get_local_cols()
-        info=psb_err_alloc_request_
-        call psb_errpush(info,name,i_err=(/2*nc2l/),&
-             & a_err='complex(psb_dpk_)')
-        goto 9999      
-      end if
+    if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-   !Compute v = Ax  
 
-    call psb_spmm(zone,a,x,zzero,v,desc_a,info)
+    ! Apply the preconditioner v=Pr
+    ! Compute w = Av  
+    call prec%apply(r,v,desc_a,info,work=aux)  
+    if (info == psb_success_) call psb_spmm(zone,a,v,zzero,w,desc_a,info)
     if (info /= psb_success_) then
       call psb_errpush(psb_err_internal_error_,name,&
-        & a_err='Error during residue')
+           & a_err='Error during residue')
       goto 9999
     end if
 
-    !Compute w = -Ax + b
 
-    call psb_geaxpby(-zone, v, zone, w, desc_a, info)   
+    vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+    vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
 
-    !Apply the preconditioner
 
-    idx=0
+    call psb_sum(ictxt, vres(1:2))
 
-    call prec%apply(w,d(idx),desc_a,info,work=aux)  
+    alpha = vres(1)
+    beta  = vres(2)
 
-    delta = psb_gedot(d(idx), w, desc_a, info)
+    ! d = v
+    call psb_geaxpby(zone, v, zzero, d, desc_a, info)   
+    ! q = w
+    call psb_geaxpby(zone, w, zzero, q, desc_a, info)   
 
+    ! compute delta=beta
+    ! then 
+    ! x = x + (alpha/delta)*d
+    ! r = r - (alpha/delta)*q
 
-    !Loop
+    delta = beta
+    theta = alpha/delta
+    call psb_geaxpby(theta, d, zone, x, desc_a, info)   
+    call psb_geaxpby(-theta, q, zone, r, desc_a, info)   
 
-    if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
+    iteration: do 
 
-    if (info /= psb_success_) Then 
-      call psb_errpush(psb_err_from_subroutine_non_,name)
-      goto 9999
-    End If
+      itx = itx + 1
 
-    iteration: do 
+      if (psb_check_conv(methdname,itx,x,r,desc_a,stopdat,info)) exit restart
 
-      call psb_spmm(zone,a,d(idx),zzero,v,desc_a,info)
+      ! Apply the preconditioner v = Pr
+      ! Compute w = Av  
+      call prec%apply(r,v,desc_a,info,work=aux)  
+      if (info == psb_success_) call psb_spmm(zone,a,v,zzero,w,desc_a,info)
       if (info /= psb_success_) then
-       call psb_errpush(psb_err_internal_error_,name,&
-          & a_err='Error during residue')
-       goto 9999
+        call psb_errpush(psb_err_internal_error_,name,&
+             & a_err='Error during residual'); goto 9999
       end if
-      tau = psb_gedot(d(idx), v, desc_a, info) 
-
 
-      alpha = delta/tau
-      !Update solution x
-      call psb_geaxpby(alpha, d(idx), zone, x, desc_a, info)   
-      !Update residual w
-      call psb_geaxpby(-alpha, v, zone, w, desc_a, info) 
+      vres(1) = psb_gedot(r, v, desc_a, info, global = .false.) 
+      vres(2) = psb_gedot(w, v, desc_a, info, global = .false.) 
+      vres(3) = psb_gedot(q, v, desc_a, info, global = .false.) 
 
-      itx = itx + 1
-      idx=mod(itx ,2)
+      call psb_sum(ictxt, vres(1:3))
 
-      call d(idx)%set(zzero)  
-      call prec%apply(w,d(idx),desc_a,info,work=aux)    
+      alpha = vres(1)
+      beta  = vres(2)
+      gamma = vres(3)
 
-      tau1= psb_gedot(d(idx), v, desc_a, info)
-      beta=tau1/tau
+      ! Compute d = v-(gamma/delta)*d
+      !         q = w-(gamma/delta)*q
+      theta= gamma/delta
+      call psb_geaxpby(zone, v, -theta, d, desc_a, info)   
+      call psb_geaxpby(zone, w, -theta, q , desc_a, info)   
 
-      if (idx == 1) then
-        call psb_geaxpby(-beta, d(idx - 1), zone, d(idx), desc_a, info)   
-      else
-        call psb_geaxpby(-beta, d(idx + 1), zone, d(idx), desc_a, info)         
-      endif
-    
-      delta = psb_gedot(w, d(idx), desc_a, info)
+      ! update delta
+      delta = beta - (gamma*gamma)/delta
 
-      if (psb_check_conv(methdname,itx ,x,w,desc_a,stopdat,info)) exit restart
-      if (info /= psb_success_) Then 
-        call psb_errpush(psb_err_from_subroutine_non_,name)
-        goto 9999
-      End If
+      ! update u and r      
+      ! u = u + (alpha/delta)*d
+      ! r = r - (alpha/delta)*q
+      theta= alpha/delta
+      call psb_geaxpby(theta, d, zone, x, desc_a, info)   
+      call psb_geaxpby(-theta, q, zone, r, desc_a, info)   
 
     end do iteration
   end do restart

util/psb_c_hbio_impl.f90

@@ -361,3 +361,338 @@ subroutine chb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
   write(psb_err_unit,*) 'Error while opening ',filename
   return
 end subroutine chb_write
+
+
+subroutine lchb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+  use psb_base_mod
+  implicit none
+  type(psb_lcspmat_type), intent(out)    :: a
+  integer(psb_ipk_), intent(out)                   :: iret
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  complex(psb_spk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+  character(len=72), optional, intent(out) :: mtitle
+
+  character  :: rhstype*3,type*3,key*8
+  character(len=72) :: mtitle_
+  character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
+  integer(psb_lpk_) :: indcrd,  ptrcrd, totcrd,&
+       & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
+  type(psb_lc_csc_sparse_mat) :: acsc
+  type(psb_lc_coo_sparse_mat) :: acoo
+  integer(psb_ipk_) :: ircode, infile, info
+  integer(psb_lpk_) :: i,nzr
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  iret = 0
+  ircode = 0
+  if (present(filename)) then
+    if (filename == '-') then 
+      infile=5
+    else
+      if (present(iunit)) then 
+        infile=iunit
+      else
+        infile=99
+      endif
+      open(infile,file=filename, status='OLD', err=901, action='READ')
+    endif
+  else 
+    if (present(iunit)) then 
+      infile=iunit
+    else
+      infile=5
+    endif
+  endif
+
+  read (infile,fmt=fmt10) mtitle_,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix
+
+  call acsc%allocate(nrow,ncol,nnzero)
+  if (ircode /= 0 ) then 
+    write(psb_err_unit,*) 'Memory allocation failed'
+    goto 993
+  end if
+
+  if (present(mtitle)) mtitle=mtitle_
+
+
+  if (psb_tolower(type(1:1)) == 'c') then 
+    if (psb_tolower(type(2:2)) == 'u') then 
+
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+      call a%mv_from(acsc)
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+    else if (psb_tolower(type(2:2)) == 's') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = acoo%val(i)
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else if (psb_tolower(type(2:2)) == 'h') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = conjg(acoo%val(i))
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else
+      write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+      iret=904
+    end if
+  else
+    write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+    iret=904
+  end if
+
+  call a%cscnv(ircode,type='csr')
+  if (infile /= 5) close(infile)
+
+  return 
+
+  ! open failed
+901 iret=901
+  write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
+  return
+902 iret=902
+  write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
+  return
+993 iret=993
+  write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
+  return
+end subroutine lchb_read
+
+subroutine lchb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+  use psb_base_mod
+  implicit none
+  type(psb_lcspmat_type), intent(in), target :: a
+  integer(psb_ipk_), intent(out)        :: iret
+  character(len=*), optional, intent(in) :: mtitle
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  character(len=*), optional, intent(in) :: key
+  complex(psb_spk_), optional             :: rhs(:), g(:), x(:)
+  integer(psb_ipk_) :: iout
+
+  character(len=*), parameter::  ptrfmt='(10I8)',indfmt='(10I8)'
+  integer(psb_ipk_), parameter :: jptr=10,jind=10
+  character(len=*), parameter::  valfmt='(4E20.12)',rhsfmt='(4E20.12)'
+  integer(psb_ipk_), parameter :: jval=2,jrhs=2
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  type(psb_lc_csc_sparse_mat)  ::  acsc
+  character(len=72) :: mtitle_
+  character(len=8)  :: key_
+
+  character  :: rhstype*3,type*3
+
+  integer(psb_lpk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
+       & nrow,ncol,nnzero, neltvl, nrhs, nrhsix
+
+  iret = 0
+
+  if (present(filename)) then 
+    if (filename == '-') then 
+      iout=6
+    else
+      if (present(iunit)) then 
+        iout = iunit
+      else
+        iout=99
+      endif
+      open(iout,file=filename, err=901, action='WRITE')
+    endif
+  else 
+    if (present(iunit)) then 
+      iout = iunit   
+    else
+      iout=6
+    endif
+  endif
+
+  if (present(mtitle)) then 
+    mtitle_ = mtitle
+  else
+    mtitle_ = 'Temporary PSBLAS title '
+  endif
+  if (present(key)) then 
+    key_ = key
+  else
+    key_ = 'PSBMAT00'
+  endif
+
+
+  call acsc%cp_from_fmt(a%a, iret)
+  if (iret /= 0) return
+
+
+  nrow   = acsc%get_nrows()
+  ncol   = acsc%get_ncols()
+  nnzero = acsc%get_nzeros()
+
+  neltvl = 0 
+
+  ptrcrd = (ncol+1)/jptr
+  if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
+  indcrd = nnzero/jind
+  if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
+  valcrd = nnzero/jval
+  if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
+  rhstype = ''
+  if (present(rhs)) then 
+    if (size(rhs)<nrow) then 
+      rhscrd = 0
+    else
+      rhscrd = nrow/jrhs
+      if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
+    endif
+    nrhs = 1
+    rhstype(1:1) = 'F'
+  else
+    rhscrd = 0
+    nrhs   = 0
+  end if
+  totcrd = ptrcrd + indcrd + valcrd + rhscrd
+
+  nrhsix  = nrhs*nrow
+
+  if (present(g)) then 
+    rhstype(2:2) = 'G'
+  end if
+  if (present(x)) then 
+    rhstype(3:3) = 'X'
+  end if
+  type    = 'CUA'
+
+  write (iout,fmt=fmt10) mtitle_,key_,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       &  type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) write (iout,fmt=fmt11) rhstype,nrhs,nrhsix
+  write (iout,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+  write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+  if (valcrd > 0) write (iout,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+  if (rhscrd > 0) write (iout,fmt=rhsfmt) (rhs(i),i=1,nrow)
+  if (present(g).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (g(i),i=1,nrow)
+  if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
+
+
+
+
+  if (iout /= 6) close(iout)
+
+
+  return
+
+901 continue 
+  iret=901
+  write(psb_err_unit,*) 'Error while opening ',filename
+  return
+end subroutine lchb_write

util/psb_d_hbio_impl.f90

@@ -314,3 +314,290 @@ subroutine dhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
   write(psb_err_unit,*) 'Error while opening ',filename
   return
 end subroutine dhb_write
+
+subroutine ldhb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+  use psb_base_mod
+  implicit none
+  type(psb_ldspmat_type), intent(out)    :: a
+  integer(psb_ipk_), intent(out)                   :: iret
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  real(psb_dpk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+  character(len=72), optional, intent(out) :: mtitle
+
+  character  :: rhstype*3,type*3,key*8
+  character(len=72) :: mtitle_
+  character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
+  integer(psb_lpk_) :: indcrd,  ptrcrd, totcrd,&
+       & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
+  type(psb_ld_csc_sparse_mat) :: acsc
+  type(psb_ld_coo_sparse_mat) :: acoo
+  integer(psb_ipk_) :: ircode, infile, info
+  integer(psb_lpk_) :: i,nzr
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  iret = 0
+  ircode = 0
+  if (present(filename)) then
+    if (filename == '-') then 
+      infile=5
+    else
+      if (present(iunit)) then 
+        infile=iunit
+      else
+        infile=99
+      endif
+      open(infile,file=filename, status='OLD', err=901, action='READ')
+    endif
+  else 
+    if (present(iunit)) then 
+      infile=iunit
+    else
+      infile=5
+    endif
+  endif
+
+  read (infile,fmt=fmt10) mtitle_,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix
+
+  call acsc%allocate(nrow,ncol,nnzero)
+  if (ircode /= 0 ) then 
+    write(psb_err_unit,*) 'Memory allocation failed'
+    goto 993
+  end if
+
+  if (present(mtitle)) mtitle=mtitle_
+
+
+  if (psb_tolower(type(1:1)) == 'r') then 
+    if (psb_tolower(type(2:2)) == 'u') then 
+
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+      call a%mv_from(acsc)
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+    else if (psb_tolower(type(2:2)) == 's') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = acoo%val(i)
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else
+      write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+      iret=904
+    end if
+  else
+    write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+    iret=904
+  end if
+
+  call a%cscnv(ircode,type='csr')
+  if (infile /= 5) close(infile)
+
+  return 
+
+  ! open failed
+901 iret=901
+  write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
+  return
+902 iret=902
+  write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
+  return
+993 iret=993
+  write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
+  return
+end subroutine ldhb_read
+
+subroutine ldhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+  use psb_base_mod
+  implicit none
+  type(psb_ldspmat_type), intent(in), target :: a
+  integer(psb_ipk_), intent(out)        :: iret
+  character(len=*), optional, intent(in) :: mtitle
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  character(len=*), optional, intent(in) :: key
+  real(psb_dpk_), optional             :: rhs(:), g(:), x(:)
+  integer(psb_ipk_) :: iout
+
+  character(len=*), parameter::  ptrfmt='(10I8)',indfmt='(10I8)'
+  integer(psb_ipk_), parameter :: jptr=10,jind=10
+  character(len=*), parameter::  valfmt='(4E20.12)',rhsfmt='(4E20.12)'
+  integer(psb_ipk_), parameter :: jval=4,jrhs=4
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  type(psb_ld_csc_sparse_mat)  ::  acsc
+  character(len=72) :: mtitle_
+  character(len=8)  :: key_
+
+  character  :: rhstype*3,type*3
+
+  integer(psb_lpk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
+       & nrow,ncol,nnzero, neltvl, nrhs, nrhsix
+
+  iret = 0
+
+  if (present(filename)) then 
+    if (filename == '-') then 
+      iout=6
+    else
+      if (present(iunit)) then 
+        iout = iunit
+      else
+        iout=99
+      endif
+      open(iout,file=filename, err=901, action='WRITE')
+    endif
+  else 
+    if (present(iunit)) then 
+      iout = iunit   
+    else
+      iout=6
+    endif
+  endif
+
+  if (present(mtitle)) then 
+    mtitle_ = mtitle
+  else
+    mtitle_ = 'Temporary PSBLAS title '
+  endif
+  if (present(key)) then 
+    key_ = key
+  else
+    key_ = 'PSBMAT00'
+  endif
+
+
+  call acsc%cp_from_fmt(a%a, iret)
+  if (iret /= 0) return
+
+
+  nrow   = acsc%get_nrows()
+  ncol   = acsc%get_ncols()
+  nnzero = acsc%get_nzeros()
+
+  neltvl = 0 
+
+  ptrcrd = (ncol+1)/jptr
+  if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
+  indcrd = nnzero/jind
+  if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
+  valcrd = nnzero/jval
+  if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
+  rhstype = ''
+  if (present(rhs)) then 
+    if (size(rhs)<nrow) then 
+      rhscrd = 0
+    else
+      rhscrd = nrow/jrhs
+      if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
+    endif
+    nrhs = 1
+    rhstype(1:1) = 'F'
+  else
+    rhscrd = 0
+    nrhs   = 0
+  end if
+  totcrd = ptrcrd + indcrd + valcrd + rhscrd
+
+  nrhsix  = nrhs*nrow
+
+  if (present(g)) then 
+    rhstype(2:2) = 'G'
+  end if
+  if (present(x)) then 
+    rhstype(3:3) = 'X'
+  end if
+  type    = 'RUA'
+
+  write (iout,fmt=fmt10) mtitle_,key_,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       &  type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) write (iout,fmt=fmt11) rhstype,nrhs,nrhsix
+  write (iout,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+  write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+  if (valcrd > 0) write (iout,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+  if (rhscrd > 0) write (iout,fmt=rhsfmt) (rhs(i),i=1,nrow)
+  if (present(g).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (g(i),i=1,nrow)
+  if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
+
+
+
+
+  if (iout /= 6) close(iout)
+
+
+  return
+
+901 continue 
+  iret=901
+  write(psb_err_unit,*) 'Error while opening ',filename
+  return
+end subroutine ldhb_write

util/psb_hbio_mod.f90

@@ -33,7 +33,10 @@ module psb_hbio_mod
 
   use psb_base_mod, only :  psb_ipk_, psb_spk_, psb_dpk_,&
        & psb_sspmat_type, psb_cspmat_type, &
-       & psb_dspmat_type, psb_zspmat_type
+       & psb_dspmat_type, psb_zspmat_type, &
+       & psb_lsspmat_type, psb_lcspmat_type, &
+       & psb_ldspmat_type, psb_lzspmat_type
+
 
 
   public hb_read, hb_write
@@ -78,6 +81,46 @@ module psb_hbio_mod
       complex(psb_dpk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
       character(len=72), optional, intent(out) :: mtitle
     end subroutine zhb_read
+    subroutine lshb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+      import :: psb_lsspmat_type, psb_spk_, psb_ipk_
+      implicit none
+      type(psb_lsspmat_type), intent(out)     :: a
+      integer(psb_ipk_), intent(out)                   :: iret
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      real(psb_spk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+      character(len=72), optional, intent(out) :: mtitle
+    end subroutine lshb_read
+    subroutine ldhb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+      import :: psb_ldspmat_type, psb_dpk_, psb_ipk_
+      implicit none
+      type(psb_ldspmat_type), intent(out)     :: a
+      integer(psb_ipk_), intent(out)                   :: iret
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      real(psb_dpk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+      character(len=72), optional, intent(out) :: mtitle
+    end subroutine ldhb_read
+    subroutine lchb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+      import :: psb_lcspmat_type, psb_spk_, psb_ipk_
+      implicit none
+      type(psb_lcspmat_type), intent(out)     :: a
+      integer(psb_ipk_), intent(out)                   :: iret
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      complex(psb_spk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+      character(len=72), optional, intent(out) :: mtitle
+    end subroutine lchb_read
+    subroutine lzhb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+      import :: psb_lzspmat_type, psb_dpk_, psb_ipk_
+      implicit none
+      type(psb_lzspmat_type), intent(out)     :: a
+      integer(psb_ipk_), intent(out)                   :: iret
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      complex(psb_dpk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+      character(len=72), optional, intent(out) :: mtitle
+    end subroutine lzhb_read
   end interface
 
   interface hb_write
@@ -125,6 +168,50 @@ module psb_hbio_mod
       character(len=*), optional, intent(in) :: key
       complex(psb_dpk_), optional             :: rhs(:), g(:), x(:)
     end subroutine zhb_write
+    subroutine lshb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+      import :: psb_lsspmat_type, psb_spk_, psb_ipk_
+      implicit none
+      type(psb_lsspmat_type), intent(inout)  :: a
+      integer(psb_ipk_), intent(out)        :: iret
+      character(len=*), optional, intent(in) :: mtitle
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      character(len=*), optional, intent(in) :: key
+      real(psb_spk_), optional             :: rhs(:), g(:), x(:)
+    end subroutine lshb_write
+    subroutine ldhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+      import :: psb_ldspmat_type, psb_dpk_, psb_ipk_
+      implicit none
+      type(psb_ldspmat_type), intent(inout)  :: a
+      integer(psb_ipk_), intent(out)        :: iret
+      character(len=*), optional, intent(in) :: mtitle
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      character(len=*), optional, intent(in) :: key
+      real(psb_dpk_), optional             :: rhs(:), g(:), x(:)
+    end subroutine ldhb_write
+    subroutine lchb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+      import :: psb_lcspmat_type, psb_spk_, psb_ipk_
+      implicit none
+      type(psb_lcspmat_type), intent(inout)  :: a
+      integer(psb_ipk_), intent(out)        :: iret
+      character(len=*), optional, intent(in) :: mtitle
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      character(len=*), optional, intent(in) :: key
+      complex(psb_spk_), optional             :: rhs(:), g(:), x(:)
+    end subroutine lchb_write
+    subroutine lzhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+      import :: psb_lzspmat_type, psb_dpk_, psb_ipk_
+      implicit none
+      type(psb_lzspmat_type), intent(inout)  :: a
+      integer(psb_ipk_), intent(out)        :: iret
+      character(len=*), optional, intent(in) :: mtitle
+      integer(psb_ipk_), optional, intent(in)          :: iunit
+      character(len=*), optional, intent(in) :: filename
+      character(len=*), optional, intent(in) :: key
+      complex(psb_dpk_), optional             :: rhs(:), g(:), x(:)
+    end subroutine lzhb_write
   end interface
 
 end module psb_hbio_mod

util/psb_s_hbio_impl.f90

@@ -314,3 +314,290 @@ subroutine shb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
   write(psb_err_unit,*) 'Error while opening ',filename
   return
 end subroutine shb_write
+
+subroutine lshb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+  use psb_base_mod
+  implicit none
+  type(psb_lsspmat_type), intent(out)    :: a
+  integer(psb_ipk_), intent(out)                   :: iret
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  real(psb_spk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+  character(len=72), optional, intent(out) :: mtitle
+
+  character  :: rhstype*3,type*3,key*8
+  character(len=72) :: mtitle_
+  character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
+  integer(psb_lpk_) :: indcrd,  ptrcrd, totcrd,&
+       & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
+  type(psb_ls_csc_sparse_mat) :: acsc
+  type(psb_ls_coo_sparse_mat) :: acoo
+  integer(psb_ipk_) :: ircode, infile, info
+  integer(psb_lpk_) :: i,nzr
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  iret = 0
+  ircode = 0
+  if (present(filename)) then
+    if (filename == '-') then 
+      infile=5
+    else
+      if (present(iunit)) then 
+        infile=iunit
+      else
+        infile=99
+      endif
+      open(infile,file=filename, status='OLD', err=901, action='READ')
+    endif
+  else 
+    if (present(iunit)) then 
+      infile=iunit
+    else
+      infile=5
+    endif
+  endif
+
+  read (infile,fmt=fmt10) mtitle_,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix
+
+  call acsc%allocate(nrow,ncol,nnzero)
+  if (ircode /= 0 ) then 
+    write(psb_err_unit,*) 'Memory allocation failed'
+    goto 993
+  end if
+
+  if (present(mtitle)) mtitle=mtitle_
+
+
+  if (psb_tolower(type(1:1)) == 'r') then 
+    if (psb_tolower(type(2:2)) == 'u') then 
+
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+      call a%mv_from(acsc)
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+    else if (psb_tolower(type(2:2)) == 's') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = acoo%val(i)
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else
+      write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+      iret=904
+    end if
+  else
+    write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+    iret=904
+  end if
+
+  call a%cscnv(ircode,type='csr')
+  if (infile /= 5) close(infile)
+
+  return 
+
+  ! open failed
+901 iret=901
+  write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
+  return
+902 iret=902
+  write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
+  return
+993 iret=993
+  write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
+  return
+end subroutine lshb_read
+
+subroutine lshb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+  use psb_base_mod
+  implicit none
+  type(psb_lsspmat_type), intent(in), target :: a
+  integer(psb_ipk_), intent(out)        :: iret
+  character(len=*), optional, intent(in) :: mtitle
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  character(len=*), optional, intent(in) :: key
+  real(psb_spk_), optional             :: rhs(:), g(:), x(:)
+  integer(psb_ipk_) :: iout
+
+  character(len=*), parameter::  ptrfmt='(10I8)',indfmt='(10I8)'
+  integer(psb_ipk_), parameter :: jptr=10,jind=10
+  character(len=*), parameter::  valfmt='(4E20.12)',rhsfmt='(4E20.12)'
+  integer(psb_ipk_), parameter :: jval=4,jrhs=4
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  type(psb_ls_csc_sparse_mat)  ::  acsc
+  character(len=72) :: mtitle_
+  character(len=8)  :: key_
+
+  character  :: rhstype*3,type*3
+
+  integer(psb_ipk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
+       & nrow,ncol,nnzero, neltvl, nrhs, nrhsix
+
+  iret = 0
+
+  if (present(filename)) then 
+    if (filename == '-') then 
+      iout=6
+    else
+      if (present(iunit)) then 
+        iout = iunit
+      else
+        iout=99
+      endif
+      open(iout,file=filename, err=901, action='WRITE')
+    endif
+  else 
+    if (present(iunit)) then 
+      iout = iunit   
+    else
+      iout=6
+    endif
+  endif
+
+  if (present(mtitle)) then 
+    mtitle_ = mtitle
+  else
+    mtitle_ = 'Temporary PSBLAS title '
+  endif
+  if (present(key)) then 
+    key_ = key
+  else
+    key_ = 'PSBMAT00'
+  endif
+
+
+  call acsc%cp_from_fmt(a%a, iret)
+  if (iret /= 0) return
+
+
+  nrow   = acsc%get_nrows()
+  ncol   = acsc%get_ncols()
+  nnzero = acsc%get_nzeros()
+
+  neltvl = 0 
+
+  ptrcrd = (ncol+1)/jptr
+  if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
+  indcrd = nnzero/jind
+  if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
+  valcrd = nnzero/jval
+  if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
+  rhstype = ''
+  if (present(rhs)) then 
+    if (size(rhs)<nrow) then 
+      rhscrd = 0
+    else
+      rhscrd = nrow/jrhs
+      if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
+    endif
+    nrhs = 1
+    rhstype(1:1) = 'F'
+  else
+    rhscrd = 0
+    nrhs   = 0
+  end if
+  totcrd = ptrcrd + indcrd + valcrd + rhscrd
+
+  nrhsix  = nrhs*nrow
+
+  if (present(g)) then 
+    rhstype(2:2) = 'G'
+  end if
+  if (present(x)) then 
+    rhstype(3:3) = 'X'
+  end if
+  type    = 'RUA'
+
+  write (iout,fmt=fmt10) mtitle_,key_,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       &  type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) write (iout,fmt=fmt11) rhstype,nrhs,nrhsix
+  write (iout,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+  write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+  if (valcrd > 0) write (iout,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+  if (rhscrd > 0) write (iout,fmt=rhsfmt) (rhs(i),i=1,nrow)
+  if (present(g).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (g(i),i=1,nrow)
+  if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
+
+
+
+
+  if (iout /= 6) close(iout)
+
+
+  return
+
+901 continue 
+  iret=901
+  write(psb_err_unit,*) 'Error while opening ',filename
+  return
+end subroutine lshb_write

util/psb_z_hbio_impl.f90

@@ -349,6 +349,336 @@ subroutine zhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
   if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
 
 
+  if (iout /= 6) close(iout)
+
+
+  return
+
+901 continue 
+  iret=901
+  write(psb_err_unit,*) 'Error while opening ',filename
+  return
+end subroutine zhb_write
+
+subroutine lzhb_read(a, iret, iunit, filename,b,g,x,mtitle)   
+  use psb_base_mod
+  implicit none
+  type(psb_lzspmat_type), intent(out)    :: a
+  integer(psb_ipk_), intent(out)                   :: iret
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  complex(psb_dpk_), optional, allocatable, intent(out)  :: b(:,:), g(:,:), x(:,:) 
+  character(len=72), optional, intent(out) :: mtitle
+
+  character  :: rhstype*3,type*3,key*8
+  character(len=72) :: mtitle_
+  character indfmt*16,ptrfmt*16,rhsfmt*20,valfmt*20
+  integer(psb_lpk_) :: indcrd,  ptrcrd, totcrd,&
+       & valcrd, rhscrd, nrow, ncol, nnzero, neltvl, nrhs, nrhsix
+  type(psb_lz_csc_sparse_mat) :: acsc
+  type(psb_lz_coo_sparse_mat) :: acoo
+  integer(psb_ipk_) :: ircode, infile, info
+  integer(psb_lpk_) :: i,nzr
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  iret = 0
+  ircode = 0
+  if (present(filename)) then
+    if (filename == '-') then 
+      infile=5
+    else
+      if (present(iunit)) then 
+        infile=iunit
+      else
+        infile=99
+      endif
+      open(infile,file=filename, status='OLD', err=901, action='READ')
+    endif
+  else 
+    if (present(iunit)) then 
+      infile=iunit
+    else
+      infile=5
+    endif
+  endif
+
+  read (infile,fmt=fmt10) mtitle_,key,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       & type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) read(infile,fmt=fmt11)rhstype,nrhs,nrhsix
+
+  call acsc%allocate(nrow,ncol,nnzero)
+  if (ircode /= 0 ) then 
+    write(psb_err_unit,*) 'Memory allocation failed'
+    goto 993
+  end if
+
+  if (present(mtitle)) mtitle=mtitle_
+
+
+  if (psb_tolower(type(1:1)) == 'c') then 
+    if (psb_tolower(type(2:2)) == 'u') then 
+
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+      call a%mv_from(acsc)
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+    else if (psb_tolower(type(2:2)) == 's') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = acoo%val(i)
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else if (psb_tolower(type(2:2)) == 'h') then 
+
+      ! we are generally working with non-symmetric matrices, so
+      ! we de-symmetrize what we are about to read
+
+      read (infile,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+      read (infile,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+      if (valcrd > 0) read (infile,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+
+
+      if (present(b)) then
+        if ((psb_toupper(rhstype(1:1)) == 'F').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,b,info)
+          read (infile,fmt=rhsfmt) (b(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(g)) then
+        if ((psb_toupper(rhstype(2:2)) == 'G').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,g,info)
+          read (infile,fmt=rhsfmt) (g(i,1),i=1,nrow)
+        endif
+      endif
+      if (present(x)) then
+        if ((psb_toupper(rhstype(3:3)) == 'X').and.(rhscrd > 0)) then 
+          call psb_realloc(nrow,1,x,info)
+          read (infile,fmt=rhsfmt) (x(i,1),i=1,nrow)
+        endif
+      endif
+
+
+      call acoo%mv_from_fmt(acsc,info)
+      call acoo%reallocate(2*nnzero)
+      ! A is now in COO format
+      nzr = nnzero
+      do i=1,nnzero
+        if (acoo%ia(i) /= acoo%ja(i)) then 
+          nzr = nzr + 1
+          acoo%val(nzr) = conjg(acoo%val(i))
+          acoo%ia(nzr) = acoo%ja(i)
+          acoo%ja(nzr) = acoo%ia(i)
+        end if
+      end do
+      call acoo%set_nzeros(nzr)
+      call acoo%fix(ircode)
+      if (ircode == 0) call a%mv_from(acoo)
+      if (ircode /= 0) goto 993
+
+    else
+      write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+      iret=904
+    end if
+  else
+    write(psb_err_unit,*) 'read_matrix: matrix type not yet supported'
+    iret=904
+  end if
+
+  call a%cscnv(ircode,type='csr')
+  if (infile /= 5) close(infile)
+
+  return 
+
+  ! open failed
+901 iret=901
+  write(psb_err_unit,*) 'read_matrix: could not open file ',filename,' for input'
+  return
+902 iret=902
+  write(psb_err_unit,*) 'HB_READ: Unexpected end of file '
+  return
+993 iret=993
+  write(psb_err_unit,*) 'HB_READ: Memory allocation failure'
+  return
+end subroutine lzhb_read
+
+subroutine lzhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
+  use psb_base_mod
+  implicit none
+  type(psb_lzspmat_type), intent(in), target :: a
+  integer(psb_ipk_), intent(out)        :: iret
+  character(len=*), optional, intent(in) :: mtitle
+  integer(psb_ipk_), optional, intent(in)          :: iunit
+  character(len=*), optional, intent(in) :: filename
+  character(len=*), optional, intent(in) :: key
+  complex(psb_dpk_), optional             :: rhs(:), g(:), x(:)
+  integer(psb_ipk_) :: iout
+
+  character(len=*), parameter::  ptrfmt='(10I8)',indfmt='(10I8)'
+  integer(psb_ipk_), parameter :: jptr=10,jind=10
+  character(len=*), parameter::  valfmt='(4E20.12)',rhsfmt='(4E20.12)'
+  integer(psb_ipk_), parameter :: jval=2,jrhs=2
+  character(len=*), parameter :: fmt10='(a72,a8,/,5i14,/,a3,11x,4i14,/,2a16,2a20)'
+  character(len=*), parameter :: fmt11='(a3,11x,2i14)'
+  character(len=*), parameter :: fmt111='(1x,a8,1x,i8,1x,a10)'
+
+  type(psb_lz_csc_sparse_mat)  ::  acsc
+  character(len=72) :: mtitle_
+  character(len=8)  :: key_
+
+  character  :: rhstype*3,type*3
+
+  integer(psb_ipk_) :: i,indcrd,ptrcrd,rhscrd,totcrd,valcrd,&
+       & nrow,ncol,nnzero, neltvl, nrhs, nrhsix
+
+  iret = 0
+
+  if (present(filename)) then 
+    if (filename == '-') then 
+      iout=6
+    else
+      if (present(iunit)) then 
+        iout = iunit
+      else
+        iout=99
+      endif
+      open(iout,file=filename, err=901, action='WRITE')
+    endif
+  else 
+    if (present(iunit)) then 
+      iout = iunit   
+    else
+      iout=6
+    endif
+  endif
+
+  if (present(mtitle)) then 
+    mtitle_ = mtitle
+  else
+    mtitle_ = 'Temporary PSBLAS title '
+  endif
+  if (present(key)) then 
+    key_ = key
+  else
+    key_ = 'PSBMAT00'
+  endif
+
+
+  call acsc%cp_from_fmt(a%a, iret)
+  if (iret /= 0) return
+
+
+  nrow   = acsc%get_nrows()
+  ncol   = acsc%get_ncols()
+  nnzero = acsc%get_nzeros()
+
+  neltvl = 0 
+
+  ptrcrd = (ncol+1)/jptr
+  if (mod(ncol+1,jptr) > 0) ptrcrd = ptrcrd + 1
+  indcrd = nnzero/jind
+  if (mod(nnzero,jind) > 0) indcrd = indcrd + 1
+  valcrd = nnzero/jval
+  if (mod(nnzero,jval) > 0) valcrd = valcrd + 1
+  rhstype = ''
+  if (present(rhs)) then 
+    if (size(rhs)<nrow) then 
+      rhscrd = 0
+    else
+      rhscrd = nrow/jrhs
+      if (mod(nrow,jrhs) > 0) rhscrd = rhscrd + 1
+    endif
+    nrhs = 1
+    rhstype(1:1) = 'F'
+  else
+    rhscrd = 0
+    nrhs   = 0
+  end if
+  totcrd = ptrcrd + indcrd + valcrd + rhscrd
+
+  nrhsix  = nrhs*nrow
+
+  if (present(g)) then 
+    rhstype(2:2) = 'G'
+  end if
+  if (present(x)) then 
+    rhstype(3:3) = 'X'
+  end if
+  type    = 'CUA'
+
+  write (iout,fmt=fmt10) mtitle_,key_,totcrd,ptrcrd,indcrd,valcrd,rhscrd,&
+       &  type,nrow,ncol,nnzero,neltvl,ptrfmt,indfmt,valfmt,rhsfmt
+  if (rhscrd > 0) write (iout,fmt=fmt11) rhstype,nrhs,nrhsix
+  write (iout,fmt=ptrfmt) (acsc%icp(i),i=1,ncol+1)
+  write (iout,fmt=indfmt) (acsc%ia(i),i=1,nnzero)
+  if (valcrd > 0) write (iout,fmt=valfmt) (acsc%val(i),i=1,nnzero)
+  if (rhscrd > 0) write (iout,fmt=rhsfmt) (rhs(i),i=1,nrow)
+  if (present(g).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (g(i),i=1,nrow)
+  if (present(x).and.(rhscrd>0)) write (iout,fmt=rhsfmt) (x(i),i=1,nrow)
 
 
   if (iout /= 6) close(iout)
@@ -360,4 +690,4 @@ subroutine zhb_write(a,iret,iunit,filename,key,rhs,g,x,mtitle)
   iret=901
   write(psb_err_unit,*) 'Error while opening ',filename
   return
-end subroutine zhb_write
+end subroutine lzhb_write