From e766e9b2fbe524d8a8acd30393190e078deaed02 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sun, 9 Jun 2019 19:49:23 +0100 Subject: [PATCH] Implementation of CSRD. To be tested. --- base/modules/serial/psb_c_csr_mat_mod.f90 | 156 ++++ base/modules/serial/psb_d_csr_mat_mod.f90 | 156 ++++ base/modules/serial/psb_s_csr_mat_mod.f90 | 156 ++++ base/modules/serial/psb_z_csr_mat_mod.f90 | 156 ++++ base/serial/impl/psb_c_base_mat_impl.F90 | 49 +- base/serial/impl/psb_c_csr_impl.f90 | 1013 ++++++++++++++++++++- base/serial/impl/psb_d_base_mat_impl.F90 | 49 +- base/serial/impl/psb_d_csr_impl.f90 | 1013 ++++++++++++++++++++- base/serial/impl/psb_s_base_mat_impl.F90 | 49 +- base/serial/impl/psb_s_csr_impl.f90 | 1013 ++++++++++++++++++++- base/serial/impl/psb_z_base_mat_impl.F90 | 49 +- base/serial/impl/psb_z_csr_impl.f90 | 1013 ++++++++++++++++++++- 12 files changed, 4812 insertions(+), 60 deletions(-) diff --git a/base/modules/serial/psb_c_csr_mat_mod.f90 b/base/modules/serial/psb_c_csr_mat_mod.f90 index 7452d5aa..0cbe5218 100644 --- a/base/modules/serial/psb_c_csr_mat_mod.f90 +++ b/base/modules/serial/psb_c_csr_mat_mod.f90 @@ -596,8 +596,120 @@ module psb_c_csr_mat_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_c_csr_scals end interface + + type, extends(psb_c_csr_sparse_mat) :: psb_c_csrd_sparse_mat + + !> Pointers to diagonal entries + integer(psb_ipk_), allocatable :: diagp(:) + + contains + procedure, nopass :: get_fmt => c_csrd_get_fmt + procedure, pass(a) :: sizeof => c_csrd_sizeof + procedure, pass(a) :: inner_cssm => psb_c_csrd_cssm + procedure, pass(a) :: inner_cssv => psb_c_csrd_cssv + procedure, pass(a) :: trmv => psb_c_csrd_trmv + !procedure, pass(a) :: reallocate_nz => psb_c_csrd_reallocate_nz + !procedure, pass(a) :: allocate_mnnz => psb_c_csrd_allocate_mnnz +!!$ procedure, pass(a) :: cp_to_coo => psb_c_cp_csrd_to_coo + procedure, pass(a) :: cp_from_coo => psb_c_cp_csrd_from_coo +!!$ procedure, pass(a) :: cp_to_fmt => psb_c_cp_csrd_to_fmt +!!$ procedure, pass(a) :: cp_from_fmt => psb_c_cp_csrd_from_fmt +!!$ procedure, pass(a) :: mv_to_coo => psb_c_mv_csrd_to_coo +!!$ procedure, pass(a) :: mv_from_coo => psb_c_mv_csrd_from_coo +!!$ procedure, pass(a) :: mv_to_fmt => psb_c_mv_csrd_to_fmt +!!$ procedure, pass(a) :: mv_from_fmt => psb_c_mv_csrd_from_fmt + procedure, pass(a) :: clean_zeros => psb_c_csrd_clean_zeros + procedure, pass(a) :: get_diag => psb_c_csrd_get_diag + !procedure, pass(a) :: reinit => psb_c_csrd_reinit + procedure, pass(a) :: trim => psb_c_csrd_trim + procedure, pass(a) :: free => c_csrd_free + procedure, pass(a) :: mold => psb_c_csrd_mold + + end type psb_c_csrd_sparse_mat + + interface + subroutine psb_c_csrd_clean_zeros(a, info) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_csrd_clean_zeros + end interface + + interface + subroutine psb_c_csrd_get_diag(a,d,info) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_csrd_get_diag + end interface + + interface + subroutine psb_c_cp_csrd_from_coo(a,b,info) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_c_coo_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + class(psb_c_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_cp_csrd_from_coo + end interface + + interface + subroutine psb_c_csrd_trim(a) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + end subroutine psb_c_csrd_trim + end interface + interface + subroutine psb_c_csrd_mold(a,b,info) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_c_base_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + class(psb_c_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_csrd_mold + end interface + + interface + subroutine psb_c_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:) + complex(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + end subroutine psb_c_csrd_trmv + end interface + interface + subroutine psb_c_csrd_cssm(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:,:) + complex(psb_spk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_c_csrd_cssm + end interface + + interface + subroutine psb_c_csrd_cssv(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_c_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:) + complex(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_c_csrd_cssv + end interface + contains @@ -717,4 +829,48 @@ contains end subroutine c_csr_free + ! == =================================== + ! + ! + ! + ! CSRD specific versions + ! + ! + ! + ! + ! + ! == =================================== + + function c_csrd_sizeof(a) result(res) + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + (2*psb_sizeof_sp) * psb_size(a%val) + res = res + psb_sizeof_int * psb_size(a%irp) + res = res + psb_sizeof_int * psb_size(a%ja) + res = res + psb_sizeof_int * psb_size(a%diagp) + + end function c_csrd_sizeof + + function c_csrd_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'CSRD' + end function c_csrd_get_fmt + + + subroutine c_csrd_free(a) + implicit none + + class(psb_c_csrd_sparse_mat), intent(inout) :: a + + if (allocated(a%diagp)) deallocate(a%diagp) + call a%psb_c_csr_sparse_mat%free() + + return + + end subroutine c_csrd_free + + end module psb_c_csr_mat_mod diff --git a/base/modules/serial/psb_d_csr_mat_mod.f90 b/base/modules/serial/psb_d_csr_mat_mod.f90 index ae5ddfd3..8b573652 100644 --- a/base/modules/serial/psb_d_csr_mat_mod.f90 +++ b/base/modules/serial/psb_d_csr_mat_mod.f90 @@ -596,8 +596,120 @@ module psb_d_csr_mat_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_d_csr_scals end interface + + type, extends(psb_d_csr_sparse_mat) :: psb_d_csrd_sparse_mat + + !> Pointers to diagonal entries + integer(psb_ipk_), allocatable :: diagp(:) + + contains + procedure, nopass :: get_fmt => d_csrd_get_fmt + procedure, pass(a) :: sizeof => d_csrd_sizeof + procedure, pass(a) :: inner_cssm => psb_d_csrd_cssm + procedure, pass(a) :: inner_cssv => psb_d_csrd_cssv + procedure, pass(a) :: trmv => psb_d_csrd_trmv + !procedure, pass(a) :: reallocate_nz => psb_d_csrd_reallocate_nz + !procedure, pass(a) :: allocate_mnnz => psb_d_csrd_allocate_mnnz +!!$ procedure, pass(a) :: cp_to_coo => psb_d_cp_csrd_to_coo + procedure, pass(a) :: cp_from_coo => psb_d_cp_csrd_from_coo +!!$ procedure, pass(a) :: cp_to_fmt => psb_d_cp_csrd_to_fmt +!!$ procedure, pass(a) :: cp_from_fmt => psb_d_cp_csrd_from_fmt +!!$ procedure, pass(a) :: mv_to_coo => psb_d_mv_csrd_to_coo +!!$ procedure, pass(a) :: mv_from_coo => psb_d_mv_csrd_from_coo +!!$ procedure, pass(a) :: mv_to_fmt => psb_d_mv_csrd_to_fmt +!!$ procedure, pass(a) :: mv_from_fmt => psb_d_mv_csrd_from_fmt + procedure, pass(a) :: clean_zeros => psb_d_csrd_clean_zeros + procedure, pass(a) :: get_diag => psb_d_csrd_get_diag + !procedure, pass(a) :: reinit => psb_d_csrd_reinit + procedure, pass(a) :: trim => psb_d_csrd_trim + procedure, pass(a) :: free => d_csrd_free + procedure, pass(a) :: mold => psb_d_csrd_mold + + end type psb_d_csrd_sparse_mat + + interface + subroutine psb_d_csrd_clean_zeros(a, info) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_csrd_clean_zeros + end interface + + interface + subroutine psb_d_csrd_get_diag(a,d,info) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_csrd_get_diag + end interface + + interface + subroutine psb_d_cp_csrd_from_coo(a,b,info) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_d_coo_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + class(psb_d_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_cp_csrd_from_coo + end interface + + interface + subroutine psb_d_csrd_trim(a) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + end subroutine psb_d_csrd_trim + end interface + interface + subroutine psb_d_csrd_mold(a,b,info) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_d_base_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + class(psb_d_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_csrd_mold + end interface + + interface + subroutine psb_d_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + end subroutine psb_d_csrd_trmv + end interface + interface + subroutine psb_d_csrd_cssm(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_d_csrd_cssm + end interface + + interface + subroutine psb_d_csrd_cssv(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_d_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_d_csrd_cssv + end interface + contains @@ -717,4 +829,48 @@ contains end subroutine d_csr_free + ! == =================================== + ! + ! + ! + ! CSRD specific versions + ! + ! + ! + ! + ! + ! == =================================== + + function d_csrd_sizeof(a) result(res) + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + psb_sizeof_dp * psb_size(a%val) + res = res + psb_sizeof_int * psb_size(a%irp) + res = res + psb_sizeof_int * psb_size(a%ja) + res = res + psb_sizeof_int * psb_size(a%diagp) + + end function d_csrd_sizeof + + function d_csrd_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'CSRD' + end function d_csrd_get_fmt + + + subroutine d_csrd_free(a) + implicit none + + class(psb_d_csrd_sparse_mat), intent(inout) :: a + + if (allocated(a%diagp)) deallocate(a%diagp) + call a%psb_d_csr_sparse_mat%free() + + return + + end subroutine d_csrd_free + + end module psb_d_csr_mat_mod diff --git a/base/modules/serial/psb_s_csr_mat_mod.f90 b/base/modules/serial/psb_s_csr_mat_mod.f90 index e48fb7c1..2432481c 100644 --- a/base/modules/serial/psb_s_csr_mat_mod.f90 +++ b/base/modules/serial/psb_s_csr_mat_mod.f90 @@ -596,8 +596,120 @@ module psb_s_csr_mat_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_s_csr_scals end interface + + type, extends(psb_s_csr_sparse_mat) :: psb_s_csrd_sparse_mat + + !> Pointers to diagonal entries + integer(psb_ipk_), allocatable :: diagp(:) + + contains + procedure, nopass :: get_fmt => s_csrd_get_fmt + procedure, pass(a) :: sizeof => s_csrd_sizeof + procedure, pass(a) :: inner_cssm => psb_s_csrd_cssm + procedure, pass(a) :: inner_cssv => psb_s_csrd_cssv + procedure, pass(a) :: trmv => psb_s_csrd_trmv + !procedure, pass(a) :: reallocate_nz => psb_s_csrd_reallocate_nz + !procedure, pass(a) :: allocate_mnnz => psb_s_csrd_allocate_mnnz +!!$ procedure, pass(a) :: cp_to_coo => psb_s_cp_csrd_to_coo + procedure, pass(a) :: cp_from_coo => psb_s_cp_csrd_from_coo +!!$ procedure, pass(a) :: cp_to_fmt => psb_s_cp_csrd_to_fmt +!!$ procedure, pass(a) :: cp_from_fmt => psb_s_cp_csrd_from_fmt +!!$ procedure, pass(a) :: mv_to_coo => psb_s_mv_csrd_to_coo +!!$ procedure, pass(a) :: mv_from_coo => psb_s_mv_csrd_from_coo +!!$ procedure, pass(a) :: mv_to_fmt => psb_s_mv_csrd_to_fmt +!!$ procedure, pass(a) :: mv_from_fmt => psb_s_mv_csrd_from_fmt + procedure, pass(a) :: clean_zeros => psb_s_csrd_clean_zeros + procedure, pass(a) :: get_diag => psb_s_csrd_get_diag + !procedure, pass(a) :: reinit => psb_s_csrd_reinit + procedure, pass(a) :: trim => psb_s_csrd_trim + procedure, pass(a) :: free => s_csrd_free + procedure, pass(a) :: mold => psb_s_csrd_mold + + end type psb_s_csrd_sparse_mat + + interface + subroutine psb_s_csrd_clean_zeros(a, info) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_csrd_clean_zeros + end interface + + interface + subroutine psb_s_csrd_get_diag(a,d,info) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_csrd_get_diag + end interface + + interface + subroutine psb_s_cp_csrd_from_coo(a,b,info) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_s_coo_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_cp_csrd_from_coo + end interface + + interface + subroutine psb_s_csrd_trim(a) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + end subroutine psb_s_csrd_trim + end interface + interface + subroutine psb_s_csrd_mold(a,b,info) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_s_base_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_csrd_mold + end interface + + interface + subroutine psb_s_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + end subroutine psb_s_csrd_trmv + end interface + interface + subroutine psb_s_csrd_cssm(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_s_csrd_cssm + end interface + + interface + subroutine psb_s_csrd_cssv(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_s_csrd_sparse_mat, psb_spk_ + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_s_csrd_cssv + end interface + contains @@ -717,4 +829,48 @@ contains end subroutine s_csr_free + ! == =================================== + ! + ! + ! + ! CSRD specific versions + ! + ! + ! + ! + ! + ! == =================================== + + function s_csrd_sizeof(a) result(res) + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + psb_sizeof_sp * psb_size(a%val) + res = res + psb_sizeof_int * psb_size(a%irp) + res = res + psb_sizeof_int * psb_size(a%ja) + res = res + psb_sizeof_int * psb_size(a%diagp) + + end function s_csrd_sizeof + + function s_csrd_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'CSRD' + end function s_csrd_get_fmt + + + subroutine s_csrd_free(a) + implicit none + + class(psb_s_csrd_sparse_mat), intent(inout) :: a + + if (allocated(a%diagp)) deallocate(a%diagp) + call a%psb_s_csr_sparse_mat%free() + + return + + end subroutine s_csrd_free + + end module psb_s_csr_mat_mod diff --git a/base/modules/serial/psb_z_csr_mat_mod.f90 b/base/modules/serial/psb_z_csr_mat_mod.f90 index 8312d960..341dc77a 100644 --- a/base/modules/serial/psb_z_csr_mat_mod.f90 +++ b/base/modules/serial/psb_z_csr_mat_mod.f90 @@ -596,8 +596,120 @@ module psb_z_csr_mat_mod integer(psb_ipk_), intent(out) :: info end subroutine psb_z_csr_scals end interface + + type, extends(psb_z_csr_sparse_mat) :: psb_z_csrd_sparse_mat + + !> Pointers to diagonal entries + integer(psb_ipk_), allocatable :: diagp(:) + + contains + procedure, nopass :: get_fmt => z_csrd_get_fmt + procedure, pass(a) :: sizeof => z_csrd_sizeof + procedure, pass(a) :: inner_cssm => psb_z_csrd_cssm + procedure, pass(a) :: inner_cssv => psb_z_csrd_cssv + procedure, pass(a) :: trmv => psb_z_csrd_trmv + !procedure, pass(a) :: reallocate_nz => psb_z_csrd_reallocate_nz + !procedure, pass(a) :: allocate_mnnz => psb_z_csrd_allocate_mnnz +!!$ procedure, pass(a) :: cp_to_coo => psb_z_cp_csrd_to_coo + procedure, pass(a) :: cp_from_coo => psb_z_cp_csrd_from_coo +!!$ procedure, pass(a) :: cp_to_fmt => psb_z_cp_csrd_to_fmt +!!$ procedure, pass(a) :: cp_from_fmt => psb_z_cp_csrd_from_fmt +!!$ procedure, pass(a) :: mv_to_coo => psb_z_mv_csrd_to_coo +!!$ procedure, pass(a) :: mv_from_coo => psb_z_mv_csrd_from_coo +!!$ procedure, pass(a) :: mv_to_fmt => psb_z_mv_csrd_to_fmt +!!$ procedure, pass(a) :: mv_from_fmt => psb_z_mv_csrd_from_fmt + procedure, pass(a) :: clean_zeros => psb_z_csrd_clean_zeros + procedure, pass(a) :: get_diag => psb_z_csrd_get_diag + !procedure, pass(a) :: reinit => psb_z_csrd_reinit + procedure, pass(a) :: trim => psb_z_csrd_trim + procedure, pass(a) :: free => z_csrd_free + procedure, pass(a) :: mold => psb_z_csrd_mold + + end type psb_z_csrd_sparse_mat + + interface + subroutine psb_z_csrd_clean_zeros(a, info) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_csrd_clean_zeros + end interface + + interface + subroutine psb_z_csrd_get_diag(a,d,info) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_csrd_get_diag + end interface + + interface + subroutine psb_z_cp_csrd_from_coo(a,b,info) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_z_coo_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + class(psb_z_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_cp_csrd_from_coo + end interface + + interface + subroutine psb_z_csrd_trim(a) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + end subroutine psb_z_csrd_trim + end interface + interface + subroutine psb_z_csrd_mold(a,b,info) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_z_base_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + class(psb_z_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_csrd_mold + end interface + + interface + subroutine psb_z_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:) + complex(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + end subroutine psb_z_csrd_trmv + end interface + interface + subroutine psb_z_csrd_cssm(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + complex(psb_dpk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_z_csrd_cssm + end interface + + interface + subroutine psb_z_csrd_cssv(alpha,a,x,beta,y,info,trans) + import :: psb_ipk_, psb_z_csrd_sparse_mat, psb_dpk_ + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:) + complex(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + end subroutine psb_z_csrd_cssv + end interface + contains @@ -717,4 +829,48 @@ contains end subroutine z_csr_free + ! == =================================== + ! + ! + ! + ! CSRD specific versions + ! + ! + ! + ! + ! + ! == =================================== + + function z_csrd_sizeof(a) result(res) + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + integer(psb_long_int_k_) :: res + res = 8 + res = res + (2*psb_sizeof_dp) * psb_size(a%val) + res = res + psb_sizeof_int * psb_size(a%irp) + res = res + psb_sizeof_int * psb_size(a%ja) + res = res + psb_sizeof_int * psb_size(a%diagp) + + end function z_csrd_sizeof + + function z_csrd_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'CSRD' + end function z_csrd_get_fmt + + + subroutine z_csrd_free(a) + implicit none + + class(psb_z_csrd_sparse_mat), intent(inout) :: a + + if (allocated(a%diagp)) deallocate(a%diagp) + call a%psb_z_csr_sparse_mat%free() + + return + + end subroutine z_csrd_free + + end module psb_z_csr_mat_mod diff --git a/base/serial/impl/psb_c_base_mat_impl.F90 b/base/serial/impl/psb_c_base_mat_impl.F90 index 0ced6dda..ac9ac82b 100644 --- a/base/serial/impl/psb_c_base_mat_impl.F90 +++ b/base/serial/impl/psb_c_base_mat_impl.F90 @@ -1184,17 +1184,28 @@ subroutine psb_c_base_csmm(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: err_act integer(psb_ipk_) :: ierr(5) - character(len=20) :: name='c_base_csmm' + integer(psb_ipk_) :: j,nc + character(len=20) :: name='c_base_csmm' logical, parameter :: debug=.false. call psb_erractionsave(err_act) - ! This is the base version. If we get here - ! it means the derived class is incomplete, - ! so we throw an error. - info = psb_err_missing_override_method_ - call psb_errpush(info,name,a_err=a%get_fmt()) + ! This is the base version. + ! It's a very inefficient implementation, + ! but it's only a fallback, if multivectors + ! are important you are supposed to implement it + ! explicitly in the derived class. + info = psb_success_ + nc = min(size(x,2),size(y,2)) + do j=1,nc + call a%spmm(alpha,x(j,:),beta,y(:,j),info,trans) + if (info /= psb_success_) goto 9999 + end do - call psb_error_handler(err_act) + call psb_erractionrestore(err_act) + + return + +9999 call psb_error_handler(err_act) end subroutine psb_c_base_csmm @@ -1266,17 +1277,29 @@ subroutine psb_c_base_inner_cssm(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: err_act integer(psb_ipk_) :: ierr(5) + integer(psb_ipk_) :: j, nc character(len=20) :: name='c_base_inner_cssm' logical, parameter :: debug=.false. call psb_erractionsave(err_act) - ! This is the base version. If we get here - ! it means the derived class is incomplete, - ! so we throw an error. - info = psb_err_missing_override_method_ - call psb_errpush(info,name,a_err=a%get_fmt()) + ! This is the base version. + ! It's a very inefficient implementation, + ! but it's only a fallback, if multivectors + ! are important you are supposed to implement it + ! explicitly in the derived class. + info = psb_success_ + nc = min(size(x,2),size(y,2)) + do j=1,nc + call a%spsm(alpha,x(j,:),beta,y(:,j),info,trans) + if (info /= psb_success_) goto 9999 + end do + + call psb_erractionrestore(err_act) + + return + +9999 call psb_error_handler(err_act) - call psb_error_handler(err_act) end subroutine psb_c_base_inner_cssm diff --git a/base/serial/impl/psb_c_csr_impl.f90 b/base/serial/impl/psb_c_csr_impl.f90 index e248c5e1..138d2963 100644 --- a/base/serial/impl/psb_c_csr_impl.f90 +++ b/base/serial/impl/psb_c_csr_impl.f90 @@ -999,8 +999,6 @@ contains end subroutine psb_c_csr_cssv - - subroutine psb_c_csr_cssm(alpha,a,x,beta,y,info,trans) use psb_error_mod use psb_string_mod @@ -3473,3 +3471,1014 @@ contains end subroutine csr_spspmm end subroutine psb_ccsrspspmm + + +! == =================================== +! +! +! +! CSRD specific versions +! +! +! +! +! +! == =================================== + +subroutine psb_c_csrd_clean_zeros(a, info) + use psb_error_mod + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_clean_zeros + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, k, nr + integer(psb_ipk_), allocatable :: ilrp(:) + + info = 0 + call a%sync() + nr = a%get_nrows() + ilrp = a%irp(:) + a%irp(1) = 1 + j = a%irp(1) + do i=1, nr + a%diagp(i) = 0 + do k = ilrp(i), ilrp(i+1) -1 + if ((a%ja(j)==i).or.(a%val(k) /= czero)) then + a%val(j) = a%val(k) + a%ja(j) = a%ja(k) + if (a%ja(j) == i) a%diagp(i) = j + j = j + 1 + end if + end do + a%irp(i+1) = j + end do + call a%trim() + call a%set_host() +end subroutine psb_c_csrd_clean_zeros + + +subroutine psb_c_csrd_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_get_diag + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + if (a%is_unit()) then + d(1:mnm) = cone + else + do i=1, mnm + if (a%diagp(i) >0) d(i) = a%val(a%diagp(i)) + end do + end if + do i=mnm+1,size(d) + d(i) = czero + end do + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_c_csrd_get_diag + +subroutine psb_c_cp_csrd_from_coo(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_c_base_mat_mod + use psb_c_csr_mat_mod, psb_protect_name => psb_c_cp_csrd_from_coo + implicit none + + class(psb_c_csrd_sparse_mat), intent(inout) :: a + class(psb_c_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + + type(psb_c_coo_sparse_mat) :: tmp + integer(psb_ipk_), allocatable :: itemp(:) + !locals + logical :: nodiag + integer(psb_ipk_) :: nza, nzt, nr, nc, i,j,k,ip,irw, err_act, ncl + integer(psb_ipk_), Parameter :: maxtry=8 + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name='c_cp_csrd_from_coo' + + info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + call a%psb_c_csr_sparse_mat%cp_from_coo(b,info) + + if (info /= 0) return + nr = a%get_nrows() + call psb_realloc(nr,a%diagp,info) + if (info /= 0) return + ! + call srch_diag(a,nodiag) + if (nodiag) then + ! + ! This means at least one of the rows does not + ! have a diagonal entry. Add explicit zeros and fix. + ! + call a%cp_to_coo(tmp,info) + nzt = tmp%get_nzeros() + call tmp%reallocate(nzt+nr) + do k=1,nr + tmp%ia(nzt+k) = k + tmp%ja(nzt+k) = k + tmp%val(nzt+k) = czero + end do + call tmp%set_nzeros(nzt+nr) + call tmp%fix(info) + call a%mv_from_coo(tmp,info) + call srch_diag(a,nodiag) + ! At this point, nodiag should not happen + if (nodiag) then + write(0,*) 'Nodiag here?' + info = -1 + end if + end if + +contains + subroutine srch_diag(a,nodiag) + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + logical, intent(out) :: nodiag + integer(psb_ipk_) :: nr, i,j,k + + nr = a%get_nrows() + nodiag = .false. + outer: do i=1, nr + a%diagp(i) = -1 + inner: do k=a%irp(i),a%irp(i+1)-1 + if (a%ja(k) == i) then + a%diagp(i) = k + exit inner + end if + end do inner + if (a%diagp(i) < 1) then + nodiag = .true. + exit outer + end if + end do outer + return + end subroutine srch_diag +end subroutine psb_c_cp_csrd_from_coo + +subroutine psb_c_csrd_trim(a) + use psb_realloc_mod + use psb_error_mod + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_trim + implicit none + class(psb_c_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + call a%psb_c_csr_sparse_mat%trim() + m = a%get_nrows() + call psb_realloc(m,a%diagp,info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_c_csrd_trim + + +subroutine psb_c_csrd_mold(a,b,info) + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_mold + use psb_error_mod + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + class(psb_c_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='csrd_mold' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + info = 0 + if (allocated(b)) then + call b%free() + deallocate(b,stat=info) + end if + if (info == 0) allocate(psb_c_csrd_sparse_mat :: b, stat=info) + + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info, name) + goto 9999 + end if + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_c_csrd_mold + +subroutine psb_c_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_trmv + use psb_error_mod + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:) + complex(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + + integer(psb_ipk_) :: i, m, n + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='c_csrd_trmv' + character :: uplo_, diag_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + if (a%is_dev()) call a%sync() + + if (present(uplo)) then + uplo_ = uplo + else + uplo_ = 'L' + end if + if (present(diag)) then + diag_ = diag + else + diag_ = 'D' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + n = a%get_ncols() + m = a%get_nrows() + + if (size(x,1) psb_c_csrd_cssm + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:,:) + complex(psb_spk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + + integer(psb_ipk_) :: jc, nc + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='c_csrd_cssm' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + + nc = min(size(x,2) , size(y,2)) + + do jc = 1, nc + call a%inner_cssv(alpha,x(:,jc),beta,y(:,jc),info,trans) + if (info /= 0) exit + end do + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='inner_trmv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_c_csrd_cssm + + +subroutine psb_c_csrd_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_c_csr_mat_mod, psb_protect_name => psb_c_csrd_cssv + implicit none + class(psb_c_csrd_sparse_mat), intent(in) :: a + complex(psb_spk_), intent(in) :: alpha, beta, x(:) + complex(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc + complex(psb_spk_) :: acc + complex(psb_spk_), allocatable :: tmp(:) + logical :: tra,ctra + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='c_csrd_cssv' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + !info = psb_err_invalid_mat_state_ + !call psb_errpush(info,name) + !goto 9999 + ! + ! The whole point of CSRD is to fake a triangle + ! "on the fly" (this is really for the sake of + ! saving memory in Gauss-Seidel lower/upper). + ! + ! + end if + + if (size(x) psb_d_csrd_clean_zeros + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, k, nr + integer(psb_ipk_), allocatable :: ilrp(:) + + info = 0 + call a%sync() + nr = a%get_nrows() + ilrp = a%irp(:) + a%irp(1) = 1 + j = a%irp(1) + do i=1, nr + a%diagp(i) = 0 + do k = ilrp(i), ilrp(i+1) -1 + if ((a%ja(j)==i).or.(a%val(k) /= dzero)) then + a%val(j) = a%val(k) + a%ja(j) = a%ja(k) + if (a%ja(j) == i) a%diagp(i) = j + j = j + 1 + end if + end do + a%irp(i+1) = j + end do + call a%trim() + call a%set_host() +end subroutine psb_d_csrd_clean_zeros + + +subroutine psb_d_csrd_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + use psb_d_csr_mat_mod, psb_protect_name => psb_d_csrd_get_diag + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + if (a%is_unit()) then + d(1:mnm) = done + else + do i=1, mnm + if (a%diagp(i) >0) d(i) = a%val(a%diagp(i)) + end do + end if + do i=mnm+1,size(d) + d(i) = dzero + end do + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_d_csrd_get_diag + +subroutine psb_d_cp_csrd_from_coo(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_d_base_mat_mod + use psb_d_csr_mat_mod, psb_protect_name => psb_d_cp_csrd_from_coo + implicit none + + class(psb_d_csrd_sparse_mat), intent(inout) :: a + class(psb_d_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + + type(psb_d_coo_sparse_mat) :: tmp + integer(psb_ipk_), allocatable :: itemp(:) + !locals + logical :: nodiag + integer(psb_ipk_) :: nza, nzt, nr, nc, i,j,k,ip,irw, err_act, ncl + integer(psb_ipk_), Parameter :: maxtry=8 + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name='d_cp_csrd_from_coo' + + info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + call a%psb_d_csr_sparse_mat%cp_from_coo(b,info) + + if (info /= 0) return + nr = a%get_nrows() + call psb_realloc(nr,a%diagp,info) + if (info /= 0) return + ! + call srch_diag(a,nodiag) + if (nodiag) then + ! + ! This means at least one of the rows does not + ! have a diagonal entry. Add explicit zeros and fix. + ! + call a%cp_to_coo(tmp,info) + nzt = tmp%get_nzeros() + call tmp%reallocate(nzt+nr) + do k=1,nr + tmp%ia(nzt+k) = k + tmp%ja(nzt+k) = k + tmp%val(nzt+k) = dzero + end do + call tmp%set_nzeros(nzt+nr) + call tmp%fix(info) + call a%mv_from_coo(tmp,info) + call srch_diag(a,nodiag) + ! At this point, nodiag should not happen + if (nodiag) then + write(0,*) 'Nodiag here?' + info = -1 + end if + end if + +contains + subroutine srch_diag(a,nodiag) + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + logical, intent(out) :: nodiag + integer(psb_ipk_) :: nr, i,j,k + + nr = a%get_nrows() + nodiag = .false. + outer: do i=1, nr + a%diagp(i) = -1 + inner: do k=a%irp(i),a%irp(i+1)-1 + if (a%ja(k) == i) then + a%diagp(i) = k + exit inner + end if + end do inner + if (a%diagp(i) < 1) then + nodiag = .true. + exit outer + end if + end do outer + return + end subroutine srch_diag +end subroutine psb_d_cp_csrd_from_coo + +subroutine psb_d_csrd_trim(a) + use psb_realloc_mod + use psb_error_mod + use psb_d_csr_mat_mod, psb_protect_name => psb_d_csrd_trim + implicit none + class(psb_d_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + call a%psb_d_csr_sparse_mat%trim() + m = a%get_nrows() + call psb_realloc(m,a%diagp,info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_d_csrd_trim + + +subroutine psb_d_csrd_mold(a,b,info) + use psb_d_csr_mat_mod, psb_protect_name => psb_d_csrd_mold + use psb_error_mod + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + class(psb_d_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='csrd_mold' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + info = 0 + if (allocated(b)) then + call b%free() + deallocate(b,stat=info) + end if + if (info == 0) allocate(psb_d_csrd_sparse_mat :: b, stat=info) + + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info, name) + goto 9999 + end if + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_d_csrd_mold + +subroutine psb_d_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + use psb_d_csr_mat_mod, psb_protect_name => psb_d_csrd_trmv + use psb_error_mod + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + + integer(psb_ipk_) :: i, m, n + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='d_csrd_trmv' + character :: uplo_, diag_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + if (a%is_dev()) call a%sync() + + if (present(uplo)) then + uplo_ = uplo + else + uplo_ = 'L' + end if + if (present(diag)) then + diag_ = diag + else + diag_ = 'D' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + n = a%get_ncols() + m = a%get_nrows() + + if (size(x,1) psb_d_csrd_cssm + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + + integer(psb_ipk_) :: jc, nc + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='d_csrd_cssm' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + + nc = min(size(x,2) , size(y,2)) + + do jc = 1, nc + call a%inner_cssv(alpha,x(:,jc),beta,y(:,jc),info,trans) + if (info /= 0) exit + end do + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='inner_trmv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_d_csrd_cssm + + +subroutine psb_d_csrd_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_d_csr_mat_mod, psb_protect_name => psb_d_csrd_cssv + implicit none + class(psb_d_csrd_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:) + logical :: tra,ctra + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='d_csrd_cssv' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + !info = psb_err_invalid_mat_state_ + !call psb_errpush(info,name) + !goto 9999 + ! + ! The whole point of CSRD is to fake a triangle + ! "on the fly" (this is really for the sake of + ! saving memory in Gauss-Seidel lower/upper). + ! + ! + end if + + if (size(x) psb_s_csrd_clean_zeros + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, k, nr + integer(psb_ipk_), allocatable :: ilrp(:) + + info = 0 + call a%sync() + nr = a%get_nrows() + ilrp = a%irp(:) + a%irp(1) = 1 + j = a%irp(1) + do i=1, nr + a%diagp(i) = 0 + do k = ilrp(i), ilrp(i+1) -1 + if ((a%ja(j)==i).or.(a%val(k) /= szero)) then + a%val(j) = a%val(k) + a%ja(j) = a%ja(k) + if (a%ja(j) == i) a%diagp(i) = j + j = j + 1 + end if + end do + a%irp(i+1) = j + end do + call a%trim() + call a%set_host() +end subroutine psb_s_csrd_clean_zeros + + +subroutine psb_s_csrd_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + use psb_s_csr_mat_mod, psb_protect_name => psb_s_csrd_get_diag + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + if (a%is_unit()) then + d(1:mnm) = sone + else + do i=1, mnm + if (a%diagp(i) >0) d(i) = a%val(a%diagp(i)) + end do + end if + do i=mnm+1,size(d) + d(i) = szero + end do + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_s_csrd_get_diag + +subroutine psb_s_cp_csrd_from_coo(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 => psb_s_cp_csrd_from_coo + implicit none + + class(psb_s_csrd_sparse_mat), intent(inout) :: a + class(psb_s_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + + type(psb_s_coo_sparse_mat) :: tmp + integer(psb_ipk_), allocatable :: itemp(:) + !locals + logical :: nodiag + integer(psb_ipk_) :: nza, nzt, nr, nc, i,j,k,ip,irw, err_act, ncl + integer(psb_ipk_), Parameter :: maxtry=8 + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name='s_cp_csrd_from_coo' + + info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + call a%psb_s_csr_sparse_mat%cp_from_coo(b,info) + + if (info /= 0) return + nr = a%get_nrows() + call psb_realloc(nr,a%diagp,info) + if (info /= 0) return + ! + call srch_diag(a,nodiag) + if (nodiag) then + ! + ! This means at least one of the rows does not + ! have a diagonal entry. Add explicit zeros and fix. + ! + call a%cp_to_coo(tmp,info) + nzt = tmp%get_nzeros() + call tmp%reallocate(nzt+nr) + do k=1,nr + tmp%ia(nzt+k) = k + tmp%ja(nzt+k) = k + tmp%val(nzt+k) = szero + end do + call tmp%set_nzeros(nzt+nr) + call tmp%fix(info) + call a%mv_from_coo(tmp,info) + call srch_diag(a,nodiag) + ! At this point, nodiag should not happen + if (nodiag) then + write(0,*) 'Nodiag here?' + info = -1 + end if + end if + +contains + subroutine srch_diag(a,nodiag) + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + logical, intent(out) :: nodiag + integer(psb_ipk_) :: nr, i,j,k + + nr = a%get_nrows() + nodiag = .false. + outer: do i=1, nr + a%diagp(i) = -1 + inner: do k=a%irp(i),a%irp(i+1)-1 + if (a%ja(k) == i) then + a%diagp(i) = k + exit inner + end if + end do inner + if (a%diagp(i) < 1) then + nodiag = .true. + exit outer + end if + end do outer + return + end subroutine srch_diag +end subroutine psb_s_cp_csrd_from_coo + +subroutine psb_s_csrd_trim(a) + use psb_realloc_mod + use psb_error_mod + use psb_s_csr_mat_mod, psb_protect_name => psb_s_csrd_trim + implicit none + class(psb_s_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + call a%psb_s_csr_sparse_mat%trim() + m = a%get_nrows() + call psb_realloc(m,a%diagp,info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_s_csrd_trim + + +subroutine psb_s_csrd_mold(a,b,info) + use psb_s_csr_mat_mod, psb_protect_name => psb_s_csrd_mold + use psb_error_mod + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + class(psb_s_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='csrd_mold' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + info = 0 + if (allocated(b)) then + call b%free() + deallocate(b,stat=info) + end if + if (info == 0) allocate(psb_s_csrd_sparse_mat :: b, stat=info) + + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info, name) + goto 9999 + end if + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_s_csrd_mold + +subroutine psb_s_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + use psb_s_csr_mat_mod, psb_protect_name => psb_s_csrd_trmv + use psb_error_mod + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + + integer(psb_ipk_) :: i, m, n + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='s_csrd_trmv' + character :: uplo_, diag_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + if (a%is_dev()) call a%sync() + + if (present(uplo)) then + uplo_ = uplo + else + uplo_ = 'L' + end if + if (present(diag)) then + diag_ = diag + else + diag_ = 'D' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + n = a%get_ncols() + m = a%get_nrows() + + if (size(x,1) psb_s_csrd_cssm + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:,:) + real(psb_spk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + + integer(psb_ipk_) :: jc, nc + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='s_csrd_cssm' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + + nc = min(size(x,2) , size(y,2)) + + do jc = 1, nc + call a%inner_cssv(alpha,x(:,jc),beta,y(:,jc),info,trans) + if (info /= 0) exit + end do + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='inner_trmv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_s_csrd_cssm + + +subroutine psb_s_csrd_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_s_csr_mat_mod, psb_protect_name => psb_s_csrd_cssv + implicit none + class(psb_s_csrd_sparse_mat), intent(in) :: a + real(psb_spk_), intent(in) :: alpha, beta, x(:) + real(psb_spk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc + real(psb_spk_) :: acc + real(psb_spk_), allocatable :: tmp(:) + logical :: tra,ctra + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='s_csrd_cssv' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + !info = psb_err_invalid_mat_state_ + !call psb_errpush(info,name) + !goto 9999 + ! + ! The whole point of CSRD is to fake a triangle + ! "on the fly" (this is really for the sake of + ! saving memory in Gauss-Seidel lower/upper). + ! + ! + end if + + if (size(x) psb_z_csrd_clean_zeros + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(out) :: info + ! + integer(psb_ipk_) :: i, j, k, nr + integer(psb_ipk_), allocatable :: ilrp(:) + + info = 0 + call a%sync() + nr = a%get_nrows() + ilrp = a%irp(:) + a%irp(1) = 1 + j = a%irp(1) + do i=1, nr + a%diagp(i) = 0 + do k = ilrp(i), ilrp(i+1) -1 + if ((a%ja(j)==i).or.(a%val(k) /= zzero)) then + a%val(j) = a%val(k) + a%ja(j) = a%ja(k) + if (a%ja(j) == i) a%diagp(i) = j + j = j + 1 + end if + end do + a%irp(i+1) = j + end do + call a%trim() + call a%set_host() +end subroutine psb_z_csrd_clean_zeros + + +subroutine psb_z_csrd_get_diag(a,d,info) + use psb_error_mod + use psb_const_mod + use psb_z_csr_mat_mod, psb_protect_name => psb_z_csrd_get_diag + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(out) :: d(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='get_diag' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + if (a%is_unit()) then + d(1:mnm) = zone + else + do i=1, mnm + if (a%diagp(i) >0) d(i) = a%val(a%diagp(i)) + end do + end if + do i=mnm+1,size(d) + d(i) = zzero + end do + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_z_csrd_get_diag + +subroutine psb_z_cp_csrd_from_coo(a,b,info) + use psb_const_mod + use psb_realloc_mod + use psb_z_base_mat_mod + use psb_z_csr_mat_mod, psb_protect_name => psb_z_cp_csrd_from_coo + implicit none + + class(psb_z_csrd_sparse_mat), intent(inout) :: a + class(psb_z_coo_sparse_mat), intent(in) :: b + integer(psb_ipk_), intent(out) :: info + + type(psb_z_coo_sparse_mat) :: tmp + integer(psb_ipk_), allocatable :: itemp(:) + !locals + logical :: nodiag + integer(psb_ipk_) :: nza, nzt, nr, nc, i,j,k,ip,irw, err_act, ncl + integer(psb_ipk_), Parameter :: maxtry=8 + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name='z_cp_csrd_from_coo' + + info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + call a%psb_z_csr_sparse_mat%cp_from_coo(b,info) + + if (info /= 0) return + nr = a%get_nrows() + call psb_realloc(nr,a%diagp,info) + if (info /= 0) return + ! + call srch_diag(a,nodiag) + if (nodiag) then + ! + ! This means at least one of the rows does not + ! have a diagonal entry. Add explicit zeros and fix. + ! + call a%cp_to_coo(tmp,info) + nzt = tmp%get_nzeros() + call tmp%reallocate(nzt+nr) + do k=1,nr + tmp%ia(nzt+k) = k + tmp%ja(nzt+k) = k + tmp%val(nzt+k) = zzero + end do + call tmp%set_nzeros(nzt+nr) + call tmp%fix(info) + call a%mv_from_coo(tmp,info) + call srch_diag(a,nodiag) + ! At this point, nodiag should not happen + if (nodiag) then + write(0,*) 'Nodiag here?' + info = -1 + end if + end if + +contains + subroutine srch_diag(a,nodiag) + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + logical, intent(out) :: nodiag + integer(psb_ipk_) :: nr, i,j,k + + nr = a%get_nrows() + nodiag = .false. + outer: do i=1, nr + a%diagp(i) = -1 + inner: do k=a%irp(i),a%irp(i+1)-1 + if (a%ja(k) == i) then + a%diagp(i) = k + exit inner + end if + end do inner + if (a%diagp(i) < 1) then + nodiag = .true. + exit outer + end if + end do outer + return + end subroutine srch_diag +end subroutine psb_z_cp_csrd_from_coo + +subroutine psb_z_csrd_trim(a) + use psb_realloc_mod + use psb_error_mod + use psb_z_csr_mat_mod, psb_protect_name => psb_z_csrd_trim + implicit none + class(psb_z_csrd_sparse_mat), intent(inout) :: a + integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='trim' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + call a%psb_z_csr_sparse_mat%trim() + m = a%get_nrows() + call psb_realloc(m,a%diagp,info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psb_z_csrd_trim + + +subroutine psb_z_csrd_mold(a,b,info) + use psb_z_csr_mat_mod, psb_protect_name => psb_z_csrd_mold + use psb_error_mod + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + class(psb_z_base_sparse_mat), intent(inout), allocatable :: b + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='csrd_mold' + logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + info = 0 + if (allocated(b)) then + call b%free() + deallocate(b,stat=info) + end if + if (info == 0) allocate(psb_z_csrd_sparse_mat :: b, stat=info) + + if (info /= 0) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info, name) + goto 9999 + end if + return + +9999 call psb_error_handler(err_act) + return + +end subroutine psb_z_csrd_mold + +subroutine psb_z_csrd_trmv(alpha,a,x,beta,y,info,uplo,diag) + use psb_z_csr_mat_mod, psb_protect_name => psb_z_csrd_trmv + use psb_error_mod + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:) + complex(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: uplo, diag + + integer(psb_ipk_) :: i, m, n + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='z_csrd_trmv' + character :: uplo_, diag_ + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + if (a%is_dev()) call a%sync() + + if (present(uplo)) then + uplo_ = uplo + else + uplo_ = 'L' + end if + if (present(diag)) then + diag_ = diag + else + diag_ = 'D' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + n = a%get_ncols() + m = a%get_nrows() + + if (size(x,1) psb_z_csrd_cssm + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + complex(psb_dpk_), intent(inout) :: y(:,:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + + integer(psb_ipk_) :: jc, nc + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='z_csrd_cssm' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + + nc = min(size(x,2) , size(y,2)) + + do jc = 1, nc + call a%inner_cssv(alpha,x(:,jc),beta,y(:,jc),info,trans) + if (info /= 0) exit + end do + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='inner_trmv') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine psb_z_csrd_cssm + + +subroutine psb_z_csrd_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psb_string_mod + use psb_z_csr_mat_mod, psb_protect_name => psb_z_csrd_cssv + implicit none + class(psb_z_csrd_sparse_mat), intent(in) :: a + complex(psb_dpk_), intent(in) :: alpha, beta, x(:) + complex(psb_dpk_), intent(inout) :: y(:) + integer(psb_ipk_), intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc + complex(psb_dpk_) :: acc + complex(psb_dpk_), allocatable :: tmp(:) + logical :: tra,ctra + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='z_csrd_cssv' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (a%is_dev()) call a%sync() + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + !info = psb_err_invalid_mat_state_ + !call psb_errpush(info,name) + !goto 9999 + ! + ! The whole point of CSRD is to fake a triangle + ! "on the fly" (this is really for the sake of + ! saving memory in Gauss-Seidel lower/upper). + ! + ! + end if + + if (size(x)