fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Implement TRIL/TRIU for CSR.

Browse Source
psblas-3.5-maint
Salvatore Filippone 7 years ago
parent 663fa5b8e5
commit 42d5a85e25
8 changed files with 1608 additions and 0 deletions
Show Stats Download Patch File Download Diff File

89
base/modules/serial/psb_c_csr_mat_mod.f90
Unescape Escape View File

@ -80,6 +80,8 @@ module psb_c_csr_mat_mod
procedure, pass(a) :: aclsum => psb_c_csr_aclsum
procedure, pass(a) :: reallocate_nz => psb_c_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_c_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_c_csr_tril
procedure, pass(a) :: triu => psb_c_csr_triu
procedure, pass(a) :: cp_to_coo => psb_c_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_c_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_c_cp_csr_to_fmt
@ -170,6 +172,93 @@ module psb_c_csr_mat_mod
integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_c_csr_print
end interface
!
!> Function tril:
!! \memberof psb_c_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param l the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param u [none] copy of the complementary triangle
!!
!
interface
subroutine psb_c_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
class(psb_c_csr_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_c_csr_tril
end interface
!
!> Function triu:
!! \memberof psb_c_csr_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!! Optionally copies the lower triangle at the same time
!!
!! \param u the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param l [none] copy of the complementary triangle
!!
!
interface
subroutine psb_c_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
class(psb_c_csr_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_c_csr_triu
end interface
!> \memberof psb_c_csr_sparse_mat
!! \see psb_c_base_mat_mod::psb_c_base_cp_to_coo

89
base/modules/serial/psb_d_csr_mat_mod.f90
Unescape Escape View File

@ -80,6 +80,8 @@ module psb_d_csr_mat_mod
procedure, pass(a) :: aclsum => psb_d_csr_aclsum
procedure, pass(a) :: reallocate_nz => psb_d_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_d_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_d_csr_tril
procedure, pass(a) :: triu => psb_d_csr_triu
procedure, pass(a) :: cp_to_coo => psb_d_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_d_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_d_cp_csr_to_fmt
@ -170,6 +172,93 @@ module psb_d_csr_mat_mod
integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_d_csr_print
end interface
!
!> Function tril:
!! \memberof psb_d_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param l the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param u [none] copy of the complementary triangle
!!
!
interface
subroutine psb_d_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import :: psb_ipk_, psb_d_csr_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
class(psb_d_csr_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_d_csr_tril
end interface
!
!> Function triu:
!! \memberof psb_d_csr_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!! Optionally copies the lower triangle at the same time
!!
!! \param u the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param l [none] copy of the complementary triangle
!!
!
interface
subroutine psb_d_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import :: psb_ipk_, psb_d_csr_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
class(psb_d_csr_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_d_csr_triu
end interface
!> \memberof psb_d_csr_sparse_mat
!! \see psb_d_base_mat_mod::psb_d_base_cp_to_coo

89
base/modules/serial/psb_s_csr_mat_mod.f90
Unescape Escape View File

@ -80,6 +80,8 @@ module psb_s_csr_mat_mod
procedure, pass(a) :: aclsum => psb_s_csr_aclsum
procedure, pass(a) :: reallocate_nz => psb_s_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_s_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_s_csr_tril
procedure, pass(a) :: triu => psb_s_csr_triu
procedure, pass(a) :: cp_to_coo => psb_s_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_s_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_s_cp_csr_to_fmt
@ -170,6 +172,93 @@ module psb_s_csr_mat_mod
integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_s_csr_print
end interface
!
!> Function tril:
!! \memberof psb_s_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param l the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param u [none] copy of the complementary triangle
!!
!
interface
subroutine psb_s_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import :: psb_ipk_, psb_s_csr_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
class(psb_s_csr_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_s_csr_tril
end interface
!
!> Function triu:
!! \memberof psb_s_csr_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!! Optionally copies the lower triangle at the same time
!!
!! \param u the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param l [none] copy of the complementary triangle
!!
!
interface
subroutine psb_s_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import :: psb_ipk_, psb_s_csr_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
class(psb_s_csr_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_s_csr_triu
end interface
!> \memberof psb_s_csr_sparse_mat
!! \see psb_s_base_mat_mod::psb_s_base_cp_to_coo

89
base/modules/serial/psb_z_csr_mat_mod.f90
Unescape Escape View File

@ -80,6 +80,8 @@ module psb_z_csr_mat_mod
procedure, pass(a) :: aclsum => psb_z_csr_aclsum
procedure, pass(a) :: reallocate_nz => psb_z_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_z_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_z_csr_tril
procedure, pass(a) :: triu => psb_z_csr_triu
procedure, pass(a) :: cp_to_coo => psb_z_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_z_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_z_cp_csr_to_fmt
@ -170,6 +172,93 @@ module psb_z_csr_mat_mod
integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_z_csr_print
end interface
!
!> Function tril:
!! \memberof psb_z_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param l the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param u [none] copy of the complementary triangle
!!
!
interface
subroutine psb_z_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import :: psb_ipk_, psb_z_csr_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
class(psb_z_csr_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_z_csr_tril
end interface
!
!> Function triu:
!! \memberof psb_z_csr_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!! Optionally copies the lower triangle at the same time
!!
!! \param u the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!! \param l [none] copy of the complementary triangle
!!
!
interface
subroutine psb_z_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import :: psb_ipk_, psb_z_csr_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
class(psb_z_csr_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_z_csr_triu
end interface
!> \memberof psb_z_csr_sparse_mat
!! \see psb_z_base_mat_mod::psb_z_base_cp_to_coo

313
base/serial/impl/psb_c_csr_impl.f90
Unescape Escape View File

@ -2233,6 +2233,319 @@ subroutine psb_c_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_c_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_c_csr_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_c_csr_mat_mod, psb_protect_name => psb_c_csr_tril
implicit none
class(psb_c_csr_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = ja(k)
l%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzin)
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_c_csr_tril
subroutine psb_c_csr_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_c_csr_mat_mod, psb_protect_name => psb_c_csr_triu
implicit none
class(psb_c_csr_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_c_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = ja(k)
u%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzin)
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_c_csr_triu
subroutine psb_c_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

313
base/serial/impl/psb_d_csr_impl.f90
Unescape Escape View File

@ -2233,6 +2233,319 @@ subroutine psb_d_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_d_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_d_csr_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_d_csr_mat_mod, psb_protect_name => psb_d_csr_tril
implicit none
class(psb_d_csr_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = ja(k)
l%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzin)
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_d_csr_tril
subroutine psb_d_csr_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_d_csr_mat_mod, psb_protect_name => psb_d_csr_triu
implicit none
class(psb_d_csr_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_d_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = ja(k)
u%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzin)
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_d_csr_triu
subroutine psb_d_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

313
base/serial/impl/psb_s_csr_impl.f90
Unescape Escape View File

@ -2233,6 +2233,319 @@ subroutine psb_s_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_s_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_s_csr_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_s_csr_mat_mod, psb_protect_name => psb_s_csr_tril
implicit none
class(psb_s_csr_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = ja(k)
l%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzin)
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_s_csr_tril
subroutine psb_s_csr_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_s_csr_mat_mod, psb_protect_name => psb_s_csr_triu
implicit none
class(psb_s_csr_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_s_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = ja(k)
u%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzin)
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_s_csr_triu
subroutine psb_s_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

313
base/serial/impl/psb_z_csr_impl.f90
Unescape Escape View File

@ -2233,6 +2233,319 @@ subroutine psb_z_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_z_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_z_csr_tril(a,l,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,u)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_z_csr_mat_mod, psb_protect_name => psb_z_csr_tril
implicit none
class(psb_z_csr_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
if (present(u)) then
nzlin = l%get_nzeros() ! At this point it should be 0
call u%allocate(mb,nb,nz)
nzuin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzlin)
call u%set_nzeros(nzuin)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
nzin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = ja(k)
l%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call l%set_nzeros(nzin)
end if
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_z_csr_tril
subroutine psb_z_csr_triu(a,u,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale,l)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_z_csr_mat_mod, psb_protect_name => psb_z_csr_triu
implicit none
class(psb_z_csr_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_z_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
if (present(l)) then
nzuin = u%get_nzeros() ! At this point it should be 0
call l%allocate(mb,nb,nz)
nzlin = l%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
else
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, irp=>a%irp)
do i=imin_,imax_
do k=irp(i),irp(i+1)-1
if ((jmin_<=j).and.(j<=jmax_)) then
if ((ja(k)-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = ja(k)
u%val(nzin) = val(k)
end if
end if
end do
end do
end associate
call u%set_nzeros(nzin)
end if
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) &
& u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
if (cscale_) &
& u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_z_csr_triu
subroutine psb_z_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

Write Preview
Loading…
Cancel
Save