Merge branch 'psblas-paraggr' into psblas-3.6-maint

psblas-3.6-maint
Salvatore Filippone 6 years ago
commit 62020199b3

@ -30,7 +30,7 @@
! !
! !
! File: psb_cspgather.f90 ! File: psb_cspgather.f90
subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc,desc_c)
use psb_desc_mod use psb_desc_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
@ -43,18 +43,21 @@ subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
#ifdef MPI_H #ifdef MPI_H
include 'mpif.h' include 'mpif.h'
#endif #endif
type(psb_cspmat_type), intent(inout) :: loca type(psb_cspmat_type), intent(inout) :: loca
type(psb_cspmat_type), intent(inout) :: globa type(psb_cspmat_type), intent(inout) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root, dupl integer(psb_ipk_), intent(in), optional :: root, dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
type(psb_c_coo_sparse_mat) :: loc_coo, glob_coo !
type(psb_c_coo_sparse_mat) :: loc_coo, glob_coo
type(psb_desc_type), pointer :: p_desc_c
integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
logical :: keepnum_, keeploc_ logical :: keepnum_, keeploc_
integer(psb_mpik_) :: ictxt,np,me integer(psb_mpik_) :: ictxt,np,me, root_
integer(psb_mpik_) :: icomm, minfo, ndx integer(psb_mpik_) :: icomm, minfo, ndx
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
@ -80,11 +83,22 @@ subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
else else
keeploc_ = .true. keeploc_ = .true.
end if end if
if (present(desc_c)) then
p_desc_c => desc_c
else
p_desc_c => desc_a
end if
if (present(root)) then
root_ = root
else
root_ = -1
end if
call globa%free() call globa%free()
if (keepnum_) then if (keepnum_) then
nrg = desc_a%get_global_rows() nrg = desc_a%get_global_rows()
ncg = desc_a%get_global_rows() ncg = p_desc_c%get_global_cols()
allocate(nzbr(np), idisp(np),stat=info) allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -102,7 +116,7 @@ subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
end if end if
nzl = loc_coo%get_nzeros() nzl = loc_coo%get_nzeros()
call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ja(1:nzl),p_desc_c,info,iact='I')
nzbr(:) = 0 nzbr(:) = 0
nzbr(me+1) = nzl nzbr(me+1) = nzl
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))

@ -30,7 +30,7 @@
! !
! !
! File: psb_dspgather.f90 ! File: psb_dspgather.f90
subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc,desc_c)
use psb_desc_mod use psb_desc_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
@ -43,18 +43,21 @@ subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
#ifdef MPI_H #ifdef MPI_H
include 'mpif.h' include 'mpif.h'
#endif #endif
type(psb_dspmat_type), intent(inout) :: loca type(psb_dspmat_type), intent(inout) :: loca
type(psb_dspmat_type), intent(inout) :: globa type(psb_dspmat_type), intent(inout) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root, dupl integer(psb_ipk_), intent(in), optional :: root, dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
type(psb_d_coo_sparse_mat) :: loc_coo, glob_coo !
type(psb_d_coo_sparse_mat) :: loc_coo, glob_coo
type(psb_desc_type), pointer :: p_desc_c
integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
logical :: keepnum_, keeploc_ logical :: keepnum_, keeploc_
integer(psb_mpik_) :: ictxt,np,me integer(psb_mpik_) :: ictxt,np,me, root_
integer(psb_mpik_) :: icomm, minfo, ndx integer(psb_mpik_) :: icomm, minfo, ndx
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
@ -80,11 +83,22 @@ subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
else else
keeploc_ = .true. keeploc_ = .true.
end if end if
if (present(desc_c)) then
p_desc_c => desc_c
else
p_desc_c => desc_a
end if
if (present(root)) then
root_ = root
else
root_ = -1
end if
call globa%free() call globa%free()
if (keepnum_) then if (keepnum_) then
nrg = desc_a%get_global_rows() nrg = desc_a%get_global_rows()
ncg = desc_a%get_global_rows() ncg = p_desc_c%get_global_cols()
allocate(nzbr(np), idisp(np),stat=info) allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -102,7 +116,7 @@ subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
end if end if
nzl = loc_coo%get_nzeros() nzl = loc_coo%get_nzeros()
call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ja(1:nzl),p_desc_c,info,iact='I')
nzbr(:) = 0 nzbr(:) = 0
nzbr(me+1) = nzl nzbr(me+1) = nzl
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))

@ -30,7 +30,7 @@
! !
! !
! File: psb_sspgather.f90 ! File: psb_sspgather.f90
subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc,desc_c)
use psb_desc_mod use psb_desc_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
@ -43,18 +43,21 @@ subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
#ifdef MPI_H #ifdef MPI_H
include 'mpif.h' include 'mpif.h'
#endif #endif
type(psb_sspmat_type), intent(inout) :: loca type(psb_sspmat_type), intent(inout) :: loca
type(psb_sspmat_type), intent(inout) :: globa type(psb_sspmat_type), intent(inout) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root, dupl integer(psb_ipk_), intent(in), optional :: root, dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
type(psb_s_coo_sparse_mat) :: loc_coo, glob_coo !
type(psb_s_coo_sparse_mat) :: loc_coo, glob_coo
type(psb_desc_type), pointer :: p_desc_c
integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
logical :: keepnum_, keeploc_ logical :: keepnum_, keeploc_
integer(psb_mpik_) :: ictxt,np,me integer(psb_mpik_) :: ictxt,np,me, root_
integer(psb_mpik_) :: icomm, minfo, ndx integer(psb_mpik_) :: icomm, minfo, ndx
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
@ -80,11 +83,22 @@ subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
else else
keeploc_ = .true. keeploc_ = .true.
end if end if
if (present(desc_c)) then
p_desc_c => desc_c
else
p_desc_c => desc_a
end if
if (present(root)) then
root_ = root
else
root_ = -1
end if
call globa%free() call globa%free()
if (keepnum_) then if (keepnum_) then
nrg = desc_a%get_global_rows() nrg = desc_a%get_global_rows()
ncg = desc_a%get_global_rows() ncg = p_desc_c%get_global_cols()
allocate(nzbr(np), idisp(np),stat=info) allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -102,7 +116,7 @@ subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
end if end if
nzl = loc_coo%get_nzeros() nzl = loc_coo%get_nzeros()
call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ja(1:nzl),p_desc_c,info,iact='I')
nzbr(:) = 0 nzbr(:) = 0
nzbr(me+1) = nzl nzbr(me+1) = nzl
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))

@ -30,7 +30,7 @@
! !
! !
! File: psb_zspgather.f90 ! File: psb_zspgather.f90
subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc,desc_c)
use psb_desc_mod use psb_desc_mod
use psb_error_mod use psb_error_mod
use psb_penv_mod use psb_penv_mod
@ -43,18 +43,21 @@ subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
#ifdef MPI_H #ifdef MPI_H
include 'mpif.h' include 'mpif.h'
#endif #endif
type(psb_zspmat_type), intent(inout) :: loca type(psb_zspmat_type), intent(inout) :: loca
type(psb_zspmat_type), intent(inout) :: globa type(psb_zspmat_type), intent(inout) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root, dupl integer(psb_ipk_), intent(in), optional :: root, dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
type(psb_z_coo_sparse_mat) :: loc_coo, glob_coo !
type(psb_z_coo_sparse_mat) :: loc_coo, glob_coo
type(psb_desc_type), pointer :: p_desc_c
integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg integer(psb_ipk_) :: err_act, dupl_, nrg, ncg, nzg
integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl integer(psb_ipk_) :: ip,naggrm1,naggrp1, i, j, k, nzl
logical :: keepnum_, keeploc_ logical :: keepnum_, keeploc_
integer(psb_mpik_) :: ictxt,np,me integer(psb_mpik_) :: ictxt,np,me, root_
integer(psb_mpik_) :: icomm, minfo, ndx integer(psb_mpik_) :: icomm, minfo, ndx
integer(psb_mpik_), allocatable :: nzbr(:), idisp(:) integer(psb_mpik_), allocatable :: nzbr(:), idisp(:)
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
@ -80,11 +83,22 @@ subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
else else
keeploc_ = .true. keeploc_ = .true.
end if end if
if (present(desc_c)) then
p_desc_c => desc_c
else
p_desc_c => desc_a
end if
if (present(root)) then
root_ = root
else
root_ = -1
end if
call globa%free() call globa%free()
if (keepnum_) then if (keepnum_) then
nrg = desc_a%get_global_rows() nrg = desc_a%get_global_rows()
ncg = desc_a%get_global_rows() ncg = p_desc_c%get_global_cols()
allocate(nzbr(np), idisp(np),stat=info) allocate(nzbr(np), idisp(np),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -102,7 +116,7 @@ subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keep
end if end if
nzl = loc_coo%get_nzeros() nzl = loc_coo%get_nzeros()
call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ia(1:nzl),desc_a,info,iact='I')
call psb_loc_to_glob(loc_coo%ja(1:nzl),desc_a,info,iact='I') call psb_loc_to_glob(loc_coo%ja(1:nzl),p_desc_c,info,iact='I')
nzbr(:) = 0 nzbr(:) = 0
nzbr(me+1) = nzl nzbr(me+1) = nzl
call psb_sum(ictxt,nzbr(1:np)) call psb_sum(ictxt,nzbr(1:np))

@ -47,13 +47,15 @@ module psb_base_linmap_mod
type(psb_desc_type), pointer :: p_desc_U=>null(), p_desc_V=>null() type(psb_desc_type), pointer :: p_desc_U=>null(), p_desc_V=>null()
type(psb_desc_type) :: desc_U, desc_V type(psb_desc_type) :: desc_U, desc_V
contains contains
procedure, pass(map) :: sizeof => base_map_sizeof procedure, pass(map) :: sizeof => base_map_sizeof
procedure, pass(map) :: is_ok => base_is_ok procedure, pass(map) :: is_ok => base_is_ok
procedure, pass(map) :: is_asb => base_is_asb procedure, pass(map) :: is_asb => base_is_asb
procedure, pass(map) :: get_kind => base_get_kind procedure, pass(map) :: get_kind => base_get_kind
procedure, pass(map) :: set_kind => base_set_kind procedure, pass(map) :: set_kind => base_set_kind
procedure, pass(map) :: free => base_free procedure, pass(map) :: is_dec_aggr => base_is_dec_aggr
procedure, pass(map) :: clone => base_clone procedure, pass(map) :: is_gen_linear => base_is_gen_linear
procedure, pass(map) :: free => base_free
procedure, pass(map) :: clone => base_clone
end type psb_base_linmap_type end type psb_base_linmap_type
@ -62,7 +64,8 @@ module psb_base_linmap_mod
end interface end interface
private :: base_map_sizeof, base_is_ok, base_is_asb,& private :: base_map_sizeof, base_is_ok, base_is_asb,&
& base_get_kind, base_set_kind, base_free, base_clone & base_get_kind, base_set_kind, base_free, base_clone,&
& base_is_dec_aggr, base_is_gen_linear
contains contains
@ -93,7 +96,7 @@ contains
res = .false. res = .false.
select case(map%get_kind()) select case(map%get_kind())
case (psb_map_aggr_) case (psb_map_dec_aggr_)
if (.not.associated(map%p_desc_U)) return if (.not.associated(map%p_desc_U)) return
if (.not.associated(map%p_desc_V)) return if (.not.associated(map%p_desc_V)) return
res = map%p_desc_U%is_ok().and.map%p_desc_V%is_ok() res = map%p_desc_U%is_ok().and.map%p_desc_V%is_ok()
@ -111,7 +114,7 @@ contains
res = .false. res = .false.
select case(map%get_kind()) select case(map%get_kind())
case (psb_map_aggr_) case (psb_map_dec_aggr_)
if (.not.associated(map%p_desc_U)) return if (.not.associated(map%p_desc_U)) return
if (.not.associated(map%p_desc_V)) return if (.not.associated(map%p_desc_V)) return
res = map%p_desc_U%is_asb().and.map%p_desc_V%is_asb() res = map%p_desc_U%is_asb().and.map%p_desc_V%is_asb()
@ -121,6 +124,24 @@ contains
end function base_is_asb end function base_is_asb
function base_is_dec_aggr(map) result(res)
use psb_desc_mod
implicit none
class(psb_base_linmap_type), intent(in) :: map
logical :: res
res = (map%get_kind() == psb_map_dec_aggr_)
end function base_is_dec_aggr
function base_is_gen_linear(map) result(res)
use psb_desc_mod
implicit none
class(psb_base_linmap_type), intent(in) :: map
logical :: res
res = (map%get_kind() == psb_map_gen_linear_)
end function base_is_gen_linear
function base_map_sizeof(map) result(val) function base_map_sizeof(map) result(val)
use psb_desc_mod use psb_desc_mod
implicit none implicit none

@ -151,15 +151,16 @@ module psb_c_comm_mod
end interface psb_scatter end interface psb_scatter
interface psb_gather interface psb_gather
subroutine psb_csp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_csp_allgather(globa, loca, desc_a, info, root,dupl,keepnum,keeploc,desc_c)
import import
implicit none implicit none
type(psb_cspmat_type), intent(inout) :: loca type(psb_cspmat_type), intent(inout) :: loca
type(psb_cspmat_type), intent(out) :: globa type(psb_cspmat_type), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root,dupl integer(psb_ipk_), intent(in), optional :: root,dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
end subroutine psb_csp_allgather end subroutine psb_csp_allgather
subroutine psb_cgatherm(globx, locx, desc_a, info, root) subroutine psb_cgatherm(globx, locx, desc_a, info, root)
import import

@ -151,15 +151,16 @@ module psb_d_comm_mod
end interface psb_scatter end interface psb_scatter
interface psb_gather interface psb_gather
subroutine psb_dsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_dsp_allgather(globa, loca, desc_a, info, root,dupl,keepnum,keeploc,desc_c)
import import
implicit none implicit none
type(psb_dspmat_type), intent(inout) :: loca type(psb_dspmat_type), intent(inout) :: loca
type(psb_dspmat_type), intent(out) :: globa type(psb_dspmat_type), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root,dupl integer(psb_ipk_), intent(in), optional :: root,dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
end subroutine psb_dsp_allgather end subroutine psb_dsp_allgather
subroutine psb_dgatherm(globx, locx, desc_a, info, root) subroutine psb_dgatherm(globx, locx, desc_a, info, root)
import import

@ -151,15 +151,16 @@ module psb_s_comm_mod
end interface psb_scatter end interface psb_scatter
interface psb_gather interface psb_gather
subroutine psb_ssp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_ssp_allgather(globa, loca, desc_a, info, root,dupl,keepnum,keeploc,desc_c)
import import
implicit none implicit none
type(psb_sspmat_type), intent(inout) :: loca type(psb_sspmat_type), intent(inout) :: loca
type(psb_sspmat_type), intent(out) :: globa type(psb_sspmat_type), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root,dupl integer(psb_ipk_), intent(in), optional :: root,dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
end subroutine psb_ssp_allgather end subroutine psb_ssp_allgather
subroutine psb_sgatherm(globx, locx, desc_a, info, root) subroutine psb_sgatherm(globx, locx, desc_a, info, root)
import import

@ -151,15 +151,16 @@ module psb_z_comm_mod
end interface psb_scatter end interface psb_scatter
interface psb_gather interface psb_gather
subroutine psb_zsp_allgather(globa, loca, desc_a, info, root, dupl,keepnum,keeploc) subroutine psb_zsp_allgather(globa, loca, desc_a, info, root,dupl,keepnum,keeploc,desc_c)
import import
implicit none implicit none
type(psb_zspmat_type), intent(inout) :: loca type(psb_zspmat_type), intent(inout) :: loca
type(psb_zspmat_type), intent(out) :: globa type(psb_zspmat_type), intent(out) :: globa
type(psb_desc_type), intent(in) :: desc_a type(psb_desc_type), intent(in), target :: desc_a
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: root,dupl integer(psb_ipk_), intent(in), optional :: root,dupl
logical, intent(in), optional :: keepnum,keeploc logical, intent(in), optional :: keepnum,keeploc
type(psb_desc_type), intent(in), optional, target :: desc_c
end subroutine psb_zsp_allgather end subroutine psb_zsp_allgather
subroutine psb_zgatherm(globx, locx, desc_a, info, root) subroutine psb_zgatherm(globx, locx, desc_a, info, root)
import import

@ -55,8 +55,8 @@ module psb_desc_const_mod
! Types of mapping between descriptors. ! Types of mapping between descriptors.
integer(psb_ipk_), parameter :: psb_map_xhal_ = 123 integer(psb_ipk_), parameter :: psb_map_xhal_ = 123
integer(psb_ipk_), parameter :: psb_map_asov_ = psb_map_xhal_+1 integer(psb_ipk_), parameter :: psb_map_asov_ = psb_map_xhal_+1
integer(psb_ipk_), parameter :: psb_map_aggr_ = psb_map_asov_+1 integer(psb_ipk_), parameter :: psb_map_dec_aggr_ = psb_map_asov_+1
integer(psb_ipk_), parameter :: psb_map_gen_linear_ = psb_map_aggr_+1 integer(psb_ipk_), parameter :: psb_map_gen_linear_ = psb_map_dec_aggr_+1
integer(psb_ipk_), parameter :: psb_ovt_xhal_ = psb_map_xhal_, psb_ovt_asov_=psb_map_asov_ integer(psb_ipk_), parameter :: psb_ovt_xhal_ = psb_map_xhal_, psb_ovt_asov_=psb_map_asov_
! !

@ -154,10 +154,10 @@ module psb_const_mod
! Duplicate coefficients handling ! Duplicate coefficients handling
! These are usually set while calling spcnv as one of its ! These are usually set while calling spcnv as one of its
! optional arugments. ! optional arugments.
integer(psb_ipk_), parameter :: psb_dupl_ovwrt_ = 0 integer(psb_ipk_), parameter :: psb_dupl_add_ = 0
integer(psb_ipk_), parameter :: psb_dupl_add_ = 1 integer(psb_ipk_), parameter :: psb_dupl_ovwrt_ = 1
integer(psb_ipk_), parameter :: psb_dupl_err_ = 2 integer(psb_ipk_), parameter :: psb_dupl_err_ = 2
integer(psb_ipk_), parameter :: psb_dupl_def_ = psb_dupl_ovwrt_ integer(psb_ipk_), parameter :: psb_dupl_def_ = psb_dupl_add_
! Matrix update mode ! Matrix update mode
integer(psb_ipk_), parameter :: psb_upd_srch_ = 98764 integer(psb_ipk_), parameter :: psb_upd_srch_ = 98764
integer(psb_ipk_), parameter :: psb_upd_perm_ = 98765 integer(psb_ipk_), parameter :: psb_upd_perm_ = 98765

@ -63,7 +63,8 @@ module psb_timers_mod
type(psb_string_item), allocatable :: timers_descr(:) type(psb_string_item), allocatable :: timers_descr(:)
logical :: wanted(timer_entries_) logical :: wanted(timer_entries_)
type(psb_string_item) :: entries_descr(timer_entries_) type(psb_string_item) :: entries_descr(timer_entries_)
save :: nsamples, timers, timers_descr, wanted, entries_descr
interface psb_realloc interface psb_realloc
module procedure psb_string_item_realloc module procedure psb_string_item_realloc
end interface psb_realloc end interface psb_realloc
@ -127,17 +128,19 @@ contains
end if end if
if (global_) then if (global_) then
allocate(ptimers(timer_entries_,size(timers,2)),stat=info) if (allocated(timers)) then
if (info /= 0) then allocate(ptimers(timer_entries_,size(timers,2)),stat=info)
write(0,*) 'Error while trying to allocate temporary ',info if (info /= 0) then
call psb_abort(ictxt) write(0,*) 'Error while trying to allocate temporary ',info
end if call psb_abort(ictxt)
ptimers = timers end if
call psb_max(ictxt,ptimers) ptimers = timers
if (me == psb_root_) then call psb_max(ictxt,ptimers)
do i=idxmin_, idxmax_ if (me == psb_root_) then
call print_timer(me, ptimers(:,i), timers_descr(i), iout) do i=idxmin_, idxmax_
end do call print_timer(me, ptimers(:,i), timers_descr(i), iout)
end do
end if
end if end if
else else
if ((proc_ == -1).or.(me==proc_)) then if ((proc_ == -1).or.(me==proc_)) then

@ -596,8 +596,121 @@ module psb_c_csr_mat_mod
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_csr_scals end subroutine psb_c_csr_scals
end interface end interface
type, extends(psb_c_csr_sparse_mat) :: psb_c_csrd_sparse_mat
!> Pointers to diagonal entries
integer(psb_ipk_), allocatable :: diagp(:)
contains
procedure, nopass :: get_fmt => c_csrd_get_fmt
procedure, pass(a) :: sizeof => c_csrd_sizeof
procedure, pass(a) :: inner_cssm => psb_c_csrd_cssm
procedure, pass(a) :: inner_cssv => psb_c_csrd_cssv
procedure, pass(a) :: trmv => psb_c_csrd_trmv
!procedure, pass(a) :: reallocate_nz => psb_c_csrd_reallocate_nz
!procedure, pass(a) :: allocate_mnnz => psb_c_csrd_allocate_mnnz
!!$ procedure, pass(a) :: cp_to_coo => psb_c_cp_csrd_to_coo
procedure, pass(a) :: cp_from_coo => psb_c_cp_csrd_from_coo
!!$ procedure, pass(a) :: cp_to_fmt => psb_c_cp_csrd_to_fmt
!!$ procedure, pass(a) :: cp_from_fmt => psb_c_cp_csrd_from_fmt
!!$ procedure, pass(a) :: mv_to_coo => psb_c_mv_csrd_to_coo
!!$ procedure, pass(a) :: mv_from_coo => psb_c_mv_csrd_from_coo
!!$ procedure, pass(a) :: mv_to_fmt => psb_c_mv_csrd_to_fmt
!!$ procedure, pass(a) :: mv_from_fmt => psb_c_mv_csrd_from_fmt
procedure, pass(a) :: clean_zeros => psb_c_csrd_clean_zeros
procedure, pass(a) :: get_diag => psb_c_csrd_get_diag
!procedure, pass(a) :: reinit => psb_c_csrd_reinit
procedure, pass(a) :: trim => psb_c_csrd_trim
procedure, pass(a) :: free => c_csrd_free
procedure, pass(a) :: mold => psb_c_csrd_mold
procedure, nopass :: has_xt_tri => c_csrd_has_xt_tri
end type psb_c_csrd_sparse_mat
interface
subroutine psb_c_csrd_clean_zeros(a, info)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_csrd_clean_zeros
end interface
interface
subroutine psb_c_csrd_get_diag(a,d,info)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
complex(psb_spk_), intent(out) :: d(:)
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_csrd_get_diag
end interface
interface
subroutine psb_c_cp_csrd_from_coo(a,b,info)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(inout) :: a
class(psb_c_coo_sparse_mat), intent(in) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_cp_csrd_from_coo
end interface
interface
subroutine psb_c_csrd_trim(a)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(inout) :: a
end subroutine psb_c_csrd_trim
end interface
interface
subroutine psb_c_csrd_mold(a,b,info)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_c_base_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
class(psb_c_base_sparse_mat), intent(inout), allocatable :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_c_csrd_mold
end interface
interface
subroutine psb_c_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
complex(psb_spk_), intent(in) :: alpha, beta, x(:)
complex(psb_spk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: uplo, diag
end subroutine psb_c_csrd_trmv
end interface
interface
subroutine psb_c_csrd_cssm(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
complex(psb_spk_), intent(in) :: alpha, beta, x(:,:)
complex(psb_spk_), intent(inout) :: y(:,:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_c_csrd_cssm
end interface
interface
subroutine psb_c_csrd_cssv(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
complex(psb_spk_), intent(in) :: alpha, beta, x(:)
complex(psb_spk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_c_csrd_cssv
end interface
contains contains
@ -717,4 +830,59 @@ contains
end subroutine c_csr_free end subroutine c_csr_free
! == ===================================
!
!
!
! CSRD specific versions
!
!
!
!
!
! == ===================================
function c_csrd_sizeof(a) result(res)
implicit none
class(psb_c_csrd_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
res = res + (2*psb_sizeof_sp) * psb_size(a%val)
res = res + psb_sizeof_int * psb_size(a%irp)
res = res + psb_sizeof_int * psb_size(a%ja)
res = res + psb_sizeof_int * psb_size(a%diagp)
end function c_csrd_sizeof
function c_csrd_get_fmt() result(res)
implicit none
character(len=5) :: res
res = 'CSRD'
end function c_csrd_get_fmt
subroutine c_csrd_free(a)
implicit none
class(psb_c_csrd_sparse_mat), intent(inout) :: a
if (allocated(a%diagp)) deallocate(a%diagp)
call a%psb_c_csr_sparse_mat%free()
return
end subroutine c_csrd_free
!
! has_xt_tri: does the current type support
! extended triangle operations?
!
function c_csrd_has_xt_tri() result(res)
implicit none
logical :: res
res = .true.
end function c_csrd_has_xt_tri
end module psb_c_csr_mat_mod end module psb_c_csr_mat_mod

@ -596,8 +596,121 @@ module psb_d_csr_mat_mod
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_csr_scals end subroutine psb_d_csr_scals
end interface end interface
type, extends(psb_d_csr_sparse_mat) :: psb_d_csrd_sparse_mat
!> Pointers to diagonal entries
integer(psb_ipk_), allocatable :: diagp(:)
contains
procedure, nopass :: get_fmt => d_csrd_get_fmt
procedure, pass(a) :: sizeof => d_csrd_sizeof
procedure, pass(a) :: inner_cssm => psb_d_csrd_cssm
procedure, pass(a) :: inner_cssv => psb_d_csrd_cssv
procedure, pass(a) :: trmv => psb_d_csrd_trmv
!procedure, pass(a) :: reallocate_nz => psb_d_csrd_reallocate_nz
!procedure, pass(a) :: allocate_mnnz => psb_d_csrd_allocate_mnnz
!!$ procedure, pass(a) :: cp_to_coo => psb_d_cp_csrd_to_coo
procedure, pass(a) :: cp_from_coo => psb_d_cp_csrd_from_coo
!!$ procedure, pass(a) :: cp_to_fmt => psb_d_cp_csrd_to_fmt
!!$ procedure, pass(a) :: cp_from_fmt => psb_d_cp_csrd_from_fmt
!!$ procedure, pass(a) :: mv_to_coo => psb_d_mv_csrd_to_coo
!!$ procedure, pass(a) :: mv_from_coo => psb_d_mv_csrd_from_coo
!!$ procedure, pass(a) :: mv_to_fmt => psb_d_mv_csrd_to_fmt
!!$ procedure, pass(a) :: mv_from_fmt => psb_d_mv_csrd_from_fmt
procedure, pass(a) :: clean_zeros => psb_d_csrd_clean_zeros
procedure, pass(a) :: get_diag => psb_d_csrd_get_diag
!procedure, pass(a) :: reinit => psb_d_csrd_reinit
procedure, pass(a) :: trim => psb_d_csrd_trim
procedure, pass(a) :: free => d_csrd_free
procedure, pass(a) :: mold => psb_d_csrd_mold
procedure, nopass :: has_xt_tri => d_csrd_has_xt_tri
end type psb_d_csrd_sparse_mat
interface
subroutine psb_d_csrd_clean_zeros(a, info)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_csrd_clean_zeros
end interface
interface
subroutine psb_d_csrd_get_diag(a,d,info)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(out) :: d(:)
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_csrd_get_diag
end interface
interface
subroutine psb_d_cp_csrd_from_coo(a,b,info)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(inout) :: a
class(psb_d_coo_sparse_mat), intent(in) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_cp_csrd_from_coo
end interface
interface
subroutine psb_d_csrd_trim(a)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(inout) :: a
end subroutine psb_d_csrd_trim
end interface
interface
subroutine psb_d_csrd_mold(a,b,info)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_d_base_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
class(psb_d_base_sparse_mat), intent(inout), allocatable :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_d_csrd_mold
end interface
interface
subroutine psb_d_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: uplo, diag
end subroutine psb_d_csrd_trmv
end interface
interface
subroutine psb_d_csrd_cssm(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
real(psb_dpk_), intent(inout) :: y(:,:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_d_csrd_cssm
end interface
interface
subroutine psb_d_csrd_cssv(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
real(psb_dpk_), intent(in) :: alpha, beta, x(:)
real(psb_dpk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_d_csrd_cssv
end interface
contains contains
@ -717,4 +830,59 @@ contains
end subroutine d_csr_free end subroutine d_csr_free
! == ===================================
!
!
!
! CSRD specific versions
!
!
!
!
!
! == ===================================
function d_csrd_sizeof(a) result(res)
implicit none
class(psb_d_csrd_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
res = res + psb_sizeof_dp * psb_size(a%val)
res = res + psb_sizeof_int * psb_size(a%irp)
res = res + psb_sizeof_int * psb_size(a%ja)
res = res + psb_sizeof_int * psb_size(a%diagp)
end function d_csrd_sizeof
function d_csrd_get_fmt() result(res)
implicit none
character(len=5) :: res
res = 'CSRD'
end function d_csrd_get_fmt
subroutine d_csrd_free(a)
implicit none
class(psb_d_csrd_sparse_mat), intent(inout) :: a
if (allocated(a%diagp)) deallocate(a%diagp)
call a%psb_d_csr_sparse_mat%free()
return
end subroutine d_csrd_free
!
! has_xt_tri: does the current type support
! extended triangle operations?
!
function d_csrd_has_xt_tri() result(res)
implicit none
logical :: res
res = .true.
end function d_csrd_has_xt_tri
end module psb_d_csr_mat_mod end module psb_d_csr_mat_mod

@ -596,8 +596,121 @@ module psb_s_csr_mat_mod
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_csr_scals end subroutine psb_s_csr_scals
end interface end interface
type, extends(psb_s_csr_sparse_mat) :: psb_s_csrd_sparse_mat
!> Pointers to diagonal entries
integer(psb_ipk_), allocatable :: diagp(:)
contains
procedure, nopass :: get_fmt => s_csrd_get_fmt
procedure, pass(a) :: sizeof => s_csrd_sizeof
procedure, pass(a) :: inner_cssm => psb_s_csrd_cssm
procedure, pass(a) :: inner_cssv => psb_s_csrd_cssv
procedure, pass(a) :: trmv => psb_s_csrd_trmv
!procedure, pass(a) :: reallocate_nz => psb_s_csrd_reallocate_nz
!procedure, pass(a) :: allocate_mnnz => psb_s_csrd_allocate_mnnz
!!$ procedure, pass(a) :: cp_to_coo => psb_s_cp_csrd_to_coo
procedure, pass(a) :: cp_from_coo => psb_s_cp_csrd_from_coo
!!$ procedure, pass(a) :: cp_to_fmt => psb_s_cp_csrd_to_fmt
!!$ procedure, pass(a) :: cp_from_fmt => psb_s_cp_csrd_from_fmt
!!$ procedure, pass(a) :: mv_to_coo => psb_s_mv_csrd_to_coo
!!$ procedure, pass(a) :: mv_from_coo => psb_s_mv_csrd_from_coo
!!$ procedure, pass(a) :: mv_to_fmt => psb_s_mv_csrd_to_fmt
!!$ procedure, pass(a) :: mv_from_fmt => psb_s_mv_csrd_from_fmt
procedure, pass(a) :: clean_zeros => psb_s_csrd_clean_zeros
procedure, pass(a) :: get_diag => psb_s_csrd_get_diag
!procedure, pass(a) :: reinit => psb_s_csrd_reinit
procedure, pass(a) :: trim => psb_s_csrd_trim
procedure, pass(a) :: free => s_csrd_free
procedure, pass(a) :: mold => psb_s_csrd_mold
procedure, nopass :: has_xt_tri => s_csrd_has_xt_tri
end type psb_s_csrd_sparse_mat
interface
subroutine psb_s_csrd_clean_zeros(a, info)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_csrd_clean_zeros
end interface
interface
subroutine psb_s_csrd_get_diag(a,d,info)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
real(psb_spk_), intent(out) :: d(:)
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_csrd_get_diag
end interface
interface
subroutine psb_s_cp_csrd_from_coo(a,b,info)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(inout) :: a
class(psb_s_coo_sparse_mat), intent(in) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_cp_csrd_from_coo
end interface
interface
subroutine psb_s_csrd_trim(a)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(inout) :: a
end subroutine psb_s_csrd_trim
end interface
interface
subroutine psb_s_csrd_mold(a,b,info)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_s_base_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
class(psb_s_base_sparse_mat), intent(inout), allocatable :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_s_csrd_mold
end interface
interface
subroutine psb_s_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
real(psb_spk_), intent(in) :: alpha, beta, x(:)
real(psb_spk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: uplo, diag
end subroutine psb_s_csrd_trmv
end interface
interface
subroutine psb_s_csrd_cssm(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
real(psb_spk_), intent(in) :: alpha, beta, x(:,:)
real(psb_spk_), intent(inout) :: y(:,:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_s_csrd_cssm
end interface
interface
subroutine psb_s_csrd_cssv(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
real(psb_spk_), intent(in) :: alpha, beta, x(:)
real(psb_spk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_s_csrd_cssv
end interface
contains contains
@ -717,4 +830,59 @@ contains
end subroutine s_csr_free end subroutine s_csr_free
! == ===================================
!
!
!
! CSRD specific versions
!
!
!
!
!
! == ===================================
function s_csrd_sizeof(a) result(res)
implicit none
class(psb_s_csrd_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
res = res + psb_sizeof_sp * psb_size(a%val)
res = res + psb_sizeof_int * psb_size(a%irp)
res = res + psb_sizeof_int * psb_size(a%ja)
res = res + psb_sizeof_int * psb_size(a%diagp)
end function s_csrd_sizeof
function s_csrd_get_fmt() result(res)
implicit none
character(len=5) :: res
res = 'CSRD'
end function s_csrd_get_fmt
subroutine s_csrd_free(a)
implicit none
class(psb_s_csrd_sparse_mat), intent(inout) :: a
if (allocated(a%diagp)) deallocate(a%diagp)
call a%psb_s_csr_sparse_mat%free()
return
end subroutine s_csrd_free
!
! has_xt_tri: does the current type support
! extended triangle operations?
!
function s_csrd_has_xt_tri() result(res)
implicit none
logical :: res
res = .true.
end function s_csrd_has_xt_tri
end module psb_s_csr_mat_mod end module psb_s_csr_mat_mod

@ -596,8 +596,121 @@ module psb_z_csr_mat_mod
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_csr_scals end subroutine psb_z_csr_scals
end interface end interface
type, extends(psb_z_csr_sparse_mat) :: psb_z_csrd_sparse_mat
!> Pointers to diagonal entries
integer(psb_ipk_), allocatable :: diagp(:)
contains
procedure, nopass :: get_fmt => z_csrd_get_fmt
procedure, pass(a) :: sizeof => z_csrd_sizeof
procedure, pass(a) :: inner_cssm => psb_z_csrd_cssm
procedure, pass(a) :: inner_cssv => psb_z_csrd_cssv
procedure, pass(a) :: trmv => psb_z_csrd_trmv
!procedure, pass(a) :: reallocate_nz => psb_z_csrd_reallocate_nz
!procedure, pass(a) :: allocate_mnnz => psb_z_csrd_allocate_mnnz
!!$ procedure, pass(a) :: cp_to_coo => psb_z_cp_csrd_to_coo
procedure, pass(a) :: cp_from_coo => psb_z_cp_csrd_from_coo
!!$ procedure, pass(a) :: cp_to_fmt => psb_z_cp_csrd_to_fmt
!!$ procedure, pass(a) :: cp_from_fmt => psb_z_cp_csrd_from_fmt
!!$ procedure, pass(a) :: mv_to_coo => psb_z_mv_csrd_to_coo
!!$ procedure, pass(a) :: mv_from_coo => psb_z_mv_csrd_from_coo
!!$ procedure, pass(a) :: mv_to_fmt => psb_z_mv_csrd_to_fmt
!!$ procedure, pass(a) :: mv_from_fmt => psb_z_mv_csrd_from_fmt
procedure, pass(a) :: clean_zeros => psb_z_csrd_clean_zeros
procedure, pass(a) :: get_diag => psb_z_csrd_get_diag
!procedure, pass(a) :: reinit => psb_z_csrd_reinit
procedure, pass(a) :: trim => psb_z_csrd_trim
procedure, pass(a) :: free => z_csrd_free
procedure, pass(a) :: mold => psb_z_csrd_mold
procedure, nopass :: has_xt_tri => z_csrd_has_xt_tri
end type psb_z_csrd_sparse_mat
interface
subroutine psb_z_csrd_clean_zeros(a, info)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(inout) :: a
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_csrd_clean_zeros
end interface
interface
subroutine psb_z_csrd_get_diag(a,d,info)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
complex(psb_dpk_), intent(out) :: d(:)
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_csrd_get_diag
end interface
interface
subroutine psb_z_cp_csrd_from_coo(a,b,info)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(inout) :: a
class(psb_z_coo_sparse_mat), intent(in) :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_cp_csrd_from_coo
end interface
interface
subroutine psb_z_csrd_trim(a)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(inout) :: a
end subroutine psb_z_csrd_trim
end interface
interface
subroutine psb_z_csrd_mold(a,b,info)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_z_base_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
class(psb_z_base_sparse_mat), intent(inout), allocatable :: b
integer(psb_ipk_), intent(out) :: info
end subroutine psb_z_csrd_mold
end interface
interface
subroutine psb_z_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
complex(psb_dpk_), intent(in) :: alpha, beta, x(:)
complex(psb_dpk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: uplo, diag
end subroutine psb_z_csrd_trmv
end interface
interface
subroutine psb_z_csrd_cssm(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
complex(psb_dpk_), intent(in) :: alpha, beta, x(:,:)
complex(psb_dpk_), intent(inout) :: y(:,:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_z_csrd_cssm
end interface
interface
subroutine psb_z_csrd_cssv(alpha,a,x,beta,y,info,trans)
import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
complex(psb_dpk_), intent(in) :: alpha, beta, x(:)
complex(psb_dpk_), intent(inout) :: y(:)
integer(psb_ipk_), intent(out) :: info
character, optional, intent(in) :: trans
end subroutine psb_z_csrd_cssv
end interface
contains contains
@ -717,4 +830,59 @@ contains
end subroutine z_csr_free end subroutine z_csr_free
! == ===================================
!
!
!
! CSRD specific versions
!
!
!
!
!
! == ===================================
function z_csrd_sizeof(a) result(res)
implicit none
class(psb_z_csrd_sparse_mat), intent(in) :: a
integer(psb_long_int_k_) :: res
res = 8
res = res + (2*psb_sizeof_dp) * psb_size(a%val)
res = res + psb_sizeof_int * psb_size(a%irp)
res = res + psb_sizeof_int * psb_size(a%ja)
res = res + psb_sizeof_int * psb_size(a%diagp)
end function z_csrd_sizeof
function z_csrd_get_fmt() result(res)
implicit none
character(len=5) :: res
res = 'CSRD'
end function z_csrd_get_fmt
subroutine z_csrd_free(a)
implicit none
class(psb_z_csrd_sparse_mat), intent(inout) :: a
if (allocated(a%diagp)) deallocate(a%diagp)
call a%psb_z_csr_sparse_mat%free()
return
end subroutine z_csrd_free
!
! has_xt_tri: does the current type support
! extended triangle operations?
!
function z_csrd_has_xt_tri() result(res)
implicit none
logical :: res
res = .true.
end function z_csrd_has_xt_tri
end module psb_z_csr_mat_mod end module psb_z_csr_mat_mod

@ -363,7 +363,6 @@ Module psb_c_tools_mod
end subroutine psb_cspins_2desc end subroutine psb_cspins_2desc
end interface end interface
interface psb_sprn interface psb_sprn
subroutine psb_csprn(a, desc_a,info,clear) subroutine psb_csprn(a, desc_a,info,clear)
import import
@ -374,5 +373,19 @@ Module psb_c_tools_mod
logical, intent(in), optional :: clear logical, intent(in), optional :: clear
end subroutine psb_csprn end subroutine psb_csprn
end interface end interface
interface psb_par_spspmm
subroutine psb_c_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
import :: psb_c_csr_sparse_mat, psb_desc_type, psb_ipk_
Implicit None
type(psb_c_csr_sparse_mat),intent(in) :: acsr
type(psb_c_csr_sparse_mat),intent(inout) :: bcsr
type(psb_c_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
End Subroutine psb_c_par_csr_spspmm
end interface psb_par_spspmm
end module psb_c_tools_mod end module psb_c_tools_mod

@ -363,7 +363,6 @@ Module psb_d_tools_mod
end subroutine psb_dspins_2desc end subroutine psb_dspins_2desc
end interface end interface
interface psb_sprn interface psb_sprn
subroutine psb_dsprn(a, desc_a,info,clear) subroutine psb_dsprn(a, desc_a,info,clear)
import import
@ -374,5 +373,19 @@ Module psb_d_tools_mod
logical, intent(in), optional :: clear logical, intent(in), optional :: clear
end subroutine psb_dsprn end subroutine psb_dsprn
end interface end interface
interface psb_par_spspmm
subroutine psb_d_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
import :: psb_d_csr_sparse_mat, psb_desc_type, psb_ipk_
Implicit None
type(psb_d_csr_sparse_mat),intent(in) :: acsr
type(psb_d_csr_sparse_mat),intent(inout) :: bcsr
type(psb_d_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
End Subroutine psb_d_par_csr_spspmm
end interface psb_par_spspmm
end module psb_d_tools_mod end module psb_d_tools_mod

@ -363,7 +363,6 @@ Module psb_s_tools_mod
end subroutine psb_sspins_2desc end subroutine psb_sspins_2desc
end interface end interface
interface psb_sprn interface psb_sprn
subroutine psb_ssprn(a, desc_a,info,clear) subroutine psb_ssprn(a, desc_a,info,clear)
import import
@ -374,5 +373,19 @@ Module psb_s_tools_mod
logical, intent(in), optional :: clear logical, intent(in), optional :: clear
end subroutine psb_ssprn end subroutine psb_ssprn
end interface end interface
interface psb_par_spspmm
subroutine psb_s_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
import :: psb_s_csr_sparse_mat, psb_desc_type, psb_ipk_
Implicit None
type(psb_s_csr_sparse_mat),intent(in) :: acsr
type(psb_s_csr_sparse_mat),intent(inout) :: bcsr
type(psb_s_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
End Subroutine psb_s_par_csr_spspmm
end interface psb_par_spspmm
end module psb_s_tools_mod end module psb_s_tools_mod

@ -363,7 +363,6 @@ Module psb_z_tools_mod
end subroutine psb_zspins_2desc end subroutine psb_zspins_2desc
end interface end interface
interface psb_sprn interface psb_sprn
subroutine psb_zsprn(a, desc_a,info,clear) subroutine psb_zsprn(a, desc_a,info,clear)
import import
@ -374,5 +373,19 @@ Module psb_z_tools_mod
logical, intent(in), optional :: clear logical, intent(in), optional :: clear
end subroutine psb_zsprn end subroutine psb_zsprn
end interface end interface
interface psb_par_spspmm
subroutine psb_z_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
import :: psb_z_csr_sparse_mat, psb_desc_type, psb_ipk_
Implicit None
type(psb_z_csr_sparse_mat),intent(in) :: acsr
type(psb_z_csr_sparse_mat),intent(inout) :: bcsr
type(psb_z_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
End Subroutine psb_z_par_csr_spspmm
end interface psb_par_spspmm
end module psb_z_tools_mod end module psb_z_tools_mod

@ -565,7 +565,7 @@ subroutine psb_cspsv_vect(alpha,a,x,beta,y,desc_a,info,&
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
complex(psb_spk_),pointer :: iwork(:), xp(:), yp(:) complex(psb_spk_) :: iwork(1)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
@ -636,32 +636,10 @@ subroutine psb_cspsv_vect(alpha,a,x,beta,y,desc_a,info,&
goto 9999 goto 9999
end if end if
iwork => null() ! With encapsulated vectors the inner routines
! check for presence/size of a work area ! do not need the work area, hence liwork is
liwork= 2*ncol ! a simple array with 1 entry to keep the interface
! happy. Might remove it entirely in the future.
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
else
aliw=.true.
endif
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
iwork => work
endif
iwork(1)=0.d0 iwork(1)=0.d0
! Perform local triangular system solve ! Perform local triangular system solve
@ -682,7 +660,6 @@ subroutine psb_cspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& cone,y%v,desc_a,iwork,info,data=psb_comm_ovr_) & cone,y%v,desc_a,iwork,info,data=psb_comm_ovr_)
if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info) if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates') call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
@ -690,8 +667,6 @@ subroutine psb_cspsv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
end if end if
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -565,7 +565,7 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_dpk_),pointer :: iwork(:), xp(:), yp(:) real(psb_dpk_) :: iwork(1)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
@ -636,32 +636,10 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
goto 9999 goto 9999
end if end if
iwork => null() ! With encapsulated vectors the inner routines
! check for presence/size of a work area ! do not need the work area, hence liwork is
liwork= 2*ncol ! a simple array with 1 entry to keep the interface
! happy. Might remove it entirely in the future.
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
else
aliw=.true.
endif
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
iwork => work
endif
iwork(1)=0.d0 iwork(1)=0.d0
! Perform local triangular system solve ! Perform local triangular system solve
@ -682,7 +660,6 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& done,y%v,desc_a,iwork,info,data=psb_comm_ovr_) & done,y%v,desc_a,iwork,info,data=psb_comm_ovr_)
if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info) if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates') call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
@ -690,8 +667,6 @@ subroutine psb_dspsv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
end if end if
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -565,7 +565,7 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
real(psb_spk_),pointer :: iwork(:), xp(:), yp(:) real(psb_spk_) :: iwork(1)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
@ -636,32 +636,10 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
goto 9999 goto 9999
end if end if
iwork => null() ! With encapsulated vectors the inner routines
! check for presence/size of a work area ! do not need the work area, hence liwork is
liwork= 2*ncol ! a simple array with 1 entry to keep the interface
! happy. Might remove it entirely in the future.
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
else
aliw=.true.
endif
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
iwork => work
endif
iwork(1)=0.d0 iwork(1)=0.d0
! Perform local triangular system solve ! Perform local triangular system solve
@ -682,7 +660,6 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& sone,y%v,desc_a,iwork,info,data=psb_comm_ovr_) & sone,y%v,desc_a,iwork,info,data=psb_comm_ovr_)
if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info) if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates') call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
@ -690,8 +667,6 @@ subroutine psb_sspsv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
end if end if
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -565,7 +565,7 @@ subroutine psb_zspsv_vect(alpha,a,x,beta,y,desc_a,info,&
character :: lscale character :: lscale
integer(psb_ipk_), parameter :: nb=4 integer(psb_ipk_), parameter :: nb=4
complex(psb_dpk_),pointer :: iwork(:), xp(:), yp(:) complex(psb_dpk_) :: iwork(1)
character :: itrans character :: itrans
character(len=20) :: name, ch_err character(len=20) :: name, ch_err
logical :: aliw logical :: aliw
@ -636,32 +636,10 @@ subroutine psb_zspsv_vect(alpha,a,x,beta,y,desc_a,info,&
goto 9999 goto 9999
end if end if
iwork => null() ! With encapsulated vectors the inner routines
! check for presence/size of a work area ! do not need the work area, hence liwork is
liwork= 2*ncol ! a simple array with 1 entry to keep the interface
! happy. Might remove it entirely in the future.
if (present(work)) then
if (size(work) >= liwork) then
aliw =.false.
else
aliw=.true.
endif
else
aliw=.true.
end if
if (aliw) then
allocate(iwork(liwork),stat=info)
if(info /= psb_success_) then
info=psb_err_from_subroutine_
ch_err='psb_realloc'
call psb_errpush(info,name,a_err=ch_err)
goto 9999
end if
else
iwork => work
endif
iwork(1)=0.d0 iwork(1)=0.d0
! Perform local triangular system solve ! Perform local triangular system solve
@ -682,7 +660,6 @@ subroutine psb_zspsv_vect(alpha,a,x,beta,y,desc_a,info,&
call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),& call psi_swapdata(ior(psb_swap_send_,psb_swap_recv_),&
& zone,y%v,desc_a,iwork,info,data=psb_comm_ovr_) & zone,y%v,desc_a,iwork,info,data=psb_comm_ovr_)
if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info) if (info == psb_success_) call psi_ovrl_upd(y%v,desc_a,choice_,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates') call psb_errpush(psb_err_from_subroutine_,name,a_err='Inner updates')
@ -690,8 +667,6 @@ subroutine psb_zspsv_vect(alpha,a,x,beta,y,desc_a,info,&
end if end if
end if end if
if (aliw) deallocate(iwork)
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return

@ -1184,17 +1184,28 @@ subroutine psb_c_base_csmm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='c_base_csmm' integer(psb_ipk_) :: j,nc
character(len=20) :: name='c_base_csmm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spmm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_error_handler(err_act) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
end subroutine psb_c_base_csmm end subroutine psb_c_base_csmm
@ -1266,17 +1277,29 @@ subroutine psb_c_base_inner_cssm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: j, nc
character(len=20) :: name='c_base_inner_cssm' character(len=20) :: name='c_base_inner_cssm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spsm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
call psb_error_handler(err_act)
end subroutine psb_c_base_inner_cssm end subroutine psb_c_base_inner_cssm

File diff suppressed because it is too large Load Diff

@ -1184,17 +1184,28 @@ subroutine psb_d_base_csmm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='d_base_csmm' integer(psb_ipk_) :: j,nc
character(len=20) :: name='d_base_csmm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spmm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_error_handler(err_act) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
end subroutine psb_d_base_csmm end subroutine psb_d_base_csmm
@ -1266,17 +1277,29 @@ subroutine psb_d_base_inner_cssm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: j, nc
character(len=20) :: name='d_base_inner_cssm' character(len=20) :: name='d_base_inner_cssm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spsm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
call psb_error_handler(err_act)
end subroutine psb_d_base_inner_cssm end subroutine psb_d_base_inner_cssm

File diff suppressed because it is too large Load Diff

@ -1184,17 +1184,28 @@ subroutine psb_s_base_csmm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='s_base_csmm' integer(psb_ipk_) :: j,nc
character(len=20) :: name='s_base_csmm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spmm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_error_handler(err_act) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
end subroutine psb_s_base_csmm end subroutine psb_s_base_csmm
@ -1266,17 +1277,29 @@ subroutine psb_s_base_inner_cssm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: j, nc
character(len=20) :: name='s_base_inner_cssm' character(len=20) :: name='s_base_inner_cssm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spsm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
call psb_error_handler(err_act)
end subroutine psb_s_base_inner_cssm end subroutine psb_s_base_inner_cssm

File diff suppressed because it is too large Load Diff

@ -1184,17 +1184,28 @@ subroutine psb_z_base_csmm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='z_base_csmm' integer(psb_ipk_) :: j,nc
character(len=20) :: name='z_base_csmm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spmm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_error_handler(err_act) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
end subroutine psb_z_base_csmm end subroutine psb_z_base_csmm
@ -1266,17 +1277,29 @@ subroutine psb_z_base_inner_cssm(alpha,a,x,beta,y,info,trans)
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
integer(psb_ipk_) :: ierr(5) integer(psb_ipk_) :: ierr(5)
integer(psb_ipk_) :: j, nc
character(len=20) :: name='z_base_inner_cssm' character(len=20) :: name='z_base_inner_cssm'
logical, parameter :: debug=.false. logical, parameter :: debug=.false.
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
! This is the base version. If we get here ! This is the base version.
! it means the derived class is incomplete, ! It's a very inefficient implementation,
! so we throw an error. ! but it's only a fallback, if multivectors
info = psb_err_missing_override_method_ ! are important you are supposed to implement it
call psb_errpush(info,name,a_err=a%get_fmt()) ! explicitly in the derived class.
info = psb_success_
nc = min(size(x,2),size(y,2))
do j=1,nc
call a%spsm(alpha,x(j,:),beta,y(:,j),info,trans)
if (info /= psb_success_) goto 9999
end do
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
call psb_error_handler(err_act)
end subroutine psb_z_base_inner_cssm end subroutine psb_z_base_inner_cssm

File diff suppressed because it is too large Load Diff

@ -19,7 +19,10 @@ FOBJS = psb_sallc.o psb_sasb.o \
psb_cspalloc.o psb_cspasb.o psb_cspfree.o\ psb_cspalloc.o psb_cspasb.o psb_cspfree.o\
psb_callc.o psb_casb.o psb_cfree.o psb_cins.o \ psb_callc.o psb_casb.o psb_cfree.o psb_cins.o \
psb_cspins.o psb_csprn.o psb_cd_set_bld.o \ psb_cspins.o psb_csprn.o psb_cd_set_bld.o \
psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o psb_s_map.o psb_d_map.o psb_c_map.o psb_z_map.o \
psb_s_par_csr_spspmm.o psb_d_par_csr_spspmm.o \
psb_c_par_csr_spspmm.o psb_z_par_csr_spspmm.o
MPFOBJS = psb_icdasb.o psb_ssphalo.o psb_dsphalo.o psb_csphalo.o psb_zsphalo.o \ MPFOBJS = psb_icdasb.o psb_ssphalo.o psb_dsphalo.o psb_csphalo.o psb_zsphalo.o \
psb_dcdbldext.o psb_zcdbldext.o psb_scdbldext.o psb_ccdbldext.o psb_dcdbldext.o psb_zcdbldext.o psb_scdbldext.o psb_ccdbldext.o

@ -64,7 +64,7 @@ subroutine psb_c_map_U2V_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
nr2 = map%p_desc_V%get_global_rows() nr2 = map%p_desc_V%get_global_rows()
@ -104,7 +104,7 @@ subroutine psb_c_map_U2V_a(alpha,x,beta,y,map,info,work)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -138,7 +138,7 @@ subroutine psb_c_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -205,7 +205,7 @@ subroutine psb_c_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -246,7 +246,7 @@ subroutine psb_c_map_V2U_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
nr2 = map%p_desc_U%get_global_rows() nr2 = map%p_desc_U%get_global_rows()
@ -319,7 +319,7 @@ subroutine psb_c_map_V2U_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -408,7 +408,7 @@ function psb_c_linmap(map_kind,desc_U, desc_V, mat_U2V, mat_V2U,iaggr,naggr) &
info = psb_success_ info = psb_success_
select case(map_kind) select case(map_kind)
case (psb_map_aggr_) case (psb_map_dec_aggr_)
! OK ! OK
if (psb_is_ok_desc(desc_U)) then if (psb_is_ok_desc(desc_U)) then

@ -0,0 +1,156 @@
!
! 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.
!
!
! File: psb_c_par_csr_spspmm.f90
!
! Subroutine: psb_c_par_csr_spspmm
! Version: complex
!
! This routine computes a parallel product of two sparse matrices
!
! C = A * B
!
! where all the matrices are stored in CSR. On input and output the matrices
! are stored with column indices in local numbering, but inermediate quantities
! are in global numbering because gathering the halo of B to multiply it
! by A implies a potential enlargement of the support.
! Also, B may have a column index space different from its row index space,
! which is obviously the same as the column space of A.
!
!
! Arguments:
! acsr - type(psb_c_csr_sparse_mat), input.
! The sparse matrix structure A
! desc_a - type(psb_desc_type), input.
! The communication descriptor of the column space of A
! bcsr - type(psb_c_csr_sparse_mat), input/output.
! The sparse matrix structure B, gets row-extended on output
! ccsr - type(psb_c_csr_sparse_mat), output
! The sparse matrix structure C
! desc_c - type(psb_desc_type), input/output.
! The communication descriptor of the column space of B
!
! info - integer, output.
! Error code.
!
Subroutine psb_c_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
use psb_base_mod, psb_protect_name => psb_c_par_csr_spspmm
Implicit None
type(psb_c_csr_sparse_mat),intent(in) :: acsr
type(psb_c_csr_sparse_mat),intent(inout) :: bcsr
type(psb_c_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
! ...local scalars....
integer(psb_ipk_) :: ictxt, np,me
integer(psb_ipk_) :: ncol, nnz
type(psb_c_csr_sparse_mat) :: tcsr1
logical :: update_desc_c
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name='psb_c_p_csr_spspmm'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
update_desc_c = desc_c%is_bld()
!
! This is a bit tricky.
! DESC_A is the descriptor of (the columns of) A, and therefore
! of the rows of B; the columns of B, in the intended usage, span
! a different space for which we have DESC_C.
! We are gathering the halo rows of B to multiply by A;
! now, the columns of B would ideally be kept in
! global numbering, so that we can call this repeatedly to accumulate
! the product of multiple operators, and convert to local numbering
! at the last possible moment. However, this would imply calling
! the serial SPSPMM with a matrix B with the GLOBAL number of columns
! and this could be very expensive in memory. The solution is to keep B
! in local numbering, so that only columns really appearing count, but to
! expand the descriptor when gathering the halo, because by performing
! the products we are extending the support of the operator; hence
! this routine is intended to be called with a temporary descriptor
! DESC_C which is in the BUILD state, to allow for such expansion
! across multiple products.
! The caller will at some later point finalize the descriptor DESC_C.
!
ncol = desc_a%get_local_cols()
call psb_sphalo(bcsr,desc_a,tcsr1,info,&
& colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data)
nnz = tcsr1%get_nzeros()
if (update_desc_c) then
call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info)
else
call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info)
end if
if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1)
if (info == psb_success_) call tcsr1%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call bcsr%set_ncols(desc_c%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols()
call psb_spspmm(acsr,bcsr,ccsr,info)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
End Subroutine psb_c_par_csr_spspmm

@ -64,7 +64,7 @@ subroutine psb_d_map_U2V_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
nr2 = map%p_desc_V%get_global_rows() nr2 = map%p_desc_V%get_global_rows()
@ -104,7 +104,7 @@ subroutine psb_d_map_U2V_a(alpha,x,beta,y,map,info,work)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -138,7 +138,7 @@ subroutine psb_d_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -205,7 +205,7 @@ subroutine psb_d_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -246,7 +246,7 @@ subroutine psb_d_map_V2U_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
nr2 = map%p_desc_U%get_global_rows() nr2 = map%p_desc_U%get_global_rows()
@ -319,7 +319,7 @@ subroutine psb_d_map_V2U_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -408,7 +408,7 @@ function psb_d_linmap(map_kind,desc_U, desc_V, mat_U2V, mat_V2U,iaggr,naggr) &
info = psb_success_ info = psb_success_
select case(map_kind) select case(map_kind)
case (psb_map_aggr_) case (psb_map_dec_aggr_)
! OK ! OK
if (psb_is_ok_desc(desc_U)) then if (psb_is_ok_desc(desc_U)) then

@ -0,0 +1,156 @@
!
! 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.
!
!
! File: psb_d_par_csr_spspmm.f90
!
! Subroutine: psb_d_par_csr_spspmm
! Version: real
!
! This routine computes a parallel product of two sparse matrices
!
! C = A * B
!
! where all the matrices are stored in CSR. On input and output the matrices
! are stored with column indices in local numbering, but inermediate quantities
! are in global numbering because gathering the halo of B to multiply it
! by A implies a potential enlargement of the support.
! Also, B may have a column index space different from its row index space,
! which is obviously the same as the column space of A.
!
!
! Arguments:
! acsr - type(psb_d_csr_sparse_mat), input.
! The sparse matrix structure A
! desc_a - type(psb_desc_type), input.
! The communication descriptor of the column space of A
! bcsr - type(psb_d_csr_sparse_mat), input/output.
! The sparse matrix structure B, gets row-extended on output
! ccsr - type(psb_d_csr_sparse_mat), output
! The sparse matrix structure C
! desc_c - type(psb_desc_type), input/output.
! The communication descriptor of the column space of B
!
! info - integer, output.
! Error code.
!
Subroutine psb_d_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
use psb_base_mod, psb_protect_name => psb_d_par_csr_spspmm
Implicit None
type(psb_d_csr_sparse_mat),intent(in) :: acsr
type(psb_d_csr_sparse_mat),intent(inout) :: bcsr
type(psb_d_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
! ...local scalars....
integer(psb_ipk_) :: ictxt, np,me
integer(psb_ipk_) :: ncol, nnz
type(psb_d_csr_sparse_mat) :: tcsr1
logical :: update_desc_c
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name='psb_d_p_csr_spspmm'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
update_desc_c = desc_c%is_bld()
!
! This is a bit tricky.
! DESC_A is the descriptor of (the columns of) A, and therefore
! of the rows of B; the columns of B, in the intended usage, span
! a different space for which we have DESC_C.
! We are gathering the halo rows of B to multiply by A;
! now, the columns of B would ideally be kept in
! global numbering, so that we can call this repeatedly to accumulate
! the product of multiple operators, and convert to local numbering
! at the last possible moment. However, this would imply calling
! the serial SPSPMM with a matrix B with the GLOBAL number of columns
! and this could be very expensive in memory. The solution is to keep B
! in local numbering, so that only columns really appearing count, but to
! expand the descriptor when gathering the halo, because by performing
! the products we are extending the support of the operator; hence
! this routine is intended to be called with a temporary descriptor
! DESC_C which is in the BUILD state, to allow for such expansion
! across multiple products.
! The caller will at some later point finalize the descriptor DESC_C.
!
ncol = desc_a%get_local_cols()
call psb_sphalo(bcsr,desc_a,tcsr1,info,&
& colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data)
nnz = tcsr1%get_nzeros()
if (update_desc_c) then
call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info)
else
call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info)
end if
if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1)
if (info == psb_success_) call tcsr1%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call bcsr%set_ncols(desc_c%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols()
call psb_spspmm(acsr,bcsr,ccsr,info)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
End Subroutine psb_d_par_csr_spspmm

@ -64,7 +64,7 @@ subroutine psb_s_map_U2V_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
nr2 = map%p_desc_V%get_global_rows() nr2 = map%p_desc_V%get_global_rows()
@ -104,7 +104,7 @@ subroutine psb_s_map_U2V_a(alpha,x,beta,y,map,info,work)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -138,7 +138,7 @@ subroutine psb_s_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -205,7 +205,7 @@ subroutine psb_s_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -246,7 +246,7 @@ subroutine psb_s_map_V2U_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
nr2 = map%p_desc_U%get_global_rows() nr2 = map%p_desc_U%get_global_rows()
@ -319,7 +319,7 @@ subroutine psb_s_map_V2U_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -408,7 +408,7 @@ function psb_s_linmap(map_kind,desc_U, desc_V, mat_U2V, mat_V2U,iaggr,naggr) &
info = psb_success_ info = psb_success_
select case(map_kind) select case(map_kind)
case (psb_map_aggr_) case (psb_map_dec_aggr_)
! OK ! OK
if (psb_is_ok_desc(desc_U)) then if (psb_is_ok_desc(desc_U)) then

@ -0,0 +1,156 @@
!
! 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.
!
!
! File: psb_s_par_csr_spspmm.f90
!
! Subroutine: psb_s_par_csr_spspmm
! Version: real
!
! This routine computes a parallel product of two sparse matrices
!
! C = A * B
!
! where all the matrices are stored in CSR. On input and output the matrices
! are stored with column indices in local numbering, but inermediate quantities
! are in global numbering because gathering the halo of B to multiply it
! by A implies a potential enlargement of the support.
! Also, B may have a column index space different from its row index space,
! which is obviously the same as the column space of A.
!
!
! Arguments:
! acsr - type(psb_s_csr_sparse_mat), input.
! The sparse matrix structure A
! desc_a - type(psb_desc_type), input.
! The communication descriptor of the column space of A
! bcsr - type(psb_s_csr_sparse_mat), input/output.
! The sparse matrix structure B, gets row-extended on output
! ccsr - type(psb_s_csr_sparse_mat), output
! The sparse matrix structure C
! desc_c - type(psb_desc_type), input/output.
! The communication descriptor of the column space of B
!
! info - integer, output.
! Error code.
!
Subroutine psb_s_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
use psb_base_mod, psb_protect_name => psb_s_par_csr_spspmm
Implicit None
type(psb_s_csr_sparse_mat),intent(in) :: acsr
type(psb_s_csr_sparse_mat),intent(inout) :: bcsr
type(psb_s_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
! ...local scalars....
integer(psb_ipk_) :: ictxt, np,me
integer(psb_ipk_) :: ncol, nnz
type(psb_s_csr_sparse_mat) :: tcsr1
logical :: update_desc_c
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name='psb_s_p_csr_spspmm'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
update_desc_c = desc_c%is_bld()
!
! This is a bit tricky.
! DESC_A is the descriptor of (the columns of) A, and therefore
! of the rows of B; the columns of B, in the intended usage, span
! a different space for which we have DESC_C.
! We are gathering the halo rows of B to multiply by A;
! now, the columns of B would ideally be kept in
! global numbering, so that we can call this repeatedly to accumulate
! the product of multiple operators, and convert to local numbering
! at the last possible moment. However, this would imply calling
! the serial SPSPMM with a matrix B with the GLOBAL number of columns
! and this could be very expensive in memory. The solution is to keep B
! in local numbering, so that only columns really appearing count, but to
! expand the descriptor when gathering the halo, because by performing
! the products we are extending the support of the operator; hence
! this routine is intended to be called with a temporary descriptor
! DESC_C which is in the BUILD state, to allow for such expansion
! across multiple products.
! The caller will at some later point finalize the descriptor DESC_C.
!
ncol = desc_a%get_local_cols()
call psb_sphalo(bcsr,desc_a,tcsr1,info,&
& colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data)
nnz = tcsr1%get_nzeros()
if (update_desc_c) then
call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info)
else
call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info)
end if
if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1)
if (info == psb_success_) call tcsr1%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call bcsr%set_ncols(desc_c%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols()
call psb_spspmm(acsr,bcsr,ccsr,info)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
End Subroutine psb_s_par_csr_spspmm

@ -64,7 +64,7 @@ subroutine psb_z_map_U2V_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
nr2 = map%p_desc_V%get_global_rows() nr2 = map%p_desc_V%get_global_rows()
@ -104,7 +104,7 @@ subroutine psb_z_map_U2V_a(alpha,x,beta,y,map,info,work)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -138,7 +138,7 @@ subroutine psb_z_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_V%get_context() ictxt = map%p_desc_V%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -205,7 +205,7 @@ subroutine psb_z_map_U2V_v(alpha,x,beta,y,map,info,work,vtx,vty)
case default case default
write(psb_err_unit,*) trim(name),' Invalid descriptor input', & write(psb_err_unit,*) trim(name),' Invalid descriptor input', &
& map_kind, psb_map_aggr_, psb_map_gen_linear_ & map_kind, psb_map_dec_aggr_, psb_map_gen_linear_
info = 1 info = 1
return return
end select end select
@ -246,7 +246,7 @@ subroutine psb_z_map_V2U_a(alpha,x,beta,y,map,info,work)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
nr2 = map%p_desc_U%get_global_rows() nr2 = map%p_desc_U%get_global_rows()
@ -319,7 +319,7 @@ subroutine psb_z_map_V2U_v(alpha,x,beta,y,map,info,work,vtx,vty)
map_kind = map%get_kind() map_kind = map%get_kind()
select case(map_kind) select case(map_kind)
case(psb_map_aggr_) case(psb_map_dec_aggr_)
ictxt = map%p_desc_U%get_context() ictxt = map%p_desc_U%get_context()
call psb_info(ictxt,iam,np) call psb_info(ictxt,iam,np)
@ -408,7 +408,7 @@ function psb_z_linmap(map_kind,desc_U, desc_V, mat_U2V, mat_V2U,iaggr,naggr) &
info = psb_success_ info = psb_success_
select case(map_kind) select case(map_kind)
case (psb_map_aggr_) case (psb_map_dec_aggr_)
! OK ! OK
if (psb_is_ok_desc(desc_U)) then if (psb_is_ok_desc(desc_U)) then

@ -0,0 +1,156 @@
!
! 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.
!
!
! File: psb_z_par_csr_spspmm.f90
!
! Subroutine: psb_z_par_csr_spspmm
! Version: complex
!
! This routine computes a parallel product of two sparse matrices
!
! C = A * B
!
! where all the matrices are stored in CSR. On input and output the matrices
! are stored with column indices in local numbering, but inermediate quantities
! are in global numbering because gathering the halo of B to multiply it
! by A implies a potential enlargement of the support.
! Also, B may have a column index space different from its row index space,
! which is obviously the same as the column space of A.
!
!
! Arguments:
! acsr - type(psb_z_csr_sparse_mat), input.
! The sparse matrix structure A
! desc_a - type(psb_desc_type), input.
! The communication descriptor of the column space of A
! bcsr - type(psb_z_csr_sparse_mat), input/output.
! The sparse matrix structure B, gets row-extended on output
! ccsr - type(psb_z_csr_sparse_mat), output
! The sparse matrix structure C
! desc_c - type(psb_desc_type), input/output.
! The communication descriptor of the column space of B
!
! info - integer, output.
! Error code.
!
Subroutine psb_z_par_csr_spspmm(acsr,desc_a,bcsr,ccsr,desc_c,info,data)
use psb_base_mod, psb_protect_name => psb_z_par_csr_spspmm
Implicit None
type(psb_z_csr_sparse_mat),intent(in) :: acsr
type(psb_z_csr_sparse_mat),intent(inout) :: bcsr
type(psb_z_csr_sparse_mat),intent(out) :: ccsr
type(psb_desc_type),intent(in) :: desc_a
type(psb_desc_type),intent(inout) :: desc_c
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_), intent(in), optional :: data
! ...local scalars....
integer(psb_ipk_) :: ictxt, np,me
integer(psb_ipk_) :: ncol, nnz
type(psb_z_csr_sparse_mat) :: tcsr1
logical :: update_desc_c
integer(psb_ipk_) :: debug_level, debug_unit, err_act
character(len=20) :: name, ch_err
if(psb_get_errstatus() /= 0) return
info=psb_success_
name='psb_z_p_csr_spspmm'
call psb_erractionsave(err_act)
if (psb_errstatus_fatal()) then
info = psb_err_internal_error_ ; goto 9999
end if
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),': Start'
update_desc_c = desc_c%is_bld()
!
! This is a bit tricky.
! DESC_A is the descriptor of (the columns of) A, and therefore
! of the rows of B; the columns of B, in the intended usage, span
! a different space for which we have DESC_C.
! We are gathering the halo rows of B to multiply by A;
! now, the columns of B would ideally be kept in
! global numbering, so that we can call this repeatedly to accumulate
! the product of multiple operators, and convert to local numbering
! at the last possible moment. However, this would imply calling
! the serial SPSPMM with a matrix B with the GLOBAL number of columns
! and this could be very expensive in memory. The solution is to keep B
! in local numbering, so that only columns really appearing count, but to
! expand the descriptor when gathering the halo, because by performing
! the products we are extending the support of the operator; hence
! this routine is intended to be called with a temporary descriptor
! DESC_C which is in the BUILD state, to allow for such expansion
! across multiple products.
! The caller will at some later point finalize the descriptor DESC_C.
!
ncol = desc_a%get_local_cols()
call psb_sphalo(bcsr,desc_a,tcsr1,info,&
& colcnv=.true.,rowscale=.true.,outcol_glob=.true.,col_desc=desc_c,data=data)
nnz = tcsr1%get_nzeros()
if (update_desc_c) then
call desc_c%indxmap%g2lip_ins(tcsr1%ja(1:nnz),info)
else
call desc_c%indxmap%g2lip(tcsr1%ja(1:nnz),info)
end if
if (info == psb_success_) call psb_rwextd(ncol,bcsr,info,b=tcsr1)
if (info == psb_success_) call tcsr1%free()
if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')
goto 9999
end if
call bcsr%set_ncols(desc_c%get_local_cols())
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm 3'
if (debug_level >= psb_debug_outer_) write(debug_unit,*) me,' ',trim(name),&
& 'starting spspmm ',acsr%get_nrows(),acsr%get_ncols(),bcsr%get_nrows(),bcsr%get_ncols()
call psb_spspmm(acsr,bcsr,ccsr,info)
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(ictxt,err_act)
return
End Subroutine psb_z_par_csr_spspmm
Loading…
Cancel
Save