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

psblas-3.99:

Browse Source 
base/modules/psb_c_base_mat_mod.f90
 base/modules/psb_c_mat_mod.f90
 base/modules/psb_d_base_mat_mod.f90
 base/modules/psb_d_mat_mod.f90
 base/modules/psb_s_base_mat_mod.f90
 base/modules/psb_s_mat_mod.f90
 base/modules/psb_z_base_mat_mod.f90
 base/modules/psb_z_mat_mod.f90
 base/serial/impl/psb_c_base_mat_impl.F90
 base/serial/impl/psb_c_mat_impl.F90
 base/serial/impl/psb_d_base_mat_impl.F90
 base/serial/impl/psb_d_mat_impl.F90
 base/serial/impl/psb_s_base_mat_impl.F90
 base/serial/impl/psb_s_mat_impl.F90
 base/serial/impl/psb_z_base_mat_impl.F90
 base/serial/impl/psb_z_mat_impl.F90

Defined TRIL and TRIU methods. To be tested.
This commit is contained in:
Salvatore Filippone
2013-10-09 16:57:56 +00:00
parent 0be3cdb6c3
commit c42aa65e47
16 changed files with 1788 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/psb_c_base_mat_mod.f90
+89 -5
View File
@@ -62,6 +62,8 @@ module psb_c_base_mat_mod
procedure, pass(a) :: csgetblk => psb_c_base_csgetblk
procedure, pass(a) :: get_diag => psb_c_base_get_diag
generic, public :: csget => csgetrow, csgetblk
procedure, pass(a) :: tril => psb_c_base_tril
procedure, pass(a) :: triu => psb_c_base_triu
procedure, pass(a) :: csclip => psb_c_base_csclip
procedure, pass(a) :: cp_to_coo => psb_c_base_cp_to_coo
procedure, pass(a) :: cp_from_coo => psb_c_base_cp_from_coo
@@ -312,7 +314,7 @@ module psb_c_base_mat_mod
!! \param info return code
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -350,7 +352,7 @@ module psb_c_base_mat_mod
!! \param imax [a%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -364,11 +366,93 @@ module psb_c_base_mat_mod
import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_c_base_csclip
end interface
!
!> Function tril:
!! \memberof psb_c_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_c_base_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_c_base_tril
end interface
!
!> Function triu:
!! \memberof psb_c_base_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_c_base_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_c_base_sparse_mat, psb_c_coo_sparse_mat, psb_spk_
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_c_base_triu
end interface
!
!> Function get_diag:
base/modules/psb_c_mat_mod.f90
+27
View File
@@ -122,6 +122,8 @@ module psb_c_mat_mod
procedure, pass(a) :: csgetrow => psb_c_csgetrow
procedure, pass(a) :: csgetblk => psb_c_csgetblk
generic, public :: csget => csgetptn, csgetrow, csgetblk
procedure, pass(a) :: tril => psb_c_tril
procedure, pass(a) :: triu => psb_c_triu
procedure, pass(a) :: m_csclip => psb_c_csclip
procedure, pass(a) :: b_csclip => psb_c_b_csclip
generic, public :: csclip => b_csclip, m_csclip
@@ -420,6 +422,31 @@ module psb_c_mat_mod
end subroutine psb_c_csgetblk
end interface
interface
subroutine psb_c_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_cspmat_type, psb_spk_
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_c_tril
end interface
interface
subroutine psb_c_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_cspmat_type, psb_spk_
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_c_triu
end interface
interface
subroutine psb_c_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
base/modules/psb_d_base_mat_mod.f90
+89 -5
View File
@@ -62,6 +62,8 @@ module psb_d_base_mat_mod
procedure, pass(a) :: csgetblk => psb_d_base_csgetblk
procedure, pass(a) :: get_diag => psb_d_base_get_diag
generic, public :: csget => csgetrow, csgetblk
procedure, pass(a) :: tril => psb_d_base_tril
procedure, pass(a) :: triu => psb_d_base_triu
procedure, pass(a) :: csclip => psb_d_base_csclip
procedure, pass(a) :: cp_to_coo => psb_d_base_cp_to_coo
procedure, pass(a) :: cp_from_coo => psb_d_base_cp_from_coo
@@ -312,7 +314,7 @@ module psb_d_base_mat_mod
!! \param info return code
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -350,7 +352,7 @@ module psb_d_base_mat_mod
!! \param imax [a%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -364,11 +366,93 @@ module psb_d_base_mat_mod
import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_d_base_csclip
end interface
!
!> Function tril:
!! \memberof psb_d_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_d_base_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_d_base_tril
end interface
!
!> Function triu:
!! \memberof psb_d_base_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_d_base_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_d_base_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_d_base_triu
end interface
!
!> Function get_diag:
base/modules/psb_d_mat_mod.f90
+27
View File
@@ -122,6 +122,8 @@ module psb_d_mat_mod
procedure, pass(a) :: csgetrow => psb_d_csgetrow
procedure, pass(a) :: csgetblk => psb_d_csgetblk
generic, public :: csget => csgetptn, csgetrow, csgetblk
procedure, pass(a) :: tril => psb_d_tril
procedure, pass(a) :: triu => psb_d_triu
procedure, pass(a) :: m_csclip => psb_d_csclip
procedure, pass(a) :: b_csclip => psb_d_b_csclip
generic, public :: csclip => b_csclip, m_csclip
@@ -420,6 +422,31 @@ module psb_d_mat_mod
end subroutine psb_d_csgetblk
end interface
interface
subroutine psb_d_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_dspmat_type, psb_dpk_
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_d_tril
end interface
interface
subroutine psb_d_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_dspmat_type, psb_dpk_
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_d_triu
end interface
interface
subroutine psb_d_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
base/modules/psb_s_base_mat_mod.f90
+89 -5
View File
@@ -62,6 +62,8 @@ module psb_s_base_mat_mod
procedure, pass(a) :: csgetblk => psb_s_base_csgetblk
procedure, pass(a) :: get_diag => psb_s_base_get_diag
generic, public :: csget => csgetrow, csgetblk
procedure, pass(a) :: tril => psb_s_base_tril
procedure, pass(a) :: triu => psb_s_base_triu
procedure, pass(a) :: csclip => psb_s_base_csclip
procedure, pass(a) :: cp_to_coo => psb_s_base_cp_to_coo
procedure, pass(a) :: cp_from_coo => psb_s_base_cp_from_coo
@@ -312,7 +314,7 @@ module psb_s_base_mat_mod
!! \param info return code
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -350,7 +352,7 @@ module psb_s_base_mat_mod
!! \param imax [a%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -364,11 +366,93 @@ module psb_s_base_mat_mod
import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_s_base_csclip
end interface
!
!> Function tril:
!! \memberof psb_s_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_s_base_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_s_base_tril
end interface
!
!> Function triu:
!! \memberof psb_s_base_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_s_base_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_s_base_sparse_mat, psb_s_coo_sparse_mat, psb_spk_
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_s_base_triu
end interface
!
!> Function get_diag:
base/modules/psb_s_mat_mod.f90
+27
View File
@@ -122,6 +122,8 @@ module psb_s_mat_mod
procedure, pass(a) :: csgetrow => psb_s_csgetrow
procedure, pass(a) :: csgetblk => psb_s_csgetblk
generic, public :: csget => csgetptn, csgetrow, csgetblk
procedure, pass(a) :: tril => psb_s_tril
procedure, pass(a) :: triu => psb_s_triu
procedure, pass(a) :: m_csclip => psb_s_csclip
procedure, pass(a) :: b_csclip => psb_s_b_csclip
generic, public :: csclip => b_csclip, m_csclip
@@ -420,6 +422,31 @@ module psb_s_mat_mod
end subroutine psb_s_csgetblk
end interface
interface
subroutine psb_s_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_sspmat_type, psb_spk_
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_s_tril
end interface
interface
subroutine psb_s_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_sspmat_type, psb_spk_
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_s_triu
end interface
interface
subroutine psb_s_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
base/modules/psb_z_base_mat_mod.f90
+89 -5
View File
@@ -62,6 +62,8 @@ module psb_z_base_mat_mod
procedure, pass(a) :: csgetblk => psb_z_base_csgetblk
procedure, pass(a) :: get_diag => psb_z_base_get_diag
generic, public :: csget => csgetrow, csgetblk
procedure, pass(a) :: tril => psb_z_base_tril
procedure, pass(a) :: triu => psb_z_base_triu
procedure, pass(a) :: csclip => psb_z_base_csclip
procedure, pass(a) :: cp_to_coo => psb_z_base_cp_to_coo
procedure, pass(a) :: cp_from_coo => psb_z_base_cp_from_coo
@@ -312,7 +314,7 @@ module psb_z_base_mat_mod
!! \param info return code
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -350,7 +352,7 @@ module psb_z_base_mat_mod
!! \param imax [a%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
@@ -364,11 +366,93 @@ module psb_z_base_mat_mod
import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_z_base_csclip
end interface
!
!> Function tril:
!! \memberof psb_z_base_sparse_mat
!! \brief Copy the lower triangle, i.e. all entries
!! A(I,J) such that J-I <= DIAG
!! default value is DIAG=0, i.e. lower triangle up to
!! the main diagonal.
!! DIAG=-1 means copy the strictly lower triangle
!! DIAG= 1 means copy the lower triangle plus the first diagonal
!! of the upper triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_z_base_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_z_base_tril
end interface
!
!> Function triu:
!! \memberof psb_z_base_sparse_mat
!! \brief Copy the upper triangle, i.e. all entries
!! A(I,J) such that DIAG <= J-I
!! default value is DIAG=0, i.e. upper triangle from
!! the main diagonal up.
!! DIAG= 1 means copy the strictly upper triangle
!! DIAG=-1 means copy the upper triangle plus the first diagonal
!! of the lower triangle.
!! Moreover, apply a clipping by copying entries A(I,J) only if
!! IMIN<=I<=IMAX
!! JMIN<=J<=JMAX
!!
!! \param b the output (sub)matrix
!! \param info return code
!! \param diag [0] the last diagonal (J-I) to be considered.
!! \param imin [1] the minimum row index we are interested in
!! \param imax [a\%get_nrows()] the minimum row index we are interested in
!! \param jmin [1] minimum col index
!! \param jmax [a\%get_ncols()] maximum col index
!! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:))
!! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1]
!! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1]
!! ( iren cannot be specified with rscale/cscale)
!! \param append [false] append to ia,ja
!! \param nzin [none] if append, then first new entry should go in entry nzin+1
!!
!
interface
subroutine psb_z_base_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_z_base_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_z_base_triu
end interface
!
!> Function get_diag:
base/modules/psb_z_mat_mod.f90
+27
View File
@@ -122,6 +122,8 @@ module psb_z_mat_mod
procedure, pass(a) :: csgetrow => psb_z_csgetrow
procedure, pass(a) :: csgetblk => psb_z_csgetblk
generic, public :: csget => csgetptn, csgetrow, csgetblk
procedure, pass(a) :: tril => psb_z_tril
procedure, pass(a) :: triu => psb_z_triu
procedure, pass(a) :: m_csclip => psb_z_csclip
procedure, pass(a) :: b_csclip => psb_z_b_csclip
generic, public :: csclip => b_csclip, m_csclip
@@ -420,6 +422,31 @@ module psb_z_mat_mod
end subroutine psb_z_csgetblk
end interface
interface
subroutine psb_z_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_zspmat_type, psb_dpk_
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_z_tril
end interface
interface
subroutine psb_z_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
import :: psb_ipk_, psb_zspmat_type, psb_dpk_
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
end subroutine psb_z_triu
end interface
interface
subroutine psb_z_csclip(a,b,info,&
& imin,imax,jmin,jmax,rscale,cscale)
base/serial/impl/psb_c_base_mat_impl.F90
+226
View File
@@ -550,6 +550,232 @@ subroutine psb_c_base_csclip(a,b,info,&
end subroutine psb_c_base_csclip
!
! Here we have the base implementation of tril and triu
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_c_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_tril
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_c_base_tril
subroutine psb_c_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_triu
implicit none
class(psb_c_base_sparse_mat), intent(in) :: a
class(psb_c_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_c_base_triu
subroutine psb_c_base_clone(a,b,info)
use psb_c_base_mat_mod, psb_protect_name => psb_c_base_clone
use psb_error_mod
base/serial/impl/psb_c_mat_impl.F90
+105
View File
@@ -942,6 +942,111 @@ subroutine psb_c_csgetblk(imin,imax,a,b,info,&
end subroutine psb_c_csgetblk
subroutine psb_c_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod
use psb_c_mat_mod, psb_protect_name => psb_c_tril
implicit none
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_c_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_c_tril
subroutine psb_c_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_c_base_mat_mod
use psb_c_mat_mod, psb_protect_name => psb_c_triu
implicit none
class(psb_cspmat_type), intent(in) :: a
class(psb_cspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_c_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_c_triu
subroutine psb_c_csclip(a,b,info,&
base/serial/impl/psb_d_base_mat_impl.F90
+226
View File
@@ -550,6 +550,232 @@ subroutine psb_d_base_csclip(a,b,info,&
end subroutine psb_d_base_csclip
!
! Here we have the base implementation of tril and triu
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_d_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_tril
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_d_base_tril
subroutine psb_d_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_triu
implicit none
class(psb_d_base_sparse_mat), intent(in) :: a
class(psb_d_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_d_base_triu
subroutine psb_d_base_clone(a,b,info)
use psb_d_base_mat_mod, psb_protect_name => psb_d_base_clone
use psb_error_mod
base/serial/impl/psb_d_mat_impl.F90
+105
View File
@@ -942,6 +942,111 @@ subroutine psb_d_csgetblk(imin,imax,a,b,info,&
end subroutine psb_d_csgetblk
subroutine psb_d_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod
use psb_d_mat_mod, psb_protect_name => psb_d_tril
implicit none
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_d_tril
subroutine psb_d_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_d_base_mat_mod
use psb_d_mat_mod, psb_protect_name => psb_d_triu
implicit none
class(psb_dspmat_type), intent(in) :: a
class(psb_dspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_d_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_d_triu
subroutine psb_d_csclip(a,b,info,&
base/serial/impl/psb_s_base_mat_impl.F90
+226
View File
@@ -550,6 +550,232 @@ subroutine psb_s_base_csclip(a,b,info,&
end subroutine psb_s_base_csclip
!
! Here we have the base implementation of tril and triu
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_s_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_tril
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_s_base_tril
subroutine psb_s_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_triu
implicit none
class(psb_s_base_sparse_mat), intent(in) :: a
class(psb_s_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_s_base_triu
subroutine psb_s_base_clone(a,b,info)
use psb_s_base_mat_mod, psb_protect_name => psb_s_base_clone
use psb_error_mod
base/serial/impl/psb_s_mat_impl.F90
+105
View File
@@ -942,6 +942,111 @@ subroutine psb_s_csgetblk(imin,imax,a,b,info,&
end subroutine psb_s_csgetblk
subroutine psb_s_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod
use psb_s_mat_mod, psb_protect_name => psb_s_tril
implicit none
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_s_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_s_tril
subroutine psb_s_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_s_base_mat_mod
use psb_s_mat_mod, psb_protect_name => psb_s_triu
implicit none
class(psb_sspmat_type), intent(in) :: a
class(psb_sspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_s_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_s_triu
subroutine psb_s_csclip(a,b,info,&
base/serial/impl/psb_z_base_mat_impl.F90
+226
View File
@@ -550,6 +550,232 @@ subroutine psb_z_base_csclip(a,b,info,&
end subroutine psb_z_base_csclip
!
! Here we have the base implementation of tril and triu
! this is just based on the getrow.
! If performance is critical it can be overridden.
!
subroutine psb_z_base_tril(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_tril
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = min(jmax_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=jmin_, jmax=k, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ <= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_lower(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_z_base_tril
subroutine psb_z_base_triu(a,b,info,&
& diag,imin,imax,jmin,jmax,rscale,cscale)
! Output is always in COO format
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_triu
implicit none
class(psb_z_base_sparse_mat), intent(in) :: a
class(psb_z_coo_sparse_mat), intent(out) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act, nzin, nzout, i, j, k
integer(psb_ipk_) :: imin_, imax_, jmin_, jmax_, mb,nb, diag_
integer(psb_ipk_) :: ierr(5)
character(len=20) :: name='tril'
logical :: rscale_, cscale_
logical, parameter :: debug=.false.
call psb_erractionsave(err_act)
info = psb_success_
if (present(diag)) then
diag_ = diag
else
diag_ = 0
end if
if (present(imin)) then
imin_ = imin
else
imin_ = 1
end if
if (present(imax)) then
imax_ = imax
else
imax_ = a%get_nrows()
end if
if (present(jmin)) then
jmin_ = jmin
else
jmin_ = 1
end if
if (present(jmax)) then
jmax_ = jmax
else
jmax_ = a%get_ncols()
end if
if (present(rscale)) then
rscale_ = rscale
else
rscale_ = .true.
end if
if (present(cscale)) then
cscale_ = cscale
else
cscale_ = .true.
end if
if (rscale_) then
mb = imax_ - imin_ +1
else
mb = imax_
endif
if (cscale_) then
nb = jmax_ - jmin_ +1
else
nb = jmax_
endif
call b%allocate(mb,nb)
nzin = b%get_nzeros() ! At this point it should be 0
do i=imin_,imax_
k = max(jmin_,i+diag_)
call a%csget(i,i,nzout,b%ia,b%ja,b%val,info,&
& jmin=k, jmax=jmax_, append=.true., &
& nzin=nzin)
if (info /= psb_success_) goto 9999
call b%set_nzeros(nzin+nzout)
end do
call b%fix(info)
nzout = b%get_nzeros()
if (rscale_) &
& b%ia(1:nzout) = b%ia(1:nzout) - imin_ + 1
if (cscale_) &
& b%ja(1:nzout) = b%ja(1:nzout) - jmin_ + 1
if ((diag_ >= 0).and.(imin_ == jmin_)) then
call b%set_triangle(.true.)
call b%set_upper(.true.)
end if
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine psb_z_base_triu
subroutine psb_z_base_clone(a,b,info)
use psb_z_base_mat_mod, psb_protect_name => psb_z_base_clone
use psb_error_mod
base/serial/impl/psb_z_mat_impl.F90
+105
View File
@@ -942,6 +942,111 @@ subroutine psb_z_csgetblk(imin,imax,a,b,info,&
end subroutine psb_z_csgetblk
subroutine psb_z_tril(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod
use psb_z_mat_mod, psb_protect_name => psb_z_tril
implicit none
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='tril'
logical, parameter :: debug=.false.
type(psb_z_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_z_tril
subroutine psb_z_triu(a,b,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
use psb_error_mod
use psb_const_mod
use psb_z_base_mat_mod
use psb_z_mat_mod, psb_protect_name => psb_z_triu
implicit none
class(psb_zspmat_type), intent(in) :: a
class(psb_zspmat_type), intent(inout) :: b
integer(psb_ipk_),intent(out) :: info
integer(psb_ipk_), intent(in), optional :: diag,imin,imax,jmin,jmax
logical, intent(in), optional :: rscale,cscale
integer(psb_ipk_) :: err_act
character(len=20) :: name='triu'
logical, parameter :: debug=.false.
type(psb_z_coo_sparse_mat), allocatable :: acoo
info = psb_success_
call psb_erractionsave(err_act)
if (a%is_null()) then
info = psb_err_invalid_mat_state_
call psb_errpush(info,name)
goto 9999
endif
allocate(acoo,stat=info)
if (info == psb_success_) then
call a%a%tril(acoo,info,diag,imin,imax,&
& jmin,jmax,rscale,cscale)
else
info = psb_err_alloc_dealloc_
end if
if (info == psb_success_) call move_alloc(acoo,b%a)
if (info == psb_success_) call b%cscnv(info,mold=a%a)
if (info /= psb_success_) goto 9999
call psb_erractionrestore(err_act)
return
9999 continue
call psb_erractionrestore(err_act)
if (err_act == psb_act_abort_) then
call psb_error()
return
end if
end subroutine psb_z_triu
subroutine psb_z_csclip(a,b,info,&