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

Merged new TRIL/TRIU from development.

Browse Source
ILmat
Salvatore Filippone 8 years ago
parent 1ca13b2a3c
commit deea842017
12 changed files with 3324 additions and 84 deletions
Show Stats Download Patch File Download Diff File

180
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
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
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
@ -518,6 +607,8 @@ module psb_c_csr_mat_mod
procedure, pass(a) :: sizeof => lc_csr_sizeof
procedure, pass(a) :: reallocate_nz => psb_lc_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_lc_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_lc_csr_tril
procedure, pass(a) :: triu => psb_lc_csr_triu
procedure, pass(a) :: cp_to_coo => psb_lc_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_lc_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_lc_cp_csr_to_fmt
@ -529,7 +620,7 @@ module psb_c_csr_mat_mod
procedure, pass(a) :: csput_a => psb_lc_csr_csput_a
procedure, pass(a) :: get_diag => psb_lc_csr_get_diag
procedure, pass(a) :: csgetptn => psb_lc_csr_csgetptn
procedure, pass(a) :: csgetrow => psb_lc_csr_csgetrow
procedure, pass(a) :: csgetrow => psb_lc_csr_csgetrow
procedure, pass(a) :: get_nz_row => lc_csr_get_nz_row
procedure, pass(a) :: reinit => psb_lc_csr_reinit
procedure, pass(a) :: trim => psb_lc_csr_trim
@ -616,6 +707,93 @@ module psb_c_csr_mat_mod
integer(psb_lpk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_lc_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_lc_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_lc_csr_sparse_mat), intent(in) :: a
class(psb_lc_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lc_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_lc_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_lc_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_lc_csr_sparse_mat), intent(in) :: a
class(psb_lc_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lc_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_lc_csr_triu
end interface
!> \memberof psb_lc_csr_sparse_mat
!! \see psb_lc_base_mat_mod::psb_lc_base_cp_to_coo

180
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
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
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
@ -518,6 +607,8 @@ module psb_d_csr_mat_mod
procedure, pass(a) :: sizeof => ld_csr_sizeof
procedure, pass(a) :: reallocate_nz => psb_ld_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_ld_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_ld_csr_tril
procedure, pass(a) :: triu => psb_ld_csr_triu
procedure, pass(a) :: cp_to_coo => psb_ld_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_ld_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_ld_cp_csr_to_fmt
@ -529,7 +620,7 @@ module psb_d_csr_mat_mod
procedure, pass(a) :: csput_a => psb_ld_csr_csput_a
procedure, pass(a) :: get_diag => psb_ld_csr_get_diag
procedure, pass(a) :: csgetptn => psb_ld_csr_csgetptn
procedure, pass(a) :: csgetrow => psb_ld_csr_csgetrow
procedure, pass(a) :: csgetrow => psb_ld_csr_csgetrow
procedure, pass(a) :: get_nz_row => ld_csr_get_nz_row
procedure, pass(a) :: reinit => psb_ld_csr_reinit
procedure, pass(a) :: trim => psb_ld_csr_trim
@ -616,6 +707,93 @@ module psb_d_csr_mat_mod
integer(psb_lpk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_ld_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_ld_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_ld_csr_sparse_mat), intent(in) :: a
class(psb_ld_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ld_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_ld_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_ld_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_ld_csr_sparse_mat), intent(in) :: a
class(psb_ld_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ld_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_ld_csr_triu
end interface
!> \memberof psb_ld_csr_sparse_mat
!! \see psb_ld_base_mat_mod::psb_ld_base_cp_to_coo

180
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
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
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
@ -518,6 +607,8 @@ module psb_s_csr_mat_mod
procedure, pass(a) :: sizeof => ls_csr_sizeof
procedure, pass(a) :: reallocate_nz => psb_ls_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_ls_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_ls_csr_tril
procedure, pass(a) :: triu => psb_ls_csr_triu
procedure, pass(a) :: cp_to_coo => psb_ls_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_ls_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_ls_cp_csr_to_fmt
@ -529,7 +620,7 @@ module psb_s_csr_mat_mod
procedure, pass(a) :: csput_a => psb_ls_csr_csput_a
procedure, pass(a) :: get_diag => psb_ls_csr_get_diag
procedure, pass(a) :: csgetptn => psb_ls_csr_csgetptn
procedure, pass(a) :: csgetrow => psb_ls_csr_csgetrow
procedure, pass(a) :: csgetrow => psb_ls_csr_csgetrow
procedure, pass(a) :: get_nz_row => ls_csr_get_nz_row
procedure, pass(a) :: reinit => psb_ls_csr_reinit
procedure, pass(a) :: trim => psb_ls_csr_trim
@ -616,6 +707,93 @@ module psb_s_csr_mat_mod
integer(psb_lpk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_ls_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_ls_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_ls_csr_sparse_mat), intent(in) :: a
class(psb_ls_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ls_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_ls_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_ls_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_ls_csr_sparse_mat), intent(in) :: a
class(psb_ls_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ls_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_ls_csr_triu
end interface
!> \memberof psb_ls_csr_sparse_mat
!! \see psb_ls_base_mat_mod::psb_ls_base_cp_to_coo

180
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
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
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
@ -518,6 +607,8 @@ module psb_z_csr_mat_mod
procedure, pass(a) :: sizeof => lz_csr_sizeof
procedure, pass(a) :: reallocate_nz => psb_lz_csr_reallocate_nz
procedure, pass(a) :: allocate_mnnz => psb_lz_csr_allocate_mnnz
procedure, pass(a) :: tril => psb_lz_csr_tril
procedure, pass(a) :: triu => psb_lz_csr_triu
procedure, pass(a) :: cp_to_coo => psb_lz_cp_csr_to_coo
procedure, pass(a) :: cp_from_coo => psb_lz_cp_csr_from_coo
procedure, pass(a) :: cp_to_fmt => psb_lz_cp_csr_to_fmt
@ -529,7 +620,7 @@ module psb_z_csr_mat_mod
procedure, pass(a) :: csput_a => psb_lz_csr_csput_a
procedure, pass(a) :: get_diag => psb_lz_csr_get_diag
procedure, pass(a) :: csgetptn => psb_lz_csr_csgetptn
procedure, pass(a) :: csgetrow => psb_lz_csr_csgetrow
procedure, pass(a) :: csgetrow => psb_lz_csr_csgetrow
procedure, pass(a) :: get_nz_row => lz_csr_get_nz_row
procedure, pass(a) :: reinit => psb_lz_csr_reinit
procedure, pass(a) :: trim => psb_lz_csr_trim
@ -616,6 +707,93 @@ module psb_z_csr_mat_mod
integer(psb_lpk_), intent(in), optional :: ivr(:), ivc(:)
end subroutine psb_lz_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_lz_csr_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_lz_csr_sparse_mat), intent(in) :: a
class(psb_lz_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lz_coo_sparse_mat), optional, intent(out) :: u
end subroutine psb_lz_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_lz_csr_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_lz_csr_sparse_mat), intent(in) :: a
class(psb_lz_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lz_coo_sparse_mat), optional, intent(out) :: l
end subroutine psb_lz_csr_triu
end interface
!> \memberof psb_lz_csr_sparse_mat
!! \see psb_lz_base_mat_mod::psb_lz_base_cp_to_coo

44
base/serial/impl/psb_c_base_mat_impl.F90
Unescape Escape View File

@ -632,13 +632,14 @@ subroutine psb_c_base_tril(a,l,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -701,12 +702,12 @@ subroutine psb_c_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -782,13 +783,14 @@ subroutine psb_c_base_triu(a,u,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -851,12 +853,12 @@ subroutine psb_c_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -2935,13 +2937,14 @@ subroutine psb_lc_base_tril(a,l,info,&
class(psb_lc_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3004,12 +3007,12 @@ subroutine psb_lc_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -3086,13 +3089,14 @@ subroutine psb_lc_base_triu(a,u,info,&
class(psb_lc_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
complex(psb_spk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3155,12 +3159,12 @@ subroutine psb_lc_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)

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

@ -2225,6 +2225,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)
@ -4139,6 +4452,321 @@ subroutine psb_lc_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_lc_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_lc_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_lc_csr_tril
implicit none
class(psb_lc_csr_sparse_mat), intent(in) :: a
class(psb_lc_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lc_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_lc_csr_tril
subroutine psb_lc_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_lc_csr_triu
implicit none
class(psb_lc_csr_sparse_mat), intent(in) :: a
class(psb_lc_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lc_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_lc_csr_triu
subroutine psb_lc_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

44
base/serial/impl/psb_d_base_mat_impl.F90
Unescape Escape View File

@ -632,13 +632,14 @@ subroutine psb_d_base_tril(a,l,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -701,12 +702,12 @@ subroutine psb_d_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -782,13 +783,14 @@ subroutine psb_d_base_triu(a,u,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -851,12 +853,12 @@ subroutine psb_d_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -2935,13 +2937,14 @@ subroutine psb_ld_base_tril(a,l,info,&
class(psb_ld_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3004,12 +3007,12 @@ subroutine psb_ld_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -3086,13 +3089,14 @@ subroutine psb_ld_base_triu(a,u,info,&
class(psb_ld_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
real(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3155,12 +3159,12 @@ subroutine psb_ld_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)

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

@ -2225,6 +2225,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)
@ -4139,6 +4452,321 @@ subroutine psb_ld_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_ld_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_ld_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_ld_csr_tril
implicit none
class(psb_ld_csr_sparse_mat), intent(in) :: a
class(psb_ld_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ld_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_ld_csr_tril
subroutine psb_ld_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_ld_csr_triu
implicit none
class(psb_ld_csr_sparse_mat), intent(in) :: a
class(psb_ld_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ld_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_ld_csr_triu
subroutine psb_ld_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

44
base/serial/impl/psb_s_base_mat_impl.F90
Unescape Escape View File

@ -632,13 +632,14 @@ subroutine psb_s_base_tril(a,l,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -701,12 +702,12 @@ subroutine psb_s_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -782,13 +783,14 @@ subroutine psb_s_base_triu(a,u,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -851,12 +853,12 @@ subroutine psb_s_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -2935,13 +2937,14 @@ subroutine psb_ls_base_tril(a,l,info,&
class(psb_ls_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3004,12 +3007,12 @@ subroutine psb_ls_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -3086,13 +3089,14 @@ subroutine psb_ls_base_triu(a,u,info,&
class(psb_ls_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
real(psb_spk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3155,12 +3159,12 @@ subroutine psb_ls_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)

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

@ -2225,6 +2225,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)
@ -4139,6 +4452,321 @@ subroutine psb_ls_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_ls_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_ls_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_ls_csr_tril
implicit none
class(psb_ls_csr_sparse_mat), intent(in) :: a
class(psb_ls_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ls_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_ls_csr_tril
subroutine psb_ls_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_ls_csr_triu
implicit none
class(psb_ls_csr_sparse_mat), intent(in) :: a
class(psb_ls_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_ls_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_ls_csr_triu
subroutine psb_ls_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

44
base/serial/impl/psb_z_base_mat_impl.F90
Unescape Escape View File

@ -632,13 +632,14 @@ subroutine psb_z_base_tril(a,l,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -701,12 +702,12 @@ subroutine psb_z_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -782,13 +783,14 @@ subroutine psb_z_base_triu(a,u,info,&
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_) :: err_act, nzin, nzout, i, j, k, ibk
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_ipk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_ipk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -851,12 +853,12 @@ subroutine psb_z_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -2935,13 +2937,14 @@ subroutine psb_lz_base_tril(a,l,info,&
class(psb_lz_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3004,12 +3007,12 @@ subroutine psb_lz_base_tril(a,l,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<=diag_) then
if ((ja(k)-ia(k))<=diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)
@ -3086,13 +3089,14 @@ subroutine psb_lz_base_triu(a,u,info,&
class(psb_lz_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: nzin, nzout, i, j, k, ibk
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), allocatable :: ia(:), ja(:)
complex(psb_dpk_), allocatable :: val(:)
character(len=20) :: name='triu'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
integer(psb_lpk_), parameter :: nbk=8
call psb_erractionsave(err_act)
info = psb_success_
@ -3155,12 +3159,12 @@ subroutine psb_lz_base_triu(a,u,info,&
call psb_realloc(max(mb,nb),ia,info)
call psb_realloc(max(mb,nb),ja,info)
call psb_realloc(max(mb,nb),val,info)
do i=imin_,imax_
call a%csget(i,i,nzout,ia,ja,val,info,&
do i=imin_,imax_, nbk
ibk = min(nbk,imax_-i+1)
call a%csget(i,i+ibk-1,nzout,ia,ja,val,info,&
& jmin=jmin_, jmax=jmax_)
do k=1, nzout
j = ja(k)
if (j-i<diag_) then
if ((ja(k)-ia(k))<diag_) then
nzlin = nzlin + 1
l%ia(nzlin) = ia(k)
l%ja(nzlin) = ja(k)

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

@ -2225,6 +2225,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)
@ -4139,6 +4452,321 @@ subroutine psb_lz_csr_csgetblk(imin,imax,a,b,info,&
end subroutine psb_lz_csr_csgetblk
!
! CSR implementation of tril/triu
!
subroutine psb_lz_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_lz_csr_tril
implicit none
class(psb_lz_csr_sparse_mat), intent(in) :: a
class(psb_lz_coo_sparse_mat), intent(out) :: l
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lz_coo_sparse_mat), optional, intent(out) :: u
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_lz_csr_tril
subroutine psb_lz_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_lz_csr_triu
implicit none
class(psb_lz_csr_sparse_mat), intent(in) :: a
class(psb_lz_coo_sparse_mat), intent(out) :: u
integer(psb_ipk_),intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
class(psb_lz_coo_sparse_mat), optional, intent(out) :: l
integer(psb_ipk_) :: err_act
integer(psb_lpk_) :: nzin, nzout, i, j, k
integer(psb_lpk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, nzlin, nzuin, nz
integer(psb_lpk_), 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_lz_csr_triu
subroutine psb_lz_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)

Write Preview
Loading…
Cancel
Save