diff --git a/base/modules/psb_c_base_mat_mod.f90 b/base/modules/psb_c_base_mat_mod.f90 index e5e5fda3..d02e6096 100644 --- a/base/modules/psb_c_base_mat_mod.f90 +++ b/base/modules/psb_c_base_mat_mod.f90 @@ -62,6 +62,8 @@ module psb_c_base_mat_mod procedure, pass(a) :: csgetblk => psb_c_base_csgetblk procedure, pass(a) :: get_diag => psb_c_base_get_diag generic, public :: csget => csgetrow, csgetblk + procedure, pass(a) :: tril => psb_c_base_tril + procedure, pass(a) :: triu => psb_c_base_triu procedure, pass(a) :: csclip => psb_c_base_csclip procedure, pass(a) :: cp_to_coo => psb_c_base_cp_to_coo procedure, pass(a) :: cp_from_coo => psb_c_base_cp_from_coo @@ -312,7 +314,7 @@ module psb_c_base_mat_mod !! \param info return code !! \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 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) @@ -350,7 +352,7 @@ module psb_c_base_mat_mod !! \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 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) @@ -364,11 +366,93 @@ module psb_c_base_mat_mod import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ class(psb_c_base_sparse_mat), intent(in) :: a class(psb_c_coo_sparse_mat), intent(out) :: b - integer(psb_ipk_),intent(out) :: info - integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax - logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale end subroutine psb_c_base_csclip 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 b 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 + !! + ! + interface + subroutine psb_c_base_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ + class(psb_c_base_sparse_mat), intent(in) :: a + class(psb_c_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_c_base_tril + end interface + + ! + !> Function triu: + !! \memberof psb_c_base_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 + !! + !! \param b 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 + !! + ! + interface + subroutine psb_c_base_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ + class(psb_c_base_sparse_mat), intent(in) :: a + class(psb_c_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_c_base_triu + end interface + ! !> Function get_diag: diff --git a/base/modules/psb_c_mat_mod.f90 b/base/modules/psb_c_mat_mod.f90 index 350d098d..0b5452c7 100644 --- a/base/modules/psb_c_mat_mod.f90 +++ b/base/modules/psb_c_mat_mod.f90 @@ -122,6 +122,8 @@ module psb_c_mat_mod procedure, pass(a) :: csgetrow => psb_c_csgetrow procedure, pass(a) :: csgetblk => psb_c_csgetblk generic, public :: csget => csgetptn, csgetrow, csgetblk + procedure, pass(a) :: tril => psb_c_tril + procedure, pass(a) :: triu => psb_c_triu procedure, pass(a) :: m_csclip => psb_c_csclip procedure, pass(a) :: b_csclip => psb_c_b_csclip generic, public :: csclip => b_csclip, m_csclip @@ -420,6 +422,31 @@ module psb_c_mat_mod end subroutine psb_c_csgetblk end interface + interface + subroutine psb_c_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_cspmat_type, psb_spk_ + class(psb_cspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_c_tril + end interface + + interface + subroutine psb_c_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_cspmat_type, psb_spk_ + class(psb_cspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_c_triu + end interface + + interface subroutine psb_c_csclip(a,b,info,& & imin,imax,jmin,jmax,rscale,cscale) diff --git a/base/modules/psb_d_base_mat_mod.f90 b/base/modules/psb_d_base_mat_mod.f90 index 90f85adf..7382eb44 100644 --- a/base/modules/psb_d_base_mat_mod.f90 +++ b/base/modules/psb_d_base_mat_mod.f90 @@ -62,6 +62,8 @@ module psb_d_base_mat_mod procedure, pass(a) :: csgetblk => psb_d_base_csgetblk procedure, pass(a) :: get_diag => psb_d_base_get_diag generic, public :: csget => csgetrow, csgetblk + procedure, pass(a) :: tril => psb_d_base_tril + procedure, pass(a) :: triu => psb_d_base_triu procedure, pass(a) :: csclip => psb_d_base_csclip procedure, pass(a) :: cp_to_coo => psb_d_base_cp_to_coo procedure, pass(a) :: cp_from_coo => psb_d_base_cp_from_coo @@ -312,7 +314,7 @@ module psb_d_base_mat_mod !! \param info return code !! \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 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) @@ -350,7 +352,7 @@ module psb_d_base_mat_mod !! \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 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) @@ -364,11 +366,93 @@ module psb_d_base_mat_mod import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ class(psb_d_base_sparse_mat), intent(in) :: a class(psb_d_coo_sparse_mat), intent(out) :: b - integer(psb_ipk_),intent(out) :: info - integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax - logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale end subroutine psb_d_base_csclip 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 b 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 + !! + ! + interface + subroutine psb_d_base_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ + class(psb_d_base_sparse_mat), intent(in) :: a + class(psb_d_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_d_base_tril + end interface + + ! + !> Function triu: + !! \memberof psb_d_base_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 + !! + !! \param b 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 + !! + ! + interface + subroutine psb_d_base_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ + class(psb_d_base_sparse_mat), intent(in) :: a + class(psb_d_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_d_base_triu + end interface + ! !> Function get_diag: diff --git a/base/modules/psb_d_mat_mod.f90 b/base/modules/psb_d_mat_mod.f90 index dfc43a5e..1b9b8e3c 100644 --- a/base/modules/psb_d_mat_mod.f90 +++ b/base/modules/psb_d_mat_mod.f90 @@ -122,6 +122,8 @@ module psb_d_mat_mod procedure, pass(a) :: csgetrow => psb_d_csgetrow procedure, pass(a) :: csgetblk => psb_d_csgetblk generic, public :: csget => csgetptn, csgetrow, csgetblk + procedure, pass(a) :: tril => psb_d_tril + procedure, pass(a) :: triu => psb_d_triu procedure, pass(a) :: m_csclip => psb_d_csclip procedure, pass(a) :: b_csclip => psb_d_b_csclip generic, public :: csclip => b_csclip, m_csclip @@ -420,6 +422,31 @@ module psb_d_mat_mod end subroutine psb_d_csgetblk end interface + interface + subroutine psb_d_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_dspmat_type, psb_dpk_ + class(psb_dspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_d_tril + end interface + + interface + subroutine psb_d_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_dspmat_type, psb_dpk_ + class(psb_dspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_d_triu + end interface + + interface subroutine psb_d_csclip(a,b,info,& & imin,imax,jmin,jmax,rscale,cscale) diff --git a/base/modules/psb_s_base_mat_mod.f90 b/base/modules/psb_s_base_mat_mod.f90 index 36251b82..dd173201 100644 --- a/base/modules/psb_s_base_mat_mod.f90 +++ b/base/modules/psb_s_base_mat_mod.f90 @@ -62,6 +62,8 @@ module psb_s_base_mat_mod procedure, pass(a) :: csgetblk => psb_s_base_csgetblk procedure, pass(a) :: get_diag => psb_s_base_get_diag generic, public :: csget => csgetrow, csgetblk + procedure, pass(a) :: tril => psb_s_base_tril + procedure, pass(a) :: triu => psb_s_base_triu procedure, pass(a) :: csclip => psb_s_base_csclip procedure, pass(a) :: cp_to_coo => psb_s_base_cp_to_coo procedure, pass(a) :: cp_from_coo => psb_s_base_cp_from_coo @@ -312,7 +314,7 @@ module psb_s_base_mat_mod !! \param info return code !! \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 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) @@ -350,7 +352,7 @@ module psb_s_base_mat_mod !! \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 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) @@ -364,11 +366,93 @@ module psb_s_base_mat_mod import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ class(psb_s_base_sparse_mat), intent(in) :: a class(psb_s_coo_sparse_mat), intent(out) :: b - integer(psb_ipk_),intent(out) :: info - integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax - logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale end subroutine psb_s_base_csclip 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 b 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 + !! + ! + interface + subroutine psb_s_base_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_s_base_tril + end interface + + ! + !> Function triu: + !! \memberof psb_s_base_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 + !! + !! \param b 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 + !! + ! + interface + subroutine psb_s_base_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_s_base_triu + end interface + ! !> Function get_diag: diff --git a/base/modules/psb_s_mat_mod.f90 b/base/modules/psb_s_mat_mod.f90 index c5dfd1ee..79fecc37 100644 --- a/base/modules/psb_s_mat_mod.f90 +++ b/base/modules/psb_s_mat_mod.f90 @@ -122,6 +122,8 @@ module psb_s_mat_mod procedure, pass(a) :: csgetrow => psb_s_csgetrow procedure, pass(a) :: csgetblk => psb_s_csgetblk generic, public :: csget => csgetptn, csgetrow, csgetblk + procedure, pass(a) :: tril => psb_s_tril + procedure, pass(a) :: triu => psb_s_triu procedure, pass(a) :: m_csclip => psb_s_csclip procedure, pass(a) :: b_csclip => psb_s_b_csclip generic, public :: csclip => b_csclip, m_csclip @@ -420,6 +422,31 @@ module psb_s_mat_mod end subroutine psb_s_csgetblk end interface + interface + subroutine psb_s_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_sspmat_type, psb_spk_ + class(psb_sspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_s_tril + end interface + + interface + subroutine psb_s_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_sspmat_type, psb_spk_ + class(psb_sspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_s_triu + end interface + + interface subroutine psb_s_csclip(a,b,info,& & imin,imax,jmin,jmax,rscale,cscale) diff --git a/base/modules/psb_z_base_mat_mod.f90 b/base/modules/psb_z_base_mat_mod.f90 index 362e0031..779f45bd 100644 --- a/base/modules/psb_z_base_mat_mod.f90 +++ b/base/modules/psb_z_base_mat_mod.f90 @@ -62,6 +62,8 @@ module psb_z_base_mat_mod procedure, pass(a) :: csgetblk => psb_z_base_csgetblk procedure, pass(a) :: get_diag => psb_z_base_get_diag generic, public :: csget => csgetrow, csgetblk + procedure, pass(a) :: tril => psb_z_base_tril + procedure, pass(a) :: triu => psb_z_base_triu procedure, pass(a) :: csclip => psb_z_base_csclip procedure, pass(a) :: cp_to_coo => psb_z_base_cp_to_coo procedure, pass(a) :: cp_from_coo => psb_z_base_cp_from_coo @@ -312,7 +314,7 @@ module psb_z_base_mat_mod !! \param info return code !! \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 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) @@ -350,7 +352,7 @@ module psb_z_base_mat_mod !! \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 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) @@ -364,11 +366,93 @@ module psb_z_base_mat_mod import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ class(psb_z_base_sparse_mat), intent(in) :: a class(psb_z_coo_sparse_mat), intent(out) :: b - integer(psb_ipk_),intent(out) :: info - integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax - logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale end subroutine psb_z_base_csclip 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 b 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 + !! + ! + interface + subroutine psb_z_base_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ + class(psb_z_base_sparse_mat), intent(in) :: a + class(psb_z_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_z_base_tril + end interface + + ! + !> Function triu: + !! \memberof psb_z_base_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 + !! + !! \param b 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 + !! + ! + interface + subroutine psb_z_base_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ + class(psb_z_base_sparse_mat), intent(in) :: a + class(psb_z_coo_sparse_mat), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_z_base_triu + end interface + ! !> Function get_diag: diff --git a/base/modules/psb_z_mat_mod.f90 b/base/modules/psb_z_mat_mod.f90 index 32ae261a..28f352f8 100644 --- a/base/modules/psb_z_mat_mod.f90 +++ b/base/modules/psb_z_mat_mod.f90 @@ -122,6 +122,8 @@ module psb_z_mat_mod procedure, pass(a) :: csgetrow => psb_z_csgetrow procedure, pass(a) :: csgetblk => psb_z_csgetblk generic, public :: csget => csgetptn, csgetrow, csgetblk + procedure, pass(a) :: tril => psb_z_tril + procedure, pass(a) :: triu => psb_z_triu procedure, pass(a) :: m_csclip => psb_z_csclip procedure, pass(a) :: b_csclip => psb_z_b_csclip generic, public :: csclip => b_csclip, m_csclip @@ -420,6 +422,31 @@ module psb_z_mat_mod end subroutine psb_z_csgetblk end interface + interface + subroutine psb_z_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_zspmat_type, psb_dpk_ + class(psb_zspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_z_tril + end interface + + interface + subroutine psb_z_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + import :: psb_ipk_, psb_zspmat_type, psb_dpk_ + class(psb_zspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + end subroutine psb_z_triu + end interface + + interface subroutine psb_z_csclip(a,b,info,& & imin,imax,jmin,jmax,rscale,cscale) diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90 index 539ad0c4..aeb62e10 100644 --- a/base/serial/impl/psb_c_base_mat_impl.F90 +++ b/base/serial/impl/psb_c_base_mat_impl.F90 @@ -550,6 +550,232 @@ subroutine psb_c_base_csclip(a,b,info,& end subroutine psb_c_base_csclip + +! +! Here we have the base implementation of tril and triu +! this is just based on the getrow. +! If performance is critical it can be overridden. +! +subroutine psb_c_base_tril(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_c_base_mat_mod, psb_protect_name => psb_c_base_tril + implicit none + + class(psb_c_base_sparse_mat), intent(in) :: a + class(psb_c_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_lower(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_c_base_tril + +subroutine psb_c_base_triu(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_c_base_mat_mod, psb_protect_name => psb_c_base_triu + implicit none + + class(psb_c_base_sparse_mat), intent(in) :: a + class(psb_c_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = max(jmin_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=k, jmax=jmax_, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_upper(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_c_base_triu + + + subroutine psb_c_base_clone(a,b,info) use psb_c_base_mat_mod, psb_protect_name => psb_c_base_clone use psb_error_mod diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index aa4231b1..2d927bfa 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -942,6 +942,111 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,& end subroutine psb_c_csgetblk +subroutine psb_c_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_c_base_mat_mod + use psb_c_mat_mod, psb_protect_name => psb_c_tril + implicit none + class(psb_cspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='tril' + logical, parameter :: debug=.false. + type(psb_c_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_c_tril + +subroutine psb_c_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_c_base_mat_mod + use psb_c_mat_mod, psb_protect_name => psb_c_triu + implicit none + class(psb_cspmat_type), intent(in) :: a + class(psb_cspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='triu' + logical, parameter :: debug=.false. + type(psb_c_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_c_triu subroutine psb_c_csclip(a,b,info,& diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90 index 70343e7d..9f8507c3 100644 --- a/base/serial/impl/psb_d_base_mat_impl.F90 +++ b/base/serial/impl/psb_d_base_mat_impl.F90 @@ -550,6 +550,232 @@ subroutine psb_d_base_csclip(a,b,info,& end subroutine psb_d_base_csclip + +! +! Here we have the base implementation of tril and triu +! this is just based on the getrow. +! If performance is critical it can be overridden. +! +subroutine psb_d_base_tril(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod, psb_protect_name => psb_d_base_tril + implicit none + + class(psb_d_base_sparse_mat), intent(in) :: a + class(psb_d_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_lower(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_d_base_tril + +subroutine psb_d_base_triu(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod, psb_protect_name => psb_d_base_triu + implicit none + + class(psb_d_base_sparse_mat), intent(in) :: a + class(psb_d_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = max(jmin_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=k, jmax=jmax_, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_upper(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_d_base_triu + + + subroutine psb_d_base_clone(a,b,info) use psb_d_base_mat_mod, psb_protect_name => psb_d_base_clone use psb_error_mod diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index e39ffbc1..848f6fa9 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -942,6 +942,111 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,& end subroutine psb_d_csgetblk +subroutine psb_d_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_d_mat_mod, psb_protect_name => psb_d_tril + implicit none + class(psb_dspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='tril' + logical, parameter :: debug=.false. + type(psb_d_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_d_tril + +subroutine psb_d_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_d_mat_mod, psb_protect_name => psb_d_triu + implicit none + class(psb_dspmat_type), intent(in) :: a + class(psb_dspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='triu' + logical, parameter :: debug=.false. + type(psb_d_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_d_triu subroutine psb_d_csclip(a,b,info,& diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90 index 48d015dc..32193ab3 100644 --- a/base/serial/impl/psb_s_base_mat_impl.F90 +++ b/base/serial/impl/psb_s_base_mat_impl.F90 @@ -550,6 +550,232 @@ subroutine psb_s_base_csclip(a,b,info,& end subroutine psb_s_base_csclip + +! +! Here we have the base implementation of tril and triu +! this is just based on the getrow. +! If performance is critical it can be overridden. +! +subroutine psb_s_base_tril(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod, psb_protect_name => psb_s_base_tril + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_lower(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_s_base_tril + +subroutine psb_s_base_triu(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod, psb_protect_name => psb_s_base_triu + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = max(jmin_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=k, jmax=jmax_, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_upper(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_s_base_triu + + + subroutine psb_s_base_clone(a,b,info) use psb_s_base_mat_mod, psb_protect_name => psb_s_base_clone use psb_error_mod diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index e8f6a541..3f8e72ac 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -942,6 +942,111 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,& end subroutine psb_s_csgetblk +subroutine psb_s_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod + use psb_s_mat_mod, psb_protect_name => psb_s_tril + implicit none + class(psb_sspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='tril' + logical, parameter :: debug=.false. + type(psb_s_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_s_tril + +subroutine psb_s_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod + use psb_s_mat_mod, psb_protect_name => psb_s_triu + implicit none + class(psb_sspmat_type), intent(in) :: a + class(psb_sspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='triu' + logical, parameter :: debug=.false. + type(psb_s_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_s_triu subroutine psb_s_csclip(a,b,info,& diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90 index fe0084fb..e5607f50 100644 --- a/base/serial/impl/psb_z_base_mat_impl.F90 +++ b/base/serial/impl/psb_z_base_mat_impl.F90 @@ -550,6 +550,232 @@ subroutine psb_z_base_csclip(a,b,info,& end subroutine psb_z_base_csclip + +! +! Here we have the base implementation of tril and triu +! this is just based on the getrow. +! If performance is critical it can be overridden. +! +subroutine psb_z_base_tril(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_z_base_mat_mod, psb_protect_name => psb_z_base_tril + implicit none + + class(psb_z_base_sparse_mat), intent(in) :: a + class(psb_z_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ <= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_lower(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_z_base_tril + +subroutine psb_z_base_triu(a,b,info,& + & diag,imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_z_base_mat_mod, psb_protect_name => psb_z_base_triu + implicit none + + class(psb_z_base_sparse_mat), intent(in) :: a + class(psb_z_coo_sparse_mat), intent(out) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k + integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_ + 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 + + call b%allocate(mb,nb) + nzin = b%get_nzeros() ! At this point it should be 0 + + do i=imin_,imax_ + k = max(jmin_,i+diag_) + call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& + & jmin=k, jmax=jmax_, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call b%set_nzeros(nzin+nzout) + end do + call b%fix(info) + nzout = b%get_nzeros() + if (rscale_) & + & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + if (cscale_) & + & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + + if ((diag_ >= 0).and.(imin_ == jmin_)) then + call b%set_triangle(.true.) + call b%set_upper(.true.) + end if + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + +end subroutine psb_z_base_triu + + + subroutine psb_z_base_clone(a,b,info) use psb_z_base_mat_mod, psb_protect_name => psb_z_base_clone use psb_error_mod diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index f313f6fd..cebd933d 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -942,6 +942,111 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,& end subroutine psb_z_csgetblk +subroutine psb_z_tril(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_z_base_mat_mod + use psb_z_mat_mod, psb_protect_name => psb_z_tril + implicit none + class(psb_zspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='tril' + logical, parameter :: debug=.false. + type(psb_z_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_z_tril + +subroutine psb_z_triu(a,b,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + use psb_error_mod + use psb_const_mod + use psb_z_base_mat_mod + use psb_z_mat_mod, psb_protect_name => psb_z_triu + implicit none + class(psb_zspmat_type), intent(in) :: a + class(psb_zspmat_type), intent(inout) :: b + integer(psb_ipk_),intent(out) :: info + integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + integer(psb_ipk_) :: err_act + character(len=20) :: name='triu' + logical, parameter :: debug=.false. + type(psb_z_coo_sparse_mat), allocatable :: acoo + + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_null()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == psb_success_) then + call a%a%tril(acoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if + if (info == psb_success_) call move_alloc(acoo,b%a) + if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + + +end subroutine psb_z_triu subroutine psb_z_csclip(a,b,info,&