diff --git a/base/modules/psb_c_base_vect_mod.f90 b/base/modules/psb_c_base_vect_mod.f90 index 0dba48b6..e3e1e818 100644 --- a/base/modules/psb_c_base_vect_mod.f90 +++ b/base/modules/psb_c_base_vect_mod.f90 @@ -1,3 +1,47 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! package: psb_c_base_vect_mod +! +! This module contains the definition of the psb_c_base_vect type which +! is a container for dense vectors. +! This is encapsulated instead of being just a simple array to allow for +! more 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. +! +! module psb_c_base_vect_mod use psb_const_mod @@ -6,14 +50,58 @@ module psb_c_base_vect_mod type psb_c_base_vect_type complex(psb_spk_), allocatable :: v(:) contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => c_base_bld_x + procedure, pass(x) :: bld_n => c_base_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => c_base_all + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => c_base_ins + procedure, pass(x) :: zero => c_base_zero + procedure, pass(x) :: asb => c_base_asb + procedure, pass(x) :: free => c_base_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 => c_base_sync + ! + ! Basic info procedure, pass(x) :: get_nrows => c_base_get_nrows procedure, pass(x) :: sizeof => c_base_sizeof + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: getCopy => c_base_getCopy + procedure, pass(x) :: cpy_vect => c_base_cpy_vect + procedure, pass(x) :: set_scal => c_base_set_scal + procedure, pass(x) :: set_vect => c_base_set_vect + generic, public :: set => set_vect, set_scal + generic, public :: assignment(=) => cpy_vect, set_scal + + ! + ! Dot product and AXPBY + ! procedure, pass(x) :: dot_v => c_base_dot_v procedure, pass(x) :: dot_a => c_base_dot_a generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => c_base_axpby_v procedure, pass(y) :: axpby_a => c_base_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 => c_base_mlt_v procedure, pass(y) :: mlt_a => c_base_mlt_a procedure, pass(z) :: mlt_a_2 => c_base_mlt_a_2 @@ -21,30 +109,22 @@ module psb_c_base_vect_mod procedure, pass(z) :: mlt_va => c_base_mlt_va procedure, pass(z) :: mlt_av => c_base_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 => c_base_scal procedure, pass(x) :: nrm2 => c_base_nrm2 procedure, pass(x) :: amax => c_base_amax procedure, pass(x) :: asum => c_base_asum - procedure, pass(x) :: all => c_base_all - procedure, pass(x) :: zero => c_base_zero - procedure, pass(x) :: asb => c_base_asb - procedure, pass(x) :: sync => c_base_sync + ! + ! Gather/scatter. These are needed for MPI interfacing. + ! May have to be reworked. + ! procedure, pass(x) :: gthab => c_base_gthab procedure, pass(x) :: gthzv => c_base_gthzv generic, public :: gth => gthab, gthzv procedure, pass(y) :: sctb => c_base_sctb generic, public :: sct => sctb - procedure, pass(x) :: free => c_base_free - procedure, pass(x) :: ins => c_base_ins - procedure, pass(x) :: bld_x => c_base_bld_x - procedure, pass(x) :: bld_n => c_base_bld_n - generic, public :: bld => bld_x, bld_n - procedure, pass(x) :: getCopy => c_base_getCopy - procedure, pass(x) :: cpy_vect => c_base_cpy_vect - generic, public :: assignment(=) => cpy_vect, set_scal - procedure, pass(x) :: set_scal => c_base_set_scal - procedure, pass(x) :: set_vect => c_base_set_vect - generic, public :: set => set_vect, set_scal end type psb_c_base_vect_type public :: psb_c_base_vect @@ -55,6 +135,33 @@ module psb_c_base_vect_mod contains + ! + ! Constructors. + ! + + function constructor(x) result(this) + complex(psb_spk_) :: x(:) + type(psb_c_base_vect_type) :: this + integer :: info + + this%v = x + call this%asb(size(x),info) + end function constructor + + + function size_const(n) result(this) + integer, intent(in) :: n + type(psb_c_base_vect_type) :: this + integer :: info + + call this%asb(n,info) + + end function size_const + + ! + ! Build from a sample + ! + subroutine c_base_bld_x(x,this) use psb_realloc_mod complex(psb_spk_), intent(in) :: this(:) @@ -70,7 +177,9 @@ contains end subroutine c_base_bld_x - + ! + ! Create with size, but no initialization + ! subroutine c_base_bld_n(x,n) use psb_realloc_mod integer, intent(in) :: n @@ -81,9 +190,172 @@ contains call x%asb(n,info) end subroutine c_base_bld_n + + subroutine c_base_all(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_c_base_vect_type), intent(out) :: x + integer, intent(out) :: info + + call psb_realloc(n,x%v,info) + + end subroutine c_base_all + + ! + ! Insert a bunch of values at specified positions. + ! + subroutine c_base_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + integer, intent(in) :: n, dupl + integer, intent(in) :: irl(:) + complex(psb_spk_), intent(in) :: val(:) + integer, intent(out) :: info + + integer :: i + + 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 + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') + return + end if + + end subroutine c_base_ins + + ! + subroutine c_base_zero(x) + use psi_serial_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=czero + + end subroutine c_base_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + + subroutine c_base_asb(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_c_base_vect_type), intent(inout) :: x + integer, intent(out) :: info + + if (x%get_nrows() < n) & + & call psb_realloc(n,x%v,info) + if (info /= 0) & + & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') + + end subroutine c_base_asb + + + subroutine c_base_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + integer, 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 c_base_free + + + + ! + ! The base version of SYNC does nothing, it's just + ! a placeholder. + ! + subroutine c_base_sync(x) + implicit none + class(psb_c_base_vect_type), intent(inout) :: x + + + end subroutine c_base_sync + + ! + ! Size info. + ! + + function c_base_get_nrows(x) result(res) + implicit none + class(psb_c_base_vect_type), intent(in) :: x + integer :: res + + res = 0 + if (allocated(x%v)) res = size(x%v) + + end function c_base_get_nrows + + function c_base_sizeof(x) result(res) + implicit none + class(psb_c_base_vect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * (2*psb_sizeof_sp)) * x%get_nrows() + + end function c_base_sizeof + + + ! + ! Two versions of extracting an array: one of them + ! overload the assignment. + ! function c_base_getCopy(x) result(res) - class(psb_c_base_vect_type), intent(in) :: x + class(psb_c_base_vect_type), intent(in) :: x complex(psb_spk_), allocatable :: res(:) integer :: info @@ -104,6 +376,9 @@ contains end subroutine c_base_cpy_vect + ! + ! Reset all values + ! subroutine c_base_set_scal(x,val) class(psb_c_base_vect_type), intent(inout) :: x complex(psb_spk_), intent(in) :: val @@ -127,46 +402,10 @@ contains end if end subroutine c_base_set_vect - - - function constructor(x) result(this) - complex(psb_spk_) :: x(:) - type(psb_c_base_vect_type) :: this - integer :: info - - this%v = x - call this%asb(size(x),info) - end function constructor - - - function size_const(n) result(this) - integer, intent(in) :: n - type(psb_c_base_vect_type) :: this - integer :: info - - call this%asb(n,info) - - end function size_const - - function c_base_get_nrows(x) result(res) - implicit none - class(psb_c_base_vect_type), intent(in) :: x - integer :: res - - res = 0 - if (allocated(x%v)) res = size(x%v) - - end function c_base_get_nrows - - function c_base_sizeof(x) result(res) - implicit none - class(psb_c_base_vect_type), intent(in) :: x - integer(psb_long_int_k_) :: res - - res = (1_psb_long_int_k_ * (2*psb_sizeof_sp)) * x%get_nrows() - - end function c_base_sizeof + ! + ! Dot products + ! function c_base_dot_v(n,x,y) result(res) implicit none class(psb_c_base_vect_type), intent(inout) :: x, y @@ -178,7 +417,9 @@ contains ! ! Note: this is the base implementation. ! When we get here, we are sure that X is of - ! TYPE psb_c_base_vect + ! TYPE psb_c_base_vect. + ! If Y is not, throw the burden on it, implicitly + ! calling dot_a ! select type(yy => y) type is (psb_c_base_vect_type) @@ -189,6 +430,9 @@ contains end function c_base_dot_v + ! + ! Base workhorse is good old BLAS1 + ! function c_base_dot_a(n,x,y) result(res) implicit none class(psb_c_base_vect_type), intent(inout) :: x @@ -201,6 +445,9 @@ contains end function c_base_dot_a + ! + ! AXPBY is invoked via Y, hence the structure below. + ! subroutine c_base_axpby_v(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -232,7 +479,16 @@ contains end subroutine c_base_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 + ! + subroutine c_base_mlt_v(x, y, info) use psi_serial_mod implicit none @@ -400,6 +656,11 @@ contains end subroutine c_base_mlt_va + + ! + ! Simple scaling + ! + subroutine c_base_scal(alpha, x) use psi_serial_mod implicit none @@ -409,7 +670,10 @@ contains if (allocated(x%v)) x%v = alpha*x%v end subroutine c_base_scal - + + ! + ! Norms 1, 2 and infinity + ! function c_base_nrm2(n,x) result(res) implicit none @@ -442,52 +706,10 @@ contains end function c_base_asum - subroutine c_base_all(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_c_base_vect_type), intent(out) :: x - integer, intent(out) :: info - - call psb_realloc(n,x%v,info) - - end subroutine c_base_all - - subroutine c_base_zero(x) - use psi_serial_mod - implicit none - class(psb_c_base_vect_type), intent(inout) :: x - - if (allocated(x%v)) x%v=czero - - end subroutine c_base_zero - - subroutine c_base_asb(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_c_base_vect_type), intent(inout) :: x - integer, intent(out) :: info - - if (x%get_nrows() < n) & - & call psb_realloc(n,x%v,info) - if (info /= 0) & - & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') - - end subroutine c_base_asb - - subroutine c_base_sync(x) - implicit none - class(psb_c_base_vect_type), intent(inout) :: x - - ! - ! The base version does nothing, it's just - ! a placeholder. - ! - - end subroutine c_base_sync + + ! + ! Gather: Y = beta * Y + alpha * X(IDX(:)) + ! subroutine c_base_gthab(n,idx,alpha,x,beta,y) use psi_serial_mod @@ -499,7 +721,9 @@ contains call psi_gth(n,idx,alpha,x%v,beta,y) end subroutine c_base_gthab - + ! + ! shortcut alpha=1 beta=0 + ! subroutine c_base_gthzv(n,idx,x,y) use psi_serial_mod integer :: n, idx(:) @@ -511,6 +735,10 @@ contains end subroutine c_base_gthzv + ! + ! Scatter: + ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) + subroutine c_base_sctb(n,idx,x,beta,y) use psi_serial_mod integer :: n, idx(:) @@ -522,76 +750,4 @@ contains end subroutine c_base_sctb - subroutine c_base_free(x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - class(psb_c_base_vect_type), intent(inout) :: x - integer, 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 c_base_free - - subroutine c_base_ins(n,irl,val,dupl,x,info) - use psi_serial_mod - implicit none - class(psb_c_base_vect_type), intent(inout) :: x - integer, intent(in) :: n, dupl - integer, intent(in) :: irl(:) - complex(psb_spk_), intent(in) :: val(:) - integer, intent(out) :: info - - integer :: i - - 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 - select case(dupl) - case(psb_dupl_ovwrt_) - do i = 1, n - !loop over all val's rows - - ! row actual block row - if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') - return - end if - - end subroutine c_base_ins - end module psb_c_base_vect_mod diff --git a/base/modules/psb_d_base_vect_mod.f90 b/base/modules/psb_d_base_vect_mod.f90 index f3422e8c..46eca791 100644 --- a/base/modules/psb_d_base_vect_mod.f90 +++ b/base/modules/psb_d_base_vect_mod.f90 @@ -1,3 +1,47 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! package: psb_d_base_vect_mod +! +! This module contains the definition of the psb_d_base_vect type which +! is a container for dense vectors. +! This is encapsulated instead of being just a simple array to allow for +! more 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. +! +! module psb_d_base_vect_mod use psb_const_mod @@ -6,14 +50,58 @@ module psb_d_base_vect_mod type psb_d_base_vect_type real(psb_dpk_), allocatable :: v(:) contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => d_base_bld_x + procedure, pass(x) :: bld_n => d_base_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => d_base_all + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => d_base_ins + procedure, pass(x) :: zero => d_base_zero + procedure, pass(x) :: asb => d_base_asb + procedure, pass(x) :: free => d_base_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_sync + ! + ! Basic info procedure, pass(x) :: get_nrows => d_base_get_nrows procedure, pass(x) :: sizeof => d_base_sizeof + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: getCopy => d_base_getCopy + procedure, pass(x) :: cpy_vect => d_base_cpy_vect + procedure, pass(x) :: set_scal => d_base_set_scal + procedure, pass(x) :: set_vect => d_base_set_vect + generic, public :: set => set_vect, set_scal + generic, public :: assignment(=) => cpy_vect, set_scal + + ! + ! Dot product and AXPBY + ! procedure, pass(x) :: dot_v => d_base_dot_v procedure, pass(x) :: dot_a => d_base_dot_a generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => d_base_axpby_v procedure, pass(y) :: axpby_a => d_base_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_mlt_v procedure, pass(y) :: mlt_a => d_base_mlt_a procedure, pass(z) :: mlt_a_2 => d_base_mlt_a_2 @@ -21,30 +109,22 @@ module psb_d_base_vect_mod procedure, pass(z) :: mlt_va => d_base_mlt_va procedure, pass(z) :: mlt_av => d_base_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_scal procedure, pass(x) :: nrm2 => d_base_nrm2 procedure, pass(x) :: amax => d_base_amax procedure, pass(x) :: asum => d_base_asum - procedure, pass(x) :: all => d_base_all - procedure, pass(x) :: zero => d_base_zero - procedure, pass(x) :: asb => d_base_asb - procedure, pass(x) :: sync => d_base_sync + ! + ! Gather/scatter. These are needed for MPI interfacing. + ! May have to be reworked. + ! procedure, pass(x) :: gthab => d_base_gthab procedure, pass(x) :: gthzv => d_base_gthzv generic, public :: gth => gthab, gthzv procedure, pass(y) :: sctb => d_base_sctb generic, public :: sct => sctb - procedure, pass(x) :: free => d_base_free - procedure, pass(x) :: ins => d_base_ins - procedure, pass(x) :: bld_x => d_base_bld_x - procedure, pass(x) :: bld_n => d_base_bld_n - generic, public :: bld => bld_x, bld_n - procedure, pass(x) :: getCopy => d_base_getCopy - procedure, pass(x) :: cpy_vect => d_base_cpy_vect - generic, public :: assignment(=) => cpy_vect, set_scal - procedure, pass(x) :: set_scal => d_base_set_scal - procedure, pass(x) :: set_vect => d_base_set_vect - generic, public :: set => set_vect, set_scal end type psb_d_base_vect_type public :: psb_d_base_vect @@ -55,6 +135,33 @@ module psb_d_base_vect_mod contains + ! + ! Constructors. + ! + + function constructor(x) result(this) + real(psb_dpk_) :: x(:) + type(psb_d_base_vect_type) :: this + integer :: info + + this%v = x + call this%asb(size(x),info) + end function constructor + + + function size_const(n) result(this) + integer, intent(in) :: n + type(psb_d_base_vect_type) :: this + integer :: info + + call this%asb(n,info) + + end function size_const + + ! + ! Build from a sample + ! + subroutine d_base_bld_x(x,this) use psb_realloc_mod real(psb_dpk_), intent(in) :: this(:) @@ -70,7 +177,9 @@ contains end subroutine d_base_bld_x - + ! + ! Create with size, but no initialization + ! subroutine d_base_bld_n(x,n) use psb_realloc_mod integer, intent(in) :: n @@ -81,9 +190,172 @@ contains call x%asb(n,info) end subroutine d_base_bld_n + + subroutine d_base_all(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_d_base_vect_type), intent(out) :: x + integer, intent(out) :: info + + call psb_realloc(n,x%v,info) + + end subroutine d_base_all + + ! + ! Insert a bunch of values at specified positions. + ! + subroutine d_base_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + integer, intent(in) :: n, dupl + integer, intent(in) :: irl(:) + real(psb_dpk_), intent(in) :: val(:) + integer, intent(out) :: info + + integer :: i + + 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 + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') + return + end if + + end subroutine d_base_ins + + ! + subroutine d_base_zero(x) + use psi_serial_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=dzero + + end subroutine d_base_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + + subroutine d_base_asb(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_d_base_vect_type), intent(inout) :: x + integer, intent(out) :: info + + if (x%get_nrows() < n) & + & call psb_realloc(n,x%v,info) + if (info /= 0) & + & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') + + end subroutine d_base_asb + + + subroutine d_base_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + integer, 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_free + + + + ! + ! The base version of SYNC does nothing, it's just + ! a placeholder. + ! + subroutine d_base_sync(x) + implicit none + class(psb_d_base_vect_type), intent(inout) :: x + + + end subroutine d_base_sync + + ! + ! Size info. + ! + + function d_base_get_nrows(x) result(res) + implicit none + class(psb_d_base_vect_type), intent(in) :: x + integer :: res + + res = 0 + if (allocated(x%v)) res = size(x%v) + + end function d_base_get_nrows + + function d_base_sizeof(x) result(res) + implicit none + class(psb_d_base_vect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * psb_sizeof_dp) * x%get_nrows() + + end function d_base_sizeof + + + ! + ! Two versions of extracting an array: one of them + ! overload the assignment. + ! function d_base_getCopy(x) result(res) - class(psb_d_base_vect_type), intent(in) :: x + class(psb_d_base_vect_type), intent(in) :: x real(psb_dpk_), allocatable :: res(:) integer :: info @@ -104,6 +376,9 @@ contains end subroutine d_base_cpy_vect + ! + ! Reset all values + ! subroutine d_base_set_scal(x,val) class(psb_d_base_vect_type), intent(inout) :: x real(psb_dpk_), intent(in) :: val @@ -127,46 +402,10 @@ contains end if end subroutine d_base_set_vect - - - function constructor(x) result(this) - real(psb_dpk_) :: x(:) - type(psb_d_base_vect_type) :: this - integer :: info - - this%v = x - call this%asb(size(x),info) - end function constructor - - - function size_const(n) result(this) - integer, intent(in) :: n - type(psb_d_base_vect_type) :: this - integer :: info - - call this%asb(n,info) - - end function size_const - - function d_base_get_nrows(x) result(res) - implicit none - class(psb_d_base_vect_type), intent(in) :: x - integer :: res - - res = 0 - if (allocated(x%v)) res = size(x%v) - - end function d_base_get_nrows - - function d_base_sizeof(x) result(res) - implicit none - class(psb_d_base_vect_type), intent(in) :: x - integer(psb_long_int_k_) :: res - - res = (1_psb_long_int_k_ * psb_sizeof_dp) * x%get_nrows() - - end function d_base_sizeof + ! + ! Dot products + ! function d_base_dot_v(n,x,y) result(res) implicit none class(psb_d_base_vect_type), intent(inout) :: x, y @@ -178,7 +417,9 @@ contains ! ! Note: this is the base implementation. ! When we get here, we are sure that X is of - ! TYPE psb_d_base_vect + ! TYPE psb_d_base_vect. + ! If Y is not, throw the burden on it, implicitly + ! calling dot_a ! select type(yy => y) type is (psb_d_base_vect_type) @@ -189,6 +430,9 @@ contains end function d_base_dot_v + ! + ! Base workhorse is good old BLAS1 + ! function d_base_dot_a(n,x,y) result(res) implicit none class(psb_d_base_vect_type), intent(inout) :: x @@ -201,6 +445,9 @@ contains end function d_base_dot_a + ! + ! AXPBY is invoked via Y, hence the structure below. + ! subroutine d_base_axpby_v(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -232,7 +479,16 @@ contains end subroutine d_base_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 + ! + subroutine d_base_mlt_v(x, y, info) use psi_serial_mod implicit none @@ -400,6 +656,11 @@ contains end subroutine d_base_mlt_va + + ! + ! Simple scaling + ! + subroutine d_base_scal(alpha, x) use psi_serial_mod implicit none @@ -409,7 +670,10 @@ contains if (allocated(x%v)) x%v = alpha*x%v end subroutine d_base_scal - + + ! + ! Norms 1, 2 and infinity + ! function d_base_nrm2(n,x) result(res) implicit none @@ -442,52 +706,10 @@ contains end function d_base_asum - subroutine d_base_all(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_d_base_vect_type), intent(out) :: x - integer, intent(out) :: info - - call psb_realloc(n,x%v,info) - - end subroutine d_base_all - - subroutine d_base_zero(x) - use psi_serial_mod - implicit none - class(psb_d_base_vect_type), intent(inout) :: x - - if (allocated(x%v)) x%v=dzero - - end subroutine d_base_zero - - subroutine d_base_asb(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_d_base_vect_type), intent(inout) :: x - integer, intent(out) :: info - - if (x%get_nrows() < n) & - & call psb_realloc(n,x%v,info) - if (info /= 0) & - & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') - - end subroutine d_base_asb - - subroutine d_base_sync(x) - implicit none - class(psb_d_base_vect_type), intent(inout) :: x - - ! - ! The base version does nothing, it's just - ! a placeholder. - ! - - end subroutine d_base_sync + + ! + ! Gather: Y = beta * Y + alpha * X(IDX(:)) + ! subroutine d_base_gthab(n,idx,alpha,x,beta,y) use psi_serial_mod @@ -499,7 +721,9 @@ contains call psi_gth(n,idx,alpha,x%v,beta,y) end subroutine d_base_gthab - + ! + ! shortcut alpha=1 beta=0 + ! subroutine d_base_gthzv(n,idx,x,y) use psi_serial_mod integer :: n, idx(:) @@ -511,6 +735,10 @@ contains end subroutine d_base_gthzv + ! + ! Scatter: + ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) + subroutine d_base_sctb(n,idx,x,beta,y) use psi_serial_mod integer :: n, idx(:) @@ -522,76 +750,4 @@ contains end subroutine d_base_sctb - subroutine d_base_free(x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - class(psb_d_base_vect_type), intent(inout) :: x - integer, 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_free - - subroutine d_base_ins(n,irl,val,dupl,x,info) - use psi_serial_mod - implicit none - class(psb_d_base_vect_type), intent(inout) :: x - integer, intent(in) :: n, dupl - integer, intent(in) :: irl(:) - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - - integer :: i - - 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 - select case(dupl) - case(psb_dupl_ovwrt_) - do i = 1, n - !loop over all val's rows - - ! row actual block row - if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') - return - end if - - end subroutine d_base_ins - end module psb_d_base_vect_mod diff --git a/base/modules/psb_s_base_vect_mod.f90 b/base/modules/psb_s_base_vect_mod.f90 index d373a644..e466b518 100644 --- a/base/modules/psb_s_base_vect_mod.f90 +++ b/base/modules/psb_s_base_vect_mod.f90 @@ -1,3 +1,47 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! package: psb_s_base_vect_mod +! +! This module contains the definition of the psb_s_base_vect type which +! is a container for dense vectors. +! This is encapsulated instead of being just a simple array to allow for +! more 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. +! +! module psb_s_base_vect_mod use psb_const_mod @@ -6,14 +50,58 @@ module psb_s_base_vect_mod type psb_s_base_vect_type real(psb_spk_), allocatable :: v(:) contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => s_base_bld_x + procedure, pass(x) :: bld_n => s_base_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => s_base_all + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => s_base_ins + procedure, pass(x) :: zero => s_base_zero + procedure, pass(x) :: asb => s_base_asb + procedure, pass(x) :: free => s_base_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 => s_base_sync + ! + ! Basic info procedure, pass(x) :: get_nrows => s_base_get_nrows procedure, pass(x) :: sizeof => s_base_sizeof + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: getCopy => s_base_getCopy + procedure, pass(x) :: cpy_vect => s_base_cpy_vect + procedure, pass(x) :: set_scal => s_base_set_scal + procedure, pass(x) :: set_vect => s_base_set_vect + generic, public :: set => set_vect, set_scal + generic, public :: assignment(=) => cpy_vect, set_scal + + ! + ! Dot product and AXPBY + ! procedure, pass(x) :: dot_v => s_base_dot_v procedure, pass(x) :: dot_a => s_base_dot_a generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => s_base_axpby_v procedure, pass(y) :: axpby_a => s_base_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 => s_base_mlt_v procedure, pass(y) :: mlt_a => s_base_mlt_a procedure, pass(z) :: mlt_a_2 => s_base_mlt_a_2 @@ -21,30 +109,22 @@ module psb_s_base_vect_mod procedure, pass(z) :: mlt_va => s_base_mlt_va procedure, pass(z) :: mlt_av => s_base_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 => s_base_scal procedure, pass(x) :: nrm2 => s_base_nrm2 procedure, pass(x) :: amax => s_base_amax procedure, pass(x) :: asum => s_base_asum - procedure, pass(x) :: all => s_base_all - procedure, pass(x) :: zero => s_base_zero - procedure, pass(x) :: asb => s_base_asb - procedure, pass(x) :: sync => s_base_sync + ! + ! Gather/scatter. These are needed for MPI interfacing. + ! May have to be reworked. + ! procedure, pass(x) :: gthab => s_base_gthab procedure, pass(x) :: gthzv => s_base_gthzv generic, public :: gth => gthab, gthzv procedure, pass(y) :: sctb => s_base_sctb generic, public :: sct => sctb - procedure, pass(x) :: free => s_base_free - procedure, pass(x) :: ins => s_base_ins - procedure, pass(x) :: bld_x => s_base_bld_x - procedure, pass(x) :: bld_n => s_base_bld_n - generic, public :: bld => bld_x, bld_n - procedure, pass(x) :: getCopy => s_base_getCopy - procedure, pass(x) :: cpy_vect => s_base_cpy_vect - generic, public :: assignment(=) => cpy_vect, set_scal - procedure, pass(x) :: set_scal => s_base_set_scal - procedure, pass(x) :: set_vect => s_base_set_vect - generic, public :: set => set_vect, set_scal end type psb_s_base_vect_type public :: psb_s_base_vect @@ -55,6 +135,33 @@ module psb_s_base_vect_mod contains + ! + ! Constructors. + ! + + function constructor(x) result(this) + real(psb_spk_) :: x(:) + type(psb_s_base_vect_type) :: this + integer :: info + + this%v = x + call this%asb(size(x),info) + end function constructor + + + function size_const(n) result(this) + integer, intent(in) :: n + type(psb_s_base_vect_type) :: this + integer :: info + + call this%asb(n,info) + + end function size_const + + ! + ! Build from a sample + ! + subroutine s_base_bld_x(x,this) use psb_realloc_mod real(psb_spk_), intent(in) :: this(:) @@ -70,7 +177,9 @@ contains end subroutine s_base_bld_x - + ! + ! Create with size, but no initialization + ! subroutine s_base_bld_n(x,n) use psb_realloc_mod integer, intent(in) :: n @@ -81,9 +190,172 @@ contains call x%asb(n,info) end subroutine s_base_bld_n + + subroutine s_base_all(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_s_base_vect_type), intent(out) :: x + integer, intent(out) :: info + + call psb_realloc(n,x%v,info) + + end subroutine s_base_all + + ! + ! Insert a bunch of values at specified positions. + ! + subroutine s_base_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + integer, intent(in) :: n, dupl + integer, intent(in) :: irl(:) + real(psb_spk_), intent(in) :: val(:) + integer, intent(out) :: info + + integer :: i + + 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 + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') + return + end if + + end subroutine s_base_ins + + ! + subroutine s_base_zero(x) + use psi_serial_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=szero + + end subroutine s_base_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + + subroutine s_base_asb(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_s_base_vect_type), intent(inout) :: x + integer, intent(out) :: info + + if (x%get_nrows() < n) & + & call psb_realloc(n,x%v,info) + if (info /= 0) & + & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') + + end subroutine s_base_asb + + + subroutine s_base_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + integer, 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 s_base_free + + + + ! + ! The base version of SYNC does nothing, it's just + ! a placeholder. + ! + subroutine s_base_sync(x) + implicit none + class(psb_s_base_vect_type), intent(inout) :: x + + + end subroutine s_base_sync + + ! + ! Size info. + ! + + function s_base_get_nrows(x) result(res) + implicit none + class(psb_s_base_vect_type), intent(in) :: x + integer :: res + + res = 0 + if (allocated(x%v)) res = size(x%v) + + end function s_base_get_nrows + + function s_base_sizeof(x) result(res) + implicit none + class(psb_s_base_vect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * psb_sizeof_sp) * x%get_nrows() + + end function s_base_sizeof + + + ! + ! Two versions of extracting an array: one of them + ! overload the assignment. + ! function s_base_getCopy(x) result(res) - class(psb_s_base_vect_type), intent(in) :: x + class(psb_s_base_vect_type), intent(in) :: x real(psb_spk_), allocatable :: res(:) integer :: info @@ -104,6 +376,9 @@ contains end subroutine s_base_cpy_vect + ! + ! Reset all values + ! subroutine s_base_set_scal(x,val) class(psb_s_base_vect_type), intent(inout) :: x real(psb_spk_), intent(in) :: val @@ -127,46 +402,10 @@ contains end if end subroutine s_base_set_vect - - - function constructor(x) result(this) - real(psb_spk_) :: x(:) - type(psb_s_base_vect_type) :: this - integer :: info - - this%v = x - call this%asb(size(x),info) - end function constructor - - - function size_const(n) result(this) - integer, intent(in) :: n - type(psb_s_base_vect_type) :: this - integer :: info - - call this%asb(n,info) - - end function size_const - - function s_base_get_nrows(x) result(res) - implicit none - class(psb_s_base_vect_type), intent(in) :: x - integer :: res - - res = 0 - if (allocated(x%v)) res = size(x%v) - - end function s_base_get_nrows - - function s_base_sizeof(x) result(res) - implicit none - class(psb_s_base_vect_type), intent(in) :: x - integer(psb_long_int_k_) :: res - - res = (1_psb_long_int_k_ * psb_sizeof_sp) * x%get_nrows() - - end function s_base_sizeof + ! + ! Dot products + ! function s_base_dot_v(n,x,y) result(res) implicit none class(psb_s_base_vect_type), intent(inout) :: x, y @@ -178,7 +417,9 @@ contains ! ! Note: this is the base implementation. ! When we get here, we are sure that X is of - ! TYPE psb_s_base_vect + ! TYPE psb_s_base_vect. + ! If Y is not, throw the burden on it, implicitly + ! calling dot_a ! select type(yy => y) type is (psb_s_base_vect_type) @@ -189,6 +430,9 @@ contains end function s_base_dot_v + ! + ! Base workhorse is good old BLAS1 + ! function s_base_dot_a(n,x,y) result(res) implicit none class(psb_s_base_vect_type), intent(inout) :: x @@ -201,6 +445,9 @@ contains end function s_base_dot_a + ! + ! AXPBY is invoked via Y, hence the structure below. + ! subroutine s_base_axpby_v(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -232,7 +479,16 @@ contains end subroutine s_base_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 + ! + subroutine s_base_mlt_v(x, y, info) use psi_serial_mod implicit none @@ -400,6 +656,11 @@ contains end subroutine s_base_mlt_va + + ! + ! Simple scaling + ! + subroutine s_base_scal(alpha, x) use psi_serial_mod implicit none @@ -409,7 +670,10 @@ contains if (allocated(x%v)) x%v = alpha*x%v end subroutine s_base_scal - + + ! + ! Norms 1, 2 and infinity + ! function s_base_nrm2(n,x) result(res) implicit none @@ -442,52 +706,10 @@ contains end function s_base_asum - subroutine s_base_all(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_s_base_vect_type), intent(out) :: x - integer, intent(out) :: info - - call psb_realloc(n,x%v,info) - - end subroutine s_base_all - - subroutine s_base_zero(x) - use psi_serial_mod - implicit none - class(psb_s_base_vect_type), intent(inout) :: x - - if (allocated(x%v)) x%v=szero - - end subroutine s_base_zero - - subroutine s_base_asb(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_s_base_vect_type), intent(inout) :: x - integer, intent(out) :: info - - if (x%get_nrows() < n) & - & call psb_realloc(n,x%v,info) - if (info /= 0) & - & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') - - end subroutine s_base_asb - - subroutine s_base_sync(x) - implicit none - class(psb_s_base_vect_type), intent(inout) :: x - - ! - ! The base version does nothing, it's just - ! a placeholder. - ! - - end subroutine s_base_sync + + ! + ! Gather: Y = beta * Y + alpha * X(IDX(:)) + ! subroutine s_base_gthab(n,idx,alpha,x,beta,y) use psi_serial_mod @@ -499,7 +721,9 @@ contains call psi_gth(n,idx,alpha,x%v,beta,y) end subroutine s_base_gthab - + ! + ! shortcut alpha=1 beta=0 + ! subroutine s_base_gthzv(n,idx,x,y) use psi_serial_mod integer :: n, idx(:) @@ -511,6 +735,10 @@ contains end subroutine s_base_gthzv + ! + ! Scatter: + ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) + subroutine s_base_sctb(n,idx,x,beta,y) use psi_serial_mod integer :: n, idx(:) @@ -522,76 +750,4 @@ contains end subroutine s_base_sctb - subroutine s_base_free(x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - class(psb_s_base_vect_type), intent(inout) :: x - integer, 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 s_base_free - - subroutine s_base_ins(n,irl,val,dupl,x,info) - use psi_serial_mod - implicit none - class(psb_s_base_vect_type), intent(inout) :: x - integer, intent(in) :: n, dupl - integer, intent(in) :: irl(:) - real(psb_spk_), intent(in) :: val(:) - integer, intent(out) :: info - - integer :: i - - 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 - select case(dupl) - case(psb_dupl_ovwrt_) - do i = 1, n - !loop over all val's rows - - ! row actual block row - if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') - return - end if - - end subroutine s_base_ins - end module psb_s_base_vect_mod diff --git a/base/modules/psb_z_base_vect_mod.f90 b/base/modules/psb_z_base_vect_mod.f90 index 990184d6..966ffd11 100644 --- a/base/modules/psb_z_base_vect_mod.f90 +++ b/base/modules/psb_z_base_vect_mod.f90 @@ -1,3 +1,47 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! package: psb_z_base_vect_mod +! +! This module contains the definition of the psb_z_base_vect type which +! is a container for dense vectors. +! This is encapsulated instead of being just a simple array to allow for +! more 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. +! +! module psb_z_base_vect_mod use psb_const_mod @@ -6,14 +50,58 @@ module psb_z_base_vect_mod type psb_z_base_vect_type complex(psb_dpk_), allocatable :: v(:) contains + ! + ! Constructors/allocators + ! + procedure, pass(x) :: bld_x => z_base_bld_x + procedure, pass(x) :: bld_n => z_base_bld_n + generic, public :: bld => bld_x, bld_n + procedure, pass(x) :: all => z_base_all + ! + ! Insert/set. Assembly and free. + ! Assembly does almost nothing here, but is important + ! in derived classes. + ! + procedure, pass(x) :: ins => z_base_ins + procedure, pass(x) :: zero => z_base_zero + procedure, pass(x) :: asb => z_base_asb + procedure, pass(x) :: free => z_base_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 => z_base_sync + ! + ! Basic info procedure, pass(x) :: get_nrows => z_base_get_nrows procedure, pass(x) :: sizeof => z_base_sizeof + ! + ! Set/get data from/to an external array; also + ! overload assignment. + ! + procedure, pass(x) :: getCopy => z_base_getCopy + procedure, pass(x) :: cpy_vect => z_base_cpy_vect + procedure, pass(x) :: set_scal => z_base_set_scal + procedure, pass(x) :: set_vect => z_base_set_vect + generic, public :: set => set_vect, set_scal + generic, public :: assignment(=) => cpy_vect, set_scal + + ! + ! Dot product and AXPBY + ! procedure, pass(x) :: dot_v => z_base_dot_v procedure, pass(x) :: dot_a => z_base_dot_a generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => z_base_axpby_v procedure, pass(y) :: axpby_a => z_base_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 => z_base_mlt_v procedure, pass(y) :: mlt_a => z_base_mlt_a procedure, pass(z) :: mlt_a_2 => z_base_mlt_a_2 @@ -21,30 +109,22 @@ module psb_z_base_vect_mod procedure, pass(z) :: mlt_va => z_base_mlt_va procedure, pass(z) :: mlt_av => z_base_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 => z_base_scal procedure, pass(x) :: nrm2 => z_base_nrm2 procedure, pass(x) :: amax => z_base_amax procedure, pass(x) :: asum => z_base_asum - procedure, pass(x) :: all => z_base_all - procedure, pass(x) :: zero => z_base_zero - procedure, pass(x) :: asb => z_base_asb - procedure, pass(x) :: sync => z_base_sync + ! + ! Gather/scatter. These are needed for MPI interfacing. + ! May have to be reworked. + ! procedure, pass(x) :: gthab => z_base_gthab procedure, pass(x) :: gthzv => z_base_gthzv generic, public :: gth => gthab, gthzv procedure, pass(y) :: sctb => z_base_sctb generic, public :: sct => sctb - procedure, pass(x) :: free => z_base_free - procedure, pass(x) :: ins => z_base_ins - procedure, pass(x) :: bld_x => z_base_bld_x - procedure, pass(x) :: bld_n => z_base_bld_n - generic, public :: bld => bld_x, bld_n - procedure, pass(x) :: getCopy => z_base_getCopy - procedure, pass(x) :: cpy_vect => z_base_cpy_vect - generic, public :: assignment(=) => cpy_vect, set_scal - procedure, pass(x) :: set_scal => z_base_set_scal - procedure, pass(x) :: set_vect => z_base_set_vect - generic, public :: set => set_vect, set_scal end type psb_z_base_vect_type public :: psb_z_base_vect @@ -55,6 +135,33 @@ module psb_z_base_vect_mod contains + ! + ! Constructors. + ! + + function constructor(x) result(this) + complex(psb_dpk_) :: x(:) + type(psb_z_base_vect_type) :: this + integer :: info + + this%v = x + call this%asb(size(x),info) + end function constructor + + + function size_const(n) result(this) + integer, intent(in) :: n + type(psb_z_base_vect_type) :: this + integer :: info + + call this%asb(n,info) + + end function size_const + + ! + ! Build from a sample + ! + subroutine z_base_bld_x(x,this) use psb_realloc_mod complex(psb_dpk_), intent(in) :: this(:) @@ -70,7 +177,9 @@ contains end subroutine z_base_bld_x - + ! + ! Create with size, but no initialization + ! subroutine z_base_bld_n(x,n) use psb_realloc_mod integer, intent(in) :: n @@ -81,9 +190,172 @@ contains call x%asb(n,info) end subroutine z_base_bld_n + + subroutine z_base_all(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_z_base_vect_type), intent(out) :: x + integer, intent(out) :: info + + call psb_realloc(n,x%v,info) + + end subroutine z_base_all + + ! + ! Insert a bunch of values at specified positions. + ! + subroutine z_base_ins(n,irl,val,dupl,x,info) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + integer, intent(in) :: n, dupl + integer, intent(in) :: irl(:) + complex(psb_dpk_), intent(in) :: val(:) + integer, intent(out) :: info + + integer :: i + + 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 + select case(dupl) + case(psb_dupl_ovwrt_) + do i = 1, n + !loop over all val's rows + + ! row actual block row + if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') + return + end if + + end subroutine z_base_ins + + ! + subroutine z_base_zero(x) + use psi_serial_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + + if (allocated(x%v)) x%v=zzero + + end subroutine z_base_zero + + + ! + ! Assembly. + ! For derived classes: after this the vector + ! storage is supposed to be in sync. + ! + + subroutine z_base_asb(n, x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + integer, intent(in) :: n + class(psb_z_base_vect_type), intent(inout) :: x + integer, intent(out) :: info + + if (x%get_nrows() < n) & + & call psb_realloc(n,x%v,info) + if (info /= 0) & + & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') + + end subroutine z_base_asb + + + subroutine z_base_free(x, info) + use psi_serial_mod + use psb_realloc_mod + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + integer, 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 z_base_free + + + + ! + ! The base version of SYNC does nothing, it's just + ! a placeholder. + ! + subroutine z_base_sync(x) + implicit none + class(psb_z_base_vect_type), intent(inout) :: x + + + end subroutine z_base_sync + + ! + ! Size info. + ! + + function z_base_get_nrows(x) result(res) + implicit none + class(psb_z_base_vect_type), intent(in) :: x + integer :: res + + res = 0 + if (allocated(x%v)) res = size(x%v) + + end function z_base_get_nrows + + function z_base_sizeof(x) result(res) + implicit none + class(psb_z_base_vect_type), intent(in) :: x + integer(psb_long_int_k_) :: res + + ! Force 8-byte integers. + res = (1_psb_long_int_k_ * (2*psb_sizeof_dp)) * x%get_nrows() + + end function z_base_sizeof + + + ! + ! Two versions of extracting an array: one of them + ! overload the assignment. + ! function z_base_getCopy(x) result(res) - class(psb_z_base_vect_type), intent(in) :: x + class(psb_z_base_vect_type), intent(in) :: x complex(psb_dpk_), allocatable :: res(:) integer :: info @@ -104,6 +376,9 @@ contains end subroutine z_base_cpy_vect + ! + ! Reset all values + ! subroutine z_base_set_scal(x,val) class(psb_z_base_vect_type), intent(inout) :: x complex(psb_dpk_), intent(in) :: val @@ -127,46 +402,10 @@ contains end if end subroutine z_base_set_vect - - - function constructor(x) result(this) - complex(psb_dpk_) :: x(:) - type(psb_z_base_vect_type) :: this - integer :: info - - this%v = x - call this%asb(size(x),info) - end function constructor - - - function size_const(n) result(this) - integer, intent(in) :: n - type(psb_z_base_vect_type) :: this - integer :: info - - call this%asb(n,info) - - end function size_const - - function z_base_get_nrows(x) result(res) - implicit none - class(psb_z_base_vect_type), intent(in) :: x - integer :: res - - res = 0 - if (allocated(x%v)) res = size(x%v) - - end function z_base_get_nrows - - function z_base_sizeof(x) result(res) - implicit none - class(psb_z_base_vect_type), intent(in) :: x - integer(psb_long_int_k_) :: res - - res = (1_psb_long_int_k_ * (2*psb_sizeof_dp)) * x%get_nrows() - - end function z_base_sizeof + ! + ! Dot products + ! function z_base_dot_v(n,x,y) result(res) implicit none class(psb_z_base_vect_type), intent(inout) :: x, y @@ -178,7 +417,9 @@ contains ! ! Note: this is the base implementation. ! When we get here, we are sure that X is of - ! TYPE psb_z_base_vect + ! TYPE psb_z_base_vect. + ! If Y is not, throw the burden on it, implicitly + ! calling dot_a ! select type(yy => y) type is (psb_z_base_vect_type) @@ -189,6 +430,9 @@ contains end function z_base_dot_v + ! + ! Base workhorse is good old BLAS1 + ! function z_base_dot_a(n,x,y) result(res) implicit none class(psb_z_base_vect_type), intent(inout) :: x @@ -201,6 +445,9 @@ contains end function z_base_dot_a + ! + ! AXPBY is invoked via Y, hence the structure below. + ! subroutine z_base_axpby_v(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -232,7 +479,16 @@ contains end subroutine z_base_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 + ! + subroutine z_base_mlt_v(x, y, info) use psi_serial_mod implicit none @@ -400,6 +656,11 @@ contains end subroutine z_base_mlt_va + + ! + ! Simple scaling + ! + subroutine z_base_scal(alpha, x) use psi_serial_mod implicit none @@ -409,7 +670,10 @@ contains if (allocated(x%v)) x%v = alpha*x%v end subroutine z_base_scal - + + ! + ! Norms 1, 2 and infinity + ! function z_base_nrm2(n,x) result(res) implicit none @@ -442,52 +706,10 @@ contains end function z_base_asum - subroutine z_base_all(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_z_base_vect_type), intent(out) :: x - integer, intent(out) :: info - - call psb_realloc(n,x%v,info) - - end subroutine z_base_all - - subroutine z_base_zero(x) - use psi_serial_mod - implicit none - class(psb_z_base_vect_type), intent(inout) :: x - - if (allocated(x%v)) x%v=zzero - - end subroutine z_base_zero - - subroutine z_base_asb(n, x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - integer, intent(in) :: n - class(psb_z_base_vect_type), intent(inout) :: x - integer, intent(out) :: info - - if (x%get_nrows() < n) & - & call psb_realloc(n,x%v,info) - if (info /= 0) & - & call psb_errpush(psb_err_alloc_dealloc_,'vect_asb') - - end subroutine z_base_asb - - subroutine z_base_sync(x) - implicit none - class(psb_z_base_vect_type), intent(inout) :: x - - ! - ! The base version does nothing, it's just - ! a placeholder. - ! - - end subroutine z_base_sync + + ! + ! Gather: Y = beta * Y + alpha * X(IDX(:)) + ! subroutine z_base_gthab(n,idx,alpha,x,beta,y) use psi_serial_mod @@ -499,7 +721,9 @@ contains call psi_gth(n,idx,alpha,x%v,beta,y) end subroutine z_base_gthab - + ! + ! shortcut alpha=1 beta=0 + ! subroutine z_base_gthzv(n,idx,x,y) use psi_serial_mod integer :: n, idx(:) @@ -511,6 +735,10 @@ contains end subroutine z_base_gthzv + ! + ! Scatter: + ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:) + subroutine z_base_sctb(n,idx,x,beta,y) use psi_serial_mod integer :: n, idx(:) @@ -522,76 +750,4 @@ contains end subroutine z_base_sctb - subroutine z_base_free(x, info) - use psi_serial_mod - use psb_realloc_mod - implicit none - class(psb_z_base_vect_type), intent(inout) :: x - integer, 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 z_base_free - - subroutine z_base_ins(n,irl,val,dupl,x,info) - use psi_serial_mod - implicit none - class(psb_z_base_vect_type), intent(inout) :: x - integer, intent(in) :: n, dupl - integer, intent(in) :: irl(:) - complex(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - - integer :: i - - 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 - select case(dupl) - case(psb_dupl_ovwrt_) - do i = 1, n - !loop over all val's rows - - ! row actual block row - if (irl(i) > 0) 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 (irl(i) > 0) 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_vect_ins') - return - end if - - end subroutine z_base_ins - end module psb_z_base_vect_mod