From daa1f91c8e33e0f34ff0820800ff1461da444672 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 10 Jul 2014 16:17:24 +0000 Subject: [PATCH] psblas3-matasb: base/modules/Makefile base/modules/psb_d_base_vect_mod.f90 base/modules/psb_d_mat_mod.f90 base/modules/psb_d_tools_mod.f90 base/modules/psb_d_vect_mod.F90 base/modules/psb_i_base_vect_mod.f90 base/modules/psb_i_vect_mod.F90 base/modules/psb_vect_mod.f90 base/serial/impl/psb_d_base_mat_impl.F90 base/serial/impl/psb_d_mat_impl.F90 base/tools/psb_dspins.f90 csputs with encapsulatd vectors. Defined encapsulated multivectors. --- base/modules/Makefile | 2 +- base/modules/psb_d_base_vect_mod.f90 | 1097 ++++++++++++++++++++++ base/modules/psb_d_mat_mod.f90 | 22 +- base/modules/psb_d_tools_mod.f90 | 14 + base/modules/psb_d_vect_mod.F90 | 637 ++++++++++++- base/modules/psb_i_base_vect_mod.f90 | 1094 +++++++++++++++++++++ base/modules/psb_i_vect_mod.F90 | 632 +++++++++++++ base/modules/psb_vect_mod.f90 | 2 + base/serial/impl/psb_d_base_mat_impl.F90 | 3 + base/serial/impl/psb_d_mat_impl.F90 | 55 +- base/tools/psb_dspins.f90 | 166 ++++ 11 files changed, 3714 insertions(+), 10 deletions(-) diff --git a/base/modules/Makefile b/base/modules/Makefile index 2866a966..265e2bbc 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -67,7 +67,7 @@ psb_z_base_mat_mod.o: psb_z_base_vect_mod.o psb_c_base_vect_mod.o psb_s_base_vect_mod.o psb_d_base_vect_mod.o psb_z_base_vect_mod.o: psb_i_base_vect_mod.o psb_i_base_vect_mod.o psb_c_base_vect_mod.o psb_s_base_vect_mod.o psb_d_base_vect_mod.o psb_z_base_vect_mod.o: psi_serial_mod.o psb_realloc_mod.o psb_s_mat_mod.o: psb_s_base_mat_mod.o psb_s_csr_mat_mod.o psb_s_csc_mat_mod.o psb_s_vect_mod.o -psb_d_mat_mod.o: psb_d_base_mat_mod.o psb_d_csr_mat_mod.o psb_d_csc_mat_mod.o psb_d_vect_mod.o +psb_d_mat_mod.o: psb_d_base_mat_mod.o psb_d_csr_mat_mod.o psb_d_csc_mat_mod.o psb_d_vect_mod.o psb_i_vect_mod.o psb_c_mat_mod.o: psb_c_base_mat_mod.o psb_c_csr_mat_mod.o psb_c_csc_mat_mod.o psb_c_vect_mod.o psb_z_mat_mod.o: psb_z_base_mat_mod.o psb_z_csr_mat_mod.o psb_z_csc_mat_mod.o psb_z_vect_mod.o psb_s_csc_mat_mod.o psb_s_csr_mat_mod.o: psb_s_base_mat_mod.o diff --git a/base/modules/psb_d_base_vect_mod.f90 b/base/modules/psb_d_base_vect_mod.f90 index 2145f070..98cdac60 100644 --- a/base/modules/psb_d_base_vect_mod.f90 +++ b/base/modules/psb_d_base_vect_mod.f90 @@ -1123,3 +1123,1100 @@ contains end subroutine d_base_sctb_x end module psb_d_base_vect_mod + + + + +module psb_d_base_multivect_mod + + use psb_const_mod + use psb_error_mod + + + !> \namespace psb_base_mod \class psb_d_base_vect_type + !! The psb_d_base_vect_type + !! defines a middle level integer(psb_ipk_) encapsulated dense vector. + !! The encapsulation is needed, in place of a simple array, to allow + !! for complicated situations, such as GPU programming, where the memory + !! area we are interested in is not easily accessible from the host/Fortran + !! side. It is also meant to be encapsulated in an outer type, to allow + !! runtime switching as per the STATE design pattern, similar to the + !! sparse matrix types. + !! + private + public :: psb_d_base_multivect, psb_d_base_multivect_type + + type psb_d_base_multivect_type + !> Values. + real(psb_dpk_), allocatable :: v(:,:) + contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => d_base_mv_bld_x + procedure, pass(x) :: bld_n => d_base_mv_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => d_base_mv_all + procedure, pass(x) :: mold => d_base_mv_mold + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => d_base_mv_ins + procedure, pass(x) :: zero => d_base_mv_zero + procedure, pass(x) :: asb => d_base_mv_asb + procedure, pass(x) :: free => d_base_mv_free + ! + ! Sync: centerpiece of handling of external storage. + ! Any derived class having extra storage upon sync + ! will guarantee that both fortran/host side and + ! external side contain the same data. The base + ! version is only a placeholder. + ! + procedure, pass(x) :: sync => d_base_mv_sync + procedure, pass(x) :: is_host => d_base_mv_is_host + procedure, pass(x) :: is_dev => d_base_mv_is_dev + procedure, pass(x) :: is_sync => d_base_mv_is_sync + procedure, pass(x) :: set_host => d_base_mv_set_host + procedure, pass(x) :: set_dev => d_base_mv_set_dev + procedure, pass(x) :: set_sync => d_base_mv_set_sync + + ! + ! Basic info + procedure, pass(x) :: get_nrows => d_base_mv_get_nrows + procedure, pass(x) :: get_ncols => d_base_mv_get_ncols + procedure, pass(x) :: sizeof => d_base_mv_sizeof + procedure, nopass :: get_fmt => d_base_mv_get_fmt + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: get_vect => d_base_mv_get_vect + procedure, pass(x) :: set_scal => d_base_mv_set_scal + procedure, pass(x) :: set_vect => d_base_mv_set_vect + generic, public :: set => set_vect, set_scal + + ! + ! Dot product and AXPBY + ! +!!$ procedure, pass(x) :: dot_v => d_base_mv_dot_v +!!$ procedure, pass(x) :: dot_a => d_base_mv_dot_a +!!$ generic, public :: dot => dot_v, dot_a +!!$ procedure, pass(y) :: axpby_v => d_base_mv_axpby_v +!!$ procedure, pass(y) :: axpby_a => d_base_mv_axpby_a +!!$ generic, public :: axpby => axpby_v, axpby_a +!!$ ! +!!$ ! Vector by vector multiplication. Need all variants +!!$ ! to handle multiple requirements from preconditioners +!!$ ! +!!$ procedure, pass(y) :: mlt_v => d_base_mv_mlt_v +!!$ procedure, pass(y) :: mlt_a => d_base_mv_mlt_a +!!$ procedure, pass(z) :: mlt_a_2 => d_base_mv_mlt_a_2 +!!$ procedure, pass(z) :: mlt_v_2 => d_base_mv_mlt_v_2 +!!$ procedure, pass(z) :: mlt_va => d_base_mv_mlt_va +!!$ procedure, pass(z) :: mlt_av => d_base_mv_mlt_av +!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va +!!$ ! +!!$ ! Scaling and norms +!!$ ! +!!$ procedure, pass(x) :: scal => d_base_mv_scal +!!$ procedure, pass(x) :: nrm2 => d_base_mv_nrm2 +!!$ procedure, pass(x) :: amax => d_base_mv_amax +!!$ procedure, pass(x) :: asum => d_base_mv_asum +!!$ ! +!!$ ! Gather/scatter. These are needed for MPI interfacing. +!!$ ! May have to be reworked. +!!$ ! +!!$ procedure, pass(x) :: gthab => d_base_mv_gthab +!!$ procedure, pass(x) :: gthzv => d_base_mv_gthzv +!!$ procedure, pass(x) :: gthzv_x => d_base_mv_gthzv_x +!!$ generic, public :: gth => gthab, gthzv, gthzv_x +!!$ procedure, pass(y) :: sctb => d_base_mv_sctb +!!$ procedure, pass(y) :: sctb_x => d_base_mv_sctb_x +!!$ generic, public :: sct => sctb, sctb_x + end type psb_d_base_multivect_type + + interface psb_d_base_multivect + module procedure constructor, size_const + end interface + +contains + + ! + ! Constructors. + ! + + !> Function constructor: + !! \brief Constructor from an array + !! \param x(:) input array to be copied + !! + function constructor(x) result(this) + real(psb_dpk_) :: x(:,:) + type(psb_d_base_multivect_type) :: this + integer(psb_ipk_) :: info + + this%v = x + call this%asb(size(x,dim=1,kind=psb_ipk_),size(x,dim=2,kind=psb_ipk_),info) + end function constructor + + + !> Function constructor: + !! \brief Constructor from size + !! \param n Size of vector to be built. + !! + function size_const(m,n) result(this) + integer(psb_ipk_), intent(in) :: m,n + type(psb_d_base_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%asb(m,n,info) + + end function size_const + + ! + ! Build from a sample + ! + + !> Function bld_x: + !! \memberof psb_d_base_multivect_type + !! \brief Build method from an array + !! \param x(:) input array to be copied + !! + subroutine d_base_mv_bld_x(x,this) + use psb_realloc_mod + real(psb_dpk_), intent(in) :: this(:,:) + class(psb_d_base_multivect_type), intent(inout) :: x + integer(psb_ipk_) :: info + + call psb_realloc(size(this,1),size(this,2),x%v,info) + if (info /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld') + return + end if + x%v(:,:) = this(:,:) + + end subroutine d_base_mv_bld_x + + ! + ! Create with size, but no initialization + ! + + !> Function bld_n: + !! \memberof psb_d_base_multivect_type + !! \brief Build method with size (uninitialized data) + !! \param n size to be allocated. + !! + subroutine d_base_mv_bld_n(x,m,n) + use psb_realloc_mod + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_base_multivect_type), intent(inout) :: x + integer(psb_ipk_) :: info + + call psb_realloc(m,n,x%v,info) + call x%asb(m,n,info) + + end subroutine d_base_mv_bld_n + + !> Function base_mv_all: + !! \memberof psb_d_base_multivect_type + !! \brief Build method with size (uninitialized data) and + !! allocation return code. + !! \param n size to be allocated. + !! \param info return code + !! + subroutine d_base_mv_all(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_base_multivect_type), intent(out) :: x + integer(psb_ipk_), intent(out) :: info + + call psb_realloc(m,n,x%v,info) + + end subroutine d_base_mv_all + + !> Function base_mv_mold: + !! \memberof psb_d_base_multivect_type + !! \brief Mold method: return a variable with the same dynamic type + !! \param y returned variable + !! \param info return code + !! + subroutine d_base_mv_mold(x, y, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + class(psb_d_base_multivect_type), intent(out), allocatable :: y + integer(psb_ipk_), intent(out) :: info + + allocate(psb_d_base_multivect_type :: y, stat=info) + + end subroutine d_base_mv_mold + + ! + ! Insert a bunch of values at specified positions. + ! + !> Function base_mv_ins: + !! \memberof psb_d_base_multivect_type + !! \brief Insert coefficients. + !! + !! + !! Given a list of N pairs + !! (IRL(i),VAL(i)) + !! record a new coefficient in X such that + !! X(IRL(1:N)) = VAL(1:N). + !! + !! - the update operation will perform either + !! X(IRL(1:n)) = VAL(1:N) + !! or + !! X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N) + !! according to the value of DUPLICATE. + !! + !! + !! \param n number of pairs in input + !! \param irl(:) the input row indices + !! \param val(:) the input coefficients + !! \param dupl how to treat duplicate entries + !! \param info return code + !! + ! + subroutine d_base_mv_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: n, dupl + integer(psb_ipk_), intent(in) :: irl(:) + real(psb_dpk_), intent(in) :: val(:,:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: i, isz + + info = 0 + if (psb_errstatus_fatal()) return + + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + else if (n > min(size(irl),size(val))) then + info = psb_err_invalid_input_ + + else + isz = size(x%v,1) + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if ((1 <= irl(i)).and.(irl(i) <= isz)) then + ! this row belongs to me + ! copy i-th row of block val in x + x%v(irl(i),:) = val(i,:) + end if + enddo + + case(psb_dupl_add_) + + do i = 1, n + !loop over all val's rows + if ((1 <= irl(i)).and.(irl(i) <= isz)) then + ! this row belongs to me + ! copy i-th row of block val in x + x%v(irl(i),:) = x%v(irl(i),:) + val(i,:) + end if + enddo + + case default + info = 321 +! !$ call psb_errpush(info,name) +! !$ goto 9999 + end select + end if + if (info /= 0) then + call psb_errpush(info,'base_mv_vect_ins') + return + end if + + end subroutine d_base_mv_ins + + ! + !> Function base_mv_zero + !! \memberof psb_d_base_multivect_type + !! \brief Zero out contents + !! + ! + subroutine d_base_mv_zero(x) + use psi_serial_mod + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=izero + + end subroutine d_base_mv_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + !> Function base_mv_asb: + !! \memberof psb_d_base_multivect_type + !! \brief Assemble vector: reallocate as necessary. + !! + !! \param n final size + !! \param info return code + !! + ! + + subroutine d_base_mv_asb(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + if ((x%get_nrows() < m).or.(x%get_ncols() Function base_mv_free: + !! \memberof psb_d_base_multivect_type + !! \brief Free vector + !! + !! \param info return code + !! + ! + subroutine d_base_mv_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v)) deallocate(x%v, stat=info) + if (info /= 0) call & + & psb_errpush(psb_err_alloc_dealloc_,'vect_free') + + end subroutine d_base_mv_free + + + + ! + ! The base version of SYNC & friends does nothing, it's just + ! a placeholder. + ! + ! + !> Function base_mv_sync: + !! \memberof psb_d_base_multivect_type + !! \brief Sync: base version is a no-op. + !! + ! + subroutine d_base_mv_sync(x) + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + + end subroutine d_base_mv_sync + + ! + !> Function base_mv_set_host: + !! \memberof psb_d_base_multivect_type + !! \brief Set_host: base version is a no-op. + !! + ! + subroutine d_base_mv_set_host(x) + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + + end subroutine d_base_mv_set_host + + ! + !> Function base_mv_set_dev: + !! \memberof psb_d_base_multivect_type + !! \brief Set_dev: base version is a no-op. + !! + ! + subroutine d_base_mv_set_dev(x) + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + + end subroutine d_base_mv_set_dev + + ! + !> Function base_mv_set_sync: + !! \memberof psb_d_base_multivect_type + !! \brief Set_sync: base version is a no-op. + !! + ! + subroutine d_base_mv_set_sync(x) + implicit none + class(psb_d_base_multivect_type), intent(inout) :: x + + end subroutine d_base_mv_set_sync + + ! + !> Function base_mv_is_dev: + !! \memberof psb_d_base_multivect_type + !! \brief Is vector on external device . + !! + ! + function d_base_mv_is_dev(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + logical :: res + + res = .false. + end function d_base_mv_is_dev + + ! + !> Function base_mv_is_host + !! \memberof psb_d_base_multivect_type + !! \brief Is vector on standard memory . + !! + ! + function d_base_mv_is_host(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + logical :: res + + res = .true. + end function d_base_mv_is_host + + ! + !> Function base_mv_is_sync + !! \memberof psb_d_base_multivect_type + !! \brief Is vector on sync . + !! + ! + function d_base_mv_is_sync(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + logical :: res + + res = .true. + end function d_base_mv_is_sync + + + ! + ! Size info. + ! + ! + !> Function base_mv_get_nrows + !! \memberof psb_d_base_multivect_type + !! \brief Number of entries + !! + ! + function d_base_mv_get_nrows(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + + res = 0 + if (allocated(x%v)) res = size(x%v,1) + + end function d_base_mv_get_nrows + + function d_base_mv_get_ncols(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + + res = 0 + if (allocated(x%v)) res = size(x%v,2) + + end function d_base_mv_get_ncols + + ! + !> Function base_mv_get_sizeof + !! \memberof psb_d_base_multivect_type + !! \brief Size in bytesa + !! + ! + function d_base_mv_sizeof(x) result(res) + implicit none + class(psb_d_base_multivect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols() + + end function d_base_mv_sizeof + + ! + !> Function base_mv_get_fmt + !! \memberof psb_d_base_multivect_type + !! \brief Format + !! + ! + function d_base_mv_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'BASE' + end function d_base_mv_get_fmt + + + ! + ! + ! + !> Function base_mv_get_vect + !! \memberof psb_d_base_multivect_type + !! \brief Extract a copy of the contents + !! + ! + function d_base_mv_get_vect(x) result(res) + class(psb_d_base_multivect_type), intent(inout) :: x + real(psb_dpk_), allocatable :: res(:,:) + integer(psb_ipk_) :: info + + if (.not.allocated(x%v)) return + call x%sync() + allocate(res(x%get_nrows(),x%get_ncols()),stat=info) + if (info /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect') + return + end if + res(:,:) = x%v(:,:) + end function d_base_mv_get_vect + + ! + ! Reset all values + ! + ! + !> Function base_mv_set_scal + !! \memberof psb_d_base_multivect_type + !! \brief Set all entries + !! \param val The value to set + !! + subroutine d_base_mv_set_scal(x,val) + class(psb_d_base_multivect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: val + + integer(psb_ipk_) :: info + x%v = val + + end subroutine d_base_mv_set_scal + + ! + !> Function base_mv_set_vect + !! \memberof psb_d_base_multivect_type + !! \brief Set all entries + !! \param val(:) The vector to be copied in + !! + subroutine d_base_mv_set_vect(x,val) + class(psb_d_base_multivect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: val(:,:) + integer(psb_ipk_) :: nr + integer(psb_ipk_) :: info + + if (allocated(x%v)) then + nr = min(size(x%v,1),size(val,1)) + nc = min(size(x%v,2),size(val,2)) + + x%v(1:nr,1:nc) = val(1:nr,1:nc) + else + x%v = val + end if + + end subroutine d_base_mv_set_vect + +!!$ ! +!!$ ! Dot products +!!$ ! +!!$ ! +!!$ !> Function base_mv_dot_v +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief Dot product by another base_mv_vector +!!$ !! \param n Number of entries to be considere +!!$ !! \param y The other (base_mv_vect) to be multiplied by +!!$ !! +!!$ function d_base_mv_dot_v(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x, y +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ real(psb_dpk_), external :: ddot +!!$ +!!$ res = izero +!!$ ! +!!$ ! Note: this is the base implementation. +!!$ ! When we get here, we are sure that X is of +!!$ ! TYPE psb_d_base_mv_vect. +!!$ ! If Y is not, throw the burden on it, implicitly +!!$ ! calling dot_a +!!$ ! +!!$ select type(yy => y) +!!$ type is (psb_d_base_multivect_type) +!!$ res = ddot(n,x%v,1,y%v,1) +!!$ class default +!!$ res = y%dot(n,x%v) +!!$ end select +!!$ +!!$ end function d_base_mv_dot_v +!!$ +!!$ ! +!!$ ! Base workhorse is good old BLAS1 +!!$ ! +!!$ ! +!!$ !> Function base_mv_dot_a +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief Dot product by a normal array +!!$ !! \param n Number of entries to be considere +!!$ !! \param y(:) The array to be multiplied by +!!$ !! +!!$ function d_base_mv_dot_a(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ integer(psb_ipk_), external :: ddot +!!$ +!!$ res = ddot(n,y,1,x%v,1) +!!$ +!!$ end function d_base_mv_dot_a + +!!$ ! +!!$ ! AXPBY is invoked via Y, hence the structure below. +!!$ ! +!!$ ! +!!$ ! +!!$ !> Function base_mv_axpby_v +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief AXPBY by a (base_mv_vect) y=alpha*x+beta*y +!!$ !! \param m Number of entries to be considere +!!$ !! \param alpha scalar alpha +!!$ !! \param x The class(base_mv_vect) to be added +!!$ !! \param beta scalar alpha +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_axpby_v(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ real(psb_dpk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ select type(xx => x) +!!$ type is (psb_d_base_multivect_type) +!!$ call psb_geaxpby(m,alpha,x%v,beta,y%v,info) +!!$ class default +!!$ call y%axpby(m,alpha,x%v,beta,info) +!!$ end select +!!$ +!!$ end subroutine d_base_mv_axpby_v +!!$ +!!$ ! +!!$ ! AXPBY is invoked via Y, hence the structure below. +!!$ ! +!!$ ! +!!$ !> Function base_mv_axpby_a +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief AXPBY by a normal array y=alpha*x+beta*y +!!$ !! \param m Number of entries to be considere +!!$ !! \param alpha scalar alpha +!!$ !! \param x(:) The array to be added +!!$ !! \param beta scalar alpha +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_axpby_a(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ real(psb_dpk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ call psb_geaxpby(m,alpha,x,beta,y%v,info) +!!$ +!!$ end subroutine d_base_mv_axpby_a + + +!!$ ! +!!$ ! Multiple variants of two operations: +!!$ ! Simple multiplication Y(:) = X(:)*Y(:) +!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:) +!!$ ! +!!$ ! Variants expanded according to the dynamic type +!!$ ! of the involved entities +!!$ ! +!!$ ! +!!$ !> Function base_mv_mlt_a +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief Vector entry-by-entry multiply by a base_mv_vect array y=x*y +!!$ !! \param x The class(base_mv_vect) to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_mlt_v(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ select type(xx => x) +!!$ type is (psb_d_base_multivect_type) +!!$ n = min(size(y%v), size(xx%v)) +!!$ do i=1, n +!!$ y%v(i) = y%v(i)*xx%v(i) +!!$ end do +!!$ class default +!!$ call y%mlt(x%v,info) +!!$ end select +!!$ +!!$ end subroutine d_base_mv_mlt_v +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_a +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y +!!$ !! \param x(:) The array to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_mlt_a(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ n = min(size(y%v), size(x)) +!!$ do i=1, n +!!$ y%v(i) = y%v(i)*x(i) +!!$ end do +!!$ +!!$ end subroutine d_base_mv_mlt_a +!!$ +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_a_2 +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays +!!$ !! z=beta*z+alpha*x*y +!!$ !! \param alpha +!!$ !! \param beta +!!$ !! \param x(:) The array to be multiplied b +!!$ !! \param y(:) The array to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_mlt_a_2(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ n = min(size(z%v), size(x), size(y)) +!!$ if (alpha == izero) then +!!$ if (beta == ione) then +!!$ return +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) +!!$ end do +!!$ end if +!!$ else +!!$ if (alpha == ione) then +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) + y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) + y(i)*x(i) +!!$ end do +!!$ end if +!!$ else if (alpha == -ione) then +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = -y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) - y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) - y(i)*x(i) +!!$ end do +!!$ end if +!!$ else +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = alpha*y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i) +!!$ end do +!!$ end if +!!$ end if +!!$ end if +!!$ end subroutine d_base_mv_mlt_a_2 +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_v_2 +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mv_vect) +!!$ !! z=beta*z+alpha*x*y +!!$ !! \param alpha +!!$ !! \param beta +!!$ !! \param x The class(base_mv_vect) to be multiplied b +!!$ !! \param y The class(base_mv_vect) to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine d_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy) +!!$ use psi_serial_mod +!!$ use psb_string_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ class(psb_d_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character(len=1), intent(in), optional :: conjgx, conjgy +!!$ integer(psb_ipk_) :: i, n +!!$ logical :: conjgx_, conjgy_ +!!$ +!!$ info = 0 +!!$ if (.not.psb_i_is_complex_) then +!!$ call z%mlt(alpha,x%v,y%v,beta,info) +!!$ else +!!$ conjgx_=.false. +!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C') +!!$ conjgy_=.false. +!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C') +!!$ if (conjgx_) x%v=(x%v) +!!$ if (conjgy_) y%v=(y%v) +!!$ call z%mlt(alpha,x%v,y%v,beta,info) +!!$ if (conjgx_) x%v=(x%v) +!!$ if (conjgy_) y%v=(y%v) +!!$ end if +!!$ end subroutine d_base_mv_mlt_v_2 +!!$ +!!$ subroutine d_base_mv_mlt_av(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_base_multivect_type), intent(inout) :: y +!!$ class(psb_d_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ call z%mlt(alpha,x,y%v,beta,info) +!!$ +!!$ end subroutine d_base_mv_mlt_av +!!$ +!!$ subroutine d_base_mv_mlt_va(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ class(psb_d_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ call z%mlt(alpha,y,x,beta,info) +!!$ +!!$ end subroutine d_base_mv_mlt_va +!!$ +!!$ +!!$ ! +!!$ ! Simple scaling +!!$ ! +!!$ !> Function base_mv_scal +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief Scale all entries x = alpha*x +!!$ !! \param alpha The multiplier +!!$ !! +!!$ subroutine d_base_mv_scal(alpha, x) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ real(psb_dpk_), intent (in) :: alpha +!!$ +!!$ if (allocated(x%v)) x%v = alpha*x%v +!!$ +!!$ end subroutine d_base_mv_scal +!!$ +!!$ ! +!!$ ! Norms 1, 2 and infinity +!!$ ! +!!$ !> Function base_mv_nrm2 +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief 2-norm |x(1:n)|_2 +!!$ !! \param n how many entries to consider +!!$ function d_base_mv_nrm2(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ integer(psb_ipk_), external :: dnrm2 +!!$ +!!$ res = dnrm2(n,x%v,1) +!!$ +!!$ end function d_base_mv_nrm2 +!!$ +!!$ ! +!!$ !> Function base_mv_amax +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief infinity-norm |x(1:n)|_\infty +!!$ !! \param n how many entries to consider +!!$ function d_base_mv_amax(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ res = maxval(abs(x%v(1:n))) +!!$ +!!$ end function d_base_mv_amax +!!$ +!!$ ! +!!$ !> Function base_mv_asum +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief 1-norm |x(1:n)|_1 +!!$ !! \param n how many entries to consider +!!$ function d_base_mv_asum(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ res = sum(abs(x%v(1:n))) +!!$ +!!$ end function d_base_mv_asum +!!$ +!!$ +!!$ ! +!!$ ! Gather: Y = beta * Y + alpha * X(IDX(:)) +!!$ ! +!!$ ! +!!$ !> Function base_mv_gthab +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief gather into an array +!!$ !! Y = beta * Y + alpha * X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ !! \param alpha +!!$ !! \param beta +!!$ subroutine d_base_mv_gthab(n,idx,alpha,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: alpha, beta, y(:) +!!$ class(psb_d_base_multivect_type) :: x +!!$ +!!$ call x%sync() +!!$ call psi_gth(n,idx,alpha,x%v,beta,y) +!!$ +!!$ end subroutine d_base_mv_gthab +!!$ ! +!!$ ! shortcut alpha=1 beta=0 +!!$ ! +!!$ !> Function base_mv_gthzv +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief gather into an array special alpha=1 beta=0 +!!$ !! Y = X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ subroutine d_base_mv_gthzv_x(i,n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: i,n +!!$ class(psb_d_base_multivect_type) :: idx +!!$ real(psb_dpk_) :: y(:) +!!$ class(psb_d_base_multivect_type) :: x +!!$ +!!$ call x%gth(n,idx%v(i:),y) +!!$ +!!$ end subroutine d_base_mv_gthzv_x +!!$ +!!$ ! +!!$ ! shortcut alpha=1 beta=0 +!!$ ! +!!$ !> Function base_mv_gthzv +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief gather into an array special alpha=1 beta=0 +!!$ !! Y = X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ subroutine d_base_mv_gthzv(n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: y(:) +!!$ class(psb_d_base_multivect_type) :: x +!!$ +!!$ call x%sync() +!!$ call psi_gth(n,idx,x%v,y) +!!$ +!!$ end subroutine d_base_mv_gthzv +!!$ +!!$ ! +!!$ ! Scatter: +!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) +!!$ ! +!!$ ! +!!$ !> Function base_mv_sctb +!!$ !! \memberof psb_d_base_multivect_type +!!$ !! \brief scatter into a class(base_mv_vect) +!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ !! \param beta +!!$ !! \param x(:) +!!$ subroutine d_base_mv_sctb(n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: beta, x(:) +!!$ class(psb_d_base_multivect_type) :: y +!!$ +!!$ call y%sync() +!!$ call psi_sct(n,idx,x,beta,y%v) +!!$ call y%set_host() +!!$ +!!$ end subroutine d_base_mv_sctb +!!$ +!!$ subroutine d_base_mv_sctb_x(i,n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: i, n +!!$ class(psb_d_base_multivect_type) :: idx +!!$ real(psb_dpk_) :: beta, x(:) +!!$ class(psb_d_base_multivect_type) :: y +!!$ +!!$ call y%sct(n,idx%v(i:),x,beta) +!!$ +!!$ end subroutine d_base_mv_sctb_x + +end module psb_d_base_multivect_mod + diff --git a/base/modules/psb_d_mat_mod.f90 b/base/modules/psb_d_mat_mod.f90 index 5369f2c9..24d8a536 100644 --- a/base/modules/psb_d_mat_mod.f90 +++ b/base/modules/psb_d_mat_mod.f90 @@ -119,7 +119,9 @@ module psb_d_mat_mod procedure, pass(a) :: csall => psb_d_csall procedure, pass(a) :: free => psb_d_free procedure, pass(a) :: trim => psb_d_trim - procedure, pass(a) :: csput => psb_d_csput + procedure, pass(a) :: csput_a => psb_d_csput_a + procedure, pass(a) :: csput_v => psb_d_csput_v + generic, public :: csput => csput_a, csput_v procedure, pass(a) :: csgetptn => psb_d_csgetptn procedure, pass(a) :: csgetrow => psb_d_csgetrow procedure, pass(a) :: csgetblk => psb_d_csgetblk @@ -379,14 +381,28 @@ module psb_d_mat_mod end interface interface - subroutine psb_d_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + subroutine psb_d_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) import :: psb_ipk_, psb_dspmat_type, psb_dpk_ class(psb_dspmat_type), intent(inout) :: a real(psb_dpk_), intent(in) :: val(:) integer(psb_ipk_), intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: gtl(:) - end subroutine psb_d_csput + end subroutine psb_d_csput_a + end interface + + interface + subroutine psb_d_csput_v(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_d_vect_mod, only : psb_d_vect_type + use psb_i_vect_mod, only : psb_i_vect_type + import :: psb_ipk_, psb_dspmat_type, psb_dpk_ + class(psb_dspmat_type), intent(inout) :: a + type(psb_d_vect_type), intent(inout) :: val + type(psb_i_vect_type), intent(inout) :: ia, ja + integer(psb_ipk_), intent(in) :: nz, imin,imax,jmin,jmax + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: gtl(:) + end subroutine psb_d_csput_v end interface interface diff --git a/base/modules/psb_d_tools_mod.f90 b/base/modules/psb_d_tools_mod.f90 index 013f9dd0..717cc56a 100644 --- a/base/modules/psb_d_tools_mod.f90 +++ b/base/modules/psb_d_tools_mod.f90 @@ -286,6 +286,20 @@ Module psb_d_tools_mod logical, intent(in), optional :: rebuild logical, intent(in), optional :: local end subroutine psb_dspins + subroutine psb_dspins_v(nz,ia,ja,val,a,desc_a,info,rebuild,local) + use psb_i_vect_mod, only : psb_i_vect_type + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_base_vect_type, psb_d_vect_type,& + & psb_dspmat_type, psb_d_base_sparse_mat + type(psb_desc_type), intent(inout) :: desc_a + type(psb_dspmat_type), intent(inout) :: a + integer(psb_ipk_), intent(in) :: nz + type(psb_i_vect_type), intent(inout) :: ia,ja + type(psb_d_vect_type), intent(inout) :: val + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: rebuild + logical, intent(in), optional :: local + end subroutine psb_dspins_v subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) import :: psb_desc_type, psb_dpk_, psb_ipk_, & & psb_d_base_vect_type, psb_d_vect_type, & diff --git a/base/modules/psb_d_vect_mod.F90 b/base/modules/psb_d_vect_mod.F90 index c7563154..7c263fd5 100644 --- a/base/modules/psb_d_vect_mod.F90 +++ b/base/modules/psb_d_vect_mod.F90 @@ -180,9 +180,9 @@ contains call mld%mold(x%v,info) #endif endif - +!!$ write(0,*) 'Size invect',size(invect) if (info == psb_success_) call x%v%bld(invect) - + end subroutine d_vect_bld_x @@ -663,3 +663,636 @@ contains end subroutine d_vect_cnv end module psb_d_vect_mod + + + + +module psb_d_multivect_mod + + use psb_d_base_multivect_mod + use psb_const_mod + + private + + type psb_d_multivect_type + class(psb_d_base_multivect_type), allocatable :: v + contains + procedure, pass(x) :: get_nrows => d_vect_get_nrows + procedure, pass(x) :: get_ncols => d_vect_get_ncols + procedure, pass(x) :: sizeof => d_vect_sizeof + procedure, pass(x) :: get_fmt => d_vect_get_fmt +!!$ procedure, pass(x) :: dot_v => d_vect_dot_v +!!$ procedure, pass(x) :: dot_a => d_vect_dot_a +!!$ generic, public :: dot => dot_v, dot_a +!!$ procedure, pass(y) :: axpby_v => d_vect_axpby_v +!!$ procedure, pass(y) :: axpby_a => d_vect_axpby_a +!!$ generic, public :: axpby => axpby_v, axpby_a +!!$ procedure, pass(y) :: mlt_v => d_vect_mlt_v +!!$ procedure, pass(y) :: mlt_a => d_vect_mlt_a +!!$ procedure, pass(z) :: mlt_a_2 => d_vect_mlt_a_2 +!!$ procedure, pass(z) :: mlt_v_2 => d_vect_mlt_v_2 +!!$ procedure, pass(z) :: mlt_va => d_vect_mlt_va +!!$ procedure, pass(z) :: mlt_av => d_vect_mlt_av +!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2,& +!!$ & mlt_v_2, mlt_av, mlt_va +!!$ procedure, pass(x) :: scal => d_vect_scal +!!$ procedure, pass(x) :: nrm2 => d_vect_nrm2 +!!$ procedure, pass(x) :: amax => d_vect_amax +!!$ procedure, pass(x) :: asum => d_vect_asum + procedure, pass(x) :: all => d_vect_all + procedure, pass(x) :: reall => d_vect_reall + procedure, pass(x) :: zero => d_vect_zero + procedure, pass(x) :: asb => d_vect_asb + procedure, pass(x) :: sync => d_vect_sync +!!$ procedure, pass(x) :: gthab => d_vect_gthab +!!$ procedure, pass(x) :: gthzv => d_vect_gthzv +!!$ generic, public :: gth => gthab, gthzv +!!$ procedure, pass(y) :: sctb => d_vect_sctb +!!$ generic, public :: sct => sctb + procedure, pass(x) :: free => d_vect_free + procedure, pass(x) :: ins => d_vect_ins + procedure, pass(x) :: bld_x => d_vect_bld_x + procedure, pass(x) :: bld_n => d_vect_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: get_vect => d_vect_get_vect + procedure, pass(x) :: cnv => d_vect_cnv + procedure, pass(x) :: set_scal => d_vect_set_scal + procedure, pass(x) :: set_vect => d_vect_set_vect + generic, public :: set => set_vect, set_scal + procedure, pass(x) :: clone => d_vect_clone + end type psb_d_multivect_type + + public :: psb_d_multivect, psb_d_multivect_type,& + & psb_set_multivect_default, psb_get_multivect_default + + interface psb_d_multivect + module procedure constructor, size_const + end interface + + class(psb_d_base_multivect_type), allocatable, target,& + & save, private :: psb_d_base_multivect_default + + interface psb_set_multivect_default + module procedure psb_d_set_multivect_default + end interface + + interface psb_get_vect_default + module procedure psb_d_get_multivect_default + end interface + + +contains + + + subroutine psb_d_set_multivect_default(v) + implicit none + class(psb_d_base_multivect_type), intent(in) :: v + + if (allocated(psb_d_base_multivect_default)) then + deallocate(psb_d_base_multivect_default) + end if + allocate(psb_d_base_multivect_default, mold=v) + + end subroutine psb_d_set_multivect_default + + function psb_d_get_multivect_default(v) result(res) + implicit none + class(psb_d_multivect_type), intent(in) :: v + class(psb_d_base_multivect_type), pointer :: res + + res => psb_d_get_base_multivect_default() + + end function psb_d_get_multivect_default + + + function psb_d_get_base_multivect_default() result(res) + implicit none + class(psb_d_base_multivect_type), pointer :: res + + if (.not.allocated(psb_d_base_multivect_default)) then + allocate(psb_d_base_multivect_type :: psb_d_base_multivect_default) + end if + + res => psb_d_base_multivect_default + + end function psb_d_get_base_multivect_default + + + subroutine d_vect_clone(x,y,info) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + class(psb_d_multivect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + call y%free(info) + if ((info==0).and.allocated(x%v)) then + call y%bld(x%get_vect(),mold=x%v) + end if + end subroutine d_vect_clone + + subroutine d_vect_bld_x(x,invect,mold) + real(psb_dpk_), intent(in) :: invect(:,:) + class(psb_d_multivect_type), intent(out) :: x + class(psb_d_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_) :: info + class(psb_d_base_multivect_type), pointer :: mld + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else +#ifdef HAVE_MOLD + allocate(x%v,stat=info, mold=psb_d_get_base_multivect_default()) +#else + mld = psb_d_get_base_multivect_default() + call mld%mold(x%v,info) +#endif + endif + + if (info == psb_success_) call x%v%bld(invect) + + end subroutine d_vect_bld_x + + + subroutine d_vect_bld_n(x,m,n,mold) + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_multivect_type), intent(out) :: x + class(psb_d_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_) :: info + class(psb_d_base_multivect_type), pointer :: mld + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else +#ifdef HAVE_MOLD + allocate(x%v,stat=info, mold=psb_d_get_base_multivect_default()) +#else + mld = psb_d_get_base_multivect_default() + call mld%mold(x%v,info) +#endif + endif + if (info == psb_success_) call x%v%bld(m,n) + + end subroutine d_vect_bld_n + + function d_vect_get_vect(x) result(res) + class(psb_d_multivect_type), intent(inout) :: x + real(psb_dpk_), allocatable :: res(:,:) + integer(psb_ipk_) :: info + + if (allocated(x%v)) then + res = x%v%get_vect() + end if + end function d_vect_get_vect + + subroutine d_vect_set_scal(x,val) + class(psb_d_multivect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: val + + integer(psb_ipk_) :: info + if (allocated(x%v)) call x%v%set(val) + + end subroutine d_vect_set_scal + + subroutine d_vect_set_vect(x,val) + class(psb_d_multivect_type), intent(inout) :: x + real(psb_dpk_), intent(in) :: val(:,:) + + integer(psb_ipk_) :: info + if (allocated(x%v)) call x%v%set(val) + + end subroutine d_vect_set_vect + + + function constructor(x) result(this) + real(psb_dpk_) :: x(:,:) + type(psb_d_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%bld(x) + call this%asb(size(x,dim=1,kind=psb_ipk_),size(x,dim=2,kind=psb_ipk_),info) + + end function constructor + + + function size_const(m,n) result(this) + integer(psb_ipk_), intent(in) :: m,n + type(psb_d_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%bld(m,n) + call this%asb(m,n,info) + + end function size_const + + function d_vect_get_nrows(x) result(res) + implicit none + class(psb_d_multivect_type), intent(in) :: x + real(psb_dpk_) :: res + res = 0 + if (allocated(x%v)) res = x%v%get_nrows() + end function d_vect_get_nrows + + function d_vect_get_ncols(x) result(res) + implicit none + class(psb_d_multivect_type), intent(in) :: x + real(psb_dpk_) :: res + res = 0 + if (allocated(x%v)) res = x%v%get_ncols() + end function d_vect_get_ncols + + function d_vect_sizeof(x) result(res) + implicit none + class(psb_d_multivect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + res = 0 + if (allocated(x%v)) res = x%v%sizeof() + end function d_vect_sizeof + + function d_vect_get_fmt(x) result(res) + implicit none + class(psb_d_multivect_type), intent(in) :: x + character(len=5) :: res + res = 'NULL' + if (allocated(x%v)) res = x%v%get_fmt() + end function d_vect_get_fmt + +!!$ function d_vect_dot_v(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x, y +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ res = izero +!!$ if (allocated(x%v).and.allocated(y%v)) & +!!$ & res = x%v%dot(n,y%v) +!!$ +!!$ end function d_vect_dot_v +!!$ +!!$ function d_vect_dot_a(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ res = izero +!!$ if (allocated(x%v)) & +!!$ & res = x%v%dot(n,y) +!!$ +!!$ end function d_vect_dot_a +!!$ +!!$ subroutine d_vect_axpby_v(m,alpha, x, beta, y, info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ real(psb_dpk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ if (allocated(x%v).and.allocated(y%v)) then +!!$ call y%v%axpby(m,alpha,x%v,beta,info) +!!$ else +!!$ info = psb_err_invalid_vect_state_ +!!$ end if +!!$ +!!$ end subroutine d_vect_axpby_v +!!$ +!!$ subroutine d_vect_axpby_a(m,alpha, x, beta, y, info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ real(psb_dpk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ if (allocated(y%v)) & +!!$ & call y%v%axpby(m,alpha,x,beta,info) +!!$ +!!$ end subroutine d_vect_axpby_a +!!$ +!!$ +!!$ subroutine d_vect_mlt_v(x, y, info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(x%v).and.allocated(y%v)) & +!!$ & call y%v%mlt(x%v,info) +!!$ +!!$ end subroutine d_vect_mlt_v +!!$ +!!$ subroutine d_vect_mlt_a(x, y, info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ +!!$ info = 0 +!!$ if (allocated(y%v)) & +!!$ & call y%v%mlt(x,info) +!!$ +!!$ end subroutine d_vect_mlt_a +!!$ +!!$ +!!$ subroutine d_vect_mlt_a_2(alpha,x,y,beta,z,info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(z%v)) & +!!$ & call z%v%mlt(alpha,x,y,beta,info) +!!$ +!!$ end subroutine d_vect_mlt_a_2 +!!$ +!!$ subroutine d_vect_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy) +!!$ use psd_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ class(psb_d_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character(len=1), intent(in), optional :: conjgx, conjgy +!!$ +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(x%v).and.allocated(y%v).and.& +!!$ & allocated(z%v)) & +!!$ & call z%v%mlt(alpha,x%v,y%v,beta,info,conjgx,conjgy) +!!$ +!!$ end subroutine d_vect_mlt_v_2 +!!$ +!!$ subroutine d_vect_mlt_av(alpha,x,y,beta,z,info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: x(:) +!!$ class(psb_d_multivect_type), intent(inout) :: y +!!$ class(psb_d_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(z%v).and.allocated(y%v)) & +!!$ & call z%v%mlt(alpha,x,y%v,beta,info) +!!$ +!!$ end subroutine d_vect_mlt_av +!!$ +!!$ subroutine d_vect_mlt_va(alpha,x,y,beta,z,info) +!!$ use psd_serial_mod +!!$ implicit none +!!$ real(psb_dpk_), intent(in) :: alpha,beta +!!$ real(psb_dpk_), intent(in) :: y(:) +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ class(psb_d_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ if (allocated(z%v).and.allocated(x%v)) & +!!$ & call z%v%mlt(alpha,x%v,y,beta,info) +!!$ +!!$ end subroutine d_vect_mlt_va +!!$ +!!$ subroutine d_vect_scal(alpha, x) +!!$ use psd_serial_mod +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ real(psb_dpk_), intent (in) :: alpha +!!$ +!!$ if (allocated(x%v)) call x%v%scal(alpha) +!!$ +!!$ end subroutine d_vect_scal +!!$ +!!$ +!!$ function d_vect_nrm2(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%nrm2(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function d_vect_nrm2 +!!$ +!!$ function d_vect_amax(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%amax(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function d_vect_amax +!!$ +!!$ function d_vect_asum(n,x) result(res) +!!$ implicit none +!!$ class(psb_d_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ real(psb_dpk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%asum(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function d_vect_asum + + subroutine d_vect_all(m,n, x, info, mold) + + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_multivect_type), intent(out) :: x + class(psb_d_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_), intent(out) :: info + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else + allocate(psb_d_base_multivect_type :: x%v,stat=info) + endif + if (info == 0) then + call x%v%all(m,n,info) + else + info = psb_err_alloc_dealloc_ + end if + + end subroutine d_vect_all + + subroutine d_vect_reall(m,n, x, info) + + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (.not.allocated(x%v)) & + & call x%all(m,n,info) + if (info == 0) & + & call x%asb(m,n,info) + + end subroutine d_vect_reall + + subroutine d_vect_zero(x) + use psi_serial_mod + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) call x%v%zero() + + end subroutine d_vect_zero + + subroutine d_vect_asb(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_d_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v)) & + & call x%v%asb(m,n,info) + + end subroutine d_vect_asb + + subroutine d_vect_sync(x) + implicit none + class(psb_d_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%sync() + + end subroutine d_vect_sync + +!!$ subroutine d_vect_gthab(n,idx,alpha,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: alpha, beta, y(:) +!!$ class(psb_d_multivect_type) :: x +!!$ +!!$ if (allocated(x%v)) & +!!$ & call x%v%gth(n,idx,alpha,beta,y) +!!$ +!!$ end subroutine d_vect_gthab +!!$ +!!$ subroutine d_vect_gthzv(n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: y(:) +!!$ class(psb_d_multivect_type) :: x +!!$ +!!$ if (allocated(x%v)) & +!!$ & call x%v%gth(n,idx,y) +!!$ +!!$ end subroutine d_vect_gthzv +!!$ +!!$ subroutine d_vect_sctb(n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ real(psb_dpk_) :: beta, x(:) +!!$ class(psb_d_multivect_type) :: y +!!$ +!!$ if (allocated(y%v)) & +!!$ & call y%v%sct(n,idx,x,beta) +!!$ +!!$ end subroutine d_vect_sctb + + subroutine d_vect_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_d_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v)) then + call x%v%free(info) + if (info == 0) deallocate(x%v,stat=info) + end if + + end subroutine d_vect_free + + subroutine d_vect_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_d_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: n, dupl + integer(psb_ipk_), intent(in) :: irl(:) + real(psb_dpk_), intent(in) :: val(:,:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: i + + info = 0 + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + return + end if + + call x%v%ins(n,irl,val,dupl,info) + + end subroutine d_vect_ins + + + subroutine d_vect_cnv(x,mold) + class(psb_d_multivect_type), intent(inout) :: x + class(psb_d_base_multivect_type), intent(in), optional :: mold + class(psb_d_base_multivect_type), allocatable :: tmp + integer(psb_ipk_) :: info + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(tmp,stat=info,mold=mold) +#else + call mold%mold(tmp,info) +#endif + if (allocated(x%v)) then + call x%v%sync() + if (info == psb_success_) call tmp%bld(x%v%v) + call x%v%free(info) + end if + call move_alloc(tmp,x%v) + end if + end subroutine d_vect_cnv + +end module psb_d_multivect_mod diff --git a/base/modules/psb_i_base_vect_mod.f90 b/base/modules/psb_i_base_vect_mod.f90 index c188e0b7..c399c84d 100644 --- a/base/modules/psb_i_base_vect_mod.f90 +++ b/base/modules/psb_i_base_vect_mod.f90 @@ -1121,3 +1121,1097 @@ contains end subroutine i_base_sctb_x end module psb_i_base_vect_mod + + +module psb_i_base_multivect_mod + + use psb_const_mod + use psb_error_mod + + + !> \namespace psb_base_mod \class psb_i_base_vect_type + !! The psb_i_base_vect_type + !! defines a middle level integer(psb_ipk_) encapsulated dense vector. + !! The encapsulation is needed, in place of a simple array, to allow + !! for complicated situations, such as GPU programming, where the memory + !! area we are interested in is not easily accessible from the host/Fortran + !! side. It is also meant to be encapsulated in an outer type, to allow + !! runtime switching as per the STATE design pattern, similar to the + !! sparse matrix types. + !! + private + public :: psb_i_base_multivect, psb_i_base_multivect_type + + type psb_i_base_multivect_type + !> Values. + integer(psb_ipk_), allocatable :: v(:,:) + contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => i_base_mv_bld_x + procedure, pass(x) :: bld_n => i_base_mv_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => i_base_mv_all + procedure, pass(x) :: mold => i_base_mv_mold + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => i_base_mv_ins + procedure, pass(x) :: zero => i_base_mv_zero + procedure, pass(x) :: asb => i_base_mv_asb + procedure, pass(x) :: free => i_base_mv_free + ! + ! Sync: centerpiece of handling of external storage. + ! Any derived class having extra storage upon sync + ! will guarantee that both fortran/host side and + ! external side contain the same data. The base + ! version is only a placeholder. + ! + procedure, pass(x) :: sync => i_base_mv_sync + procedure, pass(x) :: is_host => i_base_mv_is_host + procedure, pass(x) :: is_dev => i_base_mv_is_dev + procedure, pass(x) :: is_sync => i_base_mv_is_sync + procedure, pass(x) :: set_host => i_base_mv_set_host + procedure, pass(x) :: set_dev => i_base_mv_set_dev + procedure, pass(x) :: set_sync => i_base_mv_set_sync + + ! + ! Basic info + procedure, pass(x) :: get_nrows => i_base_mv_get_nrows + procedure, pass(x) :: get_ncols => i_base_mv_get_ncols + procedure, pass(x) :: sizeof => i_base_mv_sizeof + procedure, nopass :: get_fmt => i_base_mv_get_fmt + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: get_vect => i_base_mv_get_vect + procedure, pass(x) :: set_scal => i_base_mv_set_scal + procedure, pass(x) :: set_vect => i_base_mv_set_vect + generic, public :: set => set_vect, set_scal + + ! + ! Dot product and AXPBY + ! +!!$ procedure, pass(x) :: dot_v => i_base_mv_dot_v +!!$ procedure, pass(x) :: dot_a => i_base_mv_dot_a +!!$ generic, public :: dot => dot_v, dot_a +!!$ procedure, pass(y) :: axpby_v => i_base_mv_axpby_v +!!$ procedure, pass(y) :: axpby_a => i_base_mv_axpby_a +!!$ generic, public :: axpby => axpby_v, axpby_a +!!$ ! +!!$ ! Vector by vector multiplication. Need all variants +!!$ ! to handle multiple requirements from preconditioners +!!$ ! +!!$ procedure, pass(y) :: mlt_v => i_base_mv_mlt_v +!!$ procedure, pass(y) :: mlt_a => i_base_mv_mlt_a +!!$ procedure, pass(z) :: mlt_a_2 => i_base_mv_mlt_a_2 +!!$ procedure, pass(z) :: mlt_v_2 => i_base_mv_mlt_v_2 +!!$ procedure, pass(z) :: mlt_va => i_base_mv_mlt_va +!!$ procedure, pass(z) :: mlt_av => i_base_mv_mlt_av +!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va +!!$ ! +!!$ ! Scaling and norms +!!$ ! +!!$ procedure, pass(x) :: scal => i_base_mv_scal +!!$ procedure, pass(x) :: nrm2 => i_base_mv_nrm2 +!!$ procedure, pass(x) :: amax => i_base_mv_amax +!!$ procedure, pass(x) :: asum => i_base_mv_asum +!!$ ! +!!$ ! Gather/scatter. These are needed for MPI interfacing. +!!$ ! May have to be reworked. +!!$ ! +!!$ procedure, pass(x) :: gthab => i_base_mv_gthab +!!$ procedure, pass(x) :: gthzv => i_base_mv_gthzv +!!$ procedure, pass(x) :: gthzv_x => i_base_mv_gthzv_x +!!$ generic, public :: gth => gthab, gthzv, gthzv_x +!!$ procedure, pass(y) :: sctb => i_base_mv_sctb +!!$ procedure, pass(y) :: sctb_x => i_base_mv_sctb_x +!!$ generic, public :: sct => sctb, sctb_x + end type psb_i_base_multivect_type + + interface psb_i_base_multivect + module procedure constructor, size_const + end interface + +contains + + ! + ! Constructors. + ! + + !> Function constructor: + !! \brief Constructor from an array + !! \param x(:) input array to be copied + !! + function constructor(x) result(this) + integer(psb_ipk_) :: x(:,:) + type(psb_i_base_multivect_type) :: this + integer(psb_ipk_) :: info + + this%v = x + call this%asb(size(x,dim=1,kind=psb_ipk_),size(x,dim=2,kind=psb_ipk_),info) + end function constructor + + + !> Function constructor: + !! \brief Constructor from size + !! \param n Size of vector to be built. + !! + function size_const(m,n) result(this) + integer(psb_ipk_), intent(in) :: m,n + type(psb_i_base_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%asb(m,n,info) + + end function size_const + + ! + ! Build from a sample + ! + + !> Function bld_x: + !! \memberof psb_i_base_multivect_type + !! \brief Build method from an array + !! \param x(:) input array to be copied + !! + subroutine i_base_mv_bld_x(x,this) + use psb_realloc_mod + integer(psb_ipk_), intent(in) :: this(:,:) + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_) :: info + + call psb_realloc(size(this,1),size(this,2),x%v,info) + if (info /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,'base_mv_vect_bld') + return + end if + x%v(:,:) = this(:,:) + + end subroutine i_base_mv_bld_x + + ! + ! Create with size, but no initialization + ! + + !> Function bld_n: + !! \memberof psb_i_base_multivect_type + !! \brief Build method with size (uninitialized data) + !! \param n size to be allocated. + !! + subroutine i_base_mv_bld_n(x,m,n) + use psb_realloc_mod + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_) :: info + + call psb_realloc(m,n,x%v,info) + call x%asb(m,n,info) + + end subroutine i_base_mv_bld_n + + !> Function base_mv_all: + !! \memberof psb_i_base_multivect_type + !! \brief Build method with size (uninitialized data) and + !! allocation return code. + !! \param n size to be allocated. + !! \param info return code + !! + subroutine i_base_mv_all(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_base_multivect_type), intent(out) :: x + integer(psb_ipk_), intent(out) :: info + + call psb_realloc(m,n,x%v,info) + + end subroutine i_base_mv_all + + !> Function base_mv_mold: + !! \memberof psb_i_base_multivect_type + !! \brief Mold method: return a variable with the same dynamic type + !! \param y returned variable + !! \param info return code + !! + subroutine i_base_mv_mold(x, y, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + class(psb_i_base_multivect_type), intent(out), allocatable :: y + integer(psb_ipk_), intent(out) :: info + + allocate(psb_i_base_multivect_type :: y, stat=info) + + end subroutine i_base_mv_mold + + ! + ! Insert a bunch of values at specified positions. + ! + !> Function base_mv_ins: + !! \memberof psb_i_base_multivect_type + !! \brief Insert coefficients. + !! + !! + !! Given a list of N pairs + !! (IRL(i),VAL(i)) + !! record a new coefficient in X such that + !! X(IRL(1:N)) = VAL(1:N). + !! + !! - the update operation will perform either + !! X(IRL(1:n)) = VAL(1:N) + !! or + !! X(IRL(1:n)) = X(IRL(1:n))+VAL(1:N) + !! according to the value of DUPLICATE. + !! + !! + !! \param n number of pairs in input + !! \param irl(:) the input row indices + !! \param val(:) the input coefficients + !! \param dupl how to treat duplicate entries + !! \param info return code + !! + ! + subroutine i_base_mv_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: n, dupl + integer(psb_ipk_), intent(in) :: irl(:) + integer(psb_ipk_), intent(in) :: val(:,:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: i, isz + + info = 0 + if (psb_errstatus_fatal()) return + + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + else if (n > min(size(irl),size(val))) then + info = psb_err_invalid_input_ + + else + isz = size(x%v,1) + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if ((1 <= irl(i)).and.(irl(i) <= isz)) then + ! this row belongs to me + ! copy i-th row of block val in x + x%v(irl(i),:) = val(i,:) + end if + enddo + + case(psb_dupl_add_) + + do i = 1, n + !loop over all val's rows + if ((1 <= irl(i)).and.(irl(i) <= isz)) then + ! this row belongs to me + ! copy i-th row of block val in x + x%v(irl(i),:) = x%v(irl(i),:) + val(i,:) + end if + enddo + + case default + info = 321 +! !$ call psb_errpush(info,name) +! !$ goto 9999 + end select + end if + if (info /= 0) then + call psb_errpush(info,'base_mv_vect_ins') + return + end if + + end subroutine i_base_mv_ins + + ! + !> Function base_mv_zero + !! \memberof psb_i_base_multivect_type + !! \brief Zero out contents + !! + ! + subroutine i_base_mv_zero(x) + use psi_serial_mod + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=izero + + end subroutine i_base_mv_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + !> Function base_mv_asb: + !! \memberof psb_i_base_multivect_type + !! \brief Assemble vector: reallocate as necessary. + !! + !! \param n final size + !! \param info return code + !! + ! + + subroutine i_base_mv_asb(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + if ((x%get_nrows() < m).or.(x%get_ncols() Function base_mv_free: + !! \memberof psb_i_base_multivect_type + !! \brief Free vector + !! + !! \param info return code + !! + ! + subroutine i_base_mv_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v)) deallocate(x%v, stat=info) + if (info /= 0) call & + & psb_errpush(psb_err_alloc_dealloc_,'vect_free') + + end subroutine i_base_mv_free + + + + ! + ! The base version of SYNC & friends does nothing, it's just + ! a placeholder. + ! + ! + !> Function base_mv_sync: + !! \memberof psb_i_base_multivect_type + !! \brief Sync: base version is a no-op. + !! + ! + subroutine i_base_mv_sync(x) + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + + end subroutine i_base_mv_sync + + ! + !> Function base_mv_set_host: + !! \memberof psb_i_base_multivect_type + !! \brief Set_host: base version is a no-op. + !! + ! + subroutine i_base_mv_set_host(x) + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + + end subroutine i_base_mv_set_host + + ! + !> Function base_mv_set_dev: + !! \memberof psb_i_base_multivect_type + !! \brief Set_dev: base version is a no-op. + !! + ! + subroutine i_base_mv_set_dev(x) + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + + end subroutine i_base_mv_set_dev + + ! + !> Function base_mv_set_sync: + !! \memberof psb_i_base_multivect_type + !! \brief Set_sync: base version is a no-op. + !! + ! + subroutine i_base_mv_set_sync(x) + implicit none + class(psb_i_base_multivect_type), intent(inout) :: x + + end subroutine i_base_mv_set_sync + + ! + !> Function base_mv_is_dev: + !! \memberof psb_i_base_multivect_type + !! \brief Is vector on external device . + !! + ! + function i_base_mv_is_dev(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + logical :: res + + res = .false. + end function i_base_mv_is_dev + + ! + !> Function base_mv_is_host + !! \memberof psb_i_base_multivect_type + !! \brief Is vector on standard memory . + !! + ! + function i_base_mv_is_host(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + logical :: res + + res = .true. + end function i_base_mv_is_host + + ! + !> Function base_mv_is_sync + !! \memberof psb_i_base_multivect_type + !! \brief Is vector on sync . + !! + ! + function i_base_mv_is_sync(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + logical :: res + + res = .true. + end function i_base_mv_is_sync + + + ! + ! Size info. + ! + ! + !> Function base_mv_get_nrows + !! \memberof psb_i_base_multivect_type + !! \brief Number of entries + !! + ! + function i_base_mv_get_nrows(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + + res = 0 + if (allocated(x%v)) res = size(x%v,1) + + end function i_base_mv_get_nrows + + function i_base_mv_get_ncols(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + + res = 0 + if (allocated(x%v)) res = size(x%v,2) + + end function i_base_mv_get_ncols + + ! + !> Function base_mv_get_sizeof + !! \memberof psb_i_base_multivect_type + !! \brief Size in bytesa + !! + ! + function i_base_mv_sizeof(x) result(res) + implicit none + class(psb_i_base_multivect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * psb_sizeof_int) * x%get_nrows() * x%get_ncols() + + end function i_base_mv_sizeof + + ! + !> Function base_mv_get_fmt + !! \memberof psb_i_base_multivect_type + !! \brief Format + !! + ! + function i_base_mv_get_fmt() result(res) + implicit none + character(len=5) :: res + res = 'BASE' + end function i_base_mv_get_fmt + + + ! + ! + ! + !> Function base_mv_get_vect + !! \memberof psb_i_base_multivect_type + !! \brief Extract a copy of the contents + !! + ! + function i_base_mv_get_vect(x) result(res) + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), allocatable :: res(:,:) + integer(psb_ipk_) :: info + + if (.not.allocated(x%v)) return + call x%sync() + allocate(res(x%get_nrows(),x%get_ncols()),stat=info) + if (info /= 0) then + call psb_errpush(psb_err_alloc_dealloc_,'base_mv_get_vect') + return + end if + res(:,:) = x%v(:,:) + end function i_base_mv_get_vect + + ! + ! Reset all values + ! + ! + !> Function base_mv_set_scal + !! \memberof psb_i_base_multivect_type + !! \brief Set all entries + !! \param val The value to set + !! + subroutine i_base_mv_set_scal(x,val) + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: val + + integer(psb_ipk_) :: info + x%v = val + + end subroutine i_base_mv_set_scal + + ! + !> Function base_mv_set_vect + !! \memberof psb_i_base_multivect_type + !! \brief Set all entries + !! \param val(:) The vector to be copied in + !! + subroutine i_base_mv_set_vect(x,val) + class(psb_i_base_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: val(:,:) + integer(psb_ipk_) :: nr + integer(psb_ipk_) :: info + + if (allocated(x%v)) then + nr = min(size(x%v,1),size(val,1)) + nc = min(size(x%v,2),size(val,2)) + + x%v(1:nr,1:nc) = val(1:nr,1:nc) + else + x%v = val + end if + + end subroutine i_base_mv_set_vect + +!!$ ! +!!$ ! Dot products +!!$ ! +!!$ ! +!!$ !> Function base_mv_dot_v +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief Dot product by another base_mv_vector +!!$ !! \param n Number of entries to be considere +!!$ !! \param y The other (base_mv_vect) to be multiplied by +!!$ !! +!!$ function i_base_mv_dot_v(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x, y +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ integer(psb_ipk_), external :: idot +!!$ +!!$ res = izero +!!$ ! +!!$ ! Note: this is the base implementation. +!!$ ! When we get here, we are sure that X is of +!!$ ! TYPE psb_i_base_mv_vect. +!!$ ! If Y is not, throw the burden on it, implicitly +!!$ ! calling dot_a +!!$ ! +!!$ select type(yy => y) +!!$ type is (psb_i_base_multivect_type) +!!$ res = idot(n,x%v,1,y%v,1) +!!$ class default +!!$ res = y%dot(n,x%v) +!!$ end select +!!$ +!!$ end function i_base_mv_dot_v +!!$ +!!$ ! +!!$ ! Base workhorse is good old BLAS1 +!!$ ! +!!$ ! +!!$ !> Function base_mv_dot_a +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief Dot product by a normal array +!!$ !! \param n Number of entries to be considere +!!$ !! \param y(:) The array to be multiplied by +!!$ !! +!!$ function i_base_mv_dot_a(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ integer(psb_ipk_), external :: idot +!!$ +!!$ res = idot(n,y,1,x%v,1) +!!$ +!!$ end function i_base_mv_dot_a + +!!$ ! +!!$ ! AXPBY is invoked via Y, hence the structure below. +!!$ ! +!!$ ! +!!$ ! +!!$ !> Function base_mv_axpby_v +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief AXPBY by a (base_mv_vect) y=alpha*x+beta*y +!!$ !! \param m Number of entries to be considere +!!$ !! \param alpha scalar alpha +!!$ !! \param x The class(base_mv_vect) to be added +!!$ !! \param beta scalar alpha +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_axpby_v(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ select type(xx => x) +!!$ type is (psb_i_base_multivect_type) +!!$ call psb_geaxpby(m,alpha,x%v,beta,y%v,info) +!!$ class default +!!$ call y%axpby(m,alpha,x%v,beta,info) +!!$ end select +!!$ +!!$ end subroutine i_base_mv_axpby_v +!!$ +!!$ ! +!!$ ! AXPBY is invoked via Y, hence the structure below. +!!$ ! +!!$ ! +!!$ !> Function base_mv_axpby_a +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief AXPBY by a normal array y=alpha*x+beta*y +!!$ !! \param m Number of entries to be considere +!!$ !! \param alpha scalar alpha +!!$ !! \param x(:) The array to be added +!!$ !! \param beta scalar alpha +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_axpby_a(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ call psb_geaxpby(m,alpha,x,beta,y%v,info) +!!$ +!!$ end subroutine i_base_mv_axpby_a + + +!!$ ! +!!$ ! Multiple variants of two operations: +!!$ ! Simple multiplication Y(:) = X(:)*Y(:) +!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:) +!!$ ! +!!$ ! Variants expanded according to the dynamic type +!!$ ! of the involved entities +!!$ ! +!!$ ! +!!$ !> Function base_mv_mlt_a +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief Vector entry-by-entry multiply by a base_mv_vect array y=x*y +!!$ !! \param x The class(base_mv_vect) to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_mlt_v(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ select type(xx => x) +!!$ type is (psb_i_base_multivect_type) +!!$ n = min(size(y%v), size(xx%v)) +!!$ do i=1, n +!!$ y%v(i) = y%v(i)*xx%v(i) +!!$ end do +!!$ class default +!!$ call y%mlt(x%v,info) +!!$ end select +!!$ +!!$ end subroutine i_base_mv_mlt_v +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_a +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y +!!$ !! \param x(:) The array to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_mlt_a(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ n = min(size(y%v), size(x)) +!!$ do i=1, n +!!$ y%v(i) = y%v(i)*x(i) +!!$ end do +!!$ +!!$ end subroutine i_base_mv_mlt_a +!!$ +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_a_2 +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays +!!$ !! z=beta*z+alpha*x*y +!!$ !! \param alpha +!!$ !! \param beta +!!$ !! \param x(:) The array to be multiplied b +!!$ !! \param y(:) The array to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_mlt_a_2(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ n = min(size(z%v), size(x), size(y)) +!!$ if (alpha == izero) then +!!$ if (beta == ione) then +!!$ return +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) +!!$ end do +!!$ end if +!!$ else +!!$ if (alpha == ione) then +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) + y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) + y(i)*x(i) +!!$ end do +!!$ end if +!!$ else if (alpha == -ione) then +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = -y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) - y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) - y(i)*x(i) +!!$ end do +!!$ end if +!!$ else +!!$ if (beta == izero) then +!!$ do i=1, n +!!$ z%v(i) = alpha*y(i)*x(i) +!!$ end do +!!$ else if (beta == ione) then +!!$ do i=1, n +!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i) +!!$ end do +!!$ else +!!$ do i=1, n +!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i) +!!$ end do +!!$ end if +!!$ end if +!!$ end if +!!$ end subroutine i_base_mv_mlt_a_2 +!!$ +!!$ ! +!!$ !> Function base_mv_mlt_v_2 +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mv_vect) +!!$ !! z=beta*z+alpha*x*y +!!$ !! \param alpha +!!$ !! \param beta +!!$ !! \param x The class(base_mv_vect) to be multiplied b +!!$ !! \param y The class(base_mv_vect) to be multiplied by +!!$ !! \param info return code +!!$ !! +!!$ subroutine i_base_mv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy) +!!$ use psi_serial_mod +!!$ use psb_string_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ class(psb_i_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character(len=1), intent(in), optional :: conjgx, conjgy +!!$ integer(psb_ipk_) :: i, n +!!$ logical :: conjgx_, conjgy_ +!!$ +!!$ info = 0 +!!$ if (.not.psb_i_is_complex_) then +!!$ call z%mlt(alpha,x%v,y%v,beta,info) +!!$ else +!!$ conjgx_=.false. +!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C') +!!$ conjgy_=.false. +!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C') +!!$ if (conjgx_) x%v=(x%v) +!!$ if (conjgy_) y%v=(y%v) +!!$ call z%mlt(alpha,x%v,y%v,beta,info) +!!$ if (conjgx_) x%v=(x%v) +!!$ if (conjgy_) y%v=(y%v) +!!$ end if +!!$ end subroutine i_base_mv_mlt_v_2 +!!$ +!!$ subroutine i_base_mv_mlt_av(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_base_multivect_type), intent(inout) :: y +!!$ class(psb_i_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ call z%mlt(alpha,x,y%v,beta,info) +!!$ +!!$ end subroutine i_base_mv_mlt_av +!!$ +!!$ subroutine i_base_mv_mlt_va(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ class(psb_i_base_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ call z%mlt(alpha,y,x,beta,info) +!!$ +!!$ end subroutine i_base_mv_mlt_va +!!$ +!!$ +!!$ ! +!!$ ! Simple scaling +!!$ ! +!!$ !> Function base_mv_scal +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief Scale all entries x = alpha*x +!!$ !! \param alpha The multiplier +!!$ !! +!!$ subroutine i_base_mv_scal(alpha, x) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent (in) :: alpha +!!$ +!!$ if (allocated(x%v)) x%v = alpha*x%v +!!$ +!!$ end subroutine i_base_mv_scal +!!$ +!!$ ! +!!$ ! Norms 1, 2 and infinity +!!$ ! +!!$ !> Function base_mv_nrm2 +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief 2-norm |x(1:n)|_2 +!!$ !! \param n how many entries to consider +!!$ function i_base_mv_nrm2(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ integer(psb_ipk_), external :: inrm2 +!!$ +!!$ res = inrm2(n,x%v,1) +!!$ +!!$ end function i_base_mv_nrm2 +!!$ +!!$ ! +!!$ !> Function base_mv_amax +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief infinity-norm |x(1:n)|_\infty +!!$ !! \param n how many entries to consider +!!$ function i_base_mv_amax(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = maxval(abs(x%v(1:n))) +!!$ +!!$ end function i_base_mv_amax +!!$ +!!$ ! +!!$ !> Function base_mv_asum +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief 1-norm |x(1:n)|_1 +!!$ !! \param n how many entries to consider +!!$ function i_base_mv_asum(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_base_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = sum(abs(x%v(1:n))) +!!$ +!!$ end function i_base_mv_asum +!!$ +!!$ +!!$ ! +!!$ ! Gather: Y = beta * Y + alpha * X(IDX(:)) +!!$ ! +!!$ ! +!!$ !> Function base_mv_gthab +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief gather into an array +!!$ !! Y = beta * Y + alpha * X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ !! \param alpha +!!$ !! \param beta +!!$ subroutine i_base_mv_gthab(n,idx,alpha,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: alpha, beta, y(:) +!!$ class(psb_i_base_multivect_type) :: x +!!$ +!!$ call x%sync() +!!$ call psi_gth(n,idx,alpha,x%v,beta,y) +!!$ +!!$ end subroutine i_base_mv_gthab +!!$ ! +!!$ ! shortcut alpha=1 beta=0 +!!$ ! +!!$ !> Function base_mv_gthzv +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief gather into an array special alpha=1 beta=0 +!!$ !! Y = X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ subroutine i_base_mv_gthzv_x(i,n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: i,n +!!$ class(psb_i_base_multivect_type) :: idx +!!$ integer(psb_ipk_) :: y(:) +!!$ class(psb_i_base_multivect_type) :: x +!!$ +!!$ call x%gth(n,idx%v(i:),y) +!!$ +!!$ end subroutine i_base_mv_gthzv_x +!!$ +!!$ ! +!!$ ! shortcut alpha=1 beta=0 +!!$ ! +!!$ !> Function base_mv_gthzv +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief gather into an array special alpha=1 beta=0 +!!$ !! Y = X(IDX(:)) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ subroutine i_base_mv_gthzv(n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: y(:) +!!$ class(psb_i_base_multivect_type) :: x +!!$ +!!$ call x%sync() +!!$ call psi_gth(n,idx,x%v,y) +!!$ +!!$ end subroutine i_base_mv_gthzv +!!$ +!!$ ! +!!$ ! Scatter: +!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) +!!$ ! +!!$ ! +!!$ !> Function base_mv_sctb +!!$ !! \memberof psb_i_base_multivect_type +!!$ !! \brief scatter into a class(base_mv_vect) +!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:) +!!$ !! \param n how many entries to consider +!!$ !! \param idx(:) indices +!!$ !! \param beta +!!$ !! \param x(:) +!!$ subroutine i_base_mv_sctb(n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: beta, x(:) +!!$ class(psb_i_base_multivect_type) :: y +!!$ +!!$ call y%sync() +!!$ call psi_sct(n,idx,x,beta,y%v) +!!$ call y%set_host() +!!$ +!!$ end subroutine i_base_mv_sctb +!!$ +!!$ subroutine i_base_mv_sctb_x(i,n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: i, n +!!$ class(psb_i_base_multivect_type) :: idx +!!$ integer(psb_ipk_) :: beta, x(:) +!!$ class(psb_i_base_multivect_type) :: y +!!$ +!!$ call y%sct(n,idx%v(i:),x,beta) +!!$ +!!$ end subroutine i_base_mv_sctb_x + +end module psb_i_base_multivect_mod diff --git a/base/modules/psb_i_vect_mod.F90 b/base/modules/psb_i_vect_mod.F90 index 4b660383..ae8c2a86 100644 --- a/base/modules/psb_i_vect_mod.F90 +++ b/base/modules/psb_i_vect_mod.F90 @@ -653,3 +653,635 @@ contains end subroutine i_vect_cnv end module psb_i_vect_mod + + + +module psb_i_multivect_mod + + use psb_i_base_multivect_mod + use psb_const_mod + + private + + type psb_i_multivect_type + class(psb_i_base_multivect_type), allocatable :: v + contains + procedure, pass(x) :: get_nrows => i_vect_get_nrows + procedure, pass(x) :: get_ncols => i_vect_get_ncols + procedure, pass(x) :: sizeof => i_vect_sizeof + procedure, pass(x) :: get_fmt => i_vect_get_fmt +!!$ procedure, pass(x) :: dot_v => i_vect_dot_v +!!$ procedure, pass(x) :: dot_a => i_vect_dot_a +!!$ generic, public :: dot => dot_v, dot_a +!!$ procedure, pass(y) :: axpby_v => i_vect_axpby_v +!!$ procedure, pass(y) :: axpby_a => i_vect_axpby_a +!!$ generic, public :: axpby => axpby_v, axpby_a +!!$ procedure, pass(y) :: mlt_v => i_vect_mlt_v +!!$ procedure, pass(y) :: mlt_a => i_vect_mlt_a +!!$ procedure, pass(z) :: mlt_a_2 => i_vect_mlt_a_2 +!!$ procedure, pass(z) :: mlt_v_2 => i_vect_mlt_v_2 +!!$ procedure, pass(z) :: mlt_va => i_vect_mlt_va +!!$ procedure, pass(z) :: mlt_av => i_vect_mlt_av +!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2,& +!!$ & mlt_v_2, mlt_av, mlt_va +!!$ procedure, pass(x) :: scal => i_vect_scal +!!$ procedure, pass(x) :: nrm2 => i_vect_nrm2 +!!$ procedure, pass(x) :: amax => i_vect_amax +!!$ procedure, pass(x) :: asum => i_vect_asum + procedure, pass(x) :: all => i_vect_all + procedure, pass(x) :: reall => i_vect_reall + procedure, pass(x) :: zero => i_vect_zero + procedure, pass(x) :: asb => i_vect_asb + procedure, pass(x) :: sync => i_vect_sync +!!$ procedure, pass(x) :: gthab => i_vect_gthab +!!$ procedure, pass(x) :: gthzv => i_vect_gthzv +!!$ generic, public :: gth => gthab, gthzv +!!$ procedure, pass(y) :: sctb => i_vect_sctb +!!$ generic, public :: sct => sctb + procedure, pass(x) :: free => i_vect_free + procedure, pass(x) :: ins => i_vect_ins + procedure, pass(x) :: bld_x => i_vect_bld_x + procedure, pass(x) :: bld_n => i_vect_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: get_vect => i_vect_get_vect + procedure, pass(x) :: cnv => i_vect_cnv + procedure, pass(x) :: set_scal => i_vect_set_scal + procedure, pass(x) :: set_vect => i_vect_set_vect + generic, public :: set => set_vect, set_scal + procedure, pass(x) :: clone => i_vect_clone + end type psb_i_multivect_type + + public :: psb_i_multivect, psb_i_multivect_type,& + & psb_set_multivect_default, psb_get_multivect_default + + interface psb_i_multivect + module procedure constructor, size_const + end interface + + class(psb_i_base_multivect_type), allocatable, target,& + & save, private :: psb_i_base_multivect_default + + interface psb_set_multivect_default + module procedure psb_i_set_multivect_default + end interface + + interface psb_get_vect_default + module procedure psb_i_get_multivect_default + end interface + + +contains + + + subroutine psb_i_set_multivect_default(v) + implicit none + class(psb_i_base_multivect_type), intent(in) :: v + + if (allocated(psb_i_base_multivect_default)) then + deallocate(psb_i_base_multivect_default) + end if + allocate(psb_i_base_multivect_default, mold=v) + + end subroutine psb_i_set_multivect_default + + function psb_i_get_multivect_default(v) result(res) + implicit none + class(psb_i_multivect_type), intent(in) :: v + class(psb_i_base_multivect_type), pointer :: res + + res => psb_i_get_base_multivect_default() + + end function psb_i_get_multivect_default + + + function psb_i_get_base_multivect_default() result(res) + implicit none + class(psb_i_base_multivect_type), pointer :: res + + if (.not.allocated(psb_i_base_multivect_default)) then + allocate(psb_i_base_multivect_type :: psb_i_base_multivect_default) + end if + + res => psb_i_base_multivect_default + + end function psb_i_get_base_multivect_default + + + subroutine i_vect_clone(x,y,info) + implicit none + class(psb_i_multivect_type), intent(inout) :: x + class(psb_i_multivect_type), intent(inout) :: y + integer(psb_ipk_), intent(out) :: info + + info = psb_success_ + call y%free(info) + if ((info==0).and.allocated(x%v)) then + call y%bld(x%get_vect(),mold=x%v) + end if + end subroutine i_vect_clone + + subroutine i_vect_bld_x(x,invect,mold) + integer(psb_ipk_), intent(in) :: invect(:,:) + class(psb_i_multivect_type), intent(out) :: x + class(psb_i_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_) :: info + class(psb_i_base_multivect_type), pointer :: mld + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else +#ifdef HAVE_MOLD + allocate(x%v,stat=info, mold=psb_i_get_base_multivect_default()) +#else + mld = psb_i_get_base_multivect_default() + call mld%mold(x%v,info) +#endif + endif + + if (info == psb_success_) call x%v%bld(invect) + + end subroutine i_vect_bld_x + + + subroutine i_vect_bld_n(x,m,n,mold) + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_multivect_type), intent(out) :: x + class(psb_i_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_) :: info + class(psb_i_base_multivect_type), pointer :: mld + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else +#ifdef HAVE_MOLD + allocate(x%v,stat=info, mold=psb_i_get_base_multivect_default()) +#else + mld = psb_i_get_base_multivect_default() + call mld%mold(x%v,info) +#endif + endif + if (info == psb_success_) call x%v%bld(m,n) + + end subroutine i_vect_bld_n + + function i_vect_get_vect(x) result(res) + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), allocatable :: res(:,:) + integer(psb_ipk_) :: info + + if (allocated(x%v)) then + res = x%v%get_vect() + end if + end function i_vect_get_vect + + subroutine i_vect_set_scal(x,val) + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: val + + integer(psb_ipk_) :: info + if (allocated(x%v)) call x%v%set(val) + + end subroutine i_vect_set_scal + + subroutine i_vect_set_vect(x,val) + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: val(:,:) + + integer(psb_ipk_) :: info + if (allocated(x%v)) call x%v%set(val) + + end subroutine i_vect_set_vect + + + function constructor(x) result(this) + integer(psb_ipk_) :: x(:,:) + type(psb_i_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%bld(x) + call this%asb(size(x,dim=1,kind=psb_ipk_),size(x,dim=2,kind=psb_ipk_),info) + + end function constructor + + + function size_const(m,n) result(this) + integer(psb_ipk_), intent(in) :: m,n + type(psb_i_multivect_type) :: this + integer(psb_ipk_) :: info + + call this%bld(m,n) + call this%asb(m,n,info) + + end function size_const + + function i_vect_get_nrows(x) result(res) + implicit none + class(psb_i_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + res = 0 + if (allocated(x%v)) res = x%v%get_nrows() + end function i_vect_get_nrows + + function i_vect_get_ncols(x) result(res) + implicit none + class(psb_i_multivect_type), intent(in) :: x + integer(psb_ipk_) :: res + res = 0 + if (allocated(x%v)) res = x%v%get_ncols() + end function i_vect_get_ncols + + function i_vect_sizeof(x) result(res) + implicit none + class(psb_i_multivect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + res = 0 + if (allocated(x%v)) res = x%v%sizeof() + end function i_vect_sizeof + + function i_vect_get_fmt(x) result(res) + implicit none + class(psb_i_multivect_type), intent(in) :: x + character(len=5) :: res + res = 'NULL' + if (allocated(x%v)) res = x%v%get_fmt() + end function i_vect_get_fmt + +!!$ function i_vect_dot_v(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x, y +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = izero +!!$ if (allocated(x%v).and.allocated(y%v)) & +!!$ & res = x%v%dot(n,y%v) +!!$ +!!$ end function i_vect_dot_v +!!$ +!!$ function i_vect_dot_a(n,x,y) result(res) +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ res = izero +!!$ if (allocated(x%v)) & +!!$ & res = x%v%dot(n,y) +!!$ +!!$ end function i_vect_dot_a +!!$ +!!$ subroutine i_vect_axpby_v(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ if (allocated(x%v).and.allocated(y%v)) then +!!$ call y%v%axpby(m,alpha,x%v,beta,info) +!!$ else +!!$ info = psb_err_invalid_vect_state_ +!!$ end if +!!$ +!!$ end subroutine i_vect_axpby_v +!!$ +!!$ subroutine i_vect_axpby_a(m,alpha, x, beta, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: m +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent (in) :: alpha, beta +!!$ integer(psb_ipk_), intent(out) :: info +!!$ +!!$ if (allocated(y%v)) & +!!$ & call y%v%axpby(m,alpha,x,beta,info) +!!$ +!!$ end subroutine i_vect_axpby_a +!!$ +!!$ +!!$ subroutine i_vect_mlt_v(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(x%v).and.allocated(y%v)) & +!!$ & call y%v%mlt(x%v,info) +!!$ +!!$ end subroutine i_vect_mlt_v +!!$ +!!$ subroutine i_vect_mlt_a(x, y, info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ +!!$ info = 0 +!!$ if (allocated(y%v)) & +!!$ & call y%v%mlt(x,info) +!!$ +!!$ end subroutine i_vect_mlt_a +!!$ +!!$ +!!$ subroutine i_vect_mlt_a_2(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(z%v)) & +!!$ & call z%v%mlt(alpha,x,y,beta,info) +!!$ +!!$ end subroutine i_vect_mlt_a_2 +!!$ +!!$ subroutine i_vect_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ class(psb_i_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ character(len=1), intent(in), optional :: conjgx, conjgy +!!$ +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(x%v).and.allocated(y%v).and.& +!!$ & allocated(z%v)) & +!!$ & call z%v%mlt(alpha,x%v,y%v,beta,info,conjgx,conjgy) +!!$ +!!$ end subroutine i_vect_mlt_v_2 +!!$ +!!$ subroutine i_vect_mlt_av(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: x(:) +!!$ class(psb_i_multivect_type), intent(inout) :: y +!!$ class(psb_i_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ if (allocated(z%v).and.allocated(y%v)) & +!!$ & call z%v%mlt(alpha,x,y%v,beta,info) +!!$ +!!$ end subroutine i_vect_mlt_av +!!$ +!!$ subroutine i_vect_mlt_va(alpha,x,y,beta,z,info) +!!$ use psi_serial_mod +!!$ implicit none +!!$ integer(psb_ipk_), intent(in) :: alpha,beta +!!$ integer(psb_ipk_), intent(in) :: y(:) +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ class(psb_i_multivect_type), intent(inout) :: z +!!$ integer(psb_ipk_), intent(out) :: info +!!$ integer(psb_ipk_) :: i, n +!!$ +!!$ info = 0 +!!$ +!!$ if (allocated(z%v).and.allocated(x%v)) & +!!$ & call z%v%mlt(alpha,x%v,y,beta,info) +!!$ +!!$ end subroutine i_vect_mlt_va +!!$ +!!$ subroutine i_vect_scal(alpha, x) +!!$ use psi_serial_mod +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent (in) :: alpha +!!$ +!!$ if (allocated(x%v)) call x%v%scal(alpha) +!!$ +!!$ end subroutine i_vect_scal +!!$ +!!$ +!!$ function i_vect_nrm2(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%nrm2(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function i_vect_nrm2 +!!$ +!!$ function i_vect_amax(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%amax(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function i_vect_amax +!!$ +!!$ function i_vect_asum(n,x) result(res) +!!$ implicit none +!!$ class(psb_i_multivect_type), intent(inout) :: x +!!$ integer(psb_ipk_), intent(in) :: n +!!$ integer(psb_ipk_) :: res +!!$ +!!$ if (allocated(x%v)) then +!!$ res = x%v%asum(n) +!!$ else +!!$ res = izero +!!$ end if +!!$ +!!$ end function i_vect_asum + + subroutine i_vect_all(m,n, x, info, mold) + + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_multivect_type), intent(out) :: x + class(psb_i_base_multivect_type), intent(in), optional :: mold + integer(psb_ipk_), intent(out) :: info + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(x%v,stat=info,mold=mold) +#else + call mold%mold(x%v,info) +#endif + else + allocate(psb_i_base_multivect_type :: x%v,stat=info) + endif + if (info == 0) then + call x%v%all(m,n,info) + else + info = psb_err_alloc_dealloc_ + end if + + end subroutine i_vect_all + + subroutine i_vect_reall(m,n, x, info) + + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (.not.allocated(x%v)) & + & call x%all(m,n,info) + if (info == 0) & + & call x%asb(m,n,info) + + end subroutine i_vect_reall + + subroutine i_vect_zero(x) + use psi_serial_mod + implicit none + class(psb_i_multivect_type), intent(inout) :: x + + if (allocated(x%v)) call x%v%zero() + + end subroutine i_vect_zero + + subroutine i_vect_asb(m,n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer(psb_ipk_), intent(in) :: m,n + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v)) & + & call x%v%asb(m,n,info) + + end subroutine i_vect_asb + + subroutine i_vect_sync(x) + implicit none + class(psb_i_multivect_type), intent(inout) :: x + + if (allocated(x%v)) & + & call x%v%sync() + + end subroutine i_vect_sync + +!!$ subroutine i_vect_gthab(n,idx,alpha,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: alpha, beta, y(:) +!!$ class(psb_i_multivect_type) :: x +!!$ +!!$ if (allocated(x%v)) & +!!$ & call x%v%gth(n,idx,alpha,beta,y) +!!$ +!!$ end subroutine i_vect_gthab +!!$ +!!$ subroutine i_vect_gthzv(n,idx,x,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: y(:) +!!$ class(psb_i_multivect_type) :: x +!!$ +!!$ if (allocated(x%v)) & +!!$ & call x%v%gth(n,idx,y) +!!$ +!!$ end subroutine i_vect_gthzv +!!$ +!!$ subroutine i_vect_sctb(n,idx,x,beta,y) +!!$ use psi_serial_mod +!!$ integer(psb_ipk_) :: n, idx(:) +!!$ integer(psb_ipk_) :: beta, x(:) +!!$ class(psb_i_multivect_type) :: y +!!$ +!!$ if (allocated(y%v)) & +!!$ & call y%v%sct(n,idx,x,beta) +!!$ +!!$ end subroutine i_vect_sctb + + subroutine i_vect_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(out) :: info + + info = 0 + if (allocated(x%v)) then + call x%v%free(info) + if (info == 0) deallocate(x%v,stat=info) + end if + + end subroutine i_vect_free + + subroutine i_vect_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_i_multivect_type), intent(inout) :: x + integer(psb_ipk_), intent(in) :: n, dupl + integer(psb_ipk_), intent(in) :: irl(:) + integer(psb_ipk_), intent(in) :: val(:,:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: i + + info = 0 + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + return + end if + + call x%v%ins(n,irl,val,dupl,info) + + end subroutine i_vect_ins + + + subroutine i_vect_cnv(x,mold) + class(psb_i_multivect_type), intent(inout) :: x + class(psb_i_base_multivect_type), intent(in), optional :: mold + class(psb_i_base_multivect_type), allocatable :: tmp + integer(psb_ipk_) :: info + + if (present(mold)) then +#ifdef HAVE_MOLD + allocate(tmp,stat=info,mold=mold) +#else + call mold%mold(tmp,info) +#endif + if (allocated(x%v)) then + call x%v%sync() + if (info == psb_success_) call tmp%bld(x%v%v) + call x%v%free(info) + end if + call move_alloc(tmp,x%v) + end if + end subroutine i_vect_cnv + +end module psb_i_multivect_mod diff --git a/base/modules/psb_vect_mod.f90 b/base/modules/psb_vect_mod.f90 index 77f1f65b..2e5297bc 100644 --- a/base/modules/psb_vect_mod.f90 +++ b/base/modules/psb_vect_mod.f90 @@ -4,4 +4,6 @@ module psb_vect_mod use psb_d_vect_mod use psb_c_vect_mod use psb_z_vect_mod + use psb_i_multivect_mod + use psb_d_multivect_mod end module psb_vect_mod diff --git a/base/serial/impl/psb_d_base_mat_impl.F90 b/base/serial/impl/psb_d_base_mat_impl.F90 index 7ec570c8..f4280fff 100644 --- a/base/serial/impl/psb_d_base_mat_impl.F90 +++ b/base/serial/impl/psb_d_base_mat_impl.F90 @@ -380,6 +380,9 @@ subroutine psb_d_base_csput_v(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) info = psb_success_ if (allocated(val%v).and.allocated(ia%v).and.allocated(ja%v)) then + if (val%is_dev()) call val%sync() + if (ia%is_dev()) call ia%sync() + if (ja%is_dev()) call ja%sync() call a%csput(nz,ia%v,ja%v,val%v,imin,imax,jmin,jmax,info,gtl) else info = psb_err_invalid_mat_state_ diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index cc39bfa3..5f974803 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -733,8 +733,8 @@ end subroutine psb_d_trim -subroutine psb_d_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) - use psb_d_mat_mod, psb_protect_name => psb_d_csput +subroutine psb_d_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_d_mat_mod, psb_protect_name => psb_d_csput_a use psb_d_base_mat_mod use psb_error_mod implicit none @@ -745,7 +745,7 @@ subroutine psb_d_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) integer(psb_ipk_), intent(in), optional :: gtl(:) integer(psb_ipk_) :: err_act - character(len=20) :: name='csput' + character(len=20) :: name='csput_a' logical, parameter :: debug=.false. info = psb_success_ @@ -771,7 +771,54 @@ subroutine psb_d_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) return end if -end subroutine psb_d_csput +end subroutine psb_d_csput_a + +subroutine psb_d_csput_v(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) + use psb_d_mat_mod, psb_protect_name => psb_d_csput_v + use psb_d_base_mat_mod + use psb_d_vect_mod, only : psb_d_vect_type + use psb_i_vect_mod, only : psb_i_vect_type + use psb_error_mod + implicit none + class(psb_dspmat_type), intent(inout) :: a + type(psb_d_vect_type), intent(inout) :: val + type(psb_i_vect_type), intent(inout) :: ia, ja + integer(psb_ipk_), intent(in) :: nz, imin,imax,jmin,jmax + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: gtl(:) + + integer(psb_ipk_) :: err_act + character(len=20) :: name='csput_a' + logical, parameter :: debug=.false. + + info = psb_success_ + call psb_erractionsave(err_act) + if (.not.(a%is_bld().or.a%is_upd())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + if (allocated(val%v).and.allocated(ia%v).and.allocated(ja%v)) then + call a%a%csput(nz,ia%v,ja%v,val%v,imin,imax,jmin,jmax,info,gtl) + else + info = psb_err_invalid_mat_state_ + endif + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + +end subroutine psb_d_csput_v subroutine psb_d_csgetptn(imin,imax,a,nz,ia,ja,info,& diff --git a/base/tools/psb_dspins.f90 b/base/tools/psb_dspins.f90 index f6ee843b..7c6b4d17 100644 --- a/base/tools/psb_dspins.f90 +++ b/base/tools/psb_dspins.f90 @@ -332,3 +332,169 @@ subroutine psb_dspins_2desc(nz,ia,ja,val,a,desc_ar,desc_ac,info) end subroutine psb_dspins_2desc +subroutine psb_dspins_v(nz,ia,ja,val,a,desc_a,info,rebuild,local) + use psb_base_mod, psb_protect_name => psb_dspins_v + use psi_mod + implicit none + + !....parameters... + type(psb_desc_type), intent(inout) :: desc_a + type(psb_dspmat_type), intent(inout) :: a + integer(psb_ipk_), intent(in) :: nz + type(psb_i_vect_type), intent(inout) :: ia,ja + type(psb_d_vect_type), intent(inout) :: val + integer(psb_ipk_), intent(out) :: info + logical, intent(in), optional :: rebuild, local + !locals..... + + integer(psb_ipk_) :: nrow, err_act, ncol, spstate + integer(psb_ipk_) :: ictxt,np,me + logical, parameter :: debug=.false. + integer(psb_ipk_), parameter :: relocsz=200 + logical :: rebuild_, local_ + integer(psb_ipk_), allocatable :: ila(:),jla(:) + real(psb_dpk_) :: t1,t2,t3,tcnv,tcsput + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name + + info = psb_success_ + name = 'psb_dspins' + call psb_erractionsave(err_act) + + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + + if (nz < 0) then + info = 1111 + call psb_errpush(info,name) + goto 9999 + end if + if (ia%get_nrows() < nz) then + info = 1111 + call psb_errpush(info,name) + goto 9999 + end if + + if (ja%get_nrows() < nz) then + info = 1111 + call psb_errpush(info,name) + goto 9999 + end if + if (val%get_nrows() < nz) then + info = 1111 + call psb_errpush(info,name) + goto 9999 + end if + if (nz == 0) return + + if (present(rebuild)) then + rebuild_ = rebuild + else + rebuild_ = .false. + endif + + if (present(local)) then + local_ = local + else + local_ = .false. + endif + + if (desc_a%is_bld()) then + +!!$ if (local_) then + info = psb_err_invalid_a_and_cd_state_ + call psb_errpush(info,name) + goto 9999 +!!$ else +!!$ allocate(ila(nz),jla(nz),stat=info) +!!$ if (info /= psb_success_) then +!!$ ierr(1) = info +!!$ call psb_errpush(psb_err_from_subroutine_ai_,name,& +!!$ & a_err='allocate',i_err=ierr) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info,owned=.true.) +!!$ call desc_a%indxmap%g2l_ins(ja(1:nz),jla(1:nz),info,mask=(ila(1:nz)>0)) +!!$ +!!$ if (info /= psb_success_) then +!!$ ierr(1) = info +!!$ call psb_errpush(psb_err_from_subroutine_ai_,name,& +!!$ & a_err='psb_cdins',i_err=ierr) +!!$ goto 9999 +!!$ end if +!!$ nrow = desc_a%get_local_rows() +!!$ ncol = desc_a%get_local_cols() +!!$ +!!$ if (a%is_bld()) then +!!$ call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info) +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ call psb_errpush(info,name,a_err='a%csput') +!!$ goto 9999 +!!$ end if +!!$ else +!!$ info = psb_err_invalid_a_and_cd_state_ +!!$ call psb_errpush(info,name) +!!$ goto 9999 +!!$ end if +!!$ endif + + else if (desc_a%is_asb()) then + + nrow = desc_a%get_local_rows() + ncol = desc_a%get_local_cols() + if (local_) then + t1=psb_wtime() + call a%csput(nz,ia,ja,val,ione,nrow,ione,ncol,info) + tcsput=psb_wtime() - t1 + tcnv=0.0 + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='a%csput') + goto 9999 + end if + else + info = psb_err_invalid_cd_state_ + +!!$ t1=psb_wtime() +!!$ allocate(ila(nz),jla(nz),stat=info) +!!$ if (info /= psb_success_) then +!!$ ierr(1) = info +!!$ call psb_errpush(psb_err_from_subroutine_ai_,name,& +!!$ & a_err='allocate',i_err=ierr) +!!$ goto 9999 +!!$ end if +!!$ +!!$ call desc_a%indxmap%g2l(ia(1:nz),ila(1:nz),info) +!!$ call desc_a%indxmap%g2l(ja(1:nz),jla(1:nz),info) +!!$ t2 = psb_Wtime() +!!$ call a%csput(nz,ila,jla,val,ione,nrow,ione,ncol,info) +!!$ t3=psb_wtime() +!!$ tcnv=t2-t1 +!!$ tcsput=t3-t2 +!!$ if (info /= psb_success_) then +!!$ info=psb_err_from_subroutine_ +!!$ call psb_errpush(info,name,a_err='a%csput') +!!$ goto 9999 +!!$ end if + end if +!!$ write(0,*)'SPINS times: ',tcnv,tcsput + else + info = psb_err_invalid_cd_state_ + 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(ictxt) + return + end if + return + +end subroutine psb_dspins_v