fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-13 09:29:04 +00:00
Code Issues Packages Projects Releases Wiki Activity

Define and test GLOBTRIL/GLOBTRIU in template

Browse Source
This commit is contained in:
sfilippone
2026-09-14 14:01:31 +02:00
parent 6263059751
commit e255d401bc
15 changed files with 2247 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/serial/psb_c_base_mat_mod.F90
+87 -1
View File
@@ -375,6 +375,8 @@ module psb_c_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_lc_mv_coo_from_fmt
procedure, pass(a) :: cp_to_icoo => psb_lc_cp_coo_to_icoo
procedure, pass(a) :: cp_from_icoo => psb_lc_cp_coo_from_icoo
procedure, pass(a) :: tril => psb_lc_coo_tril
procedure, pass(a) :: triu => psb_lc_coo_triu
procedure, pass(a) :: csput_a => psb_lc_coo_csput_a
procedure, pass(a) :: get_diag => psb_lc_coo_get_diag
@@ -1798,7 +1800,6 @@ module psb_c_base_mat_mod
end subroutine psb_c_coo_mold
end interface
!
!> Function print.
!! \memberof psb_c_coo_sparse_mat
@@ -3548,7 +3549,92 @@ module psb_c_base_mat_mod
integer(psb_ipk_), intent(out) :: info
end subroutine psb_lc_coo_mold
end interface
!
!> Function tril:
!! \memberof psb_d_coo_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_coo_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_lc_coo_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_coo_tril
end interface
!
!> Function triu:
!! \memberof psb_d_coo_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_coo_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_lc_coo_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_coo_triu
end interface
!
!> Function print.
base/modules/serial/psb_d_base_mat_mod.F90
+2 -4
View File
@@ -375,8 +375,8 @@ module psb_d_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_ld_mv_coo_from_fmt
procedure, pass(a) :: cp_to_icoo => psb_ld_cp_coo_to_icoo
procedure, pass(a) :: cp_from_icoo => psb_ld_cp_coo_from_icoo
procedure, pass(a) :: tril => psb_ld_coo_tril
procedure, pass(a) :: triu => psb_ld_coo_triu
procedure, pass(a) :: tril => psb_ld_coo_tril
procedure, pass(a) :: triu => psb_ld_coo_triu
procedure, pass(a) :: csput_a => psb_ld_coo_csput_a
procedure, pass(a) :: get_diag => psb_ld_coo_get_diag
@@ -1800,7 +1800,6 @@ module psb_d_base_mat_mod
end subroutine psb_d_coo_mold
end interface
!
!> Function print.
!! \memberof psb_d_coo_sparse_mat
@@ -3637,7 +3636,6 @@ module psb_d_base_mat_mod
end subroutine psb_ld_coo_triu
end interface
!
!> Function print.
!! \memberof psb_ld_coo_sparse_mat
base/modules/serial/psb_s_base_mat_mod.F90
+87 -1
View File
@@ -375,6 +375,8 @@ module psb_s_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_ls_mv_coo_from_fmt
procedure, pass(a) :: cp_to_icoo => psb_ls_cp_coo_to_icoo
procedure, pass(a) :: cp_from_icoo => psb_ls_cp_coo_from_icoo
procedure, pass(a) :: tril => psb_ls_coo_tril
procedure, pass(a) :: triu => psb_ls_coo_triu
procedure, pass(a) :: csput_a => psb_ls_coo_csput_a
procedure, pass(a) :: get_diag => psb_ls_coo_get_diag
@@ -1798,7 +1800,6 @@ module psb_s_base_mat_mod
end subroutine psb_s_coo_mold
end interface
!
!> Function print.
!! \memberof psb_s_coo_sparse_mat
@@ -3548,7 +3549,92 @@ module psb_s_base_mat_mod
integer(psb_ipk_), intent(out) :: info
end subroutine psb_ls_coo_mold
end interface
!
!> Function tril:
!! \memberof psb_d_coo_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_coo_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_ls_coo_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_coo_tril
end interface
!
!> Function triu:
!! \memberof psb_d_coo_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_coo_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_ls_coo_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_coo_triu
end interface
!
!> Function print.
base/modules/serial/psb_z_base_mat_mod.F90
+87 -1
View File
@@ -375,6 +375,8 @@ module psb_z_base_mat_mod
procedure, pass(a) :: mv_from_fmt => psb_lz_mv_coo_from_fmt
procedure, pass(a) :: cp_to_icoo => psb_lz_cp_coo_to_icoo
procedure, pass(a) :: cp_from_icoo => psb_lz_cp_coo_from_icoo
procedure, pass(a) :: tril => psb_lz_coo_tril
procedure, pass(a) :: triu => psb_lz_coo_triu
procedure, pass(a) :: csput_a => psb_lz_coo_csput_a
procedure, pass(a) :: get_diag => psb_lz_coo_get_diag
@@ -1798,7 +1800,6 @@ module psb_z_base_mat_mod
end subroutine psb_z_coo_mold
end interface
!
!> Function print.
!! \memberof psb_z_coo_sparse_mat
@@ -3548,7 +3549,92 @@ module psb_z_base_mat_mod
integer(psb_ipk_), intent(out) :: info
end subroutine psb_lz_coo_mold
end interface
!
!> Function tril:
!! \memberof psb_d_coo_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_coo_tril(a,l,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,u)
import
class(psb_lz_coo_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_coo_tril
end interface
!
!> Function triu:
!! \memberof psb_d_coo_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_coo_triu(a,u,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale,l)
import
class(psb_lz_coo_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_coo_triu
end interface
!
!> Function print.
base/modules/tools/psb_c_tools_mod.F90
+24
View File
@@ -435,6 +435,30 @@ Module psb_c_tools_mod
end subroutine psb_c_simple_glob_transpose_ip
end interface psb_glob_transpose
interface psb_globtril
Subroutine psb_cglobtril(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_cspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_cspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_cglobtril
end interface psb_globtril
interface psb_globtriu
Subroutine psb_cglobtriu(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_cspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_cspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_cglobtriu
end interface psb_globtriu
interface psb_getelem
function psb_c_getelem(x,index,desc_a,info) result(res)
import
base/modules/tools/psb_s_tools_mod.F90
+24
View File
@@ -435,6 +435,30 @@ Module psb_s_tools_mod
end subroutine psb_s_simple_glob_transpose_ip
end interface psb_glob_transpose
interface psb_globtril
Subroutine psb_sglobtril(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_sspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_sspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_sglobtril
end interface psb_globtril
interface psb_globtriu
Subroutine psb_sglobtriu(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_sspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_sspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_sglobtriu
end interface psb_globtriu
interface psb_getelem
function psb_s_getelem(x,index,desc_a,info) result(res)
import
base/modules/tools/psb_z_tools_mod.F90
+24
View File
@@ -435,6 +435,30 @@ Module psb_z_tools_mod
end subroutine psb_z_simple_glob_transpose_ip
end interface psb_glob_transpose
interface psb_globtril
Subroutine psb_zglobtril(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_zspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_zspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_zglobtril
end interface psb_globtril
interface psb_globtriu
Subroutine psb_zglobtriu(a,desc_a,b,info,diag,imin,imax,jmin,jmax)
import
implicit none
Type(psb_zspmat_type), Intent(inout) :: a
Type(psb_desc_type), Intent(inout), target :: desc_a
Type(psb_zspmat_type), Intent(out) :: b
integer(psb_ipk_), intent(out) :: info
integer(psb_lpk_), intent(in), optional :: diag,imin,imax,jmin,jmax
end Subroutine psb_zglobtriu
end interface psb_globtriu
interface psb_getelem
function psb_z_getelem(x,index,desc_a,info) result(res)
import
base/serial/impl/psb_c_coo_impl.F90
+594
View File
@@ -5904,7 +5904,601 @@ subroutine psb_lc_coo_reinit(a,clear)
end subroutine psb_lc_coo_reinit
!
! COO implementation of tril/triu
!
subroutine psb_lc_coo_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_base_mat_mod, psb_protect_name => psb_lc_coo_tril
implicit none
class(psb_lc_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
& nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,lrws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
lnz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws) reduction(+: lnz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,lrws,info)
!$omp parallel do private(i,j,k,lpnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call l%set_nzeros(lnz)
call l%fix(info)
end if
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-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 loop1
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, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = j
l%val(nzin) = val(k)
end if
end if
end do loop2
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
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_lc_coo_tril
subroutine psb_lc_coo_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_base_mat_mod, psb_protect_name => psb_lc_coo_triu
implicit none
class(psb_lc_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
&nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <=-1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.false.)
end if
else
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,urws) reduction(+: unz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,urws,info)
!$omp parallel do private(i,j,k,upnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call u%set_nzeros(unz)
call u%fix(info)
end if
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
else
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
end if
end if
end do loop1
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = j
u%val(nzin) = val(k)
end if
end if
end do loop2
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_lower(.false.)
end if
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_lc_coo_triu
subroutine psb_lc_coo_trim(a)
use psb_c_base_mat_mod, psb_protect_name => psb_lc_coo_trim
base/serial/impl/psb_d_coo_impl.F90
+1 -2
View File
@@ -4072,7 +4072,7 @@ subroutine psb_d_coo_triu(a,u,info,&
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
@@ -6500,7 +6500,6 @@ subroutine psb_ld_coo_triu(a,u,info,&
end subroutine psb_ld_coo_triu
subroutine psb_ld_coo_trim(a)
use psb_d_base_mat_mod, psb_protect_name => psb_ld_coo_trim
use psb_realloc_mod
base/serial/impl/psb_s_coo_impl.F90
+594
View File
@@ -5904,7 +5904,601 @@ subroutine psb_ls_coo_reinit(a,clear)
end subroutine psb_ls_coo_reinit
!
! COO implementation of tril/triu
!
subroutine psb_ls_coo_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_base_mat_mod, psb_protect_name => psb_ls_coo_tril
implicit none
class(psb_ls_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
& nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,lrws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
lnz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws) reduction(+: lnz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,lrws,info)
!$omp parallel do private(i,j,k,lpnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call l%set_nzeros(lnz)
call l%fix(info)
end if
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-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 loop1
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, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = j
l%val(nzin) = val(k)
end if
end if
end do loop2
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
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_ls_coo_tril
subroutine psb_ls_coo_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_base_mat_mod, psb_protect_name => psb_ls_coo_triu
implicit none
class(psb_ls_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
&nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <=-1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.false.)
end if
else
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,urws) reduction(+: unz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,urws,info)
!$omp parallel do private(i,j,k,upnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call u%set_nzeros(unz)
call u%fix(info)
end if
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
else
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
end if
end if
end do loop1
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = j
u%val(nzin) = val(k)
end if
end if
end do loop2
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_lower(.false.)
end if
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_ls_coo_triu
subroutine psb_ls_coo_trim(a)
use psb_s_base_mat_mod, psb_protect_name => psb_ls_coo_trim
base/serial/impl/psb_z_coo_impl.F90
+594
View File
@@ -5904,7 +5904,601 @@ subroutine psb_lz_coo_reinit(a,clear)
end subroutine psb_lz_coo_reinit
!
! COO implementation of tril/triu
!
subroutine psb_lz_coo_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_base_mat_mod, psb_protect_name => psb_lz_coo_tril
implicit none
class(psb_lz_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
& nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,lrws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call l%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >=-1).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_lower(.false.)
end if
else
lnz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws) reduction(+: lnz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,lrws,info)
!$omp parallel do private(i,j,k,lpnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call l%set_nzeros(lnz)
call l%fix(info)
end if
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-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 loop1
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, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<=diag_) then
nzin = nzin + 1
l%ia(nzin) = i
l%ja(nzin) = j
l%val(nzin) = val(k)
end if
end if
end do loop2
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
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_lz_coo_tril
subroutine psb_lz_coo_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_base_mat_mod, psb_protect_name => psb_lz_coo_triu
implicit none
class(psb_lz_coo_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_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_, &
&nzlin, nzuin, nz, nzin, nzout, i, j, k
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
#if defined(PSB_OPENMP)
block
integer(psb_lpk_), allocatable :: lrws(:),urws(:)
integer(psb_lpk_) :: lpnt, upnt, lnz, unz
call psb_realloc(mb,urws,info)
!$omp workshare
urws(:) = 0
!$omp end workshare
nz = a%get_nzeros()
call u%allocate(mb,nb,nz)
!write(0,*) 'Invocation of COO%TRIL', present(u),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
if (info == 0) call psb_realloc(mb,urws,info)
!$omp workshare
lrws(:) = 0
!$omp end workshare
!write(0,*) 'omp version of COO%TRIL/TRIU'
lnz = 0
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,lrws,urws) reduction(+: lnz,unz)
loop1: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic update
lrws(i-imin_+1) = lrws(i-imin_+1) +1
!$omp end atomic
lnz = lnz + 1
else
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop1
!$omp end parallel do
call psi_exscan(mb,lrws,info)
call psi_exscan(mb,urws,info)
!write(0,*) lrws(:), urws(:)
!$omp parallel do private(i,j,k,lpnt,upnt) shared(imin_,imax_,a)
loop2: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)<diag_) then
!$omp atomic capture
lrws(i-imin_+1) = lrws(i-imin_+1) +1
lpnt = lrws(i-imin_+1)
!$omp end atomic
l%ia(lpnt) = a%ia(k)
l%ja(lpnt) = a%ja(k)
l%val(lpnt) = a%val(k)
else
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop2
!$omp end parallel do
!write(0,*) 'End of copyout',lnz,unz
call l%set_nzeros(lnz)
call l%fix(info)
call u%set_nzeros(unz)
call u%fix(info)
nzout = l%get_nzeros()
if (rscale_) then
!$omp workshare
l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ <=-1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.false.)
end if
else
unz = 0
!$omp parallel do private(i,j,k) shared(imin_,imax_,a,urws) reduction(+: unz)
loop3: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic update
urws(i-imin_+1) = urws(i-imin_+1) +1
!$omp end atomic
unz = unz + 1
end if
end if
end do loop3
!$omp end parallel do
call psi_exscan(mb,urws,info)
!$omp parallel do private(i,j,k,upnt) shared(imin_,imax_,a)
loop4: do k=1,nz
i = a%ia(k)
j = a%ja(k)
if ((i>=imin_).and.(i<=imax_).and.(jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
!$omp atomic capture
urws(i-imin_+1) = urws(i-imin_+1) +1
upnt = urws(i-imin_+1)
!$omp end atomic
u%ia(upnt) = a%ia(k)
u%ja(upnt) = a%ja(k)
u%val(upnt) = a%val(k)
end if
end if
end do loop4
!$omp end parallel do
call u%set_nzeros(unz)
call u%fix(info)
end if
nzout = u%get_nzeros()
if (rscale_) then
!$omp workshare
u%ia(1:nzout) = u%ia(1:nzout) - imin_ + 1
!$omp end workshare
end if
if (cscale_) then
!$omp workshare
u%ja(1:nzout) = u%ja(1:nzout) - jmin_ + 1
!$omp end workshare
end if
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call u%set_triangle(.true.)
call u%set_upper(.true.)
end if
end block
#else
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, ia=>a%ia)
loop1: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzuin = nzuin + 1
u%ia(nzuin) = i
u%ja(nzuin) = ja(k)
u%val(nzuin) = val(k)
else
nzlin = nzlin + 1
l%ia(nzlin) = i
l%ja(nzlin) = ja(k)
l%val(nzlin) = val(k)
end if
end if
end do loop1
end associate
call u%set_nzeros(nzuin)
call l%set_nzeros(nzlin)
call l%fix(info)
nzout = l%get_nzeros()
if (rscale_) &
& l%ia(1:nzout) = l%ia(1:nzout) - imin_ + 1
if (cscale_) &
& l%ja(1:nzout) = l%ja(1:nzout) - jmin_ + 1
if ((diag_ <=1).and.(imin_ == jmin_)) then
call l%set_triangle(.true.)
call l%set_lower(.true.)
end if
else
nzin = u%get_nzeros() ! At this point it should be 0
associate(val =>a%val, ja => a%ja, ia=>a%ia)
loop2: do k=1,nz
i = ia(k)
j = ja(k)
if ((jmin_<=j).and.(j<=jmax_)) then
if ((j-i)>=diag_) then
nzin = nzin + 1
u%ia(nzin) = i
u%ja(nzin) = j
u%val(nzin) = val(k)
end if
end if
end do loop2
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_lower(.false.)
end if
#endif
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine psb_lz_coo_triu
subroutine psb_lz_coo_trim(a)
use psb_z_base_mat_mod, psb_protect_name => psb_lz_coo_trim
test/pdegen/psb_d_pde2d.F90
+33
View File
@@ -825,6 +825,39 @@ program psb_d_pde2d
write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt()
end if
if (.false.) then
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_dspmat_type) :: atril
call psb_globtril(a,desc_a,atril,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atril-loc-",iam,".mtx"
call atril%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atril-glb-",iam,".mtx"
call atril%print(fname,iv=iv)
end block
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_dspmat_type) :: atriu
call psb_globtriu(a,desc_a,atriu,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atriu-loc-",iam,".mtx"
call atriu%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atriu-glb-",iam,".mtx"
call atriu%print(fname,iv=iv)
end block
end if
!
! cleanup storage and exit
test/pdegen/psb_d_pde3d.F90
+29 -29
View File
@@ -728,7 +728,7 @@ program psb_d_pde3d
#else
nth = 1
#endif
if (iam < 0) then
! This should not happen, but just in case
call psb_exit(ctxt)
@@ -773,33 +773,33 @@ program psb_d_pde3d
! Set the options for the BJAC preconditioner
!
if (psb_toupper(ptype) == "BJAC") then
call prec%set('sub_solve', parms%alg, info)
select case (psb_toupper(parms%alg))
case ("ILU")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('ilu_alg', parms%ilu_alg, info)
case ("ILUT")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("AINV")
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
call prec%set('ainv_alg', parms%orth_alg, info)
case ("INVK")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("INVT")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case default
! Do nothing, use default setting in the init routine
end select
call prec%set('sub_solve', parms%alg, info)
select case (psb_toupper(parms%alg))
case ("ILU")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('ilu_alg', parms%ilu_alg, info)
case ("ILUT")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("AINV")
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
call prec%set('ainv_alg', parms%orth_alg, info)
case ("INVK")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case ("INVT")
call prec%set('sub_fillin', parms%fill, info)
call prec%set('inv_fillin', parms%inv_fill, info)
call prec%set('sub_iluthrs', parms%thresh, info)
call prec%set('inv_thresh', parms%inv_thresh, info)
call prec%set('ilut_scale', parms%ilut_scale, info)
case default
! Do nothing, use default setting in the init routine
end select
else
! nothing to set for NONE or DIAG preconditioner
end if
@@ -1064,7 +1064,7 @@ contains
write(psb_out_unit,'("Orthogonalization : ",a)') parms%orth_alg
write(psb_out_unit,'("Scaling : ",a)') parms%ilut_scale
case default
write(psb_out_unit,'("Unknown diagonal solver")')
write(psb_out_unit,'("Unknown diagonal solver")')
end select
end if
write(psb_out_unit,'("Iterative method : ",a)') kmethd
test/pdegen/psb_s_pde2d.F90
+33
View File
@@ -825,6 +825,39 @@ program psb_s_pde2d
write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt()
end if
if (.false.) then
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_sspmat_type) :: atril
call psb_globtril(a,desc_a,atril,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atril-loc-",iam,".mtx"
call atril%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atril-glb-",iam,".mtx"
call atril%print(fname,iv=iv)
end block
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_sspmat_type) :: atriu
call psb_globtriu(a,desc_a,atriu,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atriu-loc-",iam,".mtx"
call atriu%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atriu-glb-",iam,".mtx"
call atriu%print(fname,iv=iv)
end block
end if
!
! cleanup storage and exit
test/pdegen/psb_s_pde3d.F90
+34
View File
@@ -881,6 +881,40 @@ program psb_s_pde3d
write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt()
end if
if (.false.) then
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_sspmat_type) :: atril
call psb_globtril(a,desc_a,atril,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atril-loc-",iam,".mtx"
call atril%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atril-glb-",iam,".mtx"
call atril%print(fname,iv=iv)
end block
block
integer(psb_lpk_), allocatable :: iv(:)
type(psb_sspmat_type) :: atriu
call psb_globtriu(a,desc_a,atriu,info)
write(fname,'(a,i0,a)') "amat-loc-",iam,".mtx"
call a%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "amat-glb-",iam,".mtx"
call a%print(fname,iv=iv)
write(fname,'(a,i0,a)') "atriu-loc-",iam,".mtx"
call atriu%print(fname)
iv = desc_a%get_global_indices(owned=.false.)
write(fname,'(a,i0,a)') "atriu-glb-",iam,".mtx"
call atriu%print(fname,iv=iv)
end block
end if
!
! cleanup storage and exit