diff --git a/base/modules/serial/psb_c_csr_mat_mod.f90 b/base/modules/serial/psb_c_csr_mat_mod.f90 index 0606c187e..222b46a48 100644 --- a/base/modules/serial/psb_c_csr_mat_mod.f90 +++ b/base/modules/serial/psb_c_csr_mat_mod.f90 @@ -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 diff --git a/base/modules/serial/psb_d_csr_mat_mod.f90 b/base/modules/serial/psb_d_csr_mat_mod.f90 index 7494bda6a..c80246b73 100644 --- a/base/modules/serial/psb_d_csr_mat_mod.f90 +++ b/base/modules/serial/psb_d_csr_mat_mod.f90 @@ -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 diff --git a/base/modules/serial/psb_s_csr_mat_mod.f90 b/base/modules/serial/psb_s_csr_mat_mod.f90 index 8d861d7eb..00dc1dfc3 100644 --- a/base/modules/serial/psb_s_csr_mat_mod.f90 +++ b/base/modules/serial/psb_s_csr_mat_mod.f90 @@ -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 diff --git a/base/modules/serial/psb_z_csr_mat_mod.f90 b/base/modules/serial/psb_z_csr_mat_mod.f90 index df09569a2..28b7d37bd 100644 --- a/base/modules/serial/psb_z_csr_mat_mod.f90 +++ b/base/modules/serial/psb_z_csr_mat_mod.f90 @@ -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 diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90 index 804e70580..55acd5500 100644 --- a/base/serial/impl/psb_c_base_mat_impl.F90 +++ b/base/serial/impl/psb_c_base_mat_impl.F90 @@ -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 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)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)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) diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90 index f359dbe85..21470d9a7 100644 --- a/base/serial/impl/psb_d_base_mat_impl.F90 +++ b/base/serial/impl/psb_d_base_mat_impl.F90 @@ -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 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)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)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) diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90 index a5ed5df05..3f52bdadc 100644 --- a/base/serial/impl/psb_s_base_mat_impl.F90 +++ b/base/serial/impl/psb_s_base_mat_impl.F90 @@ -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 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)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)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) diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90 index 17223e24a..8fdd5ba51 100644 --- a/base/serial/impl/psb_z_base_mat_impl.F90 +++ b/base/serial/impl/psb_z_base_mat_impl.F90 @@ -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 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)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)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)