From 7218ae86d6fb017fc5e1ed3bf8ecfde5da44d42f Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Wed, 9 Aug 2017 10:43:38 +0100 Subject: [PATCH] Added new optional arguments to TRIL and TRIU. Document same. --- base/modules/serial/psb_c_base_mat_mod.f90 | 21 ++- base/modules/serial/psb_c_mat_mod.f90 | 14 +- base/modules/serial/psb_d_base_mat_mod.f90 | 21 ++- base/modules/serial/psb_d_mat_mod.f90 | 14 +- base/modules/serial/psb_s_base_mat_mod.f90 | 21 ++- base/modules/serial/psb_s_mat_mod.f90 | 14 +- base/modules/serial/psb_z_base_mat_mod.f90 | 21 ++- base/modules/serial/psb_z_mat_mod.f90 | 14 +- base/serial/impl/psb_c_base_mat_impl.F90 | 179 ++++++++++++++++----- base/serial/impl/psb_c_mat_impl.F90 | 72 +++++---- base/serial/impl/psb_d_base_mat_impl.F90 | 179 ++++++++++++++++----- base/serial/impl/psb_d_mat_impl.F90 | 72 +++++---- base/serial/impl/psb_s_base_mat_impl.F90 | 179 ++++++++++++++++----- base/serial/impl/psb_s_mat_impl.F90 | 72 +++++---- base/serial/impl/psb_z_base_mat_impl.F90 | 179 ++++++++++++++++----- base/serial/impl/psb_z_mat_impl.F90 | 72 +++++---- docs/src/datastruct.tex | 22 ++- docs/src/userguide.tex | 2 +- docs/src/userhtml.tex | 2 +- 19 files changed, 816 insertions(+), 354 deletions(-) diff --git a/base/modules/serial/psb_c_base_mat_mod.f90 b/base/modules/serial/psb_c_base_mat_mod.f90 index a83ad57b..6990a210 100644 --- a/base/modules/serial/psb_c_base_mat_mod.f90 +++ b/base/modules/serial/psb_c_base_mat_mod.f90 @@ -416,7 +416,7 @@ module psb_c_base_mat_mod !! IMIN<=I<=IMAX !! JMIN<=J<=JMAX !! - !! \param b the output (sub)matrix + !! \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 @@ -429,17 +429,19 @@ module psb_c_base_mat_mod !! ( 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_base_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_c_base_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ class(psb_c_base_sparse_mat), intent(in) :: a - class(psb_c_coo_sparse_mat), intent(out) :: b + 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_base_tril end interface @@ -456,8 +458,9 @@ module psb_c_base_mat_mod !! 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 b the output (sub)matrix + !! \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 @@ -470,17 +473,19 @@ module psb_c_base_mat_mod !! ( 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_base_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_c_base_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ class(psb_c_base_sparse_mat), intent(in) :: a - class(psb_c_coo_sparse_mat), intent(out) :: b + 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_base_triu end interface diff --git a/base/modules/serial/psb_c_mat_mod.f90 b/base/modules/serial/psb_c_mat_mod.f90 index dc54771a..70fa36d7 100644 --- a/base/modules/serial/psb_c_mat_mod.f90 +++ b/base/modules/serial/psb_c_mat_mod.f90 @@ -481,26 +481,28 @@ module psb_c_mat_mod end interface interface - subroutine psb_c_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_c_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_cspmat_type, psb_spk_ class(psb_cspmat_type), intent(in) :: a - class(psb_cspmat_type), intent(inout) :: b + class(psb_cspmat_type), intent(inout) :: 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_cspmat_type), optional, intent(inout) :: u end subroutine psb_c_tril end interface interface - subroutine psb_c_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_c_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_cspmat_type, psb_spk_ class(psb_cspmat_type), intent(in) :: a - class(psb_cspmat_type), intent(inout) :: b + class(psb_cspmat_type), intent(inout) :: 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_cspmat_type), optional, intent(inout) :: l end subroutine psb_c_triu end interface diff --git a/base/modules/serial/psb_d_base_mat_mod.f90 b/base/modules/serial/psb_d_base_mat_mod.f90 index 3f1ba8d5..ac88a9d0 100644 --- a/base/modules/serial/psb_d_base_mat_mod.f90 +++ b/base/modules/serial/psb_d_base_mat_mod.f90 @@ -416,7 +416,7 @@ module psb_d_base_mat_mod !! IMIN<=I<=IMAX !! JMIN<=J<=JMAX !! - !! \param b the output (sub)matrix + !! \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 @@ -429,17 +429,19 @@ module psb_d_base_mat_mod !! ( 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_base_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_d_base_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ class(psb_d_base_sparse_mat), intent(in) :: a - class(psb_d_coo_sparse_mat), intent(out) :: b + 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_base_tril end interface @@ -456,8 +458,9 @@ module psb_d_base_mat_mod !! 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 b the output (sub)matrix + !! \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 @@ -470,17 +473,19 @@ module psb_d_base_mat_mod !! ( 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_base_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_d_base_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ class(psb_d_base_sparse_mat), intent(in) :: a - class(psb_d_coo_sparse_mat), intent(out) :: b + 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_base_triu end interface diff --git a/base/modules/serial/psb_d_mat_mod.f90 b/base/modules/serial/psb_d_mat_mod.f90 index f70ae449..2c065f70 100644 --- a/base/modules/serial/psb_d_mat_mod.f90 +++ b/base/modules/serial/psb_d_mat_mod.f90 @@ -481,26 +481,28 @@ module psb_d_mat_mod end interface interface - subroutine psb_d_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_d_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_dspmat_type, psb_dpk_ class(psb_dspmat_type), intent(in) :: a - class(psb_dspmat_type), intent(inout) :: b + class(psb_dspmat_type), intent(inout) :: 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_dspmat_type), optional, intent(inout) :: u end subroutine psb_d_tril end interface interface - subroutine psb_d_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_d_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_dspmat_type, psb_dpk_ class(psb_dspmat_type), intent(in) :: a - class(psb_dspmat_type), intent(inout) :: b + class(psb_dspmat_type), intent(inout) :: 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_dspmat_type), optional, intent(inout) :: l end subroutine psb_d_triu end interface diff --git a/base/modules/serial/psb_s_base_mat_mod.f90 b/base/modules/serial/psb_s_base_mat_mod.f90 index df05743b..5b39b607 100644 --- a/base/modules/serial/psb_s_base_mat_mod.f90 +++ b/base/modules/serial/psb_s_base_mat_mod.f90 @@ -416,7 +416,7 @@ module psb_s_base_mat_mod !! IMIN<=I<=IMAX !! JMIN<=J<=JMAX !! - !! \param b the output (sub)matrix + !! \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 @@ -429,17 +429,19 @@ module psb_s_base_mat_mod !! ( 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_base_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_s_base_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ class(psb_s_base_sparse_mat), intent(in) :: a - class(psb_s_coo_sparse_mat), intent(out) :: b + 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_base_tril end interface @@ -456,8 +458,9 @@ module psb_s_base_mat_mod !! 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 b the output (sub)matrix + !! \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 @@ -470,17 +473,19 @@ module psb_s_base_mat_mod !! ( 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_base_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_s_base_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ class(psb_s_base_sparse_mat), intent(in) :: a - class(psb_s_coo_sparse_mat), intent(out) :: b + 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_base_triu end interface diff --git a/base/modules/serial/psb_s_mat_mod.f90 b/base/modules/serial/psb_s_mat_mod.f90 index 710637d3..68cd902c 100644 --- a/base/modules/serial/psb_s_mat_mod.f90 +++ b/base/modules/serial/psb_s_mat_mod.f90 @@ -481,26 +481,28 @@ module psb_s_mat_mod end interface interface - subroutine psb_s_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_s_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_sspmat_type, psb_spk_ class(psb_sspmat_type), intent(in) :: a - class(psb_sspmat_type), intent(inout) :: b + class(psb_sspmat_type), intent(inout) :: 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_sspmat_type), optional, intent(inout) :: u end subroutine psb_s_tril end interface interface - subroutine psb_s_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_s_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_sspmat_type, psb_spk_ class(psb_sspmat_type), intent(in) :: a - class(psb_sspmat_type), intent(inout) :: b + class(psb_sspmat_type), intent(inout) :: 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_sspmat_type), optional, intent(inout) :: l end subroutine psb_s_triu end interface diff --git a/base/modules/serial/psb_z_base_mat_mod.f90 b/base/modules/serial/psb_z_base_mat_mod.f90 index 5fc7ea4b..7359aea2 100644 --- a/base/modules/serial/psb_z_base_mat_mod.f90 +++ b/base/modules/serial/psb_z_base_mat_mod.f90 @@ -416,7 +416,7 @@ module psb_z_base_mat_mod !! IMIN<=I<=IMAX !! JMIN<=J<=JMAX !! - !! \param b the output (sub)matrix + !! \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 @@ -429,17 +429,19 @@ module psb_z_base_mat_mod !! ( 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_base_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_z_base_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ class(psb_z_base_sparse_mat), intent(in) :: a - class(psb_z_coo_sparse_mat), intent(out) :: b + 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_base_tril end interface @@ -456,8 +458,9 @@ module psb_z_base_mat_mod !! 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 b the output (sub)matrix + !! \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 @@ -470,17 +473,19 @@ module psb_z_base_mat_mod !! ( 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_base_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_z_base_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ class(psb_z_base_sparse_mat), intent(in) :: a - class(psb_z_coo_sparse_mat), intent(out) :: b + 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_base_triu end interface diff --git a/base/modules/serial/psb_z_mat_mod.f90 b/base/modules/serial/psb_z_mat_mod.f90 index 02188678..0cfca6b3 100644 --- a/base/modules/serial/psb_z_mat_mod.f90 +++ b/base/modules/serial/psb_z_mat_mod.f90 @@ -481,26 +481,28 @@ module psb_z_mat_mod end interface interface - subroutine psb_z_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_z_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) import :: psb_ipk_, psb_zspmat_type, psb_dpk_ class(psb_zspmat_type), intent(in) :: a - class(psb_zspmat_type), intent(inout) :: b + class(psb_zspmat_type), intent(inout) :: 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_zspmat_type), optional, intent(inout) :: u end subroutine psb_z_tril end interface interface - subroutine psb_z_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + subroutine psb_z_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) import :: psb_ipk_, psb_zspmat_type, psb_dpk_ class(psb_zspmat_type), intent(in) :: a - class(psb_zspmat_type), intent(inout) :: b + class(psb_zspmat_type), intent(inout) :: 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_zspmat_type), optional, intent(inout) :: l end subroutine psb_z_triu end interface diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90 index e3e8ea4d..bb422cd7 100644 --- a/base/serial/impl/psb_c_base_mat_impl.F90 +++ b/base/serial/impl/psb_c_base_mat_impl.F90 @@ -630,8 +630,8 @@ end subroutine psb_c_base_csclip ! this is just based on the getrow. ! If performance is critical it can be overridden. ! -subroutine psb_c_base_tril(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_c_base_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 @@ -639,12 +639,16 @@ subroutine psb_c_base_tril(a,b,info,& implicit none class(psb_c_base_sparse_mat), intent(in) :: a - class(psb_c_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -700,28 +704,70 @@ subroutine psb_c_base_tril(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = min(jmax_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin_, jmax=k, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i<=diag_) then + nzlin = nzlin + 1 + l%ia(nzlin) = ia(k) + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + else + nzuin = nzuin + 1 + u%ia(nzuin) = ia(k) + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + end if + end do + end do + + 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 + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call l%set_nzeros(nzin+nzout) + nzin = nzin+nzout + end do + end if + call l%fix(info) + nzout = l%get_nzeros() if (rscale_) & - & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 if (cscale_) & - & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 if ((diag_ <= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_lower(.true.) + call l%set_triangle(.true.) + call l%set_lower(.true.) end if if (info /= psb_success_) goto 9999 @@ -735,8 +781,8 @@ subroutine psb_c_base_tril(a,b,info,& end subroutine psb_c_base_tril -subroutine psb_c_base_triu(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_c_base_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 @@ -744,12 +790,16 @@ subroutine psb_c_base_triu(a,b,info,& implicit none class(psb_c_base_sparse_mat), intent(in) :: a - class(psb_c_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -805,28 +855,69 @@ subroutine psb_c_base_triu(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = max(jmin_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=k, jmax=jmax_, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_upper(.true.) + call u%set_triangle(.true.) + call u%set_upper(.true.) end if if (info /= psb_success_) goto 9999 diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index 49406e19..628c52c3 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -906,24 +906,24 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,& end subroutine psb_c_csgetblk -subroutine psb_c_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_c_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) use psb_error_mod use psb_const_mod use psb_c_base_mat_mod use psb_c_mat_mod, psb_protect_name => psb_c_tril implicit none class(psb_cspmat_type), intent(in) :: a - class(psb_cspmat_type), intent(inout) :: b + class(psb_cspmat_type), intent(inout) :: 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_cspmat_type), optional, intent(inout) :: u integer(psb_ipk_) :: err_act character(len=20) :: name='tril' logical, parameter :: debug=.false. - type(psb_c_coo_sparse_mat), allocatable :: acoo - + type(psb_c_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ call psb_erractionsave(err_act) @@ -932,23 +932,30 @@ subroutine psb_c_tril(a,b,info,diag,imin,imax,& call psb_errpush(info,name) goto 9999 endif - allocate(acoo,stat=info) - call b%free() - - if (info == psb_success_) then - call a%a%tril(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + allocate(lcoo,stat=info) + call l%free() + if (present(u)) then + if (info == psb_success_) allocate(ucoo,stat=info) + call u%free() + if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,ucoo) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return @@ -956,23 +963,24 @@ subroutine psb_c_tril(a,b,info,diag,imin,imax,& end subroutine psb_c_tril -subroutine psb_c_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_c_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) use psb_error_mod use psb_const_mod use psb_c_base_mat_mod use psb_c_mat_mod, psb_protect_name => psb_c_triu implicit none class(psb_cspmat_type), intent(in) :: a - class(psb_cspmat_type), intent(inout) :: b + class(psb_cspmat_type), intent(inout) :: 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_cspmat_type), optional, intent(inout) :: l integer(psb_ipk_) :: err_act character(len=20) :: name='triu' logical, parameter :: debug=.false. - type(psb_c_coo_sparse_mat), allocatable :: acoo + type(psb_c_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ @@ -983,23 +991,31 @@ subroutine psb_c_triu(a,b,info,diag,imin,imax,& goto 9999 endif - allocate(acoo,stat=info) - call b%free() + allocate(ucoo,stat=info) + call u%free() - if (info == psb_success_) then - call a%a%triu(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + if (present(l)) then + if (info == psb_success_) allocate(lcoo,stat=info) + call l%free() + if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,lcoo) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90 index 0bfc4c81..5ab73095 100644 --- a/base/serial/impl/psb_d_base_mat_impl.F90 +++ b/base/serial/impl/psb_d_base_mat_impl.F90 @@ -630,8 +630,8 @@ end subroutine psb_d_base_csclip ! this is just based on the getrow. ! If performance is critical it can be overridden. ! -subroutine psb_d_base_tril(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_d_base_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 @@ -639,12 +639,16 @@ subroutine psb_d_base_tril(a,b,info,& implicit none class(psb_d_base_sparse_mat), intent(in) :: a - class(psb_d_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -700,28 +704,70 @@ subroutine psb_d_base_tril(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = min(jmax_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin_, jmax=k, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i<=diag_) then + nzlin = nzlin + 1 + l%ia(nzlin) = ia(k) + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + else + nzuin = nzuin + 1 + u%ia(nzuin) = ia(k) + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + end if + end do + end do + + 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 + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call l%set_nzeros(nzin+nzout) + nzin = nzin+nzout + end do + end if + call l%fix(info) + nzout = l%get_nzeros() if (rscale_) & - & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 if (cscale_) & - & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 if ((diag_ <= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_lower(.true.) + call l%set_triangle(.true.) + call l%set_lower(.true.) end if if (info /= psb_success_) goto 9999 @@ -735,8 +781,8 @@ subroutine psb_d_base_tril(a,b,info,& end subroutine psb_d_base_tril -subroutine psb_d_base_triu(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_d_base_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 @@ -744,12 +790,16 @@ subroutine psb_d_base_triu(a,b,info,& implicit none class(psb_d_base_sparse_mat), intent(in) :: a - class(psb_d_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -805,28 +855,69 @@ subroutine psb_d_base_triu(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = max(jmin_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=k, jmax=jmax_, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_upper(.true.) + call u%set_triangle(.true.) + call u%set_upper(.true.) end if if (info /= psb_success_) goto 9999 diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index 72a05a84..6e36a609 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -906,24 +906,24 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,& end subroutine psb_d_csgetblk -subroutine psb_d_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_d_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) use psb_error_mod use psb_const_mod use psb_d_base_mat_mod use psb_d_mat_mod, psb_protect_name => psb_d_tril implicit none class(psb_dspmat_type), intent(in) :: a - class(psb_dspmat_type), intent(inout) :: b + class(psb_dspmat_type), intent(inout) :: 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_dspmat_type), optional, intent(inout) :: u integer(psb_ipk_) :: err_act character(len=20) :: name='tril' logical, parameter :: debug=.false. - type(psb_d_coo_sparse_mat), allocatable :: acoo - + type(psb_d_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ call psb_erractionsave(err_act) @@ -932,23 +932,30 @@ subroutine psb_d_tril(a,b,info,diag,imin,imax,& call psb_errpush(info,name) goto 9999 endif - allocate(acoo,stat=info) - call b%free() - - if (info == psb_success_) then - call a%a%tril(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + allocate(lcoo,stat=info) + call l%free() + if (present(u)) then + if (info == psb_success_) allocate(ucoo,stat=info) + call u%free() + if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,ucoo) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return @@ -956,23 +963,24 @@ subroutine psb_d_tril(a,b,info,diag,imin,imax,& end subroutine psb_d_tril -subroutine psb_d_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_d_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) use psb_error_mod use psb_const_mod use psb_d_base_mat_mod use psb_d_mat_mod, psb_protect_name => psb_d_triu implicit none class(psb_dspmat_type), intent(in) :: a - class(psb_dspmat_type), intent(inout) :: b + class(psb_dspmat_type), intent(inout) :: 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_dspmat_type), optional, intent(inout) :: l integer(psb_ipk_) :: err_act character(len=20) :: name='triu' logical, parameter :: debug=.false. - type(psb_d_coo_sparse_mat), allocatable :: acoo + type(psb_d_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ @@ -983,23 +991,31 @@ subroutine psb_d_triu(a,b,info,diag,imin,imax,& goto 9999 endif - allocate(acoo,stat=info) - call b%free() + allocate(ucoo,stat=info) + call u%free() - if (info == psb_success_) then - call a%a%triu(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + if (present(l)) then + if (info == psb_success_) allocate(lcoo,stat=info) + call l%free() + if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,lcoo) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return diff --git a/base/serial/impl/psb_s_base_mat_impl.F90 b/base/serial/impl/psb_s_base_mat_impl.F90 index 259c6b7e..39c02e8f 100644 --- a/base/serial/impl/psb_s_base_mat_impl.F90 +++ b/base/serial/impl/psb_s_base_mat_impl.F90 @@ -630,8 +630,8 @@ end subroutine psb_s_base_csclip ! this is just based on the getrow. ! If performance is critical it can be overridden. ! -subroutine psb_s_base_tril(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_s_base_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 @@ -639,12 +639,16 @@ subroutine psb_s_base_tril(a,b,info,& implicit none class(psb_s_base_sparse_mat), intent(in) :: a - class(psb_s_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -700,28 +704,70 @@ subroutine psb_s_base_tril(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = min(jmax_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin_, jmax=k, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i<=diag_) then + nzlin = nzlin + 1 + l%ia(nzlin) = ia(k) + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + else + nzuin = nzuin + 1 + u%ia(nzuin) = ia(k) + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + end if + end do + end do + + 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 + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call l%set_nzeros(nzin+nzout) + nzin = nzin+nzout + end do + end if + call l%fix(info) + nzout = l%get_nzeros() if (rscale_) & - & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 if (cscale_) & - & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 if ((diag_ <= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_lower(.true.) + call l%set_triangle(.true.) + call l%set_lower(.true.) end if if (info /= psb_success_) goto 9999 @@ -735,8 +781,8 @@ subroutine psb_s_base_tril(a,b,info,& end subroutine psb_s_base_tril -subroutine psb_s_base_triu(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_s_base_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 @@ -744,12 +790,16 @@ subroutine psb_s_base_triu(a,b,info,& implicit none class(psb_s_base_sparse_mat), intent(in) :: a - class(psb_s_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -805,28 +855,69 @@ subroutine psb_s_base_triu(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = max(jmin_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=k, jmax=jmax_, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_upper(.true.) + call u%set_triangle(.true.) + call u%set_upper(.true.) end if if (info /= psb_success_) goto 9999 diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index 091dec6e..ed69c0f4 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -906,24 +906,24 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,& end subroutine psb_s_csgetblk -subroutine psb_s_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_s_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) use psb_error_mod use psb_const_mod use psb_s_base_mat_mod use psb_s_mat_mod, psb_protect_name => psb_s_tril implicit none class(psb_sspmat_type), intent(in) :: a - class(psb_sspmat_type), intent(inout) :: b + class(psb_sspmat_type), intent(inout) :: 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_sspmat_type), optional, intent(inout) :: u integer(psb_ipk_) :: err_act character(len=20) :: name='tril' logical, parameter :: debug=.false. - type(psb_s_coo_sparse_mat), allocatable :: acoo - + type(psb_s_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ call psb_erractionsave(err_act) @@ -932,23 +932,30 @@ subroutine psb_s_tril(a,b,info,diag,imin,imax,& call psb_errpush(info,name) goto 9999 endif - allocate(acoo,stat=info) - call b%free() - - if (info == psb_success_) then - call a%a%tril(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + allocate(lcoo,stat=info) + call l%free() + if (present(u)) then + if (info == psb_success_) allocate(ucoo,stat=info) + call u%free() + if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,ucoo) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return @@ -956,23 +963,24 @@ subroutine psb_s_tril(a,b,info,diag,imin,imax,& end subroutine psb_s_tril -subroutine psb_s_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_s_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) use psb_error_mod use psb_const_mod use psb_s_base_mat_mod use psb_s_mat_mod, psb_protect_name => psb_s_triu implicit none class(psb_sspmat_type), intent(in) :: a - class(psb_sspmat_type), intent(inout) :: b + class(psb_sspmat_type), intent(inout) :: 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_sspmat_type), optional, intent(inout) :: l integer(psb_ipk_) :: err_act character(len=20) :: name='triu' logical, parameter :: debug=.false. - type(psb_s_coo_sparse_mat), allocatable :: acoo + type(psb_s_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ @@ -983,23 +991,31 @@ subroutine psb_s_triu(a,b,info,diag,imin,imax,& goto 9999 endif - allocate(acoo,stat=info) - call b%free() + allocate(ucoo,stat=info) + call u%free() - if (info == psb_success_) then - call a%a%triu(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + if (present(l)) then + if (info == psb_success_) allocate(lcoo,stat=info) + call l%free() + if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,lcoo) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return diff --git a/base/serial/impl/psb_z_base_mat_impl.F90 b/base/serial/impl/psb_z_base_mat_impl.F90 index b8b88cb9..5f984f6c 100644 --- a/base/serial/impl/psb_z_base_mat_impl.F90 +++ b/base/serial/impl/psb_z_base_mat_impl.F90 @@ -630,8 +630,8 @@ end subroutine psb_z_base_csclip ! this is just based on the getrow. ! If performance is critical it can be overridden. ! -subroutine psb_z_base_tril(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_z_base_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 @@ -639,12 +639,16 @@ subroutine psb_z_base_tril(a,b,info,& implicit none class(psb_z_base_sparse_mat), intent(in) :: a - class(psb_z_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -700,28 +704,70 @@ subroutine psb_z_base_tril(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = min(jmax_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=jmin_, jmax=k, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i<=diag_) then + nzlin = nzlin + 1 + l%ia(nzlin) = ia(k) + l%ja(nzlin) = ja(k) + l%val(nzlin) = val(k) + else + nzuin = nzuin + 1 + u%ia(nzuin) = ia(k) + u%ja(nzuin) = ja(k) + u%val(nzuin) = val(k) + end if + end do + end do + + 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 + do i=imin_,imax_ + k = min(jmax_,i+diag_) + call a%csget(i,i,nzout,l%ia,l%ja,l%val,info,& + & jmin=jmin_, jmax=k, append=.true., & + & nzin=nzin) + if (info /= psb_success_) goto 9999 + call l%set_nzeros(nzin+nzout) + nzin = nzin+nzout + end do + end if + call l%fix(info) + nzout = l%get_nzeros() if (rscale_) & - & b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1 + & l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1 if (cscale_) & - & b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1 + & l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1 if ((diag_ <= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_lower(.true.) + call l%set_triangle(.true.) + call l%set_lower(.true.) end if if (info /= psb_success_) goto 9999 @@ -735,8 +781,8 @@ subroutine psb_z_base_tril(a,b,info,& end subroutine psb_z_base_tril -subroutine psb_z_base_triu(a,b,info,& - & diag,imin,imax,jmin,jmax,rscale,cscale) +subroutine psb_z_base_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 @@ -744,12 +790,16 @@ subroutine psb_z_base_triu(a,b,info,& implicit none class(psb_z_base_sparse_mat), intent(in) :: a - class(psb_z_coo_sparse_mat), intent(out) :: b + 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_ + 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_ @@ -805,28 +855,69 @@ subroutine psb_z_base_triu(a,b,info,& nb = jmax_ endif - call b%allocate(mb,nb) - nzin = b%get_nzeros() ! At this point it should be 0 - - do i=imin_,imax_ - k = max(jmin_,i+diag_) - call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,& - & jmin=k, jmax=jmax_, append=.true., & - & nzin=nzin) - if (info /= psb_success_) goto 9999 - call b%set_nzeros(nzin+nzout) - nzin = nzin+nzout - end do - call b%fix(info) - nzout = b%get_nzeros() + + 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 + 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,& + & jmin=jmin_, jmax=jmax_) + do k=1, nzout + j = ja(k) + if (j-i= 0).and.(imin_ == jmin_)) then - call b%set_triangle(.true.) - call b%set_upper(.true.) + call u%set_triangle(.true.) + call u%set_upper(.true.) end if if (info /= psb_success_) goto 9999 diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index 7ab5f7ba..41dd41fa 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -906,24 +906,24 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,& end subroutine psb_z_csgetblk -subroutine psb_z_tril(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_z_tril(a,l,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,u) use psb_error_mod use psb_const_mod use psb_z_base_mat_mod use psb_z_mat_mod, psb_protect_name => psb_z_tril implicit none class(psb_zspmat_type), intent(in) :: a - class(psb_zspmat_type), intent(inout) :: b + class(psb_zspmat_type), intent(inout) :: 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_zspmat_type), optional, intent(inout) :: u integer(psb_ipk_) :: err_act character(len=20) :: name='tril' logical, parameter :: debug=.false. - type(psb_z_coo_sparse_mat), allocatable :: acoo - + type(psb_z_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ call psb_erractionsave(err_act) @@ -932,23 +932,30 @@ subroutine psb_z_tril(a,b,info,diag,imin,imax,& call psb_errpush(info,name) goto 9999 endif - allocate(acoo,stat=info) - call b%free() - - if (info == psb_success_) then - call a%a%tril(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + allocate(lcoo,stat=info) + call l%free() + if (present(u)) then + if (info == psb_success_) allocate(ucoo,stat=info) + call u%free() + if (info == psb_success_) call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,ucoo) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%tril(lcoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return @@ -956,23 +963,24 @@ subroutine psb_z_tril(a,b,info,diag,imin,imax,& end subroutine psb_z_tril -subroutine psb_z_triu(a,b,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) +subroutine psb_z_triu(a,u,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,l) use psb_error_mod use psb_const_mod use psb_z_base_mat_mod use psb_z_mat_mod, psb_protect_name => psb_z_triu implicit none class(psb_zspmat_type), intent(in) :: a - class(psb_zspmat_type), intent(inout) :: b + class(psb_zspmat_type), intent(inout) :: 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_zspmat_type), optional, intent(inout) :: l integer(psb_ipk_) :: err_act character(len=20) :: name='triu' logical, parameter :: debug=.false. - type(psb_z_coo_sparse_mat), allocatable :: acoo + type(psb_z_coo_sparse_mat), allocatable :: lcoo, ucoo info = psb_success_ @@ -983,23 +991,31 @@ subroutine psb_z_triu(a,b,info,diag,imin,imax,& goto 9999 endif - allocate(acoo,stat=info) - call b%free() + allocate(ucoo,stat=info) + call u%free() - if (info == psb_success_) then - call a%a%triu(acoo,info,diag,imin,imax,& - & jmin,jmax,rscale,cscale) + if (present(l)) then + if (info == psb_success_) allocate(lcoo,stat=info) + call l%free() + if (info == psb_success_) call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale,lcoo) + if (info == psb_success_) call move_alloc(lcoo,l%a) + if (info == psb_success_) call l%cscnv(info,mold=a%a) else - info = psb_err_alloc_dealloc_ + if (info == psb_success_) then + call a%a%triu(ucoo,info,diag,imin,imax,& + & jmin,jmax,rscale,cscale) + else + info = psb_err_alloc_dealloc_ + end if end if - if (info == psb_success_) call move_alloc(acoo,b%a) - if (info == psb_success_) call b%cscnv(info,mold=a%a) + if (info == psb_success_) call move_alloc(ucoo,u%a) + if (info == psb_success_) call u%cscnv(info,mold=a%a) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) return - 9999 call psb_error_handler(err_act) return diff --git a/docs/src/datastruct.tex b/docs/src/datastruct.tex index 1ea2f816..17e7c237 100644 --- a/docs/src/datastruct.tex +++ b/docs/src/datastruct.tex @@ -826,13 +826,14 @@ A variable of type \verb|psb_Tspmat_type|. \subsubsection*{tril --- Return the lower triangle} \addcontentsline{toc}{paragraph}{tril} \begin{verbatim} - call a%tril(b,info[,& - & diag,imin,imax,jmin,jmax,rscale,cscale]) + call a%tril(l,info[,& + & diag,imin,imax,jmin,jmax,rscale,cscale,u]) \end{verbatim} Returns the lower triangular part of submatrix \verb|A(imin:imax,jmin:jmax)|, optionally rescaling row/col indices to -the range \verb|1:imax-imin+1,1:jmax-jmin+1|. +the range \verb|1:imax-imin+1,1:jmax-jmin+1| and returing the +complementary upper triangle. \begin{description} \item[Type:] Asynchronous. \item[\bf On Entry] @@ -849,7 +850,9 @@ Type: optional. \end{description} \begin{description} \item[\bf On Return] -\item[b] A copy of a subtriangle of \verb|a|.\\ +\item[l] A copy of the lower triangle of \verb|a|.\\ +A variable of type \verb|psb_Tspmat_type|. +\item[u] (optional) A copy of the upper triangle of \verb|a|.\\ A variable of type \verb|psb_Tspmat_type|. \item[info] Return code. \end{description} @@ -857,13 +860,14 @@ A variable of type \verb|psb_Tspmat_type|. \subsubsection*{triu --- Return the upper triangle} \addcontentsline{toc}{paragraph}{triu} \begin{verbatim} - call a%triu(b,info[,& - & diag,imin,imax,jmin,jmax,rscale,cscale]) + call a%triu(u,info[,& + & diag,imin,imax,jmin,jmax,rscale,cscale,l]) \end{verbatim} Returns the upper triangular part of submatrix \verb|A(imin:imax,jmin:jmax)|, optionally rescaling row/col indices to -the range \verb|1:imax-imin+1,1:jmax-jmin+1|. +the range \verb|1:imax-imin+1,1:jmax-jmin+1|, and returing the +complementary lower triangle. \begin{description} \item[Type:] Asynchronous. \item[\bf On Entry] @@ -880,7 +884,9 @@ Type: optional. \end{description} \begin{description} \item[\bf On Return] -\item[b] A copy of a subtriangle of \verb|a|.\\ +\item[u] A copy of the upper triangle of \verb|a|.\\ +A variable of type \verb|psb_Tspmat_type|. +\item[l] (optional) A copy of the lower triangle of \verb|a|.\\ A variable of type \verb|psb_Tspmat_type|. \item[info] Return code. \end{description} diff --git a/docs/src/userguide.tex b/docs/src/userguide.tex index 02b9e06e..c9095496 100644 --- a/docs/src/userguide.tex +++ b/docs/src/userguide.tex @@ -111,7 +111,7 @@ by Salvatore Filippone\\ and Alfredo Buttari}\\ University of Rome ``Tor Vergata''.\\[3ex] -May 1st, 2017 +Sep 1st, 2017 \end{minipage}} %\addtolength{\textwidth}{\centeroffset} diff --git a/docs/src/userhtml.tex b/docs/src/userhtml.tex index abf87db1..43436311 100644 --- a/docs/src/userhtml.tex +++ b/docs/src/userhtml.tex @@ -96,7 +96,7 @@ %\today Software version: 3.5.0\\ %\today -May 1st, 2017 +Sep 1st, 2017 \cleardoublepage \begingroup \renewcommand*{\thepage}{toc}