diff --git a/base/modules/Makefile b/base/modules/Makefile index 6a06611d..6a8594c2 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -7,8 +7,9 @@ UTIL_MODS = psb_string_mod.o psb_spmat_type.o \ psb_linmap_type_mod.o psb_comm_mod.o psb_psblas_mod.o \ psi_serial_mod.o psi_mod.o psb_ip_reord_mod.o\ psb_check_mod.o psb_gps_mod.o psb_linmap_mod.o psb_hash_mod.o\ - psb_base_mat_mod.o psb_d_base_mat_mod.o psb_mat_mod.o\ - psb_d_csr_mat_mod.o psb_d_mat_mod.o + psb_base_mat_mod.o psb_mat_mod.o\ + psb_d_base_mat_mod.o psb_d_csr_mat_mod.o psb_d_mat_mod.o \ + psb_s_base_mat_mod.o psb_s_csr_mat_mod.o psb_s_mat_mod.o MODULES=$(BASIC_MODS) $(UTIL_MODS) @@ -27,9 +28,10 @@ lib: $(BASIC_MODS) blacsmod $(UTIL_MODS) $(OBJS) $(LIBMOD) psb_base_mat_mod.o: psb_string_mod.o psb_sort_mod.o psb_ip_reord_mod.o psb_error_mod.o -psb_d_base_mat_mod.o: psb_base_mat_mod.o +psb_s_base_mat_mod.o psb_d_base_mat_mod.o: psb_base_mat_mod.o +psb_s_mat_mod.o: psb_s_base_mat_mod.o psb_s_csr_mat_mod.o psb_d_mat_mod.o: psb_d_base_mat_mod.o psb_d_csr_mat_mod.o -psb_mat_mod.o: psb_d_mat_mod.o +psb_mat_mod.o: psb_d_mat_mod.o psb_s_mat_mod.o psb_realloc_mod.o : psb_error_mod.o psb_spmat_type.o : psb_realloc_mod.o psb_error_mod.o psb_const_mod.o psb_string_mod.o psb_sort_mod.o psb_error_mod.o: psb_const_mod.o diff --git a/base/modules/psb_d_base_mat_mod.f03 b/base/modules/psb_d_base_mat_mod.f03 index 094c05fa..1af7ca85 100644 --- a/base/modules/psb_d_base_mat_mod.f03 +++ b/base/modules/psb_d_base_mat_mod.f03 @@ -824,8 +824,8 @@ contains use psb_error_mod implicit none class(psb_d_base_sparse_mat), intent(in) :: a - real(kind(1.d0)), intent(in) :: alpha, beta, x(:) - real(kind(1.d0)), intent(inout) :: y(:) + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) integer, intent(out) :: info character, optional, intent(in) :: trans diff --git a/base/modules/psb_d_mat_mod.f03 b/base/modules/psb_d_mat_mod.f03 index ea608e55..e4eb61ba 100644 --- a/base/modules/psb_d_mat_mod.f03 +++ b/base/modules/psb_d_mat_mod.f03 @@ -1631,8 +1631,8 @@ contains use psb_error_mod implicit none class(psb_d_sparse_mat), intent(in) :: a - real(kind(1.d0)), intent(in) :: alpha, beta, x(:,:) - real(kind(1.d0)), intent(inout) :: y(:,:) + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) integer, intent(out) :: info character, optional, intent(in) :: trans Integer :: err_act @@ -1667,8 +1667,8 @@ contains use psb_error_mod implicit none class(psb_d_sparse_mat), intent(in) :: a - real(kind(1.d0)), intent(in) :: alpha, beta, x(:) - real(kind(1.d0)), intent(inout) :: y(:) + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) integer, intent(out) :: info character, optional, intent(in) :: trans Integer :: err_act @@ -1703,8 +1703,8 @@ contains use psb_error_mod implicit none class(psb_d_sparse_mat), intent(in) :: a - real(kind(1.d0)), intent(in) :: alpha, beta, x(:,:) - real(kind(1.d0)), intent(inout) :: y(:,:) + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) integer, intent(out) :: info character, optional, intent(in) :: trans, side real(psb_dpk_), intent(in), optional :: d(:) @@ -1741,8 +1741,8 @@ contains use psb_error_mod implicit none class(psb_d_sparse_mat), intent(in) :: a - real(kind(1.d0)), intent(in) :: alpha, beta, x(:) - real(kind(1.d0)), intent(inout) :: y(:) + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) integer, intent(out) :: info character, optional, intent(in) :: trans, side real(psb_dpk_), intent(in), optional :: d(:) diff --git a/base/modules/psb_mat_mod.f03 b/base/modules/psb_mat_mod.f03 index f7d4d181..30fa1ebd 100644 --- a/base/modules/psb_mat_mod.f03 +++ b/base/modules/psb_mat_mod.f03 @@ -1,3 +1,4 @@ module psb_mat_mod + use psb_s_mat_mod use psb_d_mat_mod end module psb_mat_mod diff --git a/base/modules/psb_s_base_mat_mod.f03 b/base/modules/psb_s_base_mat_mod.f03 new file mode 100644 index 00000000..c00a290a --- /dev/null +++ b/base/modules/psb_s_base_mat_mod.f03 @@ -0,0 +1,2594 @@ +module psb_s_base_mat_mod + + use psb_base_mat_mod + + type, extends(psb_base_sparse_mat) :: psb_s_base_sparse_mat + contains + procedure, pass(a) :: s_base_csmv + procedure, pass(a) :: s_base_csmm + generic, public :: csmm => s_base_csmm, s_base_csmv + procedure, pass(a) :: s_base_cssv + procedure, pass(a) :: s_base_cssm + generic, public :: base_cssm => s_base_cssm, s_base_cssv + procedure, pass(a) :: s_cssv + procedure, pass(a) :: s_cssm + generic, public :: cssm => s_cssm, s_cssv + procedure, pass(a) :: s_scals + procedure, pass(a) :: s_scal + generic, public :: scal => s_scals, s_scal + procedure, pass(a) :: csnmi + procedure, pass(a) :: get_diag + procedure, pass(a) :: csput + + procedure, pass(a) :: s_csgetrow + procedure, pass(a) :: s_csgetblk + generic, public :: csget => s_csgetrow, s_csgetblk + procedure, pass(a) :: csclip + procedure, pass(a) :: cp_to_coo + procedure, pass(a) :: cp_from_coo + procedure, pass(a) :: cp_to_fmt + procedure, pass(a) :: cp_from_fmt + procedure, pass(a) :: mv_to_coo + procedure, pass(a) :: mv_from_coo + procedure, pass(a) :: mv_to_fmt + procedure, pass(a) :: mv_from_fmt + procedure, pass(a) :: s_base_cp_from + generic, public :: cp_from => s_base_cp_from + procedure, pass(a) :: s_base_mv_from + generic, public :: mv_from => s_base_mv_from + end type psb_s_base_sparse_mat + + private :: s_base_csmv, s_base_csmm, s_base_cssv, s_base_cssm,& + & s_scals, s_scal, csnmi, csput, s_csgetrow, s_csgetblk, & + & cp_to_coo, cp_from_coo, cp_to_fmt, cp_from_fmt, & + & mv_to_coo, mv_from_coo, mv_to_fmt, mv_from_fmt, & + & get_diag, csclip, s_cssv, s_cssm, base_cp_from, base_mv_from + + type, extends(psb_s_base_sparse_mat) :: psb_s_coo_sparse_mat + + integer :: nnz + integer, allocatable :: ia(:), ja(:) + real(psb_spk_), allocatable :: val(:) + + contains + + procedure, pass(a) :: get_size => s_coo_get_size + procedure, pass(a) :: get_nzeros => s_coo_get_nzeros + procedure, pass(a) :: set_nzeros => s_coo_set_nzeros + procedure, pass(a) :: s_base_csmm => s_coo_csmm + procedure, pass(a) :: s_base_csmv => s_coo_csmv + procedure, pass(a) :: s_base_cssm => s_coo_cssm + procedure, pass(a) :: s_base_cssv => s_coo_cssv + procedure, pass(a) :: s_scals => s_coo_scals + procedure, pass(a) :: s_scal => s_coo_scal + procedure, pass(a) :: csnmi => s_coo_csnmi + procedure, pass(a) :: csput => s_coo_csput + procedure, pass(a) :: get_diag => s_coo_get_diag + procedure, pass(a) :: reallocate_nz => s_coo_reallocate_nz + procedure, pass(a) :: allocate_mnnz => s_coo_allocate_mnnz + procedure, pass(a) :: cp_to_coo => s_cp_coo_to_coo + procedure, pass(a) :: cp_from_coo => s_cp_coo_from_coo + procedure, pass(a) :: cp_to_fmt => s_cp_coo_to_fmt + procedure, pass(a) :: cp_from_fmt => s_cp_coo_from_fmt + procedure, pass(a) :: mv_to_coo => s_mv_coo_to_coo + procedure, pass(a) :: mv_from_coo => s_mv_coo_from_coo + procedure, pass(a) :: mv_to_fmt => s_mv_coo_to_fmt + procedure, pass(a) :: mv_from_fmt => s_mv_coo_from_fmt + procedure, pass(a) :: fix => s_fix_coo + procedure, pass(a) :: free => s_coo_free + procedure, pass(a) :: trim => s_coo_trim + procedure, pass(a) :: s_csgetrow => s_coo_csgetrow + procedure, pass(a) :: csgetptn => s_coo_csgetptn + procedure, pass(a) :: print => s_coo_print + procedure, pass(a) :: get_fmt => s_coo_get_fmt + procedure, pass(a) :: get_nz_row => s_coo_get_nz_row + procedure, pass(a) :: sizeof => s_coo_sizeof + procedure, pass(a) :: reinit => s_coo_reinit + procedure, pass(a) :: s_coo_cp_from + generic, public :: cp_from => s_coo_cp_from + procedure, pass(a) :: s_coo_mv_from + generic, public :: mv_from => s_coo_mv_from + + end type psb_s_coo_sparse_mat + + private :: s_coo_get_nzeros, s_coo_set_nzeros, s_coo_get_diag, & + & s_coo_csmm, s_coo_csmv, s_coo_cssm, s_coo_cssv, s_coo_csnmi, & + & s_coo_csput, s_coo_reallocate_nz, s_coo_allocate_mnnz, & + & s_fix_coo, s_coo_free, s_coo_print, s_coo_get_fmt, & + & s_cp_coo_to_coo, s_cp_coo_from_coo, & + & s_cp_coo_to_fmt, s_cp_coo_from_fmt, & + & s_coo_scals, s_coo_scal, s_coo_csgetrow, s_coo_sizeof, & + & s_coo_csgetptn, s_coo_get_nz_row, s_coo_reinit,& + & s_coo_cp_from, s_coo_mv_from + + + interface + subroutine s_fix_coo_inner(nzin,dupl,ia,ja,val,nzout,info,idir) + use psb_const_mod + integer, intent(in) :: nzin,dupl + integer, intent(inout) :: ia(:), ja(:) + real(psb_spk_), intent(inout) :: val(:) + integer, intent(out) :: nzout, info + integer, intent(in), optional :: idir + end subroutine s_fix_coo_inner + end interface + + interface + subroutine s_fix_coo_impl(a,info,idir) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + end subroutine s_fix_coo_impl + end interface + + interface + subroutine s_cp_coo_to_coo_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_cp_coo_to_coo_impl + end interface + + interface + subroutine s_cp_coo_from_coo_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(out) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine s_cp_coo_from_coo_impl + end interface + + interface + subroutine s_cp_coo_to_fmt_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat, psb_s_base_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_cp_coo_to_fmt_impl + end interface + + interface + subroutine s_cp_coo_from_fmt_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat, psb_s_base_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine s_cp_coo_from_fmt_impl + end interface + + interface + subroutine s_mv_coo_to_coo_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_mv_coo_to_coo_impl + end interface + + interface + subroutine s_mv_coo_from_coo_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + end subroutine s_mv_coo_from_coo_impl + end interface + + interface + subroutine s_mv_coo_to_fmt_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat, psb_s_base_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_mv_coo_to_fmt_impl + end interface + + interface + subroutine s_mv_coo_from_fmt_impl(a,b,info) + use psb_const_mod + import psb_s_coo_sparse_mat, psb_s_base_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + end subroutine s_mv_coo_from_fmt_impl + end interface + + + interface + subroutine s_coo_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz,ia(:), ja(:),& + & imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + end subroutine s_coo_csput_impl + end interface + + interface + subroutine s_coo_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_s_coo_sparse_mat + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + end subroutine s_coo_csgetptn_impl + end interface + + interface + subroutine s_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_s_coo_sparse_mat + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + end subroutine s_coo_csgetrow_impl + end interface + + interface s_coo_cssm_impl + subroutine s_coo_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_coo_cssv_impl + subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_coo_cssm_impl + end interface + + interface s_coo_csmm_impl + subroutine s_coo_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_coo_csmv_impl + subroutine s_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_coo_csmm_impl + end interface + + + interface s_coo_csnmi_impl + function s_coo_csnmi_impl(a) result(res) + use psb_const_mod + import psb_s_coo_sparse_mat + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + end function s_coo_csnmi_impl + end interface + + +contains + + + !==================================== + ! + ! + ! + ! Data management + ! + ! + ! + ! + ! + !==================================== + + subroutine cp_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine cp_to_coo + + subroutine cp_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine cp_from_coo + + + subroutine cp_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_fmt' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine cp_to_fmt + + subroutine cp_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_fmt' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine cp_from_fmt + + + subroutine mv_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine mv_to_coo + + subroutine mv_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine mv_from_coo + + + subroutine mv_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_fmt' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine mv_to_fmt + + subroutine mv_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_fmt' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine mv_from_fmt + + subroutine s_base_mv_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_base_sparse_mat), intent(out) :: a + type(psb_s_base_sparse_mat), intent(inout) :: b + + + ! No new things here, very easy + call a%psb_base_sparse_mat%mv_from(b%psb_base_sparse_mat) + + return + + end subroutine s_base_mv_from + + subroutine s_base_cp_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_base_sparse_mat), intent(out) :: a + type(psb_s_base_sparse_mat), intent(in) :: b + + ! No new things here, very easy + call a%psb_base_sparse_mat%cp_from(b%psb_base_sparse_mat) + + return + + end subroutine s_base_cp_from + + + + subroutine csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_base_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + Integer :: err_act + character(len=20) :: name='csput' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine csput + + subroutine s_csgetrow(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_csgetrow + + + + subroutine s_csgetblk(imin,imax,a,b,info,& + & jmin,jmax,iren,append,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(in) :: imin,imax + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax + logical, intent(in), optional :: rscale,cscale + Integer :: err_act, nzin, nzout + character(len=20) :: name='csget' + logical :: appens_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(append)) then + appens_ = append + else + appens_ = .false. + endif + if (appens_) then + nzin = a%get_nzeros() + else + nzin = 0 + endif + + call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin, jmax=jmax, iren=iren, append=appens_, & + & nzin=nzin, rscale=rscale, cscale=cscale) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + if (info /= 0) 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 s_csgetblk + + + subroutine csclip(a,b,info,& + & imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_base_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer,intent(out) :: info + integer, intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb + character(len=20) :: name='csget' + logical :: rscale_, cscale_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + nzin = 0 + 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 = a%get_nrows() ! Should this be imax_ ?? + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = a%get_ncols() ! Should this be jmax_ ?? + endif + call b%allocate(mb,nb) + + call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=jmax_, append=.false., & + & nzin=nzin, rscale=rscale_, cscale=cscale_) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + + if (info /= 0) 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 csclip + + + + !==================================== + ! + ! + ! + ! Computational routines + ! + ! + ! + ! + ! + ! + !==================================== + + subroutine s_base_csmm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + Integer :: err_act + character(len=20) :: name='s_base_csmm' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_base_csmm + + subroutine s_base_csmv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + Integer :: err_act + character(len=20) :: name='s_base_csmv' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + + end subroutine s_base_csmv + + subroutine s_base_cssm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + Integer :: err_act + character(len=20) :: name='s_base_cssm' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_base_cssm + + subroutine s_base_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + Integer :: err_act + character(len=20) :: name='s_base_cssv' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_base_cssv + + subroutine s_cssm(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_spk_), intent(in), optional :: d(:) + + real(psb_spk_), allocatable :: tmp(:,:) + Integer :: err_act, nar,nac,nc, i + character(len=1) :: side_ + character(len=20) :: name='s_cssm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + nc = min(size(x,2), size(y,2)) + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(d)) then + if (present(side)) then + side_ = side + else + side_ = 'L' + end if + + if (psb_toupper(side_) == 'R') then + if (size(d,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nac,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nac,nc),stat=info) + if (info /= 0) info = 4000 + if (info == 0) then + do i=1, nac + tmp(i,1:nc) = d(i)*x(i,1:nc) + end do + end if + if (info == 0)& + & call a%base_cssm(alpha,tmp,beta,y,info,trans) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else if (psb_toupper(side_) == 'L') then + + if (size(d,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nar,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nar,nc),stat=info) + if (info /= 0) info = 4000 + if (info == 0)& + & call a%base_cssm(sone,x,szero,tmp,info,trans) + + if (info == 0)then + do i=1, nar + tmp(i,1:nc) = d(i)*tmp(i,1:nc) + end do + end if + if (info == 0)& + & call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else + info = 31 + call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_) + goto 9999 + end if + else + ! Side is ignored in this case + call a%base_cssm(alpha,x,beta,y,info,trans) + end if + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name, a_err='base_cssm') + goto 9999 + end if + + + return + 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 s_cssm + + subroutine s_cssv(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_spk_), intent(in), optional :: d(:) + + real(psb_spk_), allocatable :: tmp(:) + Integer :: err_act, nar,nac,nc, i + character(len=1) :: side_ + character(len=20) :: name='s_cssm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + nc = 1 + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(d)) then + if (present(side)) then + side_ = side + else + side_ = 'L' + end if + + if (psb_toupper(side_) == 'R') then + if (size(d,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nac,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nac),stat=info) + if (info /= 0) info = 4000 + if (info == 0) tmp(1:nac) = d(1:nac)*x(1:nac) + if (info == 0)& + & call a%base_cssm(alpha,tmp,beta,y,info,trans) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else if (psb_toupper(side_) == 'L') then + if (size(d,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/9,nar,0,0,0/)) + goto 9999 + end if + + allocate(tmp(nar),stat=info) + if (info /= 0) info = 4000 + if (info == 0)& + & call a%base_cssm(sone,x,szero,tmp,info,trans) + + if (info == 0) tmp(1:nar) = d(1:nar)*tmp(1:nar) + if (info == 0)& + & call daxpby(nar,nc,alpha,tmp,size(tmp,1),beta,y,size(y,1),info) + + if (info == 0) then + deallocate(tmp,stat=info) + if (info /= 0) info = 4000 + end if + + else + info = 31 + call psb_errpush(info,name,i_err=(/8,0,0,0,0/),a_err=side_) + goto 9999 + end if + else + ! Side is ignored in this case + call a%base_cssm(alpha,x,beta,y,info,trans) + end if + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name, a_err='base_cssm') + goto 9999 + end if + + + return + 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 s_cssv + + + subroutine s_scals(d,a,info) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='s_scals' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_scals + + + subroutine s_scal(d,a,info) + use psb_error_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='s_scal' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_scal + + + function csnmi(a) result(res) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + Integer :: err_act, info + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + res = -sone + + return + + end function csnmi + + subroutine get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_base_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name,a_err=a%get_fmt()) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + + return + + end subroutine get_diag + + + + + !==================================== + ! + ! + ! + ! Getters + ! + ! + ! + ! + ! + !==================================== + + + + function s_coo_sizeof(a) result(res) + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + 1 + res = res + psb_sizeof_sp * size(a%val) + res = res + psb_sizeof_int * size(a%ia) + res = res + psb_sizeof_int * size(a%ja) + + end function s_coo_sizeof + + + function s_coo_get_fmt(a) result(res) + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + character(len=5) :: res + res = 'COO' + end function s_coo_get_fmt + + + function s_coo_get_size(a) result(res) + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + integer :: res + res = -1 + + if (allocated(a%ia)) res = size(a%ia) + if (allocated(a%ja)) then + if (res >= 0) then + res = min(res,size(a%ja)) + else + res = size(a%ja) + end if + end if + if (allocated(a%val)) then + if (res >= 0) then + res = min(res,size(a%val)) + else + res = size(a%val) + end if + end if + end function s_coo_get_size + + + function s_coo_get_nzeros(a) result(res) + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + integer :: res + res = a%nnz + end function s_coo_get_nzeros + + + function s_coo_get_nz_row(idx,a) result(res) + use psb_const_mod + use psb_sort_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer :: res + integer :: nzin_, nza,ip,jp,i,k + + res = 0 + nza = a%get_nzeros() + if (a%is_sorted()) then + ! In this case we can do a binary search. + ip = psb_ibsrch(idx,nza,a%ia) + if (ip /= -1) return + jp = ip + do + if (ip < 2) exit + if (a%ia(ip-1) == idx) then + ip = ip -1 + else + exit + end if + end do + do + if (jp == nza) exit + if (a%ia(jp+1) == idx) then + jp = jp + 1 + else + exit + end if + end do + + res = jp - ip +1 + + else + + res = 0 + + do i=1, nza + if (a%ia(i) == idx) then + res = res + 1 + end if + end do + + end if + + end function s_coo_get_nz_row + + !==================================== + ! + ! + ! + ! Setters + ! + ! + ! + ! + ! + ! + !==================================== + + subroutine s_coo_set_nzeros(nz,a) + implicit none + integer, intent(in) :: nz + class(psb_s_coo_sparse_mat), intent(inout) :: a + + a%nnz = nz + + end subroutine s_coo_set_nzeros + + !==================================== + ! + ! + ! + ! Data management + ! + ! + ! + ! + ! + !==================================== + + + subroutine s_fix_coo(a,info,idir) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + Integer :: err_act + character(len=20) :: name='fix_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_fix_coo_impl(a,info,idir) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + + end subroutine s_fix_coo + + + subroutine s_cp_coo_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_coo_to_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_coo_to_coo + + subroutine s_cp_coo_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(out) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_coo_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_coo_from_coo + + + subroutine s_cp_coo_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_coo_to_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_coo_to_fmt + + subroutine s_cp_coo_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_coo_from_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_coo_from_fmt + + + + subroutine s_mv_coo_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_coo_to_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_coo_to_coo + + subroutine s_mv_coo_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_coo_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_coo_from_coo + + + + subroutine s_coo_cp_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(out) :: a + type(psb_s_coo_sparse_mat), intent(in) :: b + + + Integer :: err_act, info + character(len=20) :: name='cp_from' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_coo_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_coo_cp_from + + subroutine s_coo_mv_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(out) :: a + type(psb_s_coo_sparse_mat), intent(inout) :: b + + + Integer :: err_act, info + character(len=20) :: name='mv_from' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_coo_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_coo_mv_from + + + subroutine s_mv_coo_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_coo_to_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_coo_to_fmt + + subroutine s_mv_coo_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_coo_from_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_coo_from_fmt + + + + subroutine s_coo_reallocate_nz(nz,a) + use psb_error_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: nz + class(psb_s_coo_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='s_coo_reallocate_nz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + call psb_realloc(nz,a%ia,a%ja,a%val,info) + + if (info /= 0) then + call psb_errpush(4000,name) + goto 9999 + end if + + 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 s_coo_reallocate_nz + + + subroutine s_coo_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='s_coo_csput' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + call psb_erractionsave(err_act) + info = 0 + + if (nz <= 0) then + info = 10 + int_err(1)=1 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(ia) < nz) then + info = 35 + int_err(1)=2 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (size(ja) < nz) then + info = 35 + int_err(1)=3 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(val) < nz) then + info = 35 + int_err(1)=4 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (nz == 0) return + nza = a%get_nzeros() + call s_coo_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + if (info /= 0) 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 s_coo_csput + + + subroutine s_coo_csgetrow(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + call s_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + + if (info /= 0) 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 s_coo_csgetrow + + + subroutine s_coo_csgetptn(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + call s_coo_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + + if (info /= 0) 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 s_coo_csgetptn + + + subroutine s_coo_free(a) + implicit none + + class(psb_s_coo_sparse_mat), intent(inout) :: a + + if (allocated(a%ia)) deallocate(a%ia) + if (allocated(a%ja)) deallocate(a%ja) + if (allocated(a%val)) deallocate(a%val) + call a%set_null() + call a%set_nrows(0) + call a%set_ncols(0) + + return + + end subroutine s_coo_free + + subroutine s_coo_reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + + Integer :: err_act, info + character(len=20) :: name='reinit' + logical :: clear_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + + if (present(clear)) then + clear_ = clear + else + clear_ = .true. + end if + + if (a%is_bld() .or. a%is_upd()) then + ! do nothing + return + else if (a%is_asb()) then + if (clear_) a%val(:) = szero + call a%set_upd() + else + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + 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 s_coo_reinit + + + subroutine s_coo_trim(a) + use psb_realloc_mod + use psb_error_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + Integer :: err_act, info, nz + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + nz = a%get_nzeros() + if (info == 0) call psb_realloc(nz,a%ia,info) + if (info == 0) call psb_realloc(nz,a%ja,info) + if (info == 0) call psb_realloc(nz,a%val,info) + + if (info /= 0) 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 s_coo_trim + + subroutine s_coo_allocate_mnnz(m,n,a,nz) + use psb_error_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: m,n + class(psb_s_coo_sparse_mat), intent(inout) :: a + integer, intent(in), optional :: nz + Integer :: err_act, info, nz_ + character(len=20) :: name='allocate_mnz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + if (m < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) + goto 9999 + endif + if (n < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) + goto 9999 + endif + if (present(nz)) then + nz_ = nz + else + nz_ = max(7*m,7*n,1) + end if + if (nz_ < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) + goto 9999 + endif + if (info == 0) call psb_realloc(nz_,a%ia,info) + if (info == 0) call psb_realloc(nz_,a%ja,info) + if (info == 0) call psb_realloc(nz_,a%val,info) + if (info == 0) then + call a%set_nrows(m) + call a%set_ncols(n) + call a%set_nzeros(0) + call a%set_bld() + call a%set_triangle(.false.) + call a%set_unit(.false.) + call a%set_dupl(psb_dupl_def_) + end if + if (info /= 0) 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 s_coo_allocate_mnnz + + + subroutine s_coo_print(iout,a,iv,eirs,eics,head,ivr,ivc) + use psb_string_mod + implicit none + + integer, intent(in) :: iout + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in), optional :: iv(:) + integer, intent(in), optional :: eirs,eics + character(len=*), optional :: head + integer, intent(in), optional :: ivr(:), ivc(:) + + Integer :: err_act + character(len=20) :: name='s_coo_print' + logical, parameter :: debug=.false. + + character(len=80) :: frmtv + integer :: irs,ics,i,j, nmx, ni, nr, nc, nz + + if (present(eirs)) then + irs = eirs + else + irs = 0 + endif + if (present(eics)) then + ics = eics + else + ics = 0 + endif + + if (present(head)) then + write(iout,'(a)') '%%MatrixMarket matrix coordinate real general' + write(iout,'(a,a)') '% ',head + write(iout,'(a)') '%' + write(iout,'(a,a)') '% COO' + endif + + nr = a%get_nrows() + nc = a%get_ncols() + nz = a%get_nzeros() + nmx = max(nr,nc,1) + ni = floor(log10(1.0*nmx)) + 1 + + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))' + write(iout,*) nr, nc, nz + if(present(iv)) then + do j=1,a%get_nzeros() + write(iout,frmtv) iv(a%ia(j)),iv(a%ja(j)),a%val(j) + enddo + else + if (present(ivr).and..not.present(ivc)) then + do j=1,a%get_nzeros() + write(iout,frmtv) ivr(a%ia(j)),a%ja(j),a%val(j) + enddo + else if (present(ivr).and.present(ivc)) then + do j=1,a%get_nzeros() + write(iout,frmtv) ivr(a%ia(j)),ivc(a%ja(j)),a%val(j) + enddo + else if (.not.present(ivr).and.present(ivc)) then + do j=1,a%get_nzeros() + write(iout,frmtv) a%ia(j),ivc(a%ja(j)),a%val(j) + enddo + else if (.not.present(ivr).and..not.present(ivc)) then + do j=1,a%get_nzeros() + write(iout,frmtv) a%ia(j),a%ja(j),a%val(j) + enddo + endif + endif + + end subroutine s_coo_print + + + + + !==================================== + ! + ! + ! + ! Computational routines + ! + ! + ! + ! + ! + ! + !==================================== + + subroutine s_coo_csmv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nac, nar + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_csmv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + + call s_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 s_coo_csmv + + subroutine s_coo_csmm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac + real(psb_spk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + call s_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 s_coo_csmm + + + subroutine s_coo_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nar, nac + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_cssv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call s_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + + 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 s_coo_cssv + + + + subroutine s_coo_cssm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc, nar, nac + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + nar = a%get_nrows() + nac = a%get_ncols() + if (size(x,1) < nac) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nac,0,0,0/)) + goto 9999 + end if + if (size(y,1) < nar) then + info = 36 + call psb_errpush(info,name,i_err=(/3,nar,0,0,0/)) + goto 9999 + end if + + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call s_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + 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 s_coo_cssm + + function s_coo_csnmi(a) result(res) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + + res = s_coo_csnmi_impl(a) + + return + + end function s_coo_csnmi + + subroutine s_coo_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + d(:) = szero + + do i=1,a%get_nzeros() + j=a%ia(i) + if ((j==a%ja(i)) .and.(j <= mnm ) .and.(j>0)) then + d(j) = a%val(i) + endif + enddo + + 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 s_coo_get_diag + + subroutine s_coo_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + m = a%get_nrows() + if (size(d) < m) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + + do i=1,a%get_nzeros() + j = a%ia(i) + a%val(i) = a%val(i) * d(j) + enddo + + 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 s_coo_scal + + subroutine s_coo_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + do i=1,a%get_nzeros() + a%val(i) = a%val(i) * d + enddo + + 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 s_coo_scals + + +end module psb_s_base_mat_mod + + + diff --git a/base/modules/psb_s_csr_mat_mod.f03 b/base/modules/psb_s_csr_mat_mod.f03 new file mode 100644 index 00000000..fc546a5a --- /dev/null +++ b/base/modules/psb_s_csr_mat_mod.f03 @@ -0,0 +1,1604 @@ +module psb_s_csr_mat_mod + + use psb_s_base_mat_mod + + type, extends(psb_s_base_sparse_mat) :: psb_s_csr_sparse_mat + + integer, allocatable :: irp(:), ja(:) + real(psb_spk_), allocatable :: val(:) + + contains + procedure, pass(a) :: get_nzeros => s_csr_get_nzeros + procedure, pass(a) :: get_fmt => s_csr_get_fmt + procedure, pass(a) :: get_diag => s_csr_get_diag + procedure, pass(a) :: s_base_csmm => s_csr_csmm + procedure, pass(a) :: s_base_csmv => s_csr_csmv + procedure, pass(a) :: s_base_cssm => s_csr_cssm + procedure, pass(a) :: s_base_cssv => s_csr_cssv + procedure, pass(a) :: s_scals => s_csr_scals + procedure, pass(a) :: s_scal => s_csr_scal + procedure, pass(a) :: csnmi => s_csr_csnmi + procedure, pass(a) :: reallocate_nz => s_csr_reallocate_nz + procedure, pass(a) :: csput => s_csr_csput + procedure, pass(a) :: allocate_mnnz => s_csr_allocate_mnnz + procedure, pass(a) :: cp_to_coo => s_cp_csr_to_coo + procedure, pass(a) :: cp_from_coo => s_cp_csr_from_coo + procedure, pass(a) :: cp_to_fmt => s_cp_csr_to_fmt + procedure, pass(a) :: cp_from_fmt => s_cp_csr_from_fmt + procedure, pass(a) :: mv_to_coo => s_mv_csr_to_coo + procedure, pass(a) :: mv_from_coo => s_mv_csr_from_coo + procedure, pass(a) :: mv_to_fmt => s_mv_csr_to_fmt + procedure, pass(a) :: mv_from_fmt => s_mv_csr_from_fmt + procedure, pass(a) :: csgetptn => s_csr_csgetptn + procedure, pass(a) :: s_csgetrow => s_csr_csgetrow + procedure, pass(a) :: get_nz_row => s_csr_get_nz_row + procedure, pass(a) :: get_size => s_csr_get_size + procedure, pass(a) :: free => s_csr_free + procedure, pass(a) :: trim => s_csr_trim + procedure, pass(a) :: print => s_csr_print + procedure, pass(a) :: sizeof => s_csr_sizeof + procedure, pass(a) :: reinit => s_csr_reinit + procedure, pass(a) :: s_csr_cp_from + generic, public :: cp_from => s_csr_cp_from + procedure, pass(a) :: s_csr_mv_from + generic, public :: mv_from => s_csr_mv_from + + end type psb_s_csr_sparse_mat + + private :: s_csr_get_nzeros, s_csr_csmm, s_csr_csmv, s_csr_cssm, s_csr_cssv, & + & s_csr_csput, s_csr_reallocate_nz, s_csr_allocate_mnnz, & + & s_csr_free, s_csr_print, s_csr_get_fmt, s_csr_csnmi, get_diag, & + & s_cp_csr_to_coo, s_cp_csr_from_coo, & + & s_mv_csr_to_coo, s_mv_csr_from_coo, & + & s_cp_csr_to_fmt, s_cp_csr_from_fmt, & + & s_mv_csr_to_fmt, s_mv_csr_from_fmt, & + & s_csr_scals, s_csr_scal, s_csr_trim, s_csr_csgetrow, s_csr_get_size, & + & s_csr_sizeof, s_csr_csgetptn, s_csr_get_nz_row, s_csr_reinit +!!$, & +!!$ & s_csr_mv_from, s_csr_mv_from + + + interface + subroutine s_cp_csr_to_fmt_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_cp_csr_to_fmt_impl + end interface + + interface + subroutine s_cp_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine s_cp_csr_from_fmt_impl + end interface + + + interface + subroutine s_cp_csr_to_coo_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_cp_csr_to_coo_impl + end interface + + interface + subroutine s_cp_csr_from_coo_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine s_cp_csr_from_coo_impl + end interface + + interface + subroutine s_mv_csr_to_fmt_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_mv_csr_to_fmt_impl + end interface + + interface + subroutine s_mv_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + end subroutine s_mv_csr_from_fmt_impl + end interface + + + interface + subroutine s_mv_csr_to_coo_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine s_mv_csr_to_coo_impl + end interface + + interface + subroutine s_mv_csr_from_coo_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + end subroutine s_mv_csr_from_coo_impl + end interface + + interface + subroutine s_csr_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + end subroutine s_csr_csput_impl + end interface + + interface + subroutine s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_s_csr_sparse_mat + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + end subroutine s_csr_csgetptn_impl + end interface + + interface + subroutine s_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + use psb_const_mod + import psb_s_csr_sparse_mat + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + end subroutine s_csr_csgetrow_impl + end interface + + interface s_csr_cssm_impl + subroutine s_csr_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_csr_cssv_impl + subroutine s_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_csr_cssm_impl + end interface + + interface s_csr_csmm_impl + subroutine s_csr_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_csr_csmv_impl + subroutine s_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine s_csr_csmm_impl + end interface + + interface s_csr_csnmi_impl + function s_csr_csnmi_impl(a) result(res) + use psb_const_mod + import psb_s_csr_sparse_mat + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + end function s_csr_csnmi_impl + end interface + + + +contains + + !===================================== + ! + ! + ! + ! Getters + ! + ! + ! + ! + ! + !===================================== + + + function s_csr_sizeof(a) result(res) + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + psb_sizeof_dp * size(a%val) + res = res + psb_sizeof_int * size(a%irp) + res = res + psb_sizeof_int * size(a%ja) + + end function s_csr_sizeof + + function s_csr_get_fmt(a) result(res) + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + character(len=5) :: res + res = 'CSR' + end function s_csr_get_fmt + + function s_csr_get_nzeros(a) result(res) + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + integer :: res + res = a%irp(a%get_nrows()+1)-1 + end function s_csr_get_nzeros + + function s_csr_get_size(a) result(res) + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + integer :: res + + res = -1 + + if (allocated(a%ja)) then + if (res >= 0) then + res = min(res,size(a%ja)) + else + res = size(a%ja) + end if + end if + if (allocated(a%val)) then + if (res >= 0) then + res = min(res,size(a%val)) + else + res = size(a%val) + end if + end if + + end function s_csr_get_size + + + + function s_csr_get_nz_row(idx,a) result(res) + use psb_const_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer :: res + + res = 0 + + if ((1<=idx).and.(idx<=a%get_nrows())) then + res = a%irp(idx+1)-a%irp(idx) + end if + + end function s_csr_get_nz_row + + + + !===================================== + ! + ! + ! + ! Data management + ! + ! + ! + ! + ! + !===================================== + + + subroutine s_csr_reallocate_nz(nz,a) + use psb_error_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: nz + class(psb_s_csr_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='s_csr_reallocate_nz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + call psb_realloc(nz,a%ja,info) + if (info == 0) call psb_realloc(nz,a%val,info) + if (info == 0) call psb_realloc(& + & max(nz,a%get_nrows()+1,a%get_ncols()+1),a%irp,info) + if (info /= 0) then + call psb_errpush(4000,name) + goto 9999 + end if + + 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 s_csr_reallocate_nz + + subroutine s_csr_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='s_csr_csput' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + call psb_erractionsave(err_act) + info = 0 + + if (nz <= 0) then + info = 10 + int_err(1)=1 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(ia) < nz) then + info = 35 + int_err(1)=2 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (size(ja) < nz) then + info = 35 + int_err(1)=3 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(val) < nz) then + info = 35 + int_err(1)=4 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (nz == 0) return + + call s_csr_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + if (info /= 0) 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 s_csr_csput + + subroutine s_csr_csgetptn(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + call s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + + if (info /= 0) 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 s_csr_csgetptn + + + subroutine s_csr_csgetrow(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + call s_csr_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + + if (info /= 0) 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 s_csr_csgetrow + + + subroutine s_csr_csgetblk(imin,imax,a,b,info,& + & jmin,jmax,iren,append,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(in) :: imin,imax + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax + logical, intent(in), optional :: rscale,cscale + Integer :: err_act, nzin, nzout + character(len=20) :: name='csget' + logical :: appens_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(append)) then + appens_ = append + else + appens_ = .false. + endif + if (appens_) then + nzin = a%get_nzeros() + else + nzin = 0 + endif + + call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin, jmax=jmax, iren=iren, append=appens_, & + & nzin=nzin, rscale=rscale, cscale=cscale) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + if (info /= 0) 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 s_csr_csgetblk + + + subroutine s_csr_csclip(a,b,info,& + & imin,imax,jmin,jmax,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer,intent(out) :: info + integer, intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act, nzin, nzout, imin_, imax_, jmin_, jmax_, mb,nb + character(len=20) :: name='csget' + logical :: rscale_, cscale_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + nzin = 0 + 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 = a%get_nrows() ! Should this be imax_ ?? + endif + if (cscale_) then + nb = jmax_ - jmin_ +1 + else + nb = a%get_ncols() ! Should this be jmax_ ?? + endif + call b%allocate(mb,nb) + + call a%csget(imin_,imax_,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin_, jmax=jmax_, append=.false., & + & nzin=nzin, rscale=rscale_, cscale=cscale_) + + if (info /= 0) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + + if (info /= 0) 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 s_csr_csclip + + + subroutine s_csr_free(a) + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + + if (allocated(a%irp)) deallocate(a%irp) + if (allocated(a%ja)) deallocate(a%ja) + if (allocated(a%val)) deallocate(a%val) + call a%set_null() + call a%set_nrows(0) + call a%set_ncols(0) + + return + + end subroutine s_csr_free + + subroutine s_csr_reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + + Integer :: err_act, info + character(len=20) :: name='reinit' + logical :: clear_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + + if (present(clear)) then + clear_ = clear + else + clear_ = .true. + end if + + if (a%is_bld() .or. a%is_upd()) then + ! do nothing + return + else if (a%is_asb()) then + if (clear_) a%val(:) = szero + call a%set_upd() + else + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + 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 s_csr_reinit + + + subroutine s_csr_trim(a) + use psb_realloc_mod + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + Integer :: err_act, info, nz, m + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + m = a%get_nrows() + nz = a%get_nzeros() + if (info == 0) call psb_realloc(m+1,a%irp,info) + if (info == 0) call psb_realloc(nz,a%ja,info) + if (info == 0) call psb_realloc(nz,a%val,info) + + if (info /= 0) 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 s_csr_trim + + + subroutine s_cp_csr_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_csr_to_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_csr_to_coo + + subroutine s_cp_csr_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_csr_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_csr_from_coo + + + subroutine s_cp_csr_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_fmt' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_csr_to_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_csr_to_fmt + + subroutine s_cp_csr_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_fmt' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_cp_csr_from_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_cp_csr_from_fmt + + + subroutine s_mv_csr_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_csr_to_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_csr_to_coo + + subroutine s_mv_csr_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_csr_from_coo_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_csr_from_coo + + + subroutine s_mv_csr_to_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_fmt' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_csr_to_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_csr_to_fmt + + subroutine s_mv_csr_from_fmt(a,b,info) + use psb_error_mod + use psb_realloc_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_fmt' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call s_mv_csr_from_fmt_impl(a,b,info) + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_mv_csr_from_fmt + + + subroutine s_csr_allocate_mnnz(m,n,a,nz) + use psb_error_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: m,n + class(psb_s_csr_sparse_mat), intent(inout) :: a + integer, intent(in), optional :: nz + Integer :: err_act, info, nz_ + character(len=20) :: name='allocate_mnz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + if (m < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) + goto 9999 + endif + if (n < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) + goto 9999 + endif + if (present(nz)) then + nz_ = nz + else + nz_ = max(7*m,7*n,1) + end if + if (nz_ < 0) then + info = 10 + call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) + goto 9999 + endif + + if (info == 0) call psb_realloc(m+1,a%irp,info) + if (info == 0) call psb_realloc(nz_,a%ja,info) + if (info == 0) call psb_realloc(nz_,a%val,info) + if (info == 0) then + a%irp=0 + call a%set_nrows(m) + call a%set_ncols(n) + call a%set_bld() + call a%set_triangle(.false.) + call a%set_unit(.false.) + end if + + 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 s_csr_allocate_mnnz + + + subroutine s_csr_print(iout,a,iv,eirs,eics,head,ivr,ivc) + use psb_string_mod + implicit none + + integer, intent(in) :: iout + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in), optional :: iv(:) + integer, intent(in), optional :: eirs,eics + character(len=*), optional :: head + integer, intent(in), optional :: ivr(:), ivc(:) + + Integer :: err_act + character(len=20) :: name='s_csr_print' + logical, parameter :: debug=.false. + + character(len=80) :: frmtv + integer :: irs,ics,i,j, nmx, ni, nr, nc, nz + + if (present(eirs)) then + irs = eirs + else + irs = 0 + endif + if (present(eics)) then + ics = eics + else + ics = 0 + endif + + if (present(head)) then + write(iout,'(a)') '%%MatrixMarket matrix coordinate real general' + write(iout,'(a,a)') '% ',head + write(iout,'(a)') '%' + write(iout,'(a,a)') '% COO' + endif + + nr = a%get_nrows() + nc = a%get_ncols() + nz = a%get_nzeros() + nmx = max(nr,nc,1) + ni = floor(log10(1.0*nmx)) + 1 + + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))' + write(iout,*) nr, nc, nz + if(present(iv)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) iv(i),iv(a%ja(j)),a%val(j) + end do + enddo + else + if (present(ivr).and..not.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) ivr(i),(a%ja(j)),a%val(j) + end do + enddo + else if (present(ivr).and.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) ivr(i),ivc(a%ja(j)),a%val(j) + end do + enddo + else if (.not.present(ivr).and.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) (i),ivc(a%ja(j)),a%val(j) + end do + enddo + else if (.not.present(ivr).and..not.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) (i),(a%ja(j)),a%val(j) + end do + enddo + endif + endif + + end subroutine s_csr_print + + + subroutine s_csr_cp_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(out) :: a + type(psb_s_csr_sparse_mat), intent(in) :: b + + + Integer :: err_act, info + character(len=20) :: name='cp_from' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + info = 0 + + call a%allocate(b%get_nrows(),b%get_ncols(),b%get_nzeros()) + call a%psb_s_base_sparse_mat%cp_from(b%psb_s_base_sparse_mat) + a%irp = b%irp + a%ja = b%ja + a%val = b%val + + if (info /= 0) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_csr_cp_from + + subroutine s_csr_mv_from(a,b) + use psb_error_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(out) :: a + type(psb_s_csr_sparse_mat), intent(inout) :: b + + + Integer :: err_act, info + character(len=20) :: name='mv_from' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call a%psb_s_base_sparse_mat%mv_from(b%psb_s_base_sparse_mat) + call move_alloc(b%irp, a%irp) + call move_alloc(b%ja, a%ja) + call move_alloc(b%val, a%val) + call b%free() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine s_csr_mv_from + + + + !===================================== + ! + ! + ! + ! Computational routines + ! + ! + ! + ! + ! + ! + !===================================== + + + subroutine s_csr_csmv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_csmv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call s_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 s_csr_csmv + + subroutine s_csr_csmm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + + + call s_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 s_csr_csmm + + + subroutine s_csr_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_cssv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call s_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + + 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 s_csr_cssv + + + + subroutine s_csr_cssm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call s_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + 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 s_csr_cssm + + function s_csr_csnmi(a) result(res) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + + res = s_csr_csnmi_impl(a) + + return + + end function s_csr_csnmi + + subroutine s_csr_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act, mnm, i, j, k + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + + + do i=1, mnm + do k=a%irp(i),a%irp(i+1)-1 + j=a%ja(k) + if ((j==i) .and.(j <= mnm )) then + d(i) = a%val(k) + endif + enddo + end do + do i=mnm+1,size(d) + d(i) = szero + end do + 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 s_csr_get_diag + + + subroutine s_csr_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + m = a%get_nrows() + if (size(d) < m) then + info=35 + call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + goto 9999 + end if + + do i=1, m + do j = a%irp(i), a%irp(i+1) -1 + a%val(j) = a%val(j) * d(i) + end do + enddo + + 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 s_csr_scal + + + subroutine s_csr_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_csr_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act,mnm, i, j, m + character(len=20) :: name='scal' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + + do i=1,a%get_nzeros() + a%val(i) = a%val(i) * d + enddo + + 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 s_csr_scals + + + +end module psb_s_csr_mat_mod diff --git a/base/modules/psb_s_mat_mod.f03 b/base/modules/psb_s_mat_mod.f03 new file mode 100644 index 00000000..6d83169b --- /dev/null +++ b/base/modules/psb_s_mat_mod.f03 @@ -0,0 +1,1924 @@ +module psb_s_mat_mod + + use psb_s_base_mat_mod + use psb_s_csr_mat_mod + + type :: psb_s_sparse_mat + + class(psb_s_base_sparse_mat), allocatable :: a + + contains + ! Setters + procedure, pass(a) :: set_nrows + procedure, pass(a) :: set_ncols + procedure, pass(a) :: set_dupl + procedure, pass(a) :: set_state + procedure, pass(a) :: set_null + procedure, pass(a) :: set_bld + procedure, pass(a) :: set_upd + procedure, pass(a) :: set_asb + procedure, pass(a) :: set_sorted + procedure, pass(a) :: set_upper + procedure, pass(a) :: set_lower + procedure, pass(a) :: set_triangle + procedure, pass(a) :: set_unit + ! Getters + procedure, pass(a) :: get_nrows + procedure, pass(a) :: get_ncols + procedure, pass(a) :: get_nzeros + procedure, pass(a) :: get_nz_row + procedure, pass(a) :: get_size + procedure, pass(a) :: get_state + procedure, pass(a) :: get_dupl + procedure, pass(a) :: is_null + procedure, pass(a) :: is_bld + procedure, pass(a) :: is_upd + procedure, pass(a) :: is_asb + procedure, pass(a) :: is_sorted + procedure, pass(a) :: is_upper + procedure, pass(a) :: is_lower + procedure, pass(a) :: is_triangle + procedure, pass(a) :: is_unit + procedure, pass(a) :: get_fmt => sparse_get_fmt + procedure, pass(a) :: sizeof => s_sizeof + + + ! Memory/data management + procedure, pass(a) :: csall + procedure, pass(a) :: free + procedure, pass(a) :: trim + procedure, pass(a) :: csput + procedure, pass(a) :: s_csgetptn + procedure, pass(a) :: s_csgetrow + procedure, pass(a) :: s_csgetblk + generic, public :: csget => s_csgetptn, s_csgetrow, s_csgetblk + procedure, pass(a) :: csclip + procedure, pass(a) :: reall => reallocate_nz + procedure, pass(a) :: get_neigh + procedure, pass(a) :: s_cscnv + procedure, pass(a) :: s_cscnv_ip + generic, public :: cscnv => s_cscnv, s_cscnv_ip + procedure, pass(a) :: reinit + procedure, pass(a) :: print => sparse_print + procedure, pass(a) :: s_mv_from + generic, public :: mv_from => s_mv_from + procedure, pass(a) :: s_cp_from + generic, public :: cp_from => s_cp_from + + + ! Computational routines + procedure, pass(a) :: get_diag + procedure, pass(a) :: csnmi + procedure, pass(a) :: s_csmv + procedure, pass(a) :: s_csmm + generic, public :: csmm => s_csmm, s_csmv + procedure, pass(a) :: s_scals + procedure, pass(a) :: s_scal + generic, public :: scal => s_scals, s_scal + procedure, pass(a) :: s_cssv + procedure, pass(a) :: s_cssm + generic, public :: cssm => s_cssm, s_cssv + + end type psb_s_sparse_mat + + private :: get_nrows, get_ncols, get_nzeros, get_size, & + & get_state, get_dupl, is_null, is_bld, is_upd, & + & is_asb, is_sorted, is_upper, is_lower, is_triangle, & + & is_unit, get_neigh, csall, csput, s_csgetrow,& + & s_csgetblk, csclip, s_cscnv, s_cscnv_ip, & + & reallocate_nz, free, trim, & + & sparse_print, reinit, & + & set_nrows, set_ncols, set_dupl, & + & set_state, set_null, set_bld, & + & set_upd, set_asb, set_sorted, & + & set_upper, set_lower, set_triangle, & + & set_unit, get_diag, get_nz_row, s_csgetptn, & + & s_mv_from, s_cp_from + + interface psb_sizeof + module procedure s_sizeof + end interface + + interface psb_move_alloc + module procedure s_sparse_mat_move + end interface + + interface psb_clone + module procedure s_sparse_mat_clone + end interface + + interface psb_csmm + module procedure s_csmm, s_csmv + end interface + + interface psb_cssm + module procedure s_cssm, s_cssv + end interface + + interface psb_csnmi + module procedure csnmi + end interface + + interface psb_scal + module procedure s_scals, s_scal + end interface + +contains + + + !===================================== + ! + ! + ! + ! Getters + ! + ! + ! + ! + ! + !===================================== + + + function s_sizeof(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + + res = 0 + if (allocated(a%a)) then + res = a%a%sizeof() + end if + + end function s_sizeof + + + + function sparse_get_fmt(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + character(len=5) :: res + + if (allocated(a%a)) then + res = a%a%get_fmt() + else + res = 'NULL' + end if + + end function sparse_get_fmt + + + + function get_dupl(a) result(res) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_dupl() + else + res = psb_invalid_ + end if + end function get_dupl + + + function get_state(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_state() + else + res = psb_spmat_null_ + end if + end function get_state + + function get_nrows(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_nrows() + else + res = 0 + end if + + end function get_nrows + + function get_ncols(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_ncols() + else + res = 0 + end if + + end function get_ncols + + function is_triangle(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_triangle() + else + res = .false. + end if + + end function is_triangle + + function is_unit(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_unit() + else + res = .false. + end if + + end function is_unit + + function is_upper(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_upper() + else + res = .false. + end if + + end function is_upper + + function is_lower(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = .not. a%a%is_upper() + else + res = .false. + end if + + end function is_lower + + function is_null(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_null() + else + res = .true. + end if + + end function is_null + + function is_bld(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_bld() + else + res = .false. + end if + + end function is_bld + + function is_upd(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_upd() + else + res = .false. + end if + + end function is_upd + + function is_asb(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_asb() + else + res = .false. + end if + + end function is_asb + + function is_sorted(a) result(res) + implicit none + class(psb_s_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_sorted() + else + res = .false. + end if + + end function is_sorted + + + + function get_nzeros(a) result(res) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%get_nzeros() + + 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 function get_nzeros + + function get_size(a) result(res) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act, info + character(len=20) :: name='get_size' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%get_size() + + 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 function get_size + + + function get_nz_row(idx,a) result(res) + use psb_error_mod + implicit none + integer, intent(in) :: idx + class(psb_s_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act + + res = 0 + + if (allocated(a%a)) res = a%a%get_nz_row(idx) + + end function get_nz_row + + + + !===================================== + ! + ! + ! + ! Setters + ! + ! + ! + ! + ! + ! + !===================================== + + + subroutine set_nrows(m,a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + integer, intent(in) :: m + Integer :: err_act, info + character(len=20) :: name='set_nrows' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_nrows(m) + + 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 set_nrows + + subroutine set_ncols(n,a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + call a%a%set_ncols(n) + + 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 set_ncols + + + subroutine set_state(n,a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + call a%a%set_state(n) + + 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 set_state + + + subroutine set_dupl(n,a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_dupl(n) + + 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 set_dupl + + subroutine set_null(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_null() + + 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 set_null + + subroutine set_bld(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_bld() + + 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 set_bld + + subroutine set_upd(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_upd() + + 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 set_upd + + subroutine set_asb(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_asb() + + 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 set_asb + + subroutine set_sorted(a,val) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_sorted(val) + + 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 set_sorted + + subroutine set_triangle(a,val) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_triangle(val) + + 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 set_triangle + + subroutine set_unit(a,val) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_unit(val) + + 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 set_unit + + subroutine set_lower(a,val) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_lower(val) + + 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 set_lower + + subroutine set_upper(a,val) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + Integer :: err_act, info + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%set_upper(val) + + 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 set_upper + + + !===================================== + ! + ! + ! + ! Data management + ! + ! + ! + ! + ! + !===================================== + + + subroutine sparse_print(iout,a,iv,eirs,eics,head,ivr,ivc) + use psb_error_mod + implicit none + + integer, intent(in) :: iout + class(psb_s_sparse_mat), intent(in) :: a + integer, intent(in), optional :: iv(:) + integer, intent(in), optional :: eirs,eics + character(len=*), optional :: head + integer, intent(in), optional :: ivr(:), ivc(:) + + Integer :: err_act, info + character(len=20) :: name='sparse_print' + logical, parameter :: debug=.false. + + info = 0 + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%print(iout,iv,eirs,eics,head,ivr,ivc) + + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine sparse_print + + + + subroutine get_neigh(a,idx,neigh,n,info,lev) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer, intent(out) :: n + integer, allocatable, intent(out) :: neigh(:) + integer, intent(out) :: info + integer, optional, intent(in) :: lev + + Integer :: err_act + character(len=20) :: name='get_neigh' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%get_neigh(idx,neigh,n,info,lev) + + if (info /= 0) 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 get_neigh + + + subroutine csall(nr,nc,a,info,nz) + use psb_s_base_mat_mod + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(out) :: a + integer, intent(in) :: nr,nc + integer, intent(out) :: info + integer, intent(in), optional :: nz + + Integer :: err_act + character(len=20) :: name='csall' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + info = 0 + allocate(psb_s_coo_sparse_mat :: a%a, stat=info) + if (info /= 0) then + info = 4000 + call psb_errpush(info, name) + goto 9999 + end if + call a%a%allocate(nr,nc,nz) + call a%set_bld() + + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine csall + + subroutine reallocate_nz(nz,a) + use psb_error_mod + implicit none + integer, intent(in) :: nz + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='reallocate_nz' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%reallocate(nz) + + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine reallocate_nz + + subroutine free(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='free' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%free() + deallocate(a%a) + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine free + + subroutine trim(a) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + Integer :: err_act, info + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%trim() + + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end subroutine trim + + + subroutine csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_s_base_mat_mod + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + Integer :: err_act + character(len=20) :: name='csput' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.a%is_bld()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call a%a%csput(nz,ia,ja,val,imin,imax,jmin,jmax,info,gtl) + if (info /= 0) 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 csput + + subroutine s_csgetptn(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod + implicit none + + class(psb_s_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call a%a%csget(imin,imax,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + if (info /= 0) 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 s_csgetptn + + subroutine s_csgetrow(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod + implicit none + + class(psb_s_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call a%a%csget(imin,imax,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + if (info /= 0) 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 s_csgetrow + + + + subroutine s_csgetblk(imin,imax,a,b,info,& + & jmin,jmax,iren,append,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_s_base_mat_mod + implicit none + + class(psb_s_sparse_mat), intent(in) :: a + class(psb_s_sparse_mat), intent(out) :: b + integer, intent(in) :: imin,imax + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + type(psb_s_coo_sparse_mat), allocatable :: acoo + + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + + if (info == 0) call a%a%csget(imin,imax,acoo,info,& + & jmin,jmax,iren,append,rscale,cscale) + if (info == 0) call move_alloc(acoo,b%a) + if (info /= 0) 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 s_csgetblk + + + + subroutine csclip(a,b,info,& + & 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 + implicit none + + class(psb_s_sparse_mat), intent(in) :: a + class(psb_s_sparse_mat), intent(out) :: b + integer,intent(out) :: info + integer, intent(in), optional :: imin,imax,jmin,jmax + logical, intent(in), optional :: rscale,cscale + + Integer :: err_act + character(len=20) :: name='csclip' + logical, parameter :: debug=.false. + type(psb_s_coo_sparse_mat), allocatable :: acoo + + info = 0 + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + allocate(acoo,stat=info) + if (info == 0) call a%a%csclip(acoo,info,& + & imin,imax,jmin,jmax,rscale,cscale) + if (info == 0) call move_alloc(acoo,b%a) + if (info /= 0) 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 csclip + + + + subroutine s_cscnv(a,b,info,type,mold,upd,dupl) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + class(psb_s_sparse_mat), intent(out) :: b + integer, intent(out) :: info + integer,optional, intent(in) :: dupl, upd + character(len=*), optional, intent(in) :: type + class(psb_s_base_sparse_mat), intent(in), optional :: mold + + + class(psb_s_base_sparse_mat), allocatable :: altmp + Integer :: err_act + character(len=20) :: name='cscnv' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(dupl)) then + call b%set_dupl(dupl) + else if (a%is_bld()) then + ! Does this make sense at all?? Who knows.. + call b%set_dupl(psb_dupl_def_) + end if + + if (count( (/present(mold),present(type) /)) > 1) then + info = 583 + call psb_errpush(info,name,a_err='TYPE, MOLD') + goto 9999 + end if + + if (present(mold)) then + + allocate(altmp, source=mold,stat=info) + + else if (present(type)) then + + select case (psb_toupper(type)) + case ('CSR') + allocate(psb_s_csr_sparse_mat :: altmp, stat=info) + case ('COO') + allocate(psb_s_coo_sparse_mat :: altmp, stat=info) + case default + info = 136 + call psb_errpush(info,name,a_err=type) + goto 9999 + end select + else + allocate(psb_s_csr_sparse_mat :: altmp, stat=info) + end if + + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + + if (debug) write(0,*) 'Converting from ',& + & a%get_fmt(),' to ',altmp%get_fmt() + + call altmp%cp_from_fmt(a%a, info) + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err="mv_from") + goto 9999 + end if + + call move_alloc(altmp,b%a) + call b%set_asb() + call b%trim() + 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 s_cscnv + + + subroutine s_cscnv_ip(a,info,type,mold,dupl) + use psb_error_mod + use psb_string_mod + implicit none + + class(psb_s_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer,optional, intent(in) :: dupl + character(len=*), optional, intent(in) :: type + class(psb_s_base_sparse_mat), intent(in), optional :: mold + + + class(psb_s_base_sparse_mat), allocatable :: altmp + Integer :: err_act + character(len=20) :: name='cscnv_ip' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + if (present(dupl)) then + call a%set_dupl(dupl) + else if (a%is_bld()) then + call a%set_dupl(psb_dupl_def_) + end if + + if (count( (/present(mold),present(type) /)) > 1) then + info = 583 + call psb_errpush(info,name,a_err='TYPE, MOLD') + goto 9999 + end if + + if (present(mold)) then + + allocate(altmp, source=mold,stat=info) + + else if (present(type)) then + + select case (psb_toupper(type)) + case ('CSR') + allocate(psb_s_csr_sparse_mat :: altmp, stat=info) + case ('COO') + allocate(psb_s_coo_sparse_mat :: altmp, stat=info) + case default + info = 136 + call psb_errpush(info,name,a_err=type) + goto 9999 + end select + else + allocate(psb_s_csr_sparse_mat :: altmp, stat=info) + end if + + if (info /= 0) then + info = 4000 + call psb_errpush(info,name) + goto 9999 + end if + + if (debug) write(0,*) 'Converting in-place from ',& + & a%get_fmt(),' to ',altmp%get_fmt() + + call altmp%mv_from_fmt(a%a, info) + + if (info /= 0) then + info = 4010 + call psb_errpush(info,name,a_err="mv_from") + goto 9999 + end if + + call move_alloc(altmp,a%a) + call a%set_asb() + call a%trim() + 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 s_cscnv_ip + + subroutine s_mv_from(a,b) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_sparse_mat), intent(out) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer :: info + + allocate(a%a,source=b, stat=info) + call a%a%mv_from_fmt(b,info) + + return + end subroutine s_mv_from + + subroutine s_cp_from(a,b) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_sparse_mat), intent(out) :: a + class(psb_s_base_sparse_mat), intent(inout), allocatable :: b + Integer :: err_act, info + character(len=20) :: name='clone' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + allocate(a%a,source=b,stat=info) + if (info /= 0) info = 4000 + if (info == 0) call a%a%cp_from_fmt(b, info) + if (info /= 0) 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 s_cp_from + + subroutine s_sparse_mat_move(a,b,info) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + class(psb_s_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='move_alloc' + logical, parameter :: debug=.false. + + info = 0 + call move_alloc(a%a,b%a) + + return + end subroutine s_sparse_mat_move + + subroutine s_sparse_mat_clone(a,b,info) + use psb_error_mod + use psb_string_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + class(psb_s_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='clone' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + allocate(b%a,source=a%a,stat=info) + if (info /= 0) info = 4000 + if (info == 0) call b%a%cp_from_fmt(a%a, info) + if (info /= 0) 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 s_sparse_mat_clone + + + subroutine reinit(a,clear) + use psb_error_mod + implicit none + + class(psb_s_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: clear + Integer :: err_act, info + character(len=20) :: name='reinit' + + call psb_erractionsave(err_act) + if (a%is_null()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%reinit(clear) + + 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 reinit + + + !===================================== + ! + ! + ! + ! Computational routines + ! + ! + ! + ! + ! + ! + !===================================== + + + subroutine s_csmm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + Integer :: err_act + character(len=20) :: name='psb_csmm' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%csmm(alpha,x,beta,y,info,trans) + if (info /= 0) 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 s_csmm + + subroutine s_csmv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + Integer :: err_act + character(len=20) :: name='psb_csmv' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%csmm(alpha,x,beta,y,info,trans) + if (info /= 0) 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 s_csmv + + subroutine s_cssm(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_spk_), intent(in), optional :: d(:) + Integer :: err_act + character(len=20) :: name='psb_cssm' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%cssm(alpha,x,beta,y,info,trans,side,d) + if (info /= 0) 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 s_cssm + + subroutine s_cssv(alpha,a,x,beta,y,info,trans,side,d) + use psb_error_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans, side + real(psb_spk_), intent(in), optional :: d(:) + Integer :: err_act + character(len=20) :: name='psb_cssv' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%cssm(alpha,x,beta,y,info,trans,side,d) + + if (info /= 0) 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 s_cssv + + + function csnmi(a) result(res) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + Integer :: err_act, info + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%csnmi() + + + return + +9999 continue + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + end function csnmi + + + + subroutine get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%get_diag(d,info) + if (info /= 0) 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 get_diag + + subroutine s_scal(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d(:) + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%scal(d,info) + if (info /= 0) 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 s_scal + + subroutine s_scals(d,a,info) + use psb_error_mod + use psb_const_mod + implicit none + class(psb_s_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: d + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='csnmi' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%scal(d,info) + if (info /= 0) 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 s_scals + + +end module psb_s_mat_mod diff --git a/base/serial/f03/Makefile b/base/serial/f03/Makefile index c2797eb2..b755527a 100644 --- a/base/serial/f03/Makefile +++ b/base/serial/f03/Makefile @@ -3,7 +3,7 @@ include ../../../Make.inc # # The object files # -FOBJS = psbn_d_csr_impl.o psbn_d_coo_impl.o +FOBJS = psb_s_csr_impl.o psb_s_coo_impl.o psb_d_csr_impl.o psb_d_coo_impl.o OBJS=$(FOBJS) diff --git a/base/serial/f03/psbn_d_coo_impl.f03 b/base/serial/f03/psb_d_coo_impl.f03 similarity index 99% rename from base/serial/f03/psbn_d_coo_impl.f03 rename to base/serial/f03/psb_d_coo_impl.f03 index ba4b3963..25b58b23 100644 --- a/base/serial/f03/psbn_d_coo_impl.f03 +++ b/base/serial/f03/psb_d_coo_impl.f03 @@ -646,7 +646,7 @@ subroutine d_coo_csmv_impl(alpha,a,x,beta,y,info,trans) else if (tra) then - if (alpha.eq.done) then + if (alpha == done) then i = 1 do i=1,nnz ir = a%ja(i) @@ -654,7 +654,7 @@ subroutine d_coo_csmv_impl(alpha,a,x,beta,y,info,trans) y(ir) = y(ir) + a%val(i)*x(jc) enddo - else if (alpha.eq.-done) then + else if (alpha == -done) then do i=1,nnz ir = a%ja(i) @@ -811,7 +811,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) end if else if (tra) then - if (alpha.eq.done) then + if (alpha == done) then i = 1 do i=1,nnz ir = a%ja(i) @@ -819,7 +819,7 @@ subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) y(ir,1:nc) = y(ir,1:nc) + a%val(i)*x(jc,1:nc) enddo - else if (alpha.eq.-done) then + else if (alpha == -done) then do i=1,nnz ir = a%ja(i) diff --git a/base/serial/f03/psbn_d_csr_impl.f03 b/base/serial/f03/psb_d_csr_impl.f03 similarity index 99% rename from base/serial/f03/psbn_d_csr_impl.f03 rename to base/serial/f03/psb_d_csr_impl.f03 index 28413f4f..41de3bbd 100644 --- a/base/serial/f03/psbn_d_csr_impl.f03 +++ b/base/serial/f03/psb_d_csr_impl.f03 @@ -254,7 +254,7 @@ contains end do end if - if (alpha.eq.done) then + if (alpha == done) then do i=1,n do j=irp(i), irp(i+1)-1 @@ -263,7 +263,7 @@ contains end do enddo - else if (alpha.eq.-done) then + else if (alpha == -done) then do i=1,n do j=irp(i), irp(i+1)-1 @@ -545,7 +545,7 @@ contains end do end if - if (alpha.eq.done) then + if (alpha == done) then do i=1,n do j=irp(i), irp(i+1)-1 @@ -554,7 +554,7 @@ contains end do enddo - else if (alpha.eq.-done) then + else if (alpha == -done) then do i=1,n do j=irp(i), irp(i+1)-1 diff --git a/base/serial/f03/psb_s_coo_impl.f03 b/base/serial/f03/psb_s_coo_impl.f03 new file mode 100644 index 00000000..5e6a6725 --- /dev/null +++ b/base/serial/f03/psb_s_coo_impl.f03 @@ -0,0 +1,2463 @@ + +subroutine s_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psb_string_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_cssm_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_base_csmm' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') + m = a%get_nrows() + nc = min(size(x,2) , size(y,2)) + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i,1:nc) = szero + enddo + else + do i = 1, m + y(i,1:nc) = beta*y(i,1:nc) + end do + endif + return + end if + + if (beta == szero) then + call inner_coosm(tra,a,x(:,1:nc),y(:,1:nc),info) + do i = 1, m + y(i,1:nc) = alpha*y(i,1:nc) + end do + else + allocate(tmp(m,nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m,1:nc) = x(1:m,1:nc) + call inner_coosm(tra,a,tmp(:,1:nc),y(:,1:nc),info) + do i = 1, m + y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) + end do + end if + + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='inner_coosm') + goto 9999 + end if + + 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 + + +contains + + subroutine inner_coosm(tra,a,x,y,info) + implicit none + logical, intent(in) :: tra + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: x(:,:) + real(psb_spk_), intent(out) :: y(:,:) + integer, intent(out) :: info + integer :: i,j,k,m, ir, jc + real(psb_spk_), allocatable :: acc(:) + + info = 0 + allocate(acc(size(x,2)), stat=info) + if(info /= 0) then + info=4010 + return + end if + + + if (.not.a%is_sorted()) then + info = 1121 + return + end if + + nnz = a%get_nzeros() + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = szero + do + if (j > nnz) exit + if (a%ia(j) > i) exit + acc = acc + a%val(j)*y(a%ja(j),:) + j = j + 1 + end do + y(i,:) = x(i,:) - acc + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = szero + do + if (j > nnz) exit + if (a%ia(j) > i) exit + if (a%ja(j) == i) then + y(i,:) = (x(i,:) - acc)/a%val(j) + j = j + 1 + exit + end if + acc = acc + a%val(j)*y(a%ja(j),:) + j = j + 1 + end do + end do + end if + + else if (a%is_upper()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = szero + do + if (j < 1) exit + if (a%ia(j) < i) exit + acc = acc + a%val(j)*x(a%ja(j),:) + j = j - 1 + end do + y(i,:) = x(i,:) - acc + end do + + else if (.not.a%is_unit()) then + + j = nnz + do i=a%get_nrows(), 1, -1 + acc = szero + do + if (j < 1) exit + if (a%ia(j) < i) exit + if (a%ja(j) == i) then + y(i,:) = (x(i,:) - acc)/a%val(j) + j = j - 1 + exit + end if + acc = acc + a%val(j)*y(a%ja(j),:) + j = j - 1 + end do + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i,:) = x(i,:) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = y(i,:) + do + if (j < 1) exit + if (a%ia(j) < i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j - 1 + end do + end do + else if (.not.a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + if (a%ja(j) == i) then + y(i,:) = y(i,:) /a%val(j) + j = j - 1 + end if + acc = y(i,:) + do + if (j < 1) exit + if (a%ia(j) < i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j - 1 + end do + end do + + else if (a%is_upper()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = y(i,:) + do + if (j > nnz) exit + if (a%ia(j) > i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j + 1 + end do + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + if (a%ja(j) == i) then + y(i,:) = y(i,:) /a%val(j) + j = j + 1 + end if + acc = y(i,:) + do + if (j > nnz) exit + if (a%ia(j) > i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j + 1 + end do + end do + end if + end if + end if + end if + end subroutine inner_coosm + +end subroutine s_coo_cssm_impl + + + +subroutine s_coo_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psb_string_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_cssv_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_cssv_impl' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i) = szero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + end if + + if (beta == szero) then + call inner_coosv(tra,a,x,y,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + do i = 1, m + y(i) = alpha*y(i) + end do + else + allocate(tmp(m), stat=info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m) = x(1:m) + call inner_coosv(tra,a,tmp,y,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + do i = 1, m + y(i) = alpha*tmp(i) + beta*y(i) + end do + end if + + 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 + +contains + + subroutine inner_coosv(tra,a,x,y,info) + implicit none + logical, intent(in) :: tra + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(out) :: y(:) + integer, intent(out) :: info + + integer :: i,j,k,m, ir, jc, nnz + real(psb_spk_) :: acc + + info = 0 + if (.not.a%is_sorted()) then + info = 1121 + return + end if + + nnz = a%get_nzeros() + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = szero + do + if (j > nnz) exit + if (a%ia(j) > i) exit + acc = acc + a%val(j)*y(a%ja(j)) + j = j + 1 + end do + y(i) = x(i) - acc + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = szero + do + if (j > nnz) exit + if (a%ia(j) > i) exit + if (a%ja(j) == i) then + y(i) = (x(i) - acc)/a%val(j) + j = j + 1 + exit + end if + acc = acc + a%val(j)*y(a%ja(j)) + j = j + 1 + end do + end do + end if + + else if (a%is_upper()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = szero + do + if (j < 1) exit + if (a%ia(j) < i) exit + acc = acc + a%val(j)*y(a%ja(j)) + j = j - 1 + end do + y(i) = x(i) - acc + end do + + else if (.not.a%is_unit()) then + + j = nnz + do i=a%get_nrows(), 1, -1 + acc = szero + do + if (j < 1) exit + if (a%ia(j) < i) exit + if (a%ja(j) == i) then + y(i) = (x(i) - acc)/a%val(j) + j = j - 1 + exit + end if + acc = acc + a%val(j)*y(a%ja(j)) + j = j - 1 + end do + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i) = x(i) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = y(i) + do + if (j < 1) exit + if (a%ia(j) < i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j - 1 + end do + end do + else if (.not.a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + if (a%ja(j) == i) then + y(i) = y(i) /a%val(j) + j = j - 1 + end if + acc = y(i) + do + if (j < 1) exit + if (a%ia(j) < i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j - 1 + end do + end do + + else if (a%is_upper()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = y(i) + do + if (j > nnz) exit + if (a%ia(j) > i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j + 1 + end do + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + if (a%ja(j) == i) then + y(i) = y(i) /a%val(j) + j = j + 1 + end if + acc = y(i) + do + if (j > nnz) exit + if (a%ia(j) > i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j + 1 + end do + end do + end if + end if + end if + end if + + end subroutine inner_coosv + + +end subroutine s_coo_cssv_impl + +subroutine s_coo_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csMv_impl + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_csmv_impl' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + tra = ((trans_=='T').or.(trans_=='t')) + + + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + nnz = a%get_nzeros() + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i) = szero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + else + if (a%is_triangle().and.a%is_unit()) then + if (beta == szero) then + do i = 1, min(m,n) + y(i) = alpha*x(i) + enddo + do i = min(m,n)+1, m + y(i) = szero + enddo + else + do i = 1, min(m,n) + y(i) = beta*y(i) + alpha*x(i) + end do + do i = min(m,n)+1, m + y(i) = beta*y(i) + enddo + endif + else + if (beta == szero) then + do i = 1, m + y(i) = szero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + + endif + + end if + + if (.not.tra) then + i = 1 + j = i + if (nnz > 0) then + ir = a%ia(1) + acc = szero + do + if (i>nnz) then + y(ir) = y(ir) + alpha * acc + exit + endif + if (a%ia(i) /= ir) then + y(ir) = y(ir) + alpha * acc + ir = a%ia(i) + acc = szero + endif + acc = acc + a%val(i) * x(a%ja(i)) + i = i + 1 + enddo + end if + + else if (tra) then + + if (alpha == sone) then + i = 1 + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir) = y(ir) + a%val(i)*x(jc) + enddo + + else if (alpha == -sone) then + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir) = y(ir) - a%val(i)*x(jc) + enddo + + else + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir) = y(ir) + alpha*a%val(i)*x(jc) + enddo + + end if !.....end testing on alpha + + endif + + 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 s_coo_csmv_impl + + +subroutine s_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csmm_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_coo_csmm_impl' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + + tra = ((trans_=='T').or.(trans_=='t')) + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + nnz = a%get_nzeros() + + nc = min(size(x,2), size(y,2)) + allocate(acc(nc),stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i,1:nc) = szero + enddo + else + do i = 1, m + y(i,1:nc) = beta*y(i,1:nc) + end do + endif + return + else + if (a%is_triangle().and.a%is_unit()) then + if (beta == szero) then + do i = 1, min(m,n) + y(i,1:nc) = alpha*x(i,1:nc) + enddo + do i = min(m,n)+1, m + y(i,1:nc) = szero + enddo + else + do i = 1, min(m,n) + y(i,1:nc) = beta*y(i,1:nc) + alpha*x(i,1:nc) + end do + do i = min(m,n)+1, m + y(i,1:nc) = beta*y(i,1:nc) + enddo + endif + else + if (beta == szero) then + do i = 1, m + y(i,1:nc) = szero + enddo + else + do i = 1, m + y(i,1:nc) = beta*y(i,1:nc) + end do + endif + + endif + + end if + + if (.not.tra) then + i = 1 + j = i + if (nnz > 0) then + ir = a%ia(1) + acc = szero + do + if (i>nnz) then + y(ir,1:nc) = y(ir,1:nc) + alpha * acc + exit + endif + if (a%ia(i) /= ir) then + y(ir,1:nc) = y(ir,1:nc) + alpha * acc + ir = a%ia(i) + acc = szero + endif + acc = acc + a%val(i) * x(a%ja(i),1:nc) + i = i + 1 + enddo + end if + + else if (tra) then + if (alpha == sone) then + i = 1 + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir,1:nc) = y(ir,1:nc) + a%val(i)*x(jc,1:nc) + enddo + + else if (alpha == -sone) then + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir,1:nc) = y(ir,1:nc) - a%val(i)*x(jc,1:nc) + enddo + + else + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir,1:nc) = y(ir,1:nc) + alpha*a%val(i)*x(jc,1:nc) + enddo + + end if !.....end testing on alpha + + endif + + 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 s_coo_csmm_impl + +function s_coo_csnmi_impl(a) result(res) + use psb_error_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csnmi_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_base_csnmi' + logical, parameter :: debug=.false. + + + res = szero + nnz = a%get_nzeros() + i = 1 + j = i + do while (i<=nnz) + do while ((a%ia(j) == a%ia(i)).and. (j <= nnz)) + j = j+1 + enddo + acc = szero + do k=i, j-1 + acc = acc + abs(a%val(k)) + end do + res = max(res,acc) + i = j + end do + +end function s_coo_csnmi_impl + + + +!==================================== +! +! +! +! Data management +! +! +! +! +! +!==================================== + + + +subroutine s_coo_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csgetptn_impl + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + logical :: appens_, rscale_, cscale_ + integer :: nzin_, jmin_, jmax_, err_act, i + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax= psb_debug_serial_)& + & write(debug_unit,*) trim(name), ': srtdcoo ' + do + ip = psb_ibsrch(irw,nza,a%ia) + if (ip /= -1) exit + irw = irw + 1 + if (irw > imax) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error? ',& + & irw,lrw,imin + exit + end if + end do + + if (ip /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (ip < 2) exit + if (a%ia(ip-1) == irw) then + ip = ip -1 + else + exit + end if + end do + + end if + + do + jp = psb_ibsrch(lrw,nza,a%ia) + if (jp /= -1) exit + lrw = lrw - 1 + if (irw > lrw) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error?' + exit + end if + end do + + if (jp /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (jp == nza) exit + if (a%ia(jp+1) == lrw) then + jp = jp + 1 + else + exit + end if + end do + end if + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': ip jp',ip,jp,nza + if ((ip /= -1) .and.(jp /= -1)) then + ! Now do the copy. + nzt = jp - ip +1 + nz = 0 + + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info /= 0) return + + if (present(iren)) then + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + ia(nzin_) = iren(a%ia(i)) + ja(nzin_) = iren(a%ja(i)) + end if + enddo + else + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + ia(nzin_) = a%ia(i) + ja(nzin_) = a%ja(i) + end if + enddo + end if + else + nz = 0 + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': unsorted ' + + nzt = (nza*(lrw-irw+1))/max(a%get_nrows(),1) + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info /= 0) return + + if (present(iren)) then + k = 0 + do i=1, a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info /= 0) return + end if + ia(nzin_+k) = iren(a%ia(i)) + ja(nzin_+k) = iren(a%ja(i)) + endif + enddo + else + k = 0 + do i=1,a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info /= 0) return + + end if + ia(nzin_+k) = (a%ia(i)) + ja(nzin_+k) = (a%ja(i)) + endif + enddo + nzin_=nzin_+k + end if + nz = k + end if + + end subroutine coo_getptn + +end subroutine s_coo_csgetptn_impl + + +subroutine s_coo_csgetrow_impl(imin,imax,a,nz,ia,ja,val,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csgetrow_impl + implicit none + + class(psb_s_coo_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + logical :: appens_, rscale_, cscale_ + integer :: nzin_, jmin_, jmax_, err_act, i + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax= psb_debug_serial_)& + & write(debug_unit,*) trim(name), ': srtdcoo ' + do + ip = psb_ibsrch(irw,nza,a%ia) + if (ip /= -1) exit + irw = irw + 1 + if (irw > imax) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error? ',& + & irw,lrw,imin + exit + end if + end do + + if (ip /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (ip < 2) exit + if (a%ia(ip-1) == irw) then + ip = ip -1 + else + exit + end if + end do + + end if + + do + jp = psb_ibsrch(lrw,nza,a%ia) + if (jp /= -1) exit + lrw = lrw - 1 + if (irw > lrw) then + write(debug_unit,*) trim(name),& + & 'Warning : did not find any rows. Is this an error?' + exit + end if + end do + + if (jp /= -1) then + ! expand [ip,jp] to contain all row entries. + do + if (jp == nza) exit + if (a%ia(jp+1) == lrw) then + jp = jp + 1 + else + exit + end if + end do + end if + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': ip jp',ip,jp,nza + if ((ip /= -1) .and.(jp /= -1)) then + ! Now do the copy. + nzt = jp - ip +1 + nz = 0 + + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + if (present(iren)) then + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + val(nzin_) = a%val(i) + ia(nzin_) = iren(a%ia(i)) + ja(nzin_) = iren(a%ja(i)) + end if + enddo + else + do i=ip,jp + if ((jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + nzin_ = nzin_ + 1 + nz = nz + 1 + val(nzin_) = a%val(i) + ia(nzin_) = a%ia(i) + ja(nzin_) = a%ja(i) + end if + enddo + end if + else + nz = 0 + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': unsorted ' + + nzt = (nza*(lrw-irw+1))/max(a%get_nrows(),1) + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + if (present(iren)) then + k = 0 + do i=1, a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + end if + val(nzin_+k) = a%val(i) + ia(nzin_+k) = iren(a%ia(i)) + ja(nzin_+k) = iren(a%ja(i)) + endif + enddo + else + k = 0 + do i=1,a%get_nzeros() + if ((a%ia(i)>=irw).and.(a%ia(i)<=lrw).and.& + & (jmin <= a%ja(i)).and.(a%ja(i)<=jmax)) then + k = k + 1 + if (k > nzt) then + nzt = k + call psb_ensure_size(nzin_+nzt,ia,info) + if (info==0) call psb_ensure_size(nzin_+nzt,ja,info) + if (info==0) call psb_ensure_size(nzin_+nzt,val,info) + if (info /= 0) return + + end if + val(nzin_+k) = a%val(i) + ia(nzin_+k) = (a%ia(i)) + ja(nzin_+k) = (a%ja(i)) + endif + enddo + nzin_=nzin_+k + end if + nz = k + end if + + end subroutine coo_getrow + +end subroutine s_coo_csgetrow_impl + + +subroutine s_coo_csput_impl(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + use psb_sort_mod + use psb_s_base_mat_mod, psb_protect_name => s_coo_csput_impl + implicit none + + class(psb_s_coo_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='s_coo_csput_impl' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + info = 0 + call psb_erractionsave(err_act) + + if (nz <= 0) then + info = 10 + int_err(1)=1 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(ia) < nz) then + info = 35 + int_err(1)=2 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (size(ja) < nz) then + info = 35 + int_err(1)=3 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(val) < nz) then + info = 35 + int_err(1)=4 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (nz == 0) return + + + nza = a%get_nzeros() + isza = a%get_size() + if (a%is_bld()) then + ! Build phase. Must handle reallocations in a sensible way. + if (isza < (nza+nz)) then + call a%reallocate(max(nza+nz,int(1.5*isza))) + isza = a%get_size() + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,& + & imin,imax,jmin,jmax,info,gtl) + call a%set_nzeros(nza) + call a%set_sorted(.false.) + + + else if (a%is_upd()) then + + call s_coo_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + if (info /= 0) then + info = 1121 + end if + + else + ! State is wrong. + info = 1121 + end if + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + 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 + + +contains + + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,maxsz,& + & imin,imax,jmin,jmax,info,gtl) + implicit none + + integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz + integer, intent(in) :: ia(:),ja(:) + integer, intent(inout) :: nza,ia1(:),ia2(:) + real(psb_spk_), intent(in) :: val(:) + real(psb_spk_), intent(inout) :: aspk(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + integer :: i,ir,ic,ng + + info = 0 + if (present(gtl)) then + ng = size(gtl) + + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + if (nza > maxsz) then + info = -91 + return + endif + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end if + end do + else + + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + if (nza > maxsz) then + info = -92 + return + endif + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end do + end if + + end subroutine psb_inner_ins + + + subroutine s_coo_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + + use psb_const_mod + use psb_realloc_mod + use psb_string_mod + implicit none + + class(psb_s_coo_sparse_mat), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: ia(:),ja(:) + real(psb_spk_), intent(in) :: val(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + integer :: i,ir,ic, ilr, ilc, ip, & + & i1,i2,nc,nnz,dupl,ng, nr + integer :: debug_level, debug_unit + character(len=20) :: name='s_coo_srch_upd' + + info = 0 + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + dupl = a%get_dupl() + + if (.not.a%is_sorted()) then + info = -4 + return + end if + + ilr = -1 + ilc = -1 + nnz = a%get_nzeros() + nr = a%get_nrows() + nc = a%get_ncols() + + + if (present(gtl)) then + ng = size(gtl) + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + if ((ir > 0).and.(ir <= nr)) then + ic = gtl(ic) + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nnz,a%ia) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia(i2+1) /= a%ia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia(i1-1) /= a%ia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,a%ja(i1:i2)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + endif + end if + end do + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= nr)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nnz,a%ia) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia(i2+1) /= a%ia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia(i1-1) /= a%ia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,a%ja(i1:i2)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nr)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nnz,a%ia) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia(i2+1) /= a%ia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia(i1-1) /= a%ia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,a%ja(i1:i2)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + info = i + return + end if + end if + end do + + case(psb_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nr)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nnz,a%ia) + i2 = i1 + do + if (i2+1 > nnz) exit + if (a%ia(i2+1) /= a%ia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (a%ia(i1-1) /= a%ia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,a%ja(i1:i2)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + info = i + return + end if + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + end if + + end subroutine s_coo_srch_upd + +end subroutine s_coo_csput_impl + + +subroutine s_cp_coo_to_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_cp_coo_to_coo_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + info = 0 + call b%psb_s_base_sparse_mat%cp_from(a%psb_s_base_sparse_mat) + + call b%set_nzeros(a%get_nzeros()) + call b%reallocate(a%get_nzeros()) + + b%ia(:) = a%ia(:) + b%ja(:) = a%ja(:) + b%val(:) = a%val(:) + + call b%fix(info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_cp_coo_to_coo_impl + +subroutine s_cp_coo_from_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_cp_coo_from_coo_impl + implicit none + class(psb_s_coo_sparse_mat), intent(out) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + integer :: m,n,nz + + + call psb_erractionsave(err_act) + info = 0 + call a%psb_s_base_sparse_mat%cp_from(b%psb_s_base_sparse_mat) + call a%set_nzeros(b%get_nzeros()) + call a%reallocate(b%get_nzeros()) + + a%ia(:) = b%ia(:) + a%ja(:) = b%ja(:) + a%val(:) = b%val(:) + + call a%fix(info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_cp_coo_from_coo_impl + + +subroutine s_cp_coo_to_fmt_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_cp_coo_to_fmt_impl + implicit none + class(psb_s_coo_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + info = 0 + + call b%cp_from_coo(a,info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_cp_coo_to_fmt_impl + +subroutine s_cp_coo_from_fmt_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_cp_coo_from_fmt_impl + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + integer :: m,n,nz + + + call psb_erractionsave(err_act) + info = 0 + + call b%cp_to_coo(a,info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_cp_coo_from_fmt_impl + + +subroutine s_fix_coo_impl(a,info,idir) + use psb_const_mod + use psb_error_mod + use psb_realloc_mod + use psb_string_mod + use psb_ip_reord_mod + use psb_s_base_mat_mod, psb_protect_name => s_fix_coo_impl + implicit none + + class(psb_s_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + integer, allocatable :: iaux(:) + !locals + Integer :: nza, nzl,iret,idir_, dupl_ + integer :: i,j, irw, icl, err_act + integer :: debug_level, debug_unit + character(len=20) :: name = 'psb_fixcoo' + + info = 0 + + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if(debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': start ',& + & size(a%ia),size(a%ja) + if (present(idir)) then + idir_ = idir + else + idir_ = 0 + endif + + nza = a%get_nzeros() + if (nza < 2) return + + dupl_ = a%get_dupl() + + call s_fix_coo_inner(nza,dupl_,a%ia,a%ja,a%val,i,info,idir_) + + call a%set_sorted() + call a%set_nzeros(i) + call a%set_asb() + + + 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 s_fix_coo_impl + + + +subroutine s_fix_coo_inner(nzin,dupl,ia,ja,val,nzout,info,idir) + use psb_const_mod + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_fix_coo_inner + use psb_string_mod + use psb_ip_reord_mod + implicit none + + integer, intent(in) :: nzin, dupl + integer, intent(inout) :: ia(:), ja(:) + real(psb_spk_), intent(inout) :: val(:) + integer, intent(out) :: nzout, info + integer, intent(in), optional :: idir + !locals + integer, allocatable :: iaux(:) + Integer :: nza, nzl,iret,idir_, dupl_ + integer :: i,j, irw, icl, err_act + integer :: debug_level, debug_unit + character(len=20) :: name = 'psb_fixcoo' + + info = 0 + + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if(debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': start ',& + & size(ia),size(ja) + if (present(idir)) then + idir_ = idir + else + idir_ = 0 + endif + + + if (nzin < 2) return + + dupl_ = dupl + + allocate(iaux(nzin+2),stat=info) + if (info /= 0) return + + + select case(idir_) + + case(0) ! Row major order + + call msort_up(nzin,ia(1),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzin,val,ia,ja,iaux) + i = 1 + j = i + do while (i <= nzin) + do while ((ia(j) == ia(i))) + j = j+1 + if (j > nzin) exit + enddo + nzl = j - i + call msort_up(nzl,ja(i),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzl,val(i:i+nzl-1),& + & ia(i:i+nzl-1),ja(i:i+nzl-1),iaux) + i = j + enddo + + i = 1 + irw = ia(i) + icl = ja(i) + j = 1 + + select case(dupl_) + case(psb_dupl_ovwrt_) + + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + val(i) = val(j) + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + + case(psb_dupl_add_) + + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + val(i) = val(i) + val(j) + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + + case(psb_dupl_err_) + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + call psb_errpush(130,name) + goto 9999 + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + case default + write(0,*) 'Error in fix_coo: unsafe dupl',dupl_ + + end select + + + if(debug_level >= psb_debug_serial_)& + & write(debug_unit,*) trim(name),': end second loop' + + case(1) ! Col major order + + call msort_up(nzin,ja(1),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzin,val,ia,ja,iaux) + i = 1 + j = i + do while (i <= nzin) + do while ((ja(j) == ja(i))) + j = j+1 + if (j > nzin) exit + enddo + nzl = j - i + call msort_up(nzl,ia(i),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzl,val(i:i+nzl-1),& + & ia(i:i+nzl-1),ja(i:i+nzl-1),iaux) + i = j + enddo + + i = 1 + irw = ia(i) + icl = ja(i) + j = 1 + + + select case(dupl_) + case(psb_dupl_ovwrt_) + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + val(i) = val(j) + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + + case(psb_dupl_add_) + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + val(i) = val(i) + val(j) + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + + case(psb_dupl_err_) + do + j = j + 1 + if (j > nzin) exit + if ((ia(j) == irw).and.(ja(j) == icl)) then + call psb_errpush(130,name) + goto 9999 + else + i = i+1 + val(i) = val(j) + ia(i) = ia(j) + ja(i) = ja(j) + irw = ia(i) + icl = ja(i) + endif + enddo + case default + write(0,*) 'Error in fix_coo: unsafe dupl',dupl_ + end select + if (debug_level >= psb_debug_serial_)& + & write(debug_unit,*) trim(name),': end second loop' + case default + write(debug_unit,*) trim(name),': unknown direction ',idir_ + end select + + nzout = i + + deallocate(iaux) + + 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 s_fix_coo_inner + + + + +subroutine s_mv_coo_to_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_mv_coo_to_coo_impl + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + info = 0 + call b%psb_s_base_sparse_mat%mv_from(a%psb_s_base_sparse_mat) + call b%set_nzeros(a%get_nzeros()) + call b%reallocate(a%get_nzeros()) + + call move_alloc(a%ia, b%ia) + call move_alloc(a%ja, b%ja) + call move_alloc(a%val, b%val) + call a%free() + + call b%fix(info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_mv_coo_to_coo_impl + +subroutine s_mv_coo_from_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_mv_coo_from_coo_impl + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + integer :: m,n,nz + + + call psb_erractionsave(err_act) + info = 0 + call a%psb_s_base_sparse_mat%mv_from(b%psb_s_base_sparse_mat) + call a%set_nzeros(b%get_nzeros()) + call a%reallocate(b%get_nzeros()) + + call move_alloc(b%ia , a%ia ) + call move_alloc(b%ja , a%ja ) + call move_alloc(b%val, a%val ) + call b%free() + + a%ia(:) = b%ia(:) + a%ja(:) = b%ja(:) + a%val(:) = b%val(:) + + call a%fix(info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_mv_coo_from_coo_impl + + +subroutine s_mv_coo_to_fmt_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_mv_coo_to_fmt_impl + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + info = 0 + + call b%mv_from_coo(a,info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_mv_coo_to_fmt_impl + +subroutine s_mv_coo_from_fmt_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psb_s_base_mat_mod, psb_protect_name => s_mv_coo_from_fmt_impl + implicit none + class(psb_s_coo_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + integer :: m,n,nz + + + call psb_erractionsave(err_act) + info = 0 + + call b%mv_to_coo(a,info) + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + +end subroutine s_mv_coo_from_fmt_impl diff --git a/base/serial/f03/psb_s_csr_impl.f03 b/base/serial/f03/psb_s_csr_impl.f03 new file mode 100644 index 00000000..12a675a5 --- /dev/null +++ b/base/serial/f03/psb_s_csr_impl.f03 @@ -0,0 +1,2015 @@ + +!===================================== +! +! +! +! Computational routines +! +! +! +! +! +! +!===================================== + +subroutine s_csr_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_csmv_impl + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_csmv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + call s_csr_csmv_inner(m,n,alpha,a%irp,a%ja,a%val,& + & a%is_triangle(),a%is_unit(),& + & x,beta,y,tra) + + 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 + +contains + subroutine s_csr_csmv_inner(m,n,alpha,irp,ja,val,is_triangle,is_unit,& + & x,beta,y,tra) + integer, intent(in) :: m,n,irp(*),ja(*) + real(psb_spk_), intent(in) :: alpha, beta, x(*),val(*) + real(psb_spk_), intent(inout) :: y(*) + logical, intent(in) :: is_triangle,is_unit,tra + + + integer :: i,j,k, ir, jc + real(psb_spk_) :: acc + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i) = szero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + end if + + + if (.not.tra) then + + if (beta == szero) then + + if (alpha == sone) then + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = acc + end do + + else if (alpha == -sone) then + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = -acc + end do + + else + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = alpha*acc + end do + + end if + + + else if (beta == sone) then + + if (alpha == sone) then + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = y(i) + acc + end do + + else if (alpha == -sone) then + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = y(i) -acc + end do + + else + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = y(i) + alpha*acc + end do + + end if + + else if (beta == -sone) then + + if (alpha == sone) then + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = -y(i) + acc + end do + + else if (alpha == -sone) then + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = -y(i) -acc + end do + + else + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = -y(i) + alpha*acc + end do + + end if + + else + + if (alpha == sone) then + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = beta*y(i) + acc + end do + + else if (alpha == -sone) then + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = beta*y(i) - acc + end do + + else + + do i=1,m + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j) * x(ja(j)) + enddo + y(i) = beta*y(i) + alpha*acc + end do + + end if + + end if + + else if (tra) then + + if (beta == szero) then + do i=1, m + y(i) = szero + end do + else if (beta == sone) then + ! Do nothing + else if (beta == -sone) then + do i=1, m + y(i) = -y(i) + end do + else + do i=1, m + y(i) = beta*y(i) + end do + end if + + if (alpha == sone) then + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir) = y(ir) + val(j)*x(i) + end do + enddo + + else if (alpha == -sone) then + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir) = y(ir) - val(j)*x(i) + end do + enddo + + else + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir) = y(ir) + alpha*val(j)*x(i) + end do + enddo + + end if + + endif + + if (is_triangle.and.is_unit) then + do i=1, min(m,n) + y(i) = y(i) + alpha*x(i) + end do + end if + + + end subroutine s_csr_csmv_inner + + +end subroutine s_csr_csmv_impl + +subroutine s_csr_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_csmm_impl + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_csmm' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + tra = (psb_toupper(trans_)=='T').or.(psb_toupper(trans_)=='C') + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + nc = min(size(x,2) , size(y,2) ) + + allocate(acc(nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + call s_csr_csmm_inner(m,n,nc,alpha,a%irp,a%ja,a%val, & + & a%is_triangle(),a%is_unit(),x,size(x,1), & + & beta,y,size(y,1),tra,acc) + + + 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 + +contains + subroutine s_csr_csmm_inner(m,n,nc,alpha,irp,ja,val,& + & is_triangle,is_unit,x,ldx,beta,y,ldy,tra,acc) + integer, intent(in) :: m,n,ldx,ldy,nc,irp(*),ja(*) + real(psb_spk_), intent(in) :: alpha, beta, x(ldx,*),val(*) + real(psb_spk_), intent(inout) :: y(ldy,*) + logical, intent(in) :: is_triangle,is_unit,tra + + real(psb_spk_), intent(inout) :: acc(*) + integer :: i,j,k, ir, jc + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i,1:nc) = szero + enddo + else + do i = 1, m + y(i,1:nc) = beta*y(i,1:nc) + end do + endif + return + end if + + if (.not.tra) then + + if (beta == szero) then + + if (alpha == sone) then + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = acc(1:nc) + end do + + else if (alpha == -sone) then + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = -acc(1:nc) + end do + + else + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = alpha*acc(1:nc) + end do + + end if + + + else if (beta == sone) then + + if (alpha == sone) then + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = y(i,1:nc) + acc(1:nc) + end do + + else if (alpha == -sone) then + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = y(i,1:nc) -acc(1:nc) + end do + + else + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = y(i,1:nc) + alpha*acc(1:nc) + end do + + end if + + else if (beta == -sone) then + + if (alpha == sone) then + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = -y(i,1:nc) + acc(1:nc) + end do + + else if (alpha == -sone) then + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = -y(i,1:nc) -acc(1:nc) + end do + + else + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = -y(i,1:nc) + alpha*acc(1:nc) + end do + + end if + + else + + if (alpha == sone) then + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = beta*y(i,1:nc) + acc(1:nc) + end do + + else if (alpha == -sone) then + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = beta*y(i,1:nc) - acc(1:nc) + end do + + else + + do i=1,m + acc(1:nc) = szero + do j=irp(i), irp(i+1)-1 + acc(1:nc) = acc(1:nc) + val(j) * x(ja(j),1:nc) + enddo + y(i,1:nc) = beta*y(i,1:nc) + alpha*acc(1:nc) + end do + + end if + + end if + + else if (tra) then + + if (beta == szero) then + do i=1, m + y(i,1:nc) = szero + end do + else if (beta == sone) then + ! Do nothing + else if (beta == -sone) then + do i=1, m + y(i,1:nc) = -y(i,1:nc) + end do + else + do i=1, m + y(i,1:nc) = beta*y(i,1:nc) + end do + end if + + if (alpha == sone) then + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir,1:nc) = y(ir,1:nc) + val(j)*x(i,1:nc) + end do + enddo + + else if (alpha == -sone) then + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir,1:nc) = y(ir,1:nc) - val(j)*x(i,1:nc) + end do + enddo + + else + + do i=1,n + do j=irp(i), irp(i+1)-1 + ir = ja(j) + y(ir,1:nc) = y(ir,1:nc) + alpha*val(j)*x(i,1:nc) + end do + enddo + + end if + + endif + + if (is_triangle.and.is_unit) then + do i=1, min(m,n) + y(i,1:nc) = y(i,1:nc) + alpha*x(i,1:nc) + end do + end if + +end subroutine s_csr_csmm_inner + +end subroutine s_csr_csmm_impl + + +subroutine s_csr_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_cssv_impl + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csr_cssv' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i) = szero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + end if + + if (beta == szero) then + + call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& + & a%irp,a%ja,a%val,x,y) + if (alpha == sone) then + ! do nothing + else if (alpha == -sone) then + do i = 1, m + y(i) = -y(i) + end do + else + do i = 1, m + y(i) = alpha*y(i) + end do + end if + else + allocate(tmp(m), stat=info) + if (info /= 0) then + return + end if + tmp(1:m) = x(1:m) + call inner_csrsv(tra,a%is_lower(),a%is_unit(),a%get_nrows(),& + & a%irp,a%ja,a%val,tmp,y) + do i = 1, m + y(i) = alpha*tmp(i) + beta*y(i) + end do + end if + + 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 + +contains + + subroutine inner_csrsv(tra,lower,unit,n,irp,ja,val,x,y) + implicit none + logical, intent(in) :: tra,lower,unit + integer, intent(in) :: irp(*), ja(*),n + real(psb_spk_), intent(in) :: val(*) + real(psb_spk_), intent(in) :: x(*) + real(psb_spk_), intent(out) :: y(*) + + integer :: i,j,k,m, ir, jc + real(psb_spk_) :: acc + + if (.not.tra) then + + if (lower) then + if (unit) then + do i=1, n + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) + end do + y(i) = x(i) - acc + end do + else if (.not.unit) then + do i=1, n + acc = szero + do j=irp(i), irp(i+1)-2 + acc = acc + val(j)*y(ja(j)) + end do + y(i) = (x(i) - acc)/val(irp(i+1)-1) + end do + end if + else if (.not.lower) then + + if (unit) then + do i=n, 1, -1 + acc = szero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) + end do + y(i) = x(i) - acc + end do + else if (.not.unit) then + do i=n, 1, -1 + acc = szero + do j=irp(i)+1, irp(i+1)-1 + acc = acc + val(j)*y(ja(j)) + end do + y(i) = (x(i) - acc)/val(irp(i)) + end do + end if + + end if + + else if (tra) then + + do i=1, n + y(i) = x(i) + end do + + if (lower) then + if (unit) then + do i=n, 1, -1 + acc = y(i) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc + end do + end do + else if (.not.unit) then + do i=n, 1, -1 + y(i) = y(i)/val(irp(i+1)-1) + acc = y(i) + do j=irp(i), irp(i+1)-2 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc + end do + end do + end if + else if (.not.lower) then + + if (unit) then + do i=1, n + acc = y(i) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc + end do + end do + else if (.not.unit) then + do i=1, n + y(i) = y(i)/val(irp(i)) + acc = y(i) + do j=irp(i)+1, irp(i+1)-1 + jc = ja(j) + y(jc) = y(jc) - val(j)*acc + end do + end do + end if + + end if + end if + end subroutine inner_csrsv + +end subroutine s_csr_cssv_impl + + + +subroutine s_csr_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_cssm_impl + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='s_base_csmm' + logical, parameter :: debug=.false. + + info = 0 + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + nc = min(size(x,2) , size(y,2)) + + if (.not. (a%is_triangle())) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + + if (alpha == szero) then + if (beta == szero) then + do i = 1, m + y(i,:) = szero + enddo + else + do i = 1, m + y(i,:) = beta*y(i,:) + end do + endif + return + end if + + if (beta == szero) then + call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1),y,size(y,1),info) + do i = 1, m + y(i,1:nc) = alpha*y(i,1:nc) + end do + else + allocate(tmp(m,nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m,:) = x(1:m,1:nc) + call inner_csrsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,tmp,size(tmp,1),y,size(y,1),info) + do i = 1, m + y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) + end do + end if + + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='inner_csrsm') + goto 9999 + end if + + 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 + + +contains + + subroutine inner_csrsm(tra,lower,unit,nr,nc,& + & irp,ja,val,x,ldx,y,ldy,info) + implicit none + logical, intent(in) :: tra,lower,unit + integer, intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*) + real(psb_spk_), intent(in) :: val(*), x(ldx,*) + real(psb_spk_), intent(out) :: y(ldy,*) + integer, intent(out) :: info + integer :: i,j,k,m, ir, jc + real(psb_spk_), allocatable :: acc(:) + + info = 0 + allocate(acc(nc), stat=info) + if(info /= 0) then + info=4010 + return + end if + + + if (.not.tra) then + + if (lower) then + if (unit) then + do i=1, nr + acc = szero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*y(a%ja(j),1:nc) + end do + y(i,1:nc) = x(i,1:nc) - acc + end do + else if (.not.unit) then + do i=1, nr + acc = szero + do j=a%irp(i), a%irp(i+1)-2 + acc = acc + a%val(j)*y(a%ja(j),1:nc) + end do + y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i+1)-1) + end do + end if + else if (.not.lower) then + + if (unit) then + do i=nr, 1, -1 + acc = szero + do j=a%irp(i), a%irp(i+1)-1 + acc = acc + a%val(j)*y(a%ja(j),1:nc) + end do + y(i,1:nc) = x(i,1:nc) - acc + end do + else if (.not.unit) then + do i=nr, 1, -1 + acc = szero + do j=a%irp(i)+1, a%irp(i+1)-1 + acc = acc + a%val(j)*y(a%ja(j),1:nc) + end do + y(i,1:nc) = (x(i,1:nc) - acc)/a%val(a%irp(i)) + end do + end if + + end if + + else if (tra) then + + do i=1, nr + y(i,1:nc) = x(i,1:nc) + end do + + if (lower) then + if (unit) then + do i=nr, 1, -1 + acc = y(i,1:nc) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc + end do + end do + else if (.not.unit) then + do i=nr, 1, -1 + y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i+1)-1) + acc = y(i,1:nc) + do j=a%irp(i), a%irp(i+1)-2 + jc = a%ja(j) + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc + end do + end do + end if + else if (.not.lower) then + + if (unit) then + do i=1, nr + acc = y(i,1:nc) + do j=a%irp(i), a%irp(i+1)-1 + jc = a%ja(j) + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc + end do + end do + else if (.not.unit) then + do i=1, nr + y(i,1:nc) = y(i,1:nc)/a%val(a%irp(i)) + acc = y(i,1:nc) + do j=a%irp(i)+1, a%irp(i+1)-1 + jc = a%ja(j) + y(jc,1:nc) = y(jc,1:nc) - a%val(j)*acc + end do + end do + end if + + end if + end if + end subroutine inner_csrsm + +end subroutine s_csr_cssm_impl + +function s_csr_csnmi_impl(a) result(res) + use psb_error_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_csnmi_impl + implicit none + class(psb_s_csr_sparse_mat), intent(in) :: a + real(psb_spk_) :: res + + integer :: i,j,k,m,n, nr, ir, jc, nc + real(psb_spk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='s_csnmi' + logical, parameter :: debug=.false. + + + res = szero + + do i = 1, a%get_nrows() + acc = szero + do j=a%irp(i),a%irp(i+1)-1 + acc = acc + abs(a%val(j)) + end do + res = max(res,acc) + end do + +end function s_csr_csnmi_impl + +!===================================== +! +! +! +! Data management +! +! +! +! +! +!===================================== + + +subroutine s_csr_csgetptn_impl(imin,imax,a,nz,ia,ja,info,& + & jmin,jmax,iren,append,nzin,rscale,cscale) + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_csr_csgetptn_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + logical :: appens_, rscale_, cscale_ + integer :: nzin_, jmin_, jmax_, err_act, i + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax s_csr_csgetrow_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + integer, intent(in) :: imin,imax + integer, intent(out) :: nz + integer, allocatable, intent(inout) :: ia(:), ja(:) + real(psb_spk_), allocatable, intent(inout) :: val(:) + integer,intent(out) :: info + logical, intent(in), optional :: append + integer, intent(in), optional :: iren(:) + integer, intent(in), optional :: jmin,jmax, nzin + logical, intent(in), optional :: rscale,cscale + + logical :: appens_, rscale_, cscale_ + integer :: nzin_, jmin_, jmax_, err_act, i + character(len=20) :: name='csget' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax s_csr_csput_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + real(psb_spk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='s_csr_csput' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + info = 0 + nza = a%get_nzeros() + + if (a%is_bld()) then + ! Build phase should only ever be in COO + info = 1121 + + else if (a%is_upd()) then + call s_csr_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + + if (info /= 0) then + + info = 1121 + end if + + else + ! State is wrong. + info = 1121 + end if + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + + 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 + + +contains + + subroutine s_csr_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + + use psb_const_mod + use psb_realloc_mod + use psb_string_mod + use psb_sort_mod + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + integer, intent(in) :: nz, imin,imax,jmin,jmax + integer, intent(in) :: ia(:),ja(:) + real(psb_spk_), intent(in) :: val(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + integer :: i,ir,ic, ilr, ilc, ip, & + & i1,i2,nr,nc,nnz,dupl,ng + integer :: debug_level, debug_unit + character(len=20) :: name='s_csr_srch_upd' + + info = 0 + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + dupl = a%get_dupl() + + if (.not.a%is_sorted()) then + info = -4 + return + end if + + ilr = -1 + ilc = -1 + nnz = a%get_nzeros() + nr = a%get_nrows() + nc = a%get_ncols() + + if (present(gtl)) then + ng = size(gtl) + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + + else + + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + + if ((ir > 0).and.(ir <= nr)) then + + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + end if + + end subroutine s_csr_srch_upd + +end subroutine s_csr_csput_impl + + + +subroutine s_cp_csr_from_coo_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_cp_csr_from_coo_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + type(psb_s_coo_sparse_mat) :: tmp + integer, allocatable :: itemp(:) + !locals + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + ! This is to have fix_coo called behind the scenes + call tmp%cp_from_coo(b,info) + if (info ==0) call a%mv_from_coo(tmp,info) + +end subroutine s_cp_csr_from_coo_impl + + + +subroutine s_cp_csr_to_coo_impl(a,b,info) + use psb_const_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_cp_csr_to_coo_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + integer, allocatable :: itemp(:) + !locals + logical :: rwshr_ + Integer :: nza, nr, nc,i,j,irw, idl,err_act + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + nr = a%get_nrows() + nc = a%get_ncols() + nza = a%get_nzeros() + + call b%allocate(nr,nc,nza) + call b%psb_s_base_sparse_mat%cp_from(a%psb_s_base_sparse_mat) + + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + b%ia(j) = i + b%ja(j) = a%ja(j) + b%val(j) = a%val(j) + end do + end do + call b%set_nzeros(a%get_nzeros()) + call b%fix(info) + + +end subroutine s_cp_csr_to_coo_impl + + +subroutine s_mv_csr_to_coo_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_mv_csr_to_coo_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + integer, allocatable :: itemp(:) + !locals + logical :: rwshr_ + Integer :: nza, nr, nc,i,j,irw, idl,err_act + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + nr = a%get_nrows() + nc = a%get_ncols() + nza = a%get_nzeros() + + call b%psb_s_base_sparse_mat%mv_from(a%psb_s_base_sparse_mat) + call b%set_nzeros(a%get_nzeros()) + call move_alloc(a%ja,b%ja) + call move_alloc(a%val,b%val) + call psb_realloc(nza,b%ia,info) + if (info /= 0) return + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + b%ia(j) = i + end do + end do + call a%free() + call b%fix(info) + + +end subroutine s_mv_csr_to_coo_impl + + + +subroutine s_mv_csr_from_coo_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_mv_csr_from_coo_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + integer, allocatable :: itemp(:) + !locals + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + call b%fix(info) + if (info /= 0) return + + nr = b%get_nrows() + nc = b%get_ncols() + nza = b%get_nzeros() + + call a%psb_s_base_sparse_mat%mv_from(b%psb_s_base_sparse_mat) + + ! Dirty trick: call move_alloc to have the new data allocated just once. + call move_alloc(b%ia,itemp) + call move_alloc(b%ja,a%ja) + call move_alloc(b%val,a%val) + call psb_realloc(max(nr+1,nc+1),a%irp,info) + call b%free() + + if (nza <= 0) then + a%irp(:) = 1 + else + a%irp(1) = 1 + if (nr < itemp(nza)) then + write(debug_unit,*) trim(name),': RWSHR=.false. : ',& + &nr,itemp(nza),' Expect trouble!' + info = 12 + end if + + j = 1 + i = 1 + irw = itemp(j) + + outer: do + inner: do + if (i >= irw) exit inner + if (i>nr) then + write(debug_unit,*) trim(name),& + & 'Strange situation: i>nr ',i,nr,j,nza,irw,idl + exit outer + end if + a%irp(i+1) = a%irp(i) + i = i + 1 + end do inner + j = j + 1 + if (j > nza) exit + if (itemp(j) /= irw) then + a%irp(i+1) = j + irw = itemp(j) + i = i + 1 + endif + if (i>nr) exit + enddo outer + ! + ! Cleanup empty rows at the end + ! + if (j /= (nza+1)) then + write(debug_unit,*) trim(name),': Problem from loop :',j,nza + info = 13 + endif + do + if (i>nr) exit + a%irp(i+1) = j + i = i + 1 + end do + + endif + + +end subroutine s_mv_csr_from_coo_impl + + +subroutine s_mv_csr_to_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_mv_csr_to_fmt_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + !locals + type(psb_s_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + select type (b) + type is (psb_s_coo_sparse_mat) + call a%mv_to_coo(b,info) + ! Need to fix trivial copies! + type is (psb_s_csr_sparse_mat) + call b%psb_s_base_sparse_mat%mv_from(a%psb_s_base_sparse_mat) + call move_alloc(a%irp, b%irp) + call move_alloc(a%ja, b%ja) + call move_alloc(a%val, b%val) + call a%free() + + class default + call tmp%mv_from_fmt(a,info) + if (info == 0) call b%mv_from_coo(tmp,info) + end select + +end subroutine s_mv_csr_to_fmt_impl + + +subroutine s_cp_csr_to_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_cp_csr_to_fmt_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + !locals + type(psb_s_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + + select type (b) + type is (psb_s_coo_sparse_mat) + call a%cp_to_coo(b,info) + + type is (psb_s_csr_sparse_mat) + b = a + + class default + call tmp%cp_from_fmt(a,info) + if (info == 0) call b%mv_from_coo(tmp,info) + end select + +end subroutine s_cp_csr_to_fmt_impl + + +subroutine s_mv_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_mv_csr_from_fmt_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(inout) :: b + integer, intent(out) :: info + + !locals + type(psb_s_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nza, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + select type (b) + type is (psb_s_coo_sparse_mat) + call a%mv_from_coo(b,info) + + type is (psb_s_csr_sparse_mat) + call a%psb_s_base_sparse_mat%mv_from(b%psb_s_base_sparse_mat) + call move_alloc(b%irp, a%irp) + call move_alloc(b%ja, a%ja) + call move_alloc(b%val, a%val) + call b%free() + + class default + call tmp%mv_from_fmt(b,info) + if (info == 0) call a%mv_from_coo(tmp,info) + end select + +end subroutine s_mv_csr_from_fmt_impl + + + +subroutine s_cp_csr_from_fmt_impl(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_s_base_mat_mod + use psb_s_csr_mat_mod, psb_protect_name => s_cp_csr_from_fmt_impl + implicit none + + class(psb_s_csr_sparse_mat), intent(inout) :: a + class(psb_s_base_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + !locals + type(psb_s_coo_sparse_mat) :: tmp + logical :: rwshr_ + Integer :: nz, nr, i,j,irw, idl,err_act, nc + Integer, Parameter :: maxtry=8 + integer :: debug_level, debug_unit + character(len=20) :: name + + info = 0 + + select type (b) + type is (psb_s_coo_sparse_mat) + call a%cp_from_coo(b,info) + + type is (psb_s_csr_sparse_mat) + call a%psb_s_base_sparse_mat%cp_from(b%psb_s_base_sparse_mat) + a%irp = b%irp + a%ja = b%ja + a%val = b%val + + class default + call tmp%cp_from_fmt(b,info) + if (info == 0) call a%mv_from_coo(tmp,info) + end select +end subroutine s_cp_csr_from_fmt_impl +