fdurastante
/
psblas3
mirror of https://github.com/sfilippone/psblas3.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

psblas3:

Browse Source 
base/modules/psb_c_base_vect_mod.f90
 base/modules/psb_d_base_vect_mod.f90
 base/modules/psb_s_base_vect_mod.f90
 base/modules/psb_z_base_vect_mod.f90

Internal docs for base_vect.
psblas3-type-indexed
Salvatore Filippone 13 years ago
parent 5087a02aaf
commit f3f670a15e
4 changed files with 1336 additions and 712 deletions
Show Stats Download Patch File Download Diff File

500
base/modules/psb_c_base_vect_mod.f90
Unescape Escape View File

@ -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 module psb_c_base_vect_mod
use psb_const_mod use psb_const_mod
@ -6,14 +50,58 @@ module psb_c_base_vect_mod
type psb_c_base_vect_type type psb_c_base_vect_type
complex(psb_spk_), allocatable :: v(:) complex(psb_spk_), allocatable :: v(:)
contains 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) :: get_nrows => c_base_get_nrows
procedure, pass(x) :: sizeof => c_base_sizeof 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_v => c_base_dot_v
procedure, pass(x) :: dot_a => c_base_dot_a procedure, pass(x) :: dot_a => c_base_dot_a
generic, public :: dot => dot_v, dot_a generic, public :: dot => dot_v, dot_a
procedure, pass(y) :: axpby_v => c_base_axpby_v procedure, pass(y) :: axpby_v => c_base_axpby_v
procedure, pass(y) :: axpby_a => c_base_axpby_a procedure, pass(y) :: axpby_a => c_base_axpby_a
generic, public :: axpby => axpby_v, 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_v => c_base_mlt_v
procedure, pass(y) :: mlt_a => c_base_mlt_a procedure, pass(y) :: mlt_a => c_base_mlt_a
procedure, pass(z) :: mlt_a_2 => c_base_mlt_a_2 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_va => c_base_mlt_va
procedure, pass(z) :: mlt_av => c_base_mlt_av 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 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) :: scal => c_base_scal
procedure, pass(x) :: nrm2 => c_base_nrm2 procedure, pass(x) :: nrm2 => c_base_nrm2
procedure, pass(x) :: amax => c_base_amax procedure, pass(x) :: amax => c_base_amax
procedure, pass(x) :: asum => c_base_asum procedure, pass(x) :: asum => c_base_asum
procedure, pass(x) :: all => c_base_all !
procedure, pass(x) :: zero => c_base_zero ! Gather/scatter. These are needed for MPI interfacing.
procedure, pass(x) :: asb => c_base_asb ! May have to be reworked.
procedure, pass(x) :: sync => c_base_sync !
procedure, pass(x) :: gthab => c_base_gthab procedure, pass(x) :: gthab => c_base_gthab
procedure, pass(x) :: gthzv => c_base_gthzv procedure, pass(x) :: gthzv => c_base_gthzv
generic, public :: gth => gthab, gthzv generic, public :: gth => gthab, gthzv
procedure, pass(y) :: sctb => c_base_sctb procedure, pass(y) :: sctb => c_base_sctb
generic, public :: sct => 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 end type psb_c_base_vect_type
public :: psb_c_base_vect public :: psb_c_base_vect
@ -55,6 +135,33 @@ module psb_c_base_vect_mod
contains 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) subroutine c_base_bld_x(x,this)
use psb_realloc_mod use psb_realloc_mod
complex(psb_spk_), intent(in) :: this(:) complex(psb_spk_), intent(in) :: this(:)
@ -70,7 +177,9 @@ contains
end subroutine c_base_bld_x end subroutine c_base_bld_x
!
! Create with size, but no initialization
!
subroutine c_base_bld_n(x,n) subroutine c_base_bld_n(x,n)
use psb_realloc_mod use psb_realloc_mod
integer, intent(in) :: n integer, intent(in) :: n
@ -82,6 +191,169 @@ contains
end subroutine c_base_bld_n 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) 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(:) complex(psb_spk_), allocatable :: res(:)
@ -104,6 +376,9 @@ contains
end subroutine c_base_cpy_vect end subroutine c_base_cpy_vect
!
! Reset all values
!
subroutine c_base_set_scal(x,val) subroutine c_base_set_scal(x,val)
class(psb_c_base_vect_type), intent(inout) :: x class(psb_c_base_vect_type), intent(inout) :: x
complex(psb_spk_), intent(in) :: val complex(psb_spk_), intent(in) :: val
@ -128,45 +403,9 @@ contains
end subroutine c_base_set_vect end subroutine c_base_set_vect
!
function constructor(x) result(this) ! Dot products
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
function c_base_dot_v(n,x,y) result(res) function c_base_dot_v(n,x,y) result(res)
implicit none implicit none
class(psb_c_base_vect_type), intent(inout) :: x, y class(psb_c_base_vect_type), intent(inout) :: x, y
@ -178,7 +417,9 @@ contains
! !
! Note: this is the base implementation. ! Note: this is the base implementation.
! When we get here, we are sure that X is of ! 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) select type(yy => y)
type is (psb_c_base_vect_type) type is (psb_c_base_vect_type)
@ -189,6 +430,9 @@ contains
end function c_base_dot_v end function c_base_dot_v
!
! Base workhorse is good old BLAS1
!
function c_base_dot_a(n,x,y) result(res) function c_base_dot_a(n,x,y) result(res)
implicit none implicit none
class(psb_c_base_vect_type), intent(inout) :: x class(psb_c_base_vect_type), intent(inout) :: x
@ -201,6 +445,9 @@ contains
end function c_base_dot_a 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) subroutine c_base_axpby_v(m,alpha, x, beta, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -233,6 +480,15 @@ contains
end subroutine c_base_axpby_a 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) subroutine c_base_mlt_v(x, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -400,6 +656,11 @@ contains
end subroutine c_base_mlt_va end subroutine c_base_mlt_va
!
! Simple scaling
!
subroutine c_base_scal(alpha, x) subroutine c_base_scal(alpha, x)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -410,6 +671,9 @@ contains
end subroutine c_base_scal end subroutine c_base_scal
!
! Norms 1, 2 and infinity
!
function c_base_nrm2(n,x) result(res) function c_base_nrm2(n,x) result(res)
implicit none implicit none
@ -442,53 +706,11 @@ contains
end function c_base_asum 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 ! Gather: Y = beta * Y + alpha * X(IDX(:))
! a placeholder.
! !
end subroutine c_base_sync
subroutine c_base_gthab(n,idx,alpha,x,beta,y) subroutine c_base_gthab(n,idx,alpha,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -499,7 +721,9 @@ contains
call psi_gth(n,idx,alpha,x%v,beta,y) call psi_gth(n,idx,alpha,x%v,beta,y)
end subroutine c_base_gthab end subroutine c_base_gthab
!
! shortcut alpha=1 beta=0
!
subroutine c_base_gthzv(n,idx,x,y) subroutine c_base_gthzv(n,idx,x,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -511,6 +735,10 @@ contains
end subroutine c_base_gthzv end subroutine c_base_gthzv
!
! Scatter:
! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
subroutine c_base_sctb(n,idx,x,beta,y) subroutine c_base_sctb(n,idx,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -522,76 +750,4 @@ contains
end subroutine c_base_sctb 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 end module psb_c_base_vect_mod

500
base/modules/psb_d_base_vect_mod.f90
Unescape Escape View File

@ -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 module psb_d_base_vect_mod
use psb_const_mod use psb_const_mod
@ -6,14 +50,58 @@ module psb_d_base_vect_mod
type psb_d_base_vect_type type psb_d_base_vect_type
real(psb_dpk_), allocatable :: v(:) real(psb_dpk_), allocatable :: v(:)
contains 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) :: get_nrows => d_base_get_nrows
procedure, pass(x) :: sizeof => d_base_sizeof 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_v => d_base_dot_v
procedure, pass(x) :: dot_a => d_base_dot_a procedure, pass(x) :: dot_a => d_base_dot_a
generic, public :: dot => dot_v, dot_a generic, public :: dot => dot_v, dot_a
procedure, pass(y) :: axpby_v => d_base_axpby_v procedure, pass(y) :: axpby_v => d_base_axpby_v
procedure, pass(y) :: axpby_a => d_base_axpby_a procedure, pass(y) :: axpby_a => d_base_axpby_a
generic, public :: axpby => axpby_v, 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_v => d_base_mlt_v
procedure, pass(y) :: mlt_a => d_base_mlt_a procedure, pass(y) :: mlt_a => d_base_mlt_a
procedure, pass(z) :: mlt_a_2 => d_base_mlt_a_2 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_va => d_base_mlt_va
procedure, pass(z) :: mlt_av => d_base_mlt_av 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 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) :: scal => d_base_scal
procedure, pass(x) :: nrm2 => d_base_nrm2 procedure, pass(x) :: nrm2 => d_base_nrm2
procedure, pass(x) :: amax => d_base_amax procedure, pass(x) :: amax => d_base_amax
procedure, pass(x) :: asum => d_base_asum procedure, pass(x) :: asum => d_base_asum
procedure, pass(x) :: all => d_base_all !
procedure, pass(x) :: zero => d_base_zero ! Gather/scatter. These are needed for MPI interfacing.
procedure, pass(x) :: asb => d_base_asb ! May have to be reworked.
procedure, pass(x) :: sync => d_base_sync !
procedure, pass(x) :: gthab => d_base_gthab procedure, pass(x) :: gthab => d_base_gthab
procedure, pass(x) :: gthzv => d_base_gthzv procedure, pass(x) :: gthzv => d_base_gthzv
generic, public :: gth => gthab, gthzv generic, public :: gth => gthab, gthzv
procedure, pass(y) :: sctb => d_base_sctb procedure, pass(y) :: sctb => d_base_sctb
generic, public :: sct => 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 end type psb_d_base_vect_type
public :: psb_d_base_vect public :: psb_d_base_vect
@ -55,6 +135,33 @@ module psb_d_base_vect_mod
contains 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) subroutine d_base_bld_x(x,this)
use psb_realloc_mod use psb_realloc_mod
real(psb_dpk_), intent(in) :: this(:) real(psb_dpk_), intent(in) :: this(:)
@ -70,7 +177,9 @@ contains
end subroutine d_base_bld_x end subroutine d_base_bld_x
!
! Create with size, but no initialization
!
subroutine d_base_bld_n(x,n) subroutine d_base_bld_n(x,n)
use psb_realloc_mod use psb_realloc_mod
integer, intent(in) :: n integer, intent(in) :: n
@ -82,6 +191,169 @@ contains
end subroutine d_base_bld_n 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) 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(:) real(psb_dpk_), allocatable :: res(:)
@ -104,6 +376,9 @@ contains
end subroutine d_base_cpy_vect end subroutine d_base_cpy_vect
!
! Reset all values
!
subroutine d_base_set_scal(x,val) subroutine d_base_set_scal(x,val)
class(psb_d_base_vect_type), intent(inout) :: x class(psb_d_base_vect_type), intent(inout) :: x
real(psb_dpk_), intent(in) :: val real(psb_dpk_), intent(in) :: val
@ -128,45 +403,9 @@ contains
end subroutine d_base_set_vect end subroutine d_base_set_vect
!
function constructor(x) result(this) ! Dot products
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
function d_base_dot_v(n,x,y) result(res) function d_base_dot_v(n,x,y) result(res)
implicit none implicit none
class(psb_d_base_vect_type), intent(inout) :: x, y class(psb_d_base_vect_type), intent(inout) :: x, y
@ -178,7 +417,9 @@ contains
! !
! Note: this is the base implementation. ! Note: this is the base implementation.
! When we get here, we are sure that X is of ! 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) select type(yy => y)
type is (psb_d_base_vect_type) type is (psb_d_base_vect_type)
@ -189,6 +430,9 @@ contains
end function d_base_dot_v end function d_base_dot_v
!
! Base workhorse is good old BLAS1
!
function d_base_dot_a(n,x,y) result(res) function d_base_dot_a(n,x,y) result(res)
implicit none implicit none
class(psb_d_base_vect_type), intent(inout) :: x class(psb_d_base_vect_type), intent(inout) :: x
@ -201,6 +445,9 @@ contains
end function d_base_dot_a 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) subroutine d_base_axpby_v(m,alpha, x, beta, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -233,6 +480,15 @@ contains
end subroutine d_base_axpby_a 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) subroutine d_base_mlt_v(x, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -400,6 +656,11 @@ contains
end subroutine d_base_mlt_va end subroutine d_base_mlt_va
!
! Simple scaling
!
subroutine d_base_scal(alpha, x) subroutine d_base_scal(alpha, x)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -410,6 +671,9 @@ contains
end subroutine d_base_scal end subroutine d_base_scal
!
! Norms 1, 2 and infinity
!
function d_base_nrm2(n,x) result(res) function d_base_nrm2(n,x) result(res)
implicit none implicit none
@ -442,53 +706,11 @@ contains
end function d_base_asum 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 ! Gather: Y = beta * Y + alpha * X(IDX(:))
! a placeholder.
! !
end subroutine d_base_sync
subroutine d_base_gthab(n,idx,alpha,x,beta,y) subroutine d_base_gthab(n,idx,alpha,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -499,7 +721,9 @@ contains
call psi_gth(n,idx,alpha,x%v,beta,y) call psi_gth(n,idx,alpha,x%v,beta,y)
end subroutine d_base_gthab end subroutine d_base_gthab
!
! shortcut alpha=1 beta=0
!
subroutine d_base_gthzv(n,idx,x,y) subroutine d_base_gthzv(n,idx,x,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -511,6 +735,10 @@ contains
end subroutine d_base_gthzv end subroutine d_base_gthzv
!
! Scatter:
! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
subroutine d_base_sctb(n,idx,x,beta,y) subroutine d_base_sctb(n,idx,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -522,76 +750,4 @@ contains
end subroutine d_base_sctb 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 end module psb_d_base_vect_mod

500
base/modules/psb_s_base_vect_mod.f90
Unescape Escape View File

@ -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 module psb_s_base_vect_mod
use psb_const_mod use psb_const_mod
@ -6,14 +50,58 @@ module psb_s_base_vect_mod
type psb_s_base_vect_type type psb_s_base_vect_type
real(psb_spk_), allocatable :: v(:) real(psb_spk_), allocatable :: v(:)
contains 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) :: get_nrows => s_base_get_nrows
procedure, pass(x) :: sizeof => s_base_sizeof 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_v => s_base_dot_v
procedure, pass(x) :: dot_a => s_base_dot_a procedure, pass(x) :: dot_a => s_base_dot_a
generic, public :: dot => dot_v, dot_a generic, public :: dot => dot_v, dot_a
procedure, pass(y) :: axpby_v => s_base_axpby_v procedure, pass(y) :: axpby_v => s_base_axpby_v
procedure, pass(y) :: axpby_a => s_base_axpby_a procedure, pass(y) :: axpby_a => s_base_axpby_a
generic, public :: axpby => axpby_v, 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_v => s_base_mlt_v
procedure, pass(y) :: mlt_a => s_base_mlt_a procedure, pass(y) :: mlt_a => s_base_mlt_a
procedure, pass(z) :: mlt_a_2 => s_base_mlt_a_2 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_va => s_base_mlt_va
procedure, pass(z) :: mlt_av => s_base_mlt_av 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 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) :: scal => s_base_scal
procedure, pass(x) :: nrm2 => s_base_nrm2 procedure, pass(x) :: nrm2 => s_base_nrm2
procedure, pass(x) :: amax => s_base_amax procedure, pass(x) :: amax => s_base_amax
procedure, pass(x) :: asum => s_base_asum procedure, pass(x) :: asum => s_base_asum
procedure, pass(x) :: all => s_base_all !
procedure, pass(x) :: zero => s_base_zero ! Gather/scatter. These are needed for MPI interfacing.
procedure, pass(x) :: asb => s_base_asb ! May have to be reworked.
procedure, pass(x) :: sync => s_base_sync !
procedure, pass(x) :: gthab => s_base_gthab procedure, pass(x) :: gthab => s_base_gthab
procedure, pass(x) :: gthzv => s_base_gthzv procedure, pass(x) :: gthzv => s_base_gthzv
generic, public :: gth => gthab, gthzv generic, public :: gth => gthab, gthzv
procedure, pass(y) :: sctb => s_base_sctb procedure, pass(y) :: sctb => s_base_sctb
generic, public :: sct => 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 end type psb_s_base_vect_type
public :: psb_s_base_vect public :: psb_s_base_vect
@ -55,6 +135,33 @@ module psb_s_base_vect_mod
contains 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) subroutine s_base_bld_x(x,this)
use psb_realloc_mod use psb_realloc_mod
real(psb_spk_), intent(in) :: this(:) real(psb_spk_), intent(in) :: this(:)
@ -70,7 +177,9 @@ contains
end subroutine s_base_bld_x end subroutine s_base_bld_x
!
! Create with size, but no initialization
!
subroutine s_base_bld_n(x,n) subroutine s_base_bld_n(x,n)
use psb_realloc_mod use psb_realloc_mod
integer, intent(in) :: n integer, intent(in) :: n
@ -82,6 +191,169 @@ contains
end subroutine s_base_bld_n 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) 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(:) real(psb_spk_), allocatable :: res(:)
@ -104,6 +376,9 @@ contains
end subroutine s_base_cpy_vect end subroutine s_base_cpy_vect
!
! Reset all values
!
subroutine s_base_set_scal(x,val) subroutine s_base_set_scal(x,val)
class(psb_s_base_vect_type), intent(inout) :: x class(psb_s_base_vect_type), intent(inout) :: x
real(psb_spk_), intent(in) :: val real(psb_spk_), intent(in) :: val
@ -128,45 +403,9 @@ contains
end subroutine s_base_set_vect end subroutine s_base_set_vect
!
function constructor(x) result(this) ! Dot products
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
function s_base_dot_v(n,x,y) result(res) function s_base_dot_v(n,x,y) result(res)
implicit none implicit none
class(psb_s_base_vect_type), intent(inout) :: x, y class(psb_s_base_vect_type), intent(inout) :: x, y
@ -178,7 +417,9 @@ contains
! !
! Note: this is the base implementation. ! Note: this is the base implementation.
! When we get here, we are sure that X is of ! 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) select type(yy => y)
type is (psb_s_base_vect_type) type is (psb_s_base_vect_type)
@ -189,6 +430,9 @@ contains
end function s_base_dot_v end function s_base_dot_v
!
! Base workhorse is good old BLAS1
!
function s_base_dot_a(n,x,y) result(res) function s_base_dot_a(n,x,y) result(res)
implicit none implicit none
class(psb_s_base_vect_type), intent(inout) :: x class(psb_s_base_vect_type), intent(inout) :: x
@ -201,6 +445,9 @@ contains
end function s_base_dot_a 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) subroutine s_base_axpby_v(m,alpha, x, beta, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -233,6 +480,15 @@ contains
end subroutine s_base_axpby_a 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) subroutine s_base_mlt_v(x, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -400,6 +656,11 @@ contains
end subroutine s_base_mlt_va end subroutine s_base_mlt_va
!
! Simple scaling
!
subroutine s_base_scal(alpha, x) subroutine s_base_scal(alpha, x)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -410,6 +671,9 @@ contains
end subroutine s_base_scal end subroutine s_base_scal
!
! Norms 1, 2 and infinity
!
function s_base_nrm2(n,x) result(res) function s_base_nrm2(n,x) result(res)
implicit none implicit none
@ -442,53 +706,11 @@ contains
end function s_base_asum 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 ! Gather: Y = beta * Y + alpha * X(IDX(:))
! a placeholder.
! !
end subroutine s_base_sync
subroutine s_base_gthab(n,idx,alpha,x,beta,y) subroutine s_base_gthab(n,idx,alpha,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -499,7 +721,9 @@ contains
call psi_gth(n,idx,alpha,x%v,beta,y) call psi_gth(n,idx,alpha,x%v,beta,y)
end subroutine s_base_gthab end subroutine s_base_gthab
!
! shortcut alpha=1 beta=0
!
subroutine s_base_gthzv(n,idx,x,y) subroutine s_base_gthzv(n,idx,x,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -511,6 +735,10 @@ contains
end subroutine s_base_gthzv end subroutine s_base_gthzv
!
! Scatter:
! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
subroutine s_base_sctb(n,idx,x,beta,y) subroutine s_base_sctb(n,idx,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -522,76 +750,4 @@ contains
end subroutine s_base_sctb 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 end module psb_s_base_vect_mod

500
base/modules/psb_z_base_vect_mod.f90
Unescape Escape View File

@ -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 module psb_z_base_vect_mod
use psb_const_mod use psb_const_mod
@ -6,14 +50,58 @@ module psb_z_base_vect_mod
type psb_z_base_vect_type type psb_z_base_vect_type
complex(psb_dpk_), allocatable :: v(:) complex(psb_dpk_), allocatable :: v(:)
contains 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) :: get_nrows => z_base_get_nrows
procedure, pass(x) :: sizeof => z_base_sizeof 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_v => z_base_dot_v
procedure, pass(x) :: dot_a => z_base_dot_a procedure, pass(x) :: dot_a => z_base_dot_a
generic, public :: dot => dot_v, dot_a generic, public :: dot => dot_v, dot_a
procedure, pass(y) :: axpby_v => z_base_axpby_v procedure, pass(y) :: axpby_v => z_base_axpby_v
procedure, pass(y) :: axpby_a => z_base_axpby_a procedure, pass(y) :: axpby_a => z_base_axpby_a
generic, public :: axpby => axpby_v, 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_v => z_base_mlt_v
procedure, pass(y) :: mlt_a => z_base_mlt_a procedure, pass(y) :: mlt_a => z_base_mlt_a
procedure, pass(z) :: mlt_a_2 => z_base_mlt_a_2 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_va => z_base_mlt_va
procedure, pass(z) :: mlt_av => z_base_mlt_av 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 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) :: scal => z_base_scal
procedure, pass(x) :: nrm2 => z_base_nrm2 procedure, pass(x) :: nrm2 => z_base_nrm2
procedure, pass(x) :: amax => z_base_amax procedure, pass(x) :: amax => z_base_amax
procedure, pass(x) :: asum => z_base_asum procedure, pass(x) :: asum => z_base_asum
procedure, pass(x) :: all => z_base_all !
procedure, pass(x) :: zero => z_base_zero ! Gather/scatter. These are needed for MPI interfacing.
procedure, pass(x) :: asb => z_base_asb ! May have to be reworked.
procedure, pass(x) :: sync => z_base_sync !
procedure, pass(x) :: gthab => z_base_gthab procedure, pass(x) :: gthab => z_base_gthab
procedure, pass(x) :: gthzv => z_base_gthzv procedure, pass(x) :: gthzv => z_base_gthzv
generic, public :: gth => gthab, gthzv generic, public :: gth => gthab, gthzv
procedure, pass(y) :: sctb => z_base_sctb procedure, pass(y) :: sctb => z_base_sctb
generic, public :: sct => 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 end type psb_z_base_vect_type
public :: psb_z_base_vect public :: psb_z_base_vect
@ -55,6 +135,33 @@ module psb_z_base_vect_mod
contains 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) subroutine z_base_bld_x(x,this)
use psb_realloc_mod use psb_realloc_mod
complex(psb_dpk_), intent(in) :: this(:) complex(psb_dpk_), intent(in) :: this(:)
@ -70,7 +177,9 @@ contains
end subroutine z_base_bld_x end subroutine z_base_bld_x
!
! Create with size, but no initialization
!
subroutine z_base_bld_n(x,n) subroutine z_base_bld_n(x,n)
use psb_realloc_mod use psb_realloc_mod
integer, intent(in) :: n integer, intent(in) :: n
@ -82,6 +191,169 @@ contains
end subroutine z_base_bld_n 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) 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(:) complex(psb_dpk_), allocatable :: res(:)
@ -104,6 +376,9 @@ contains
end subroutine z_base_cpy_vect end subroutine z_base_cpy_vect
!
! Reset all values
!
subroutine z_base_set_scal(x,val) subroutine z_base_set_scal(x,val)
class(psb_z_base_vect_type), intent(inout) :: x class(psb_z_base_vect_type), intent(inout) :: x
complex(psb_dpk_), intent(in) :: val complex(psb_dpk_), intent(in) :: val
@ -128,45 +403,9 @@ contains
end subroutine z_base_set_vect end subroutine z_base_set_vect
!
function constructor(x) result(this) ! Dot products
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
function z_base_dot_v(n,x,y) result(res) function z_base_dot_v(n,x,y) result(res)
implicit none implicit none
class(psb_z_base_vect_type), intent(inout) :: x, y class(psb_z_base_vect_type), intent(inout) :: x, y
@ -178,7 +417,9 @@ contains
! !
! Note: this is the base implementation. ! Note: this is the base implementation.
! When we get here, we are sure that X is of ! 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) select type(yy => y)
type is (psb_z_base_vect_type) type is (psb_z_base_vect_type)
@ -189,6 +430,9 @@ contains
end function z_base_dot_v end function z_base_dot_v
!
! Base workhorse is good old BLAS1
!
function z_base_dot_a(n,x,y) result(res) function z_base_dot_a(n,x,y) result(res)
implicit none implicit none
class(psb_z_base_vect_type), intent(inout) :: x class(psb_z_base_vect_type), intent(inout) :: x
@ -201,6 +445,9 @@ contains
end function z_base_dot_a 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) subroutine z_base_axpby_v(m,alpha, x, beta, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -233,6 +480,15 @@ contains
end subroutine z_base_axpby_a 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) subroutine z_base_mlt_v(x, y, info)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -400,6 +656,11 @@ contains
end subroutine z_base_mlt_va end subroutine z_base_mlt_va
!
! Simple scaling
!
subroutine z_base_scal(alpha, x) subroutine z_base_scal(alpha, x)
use psi_serial_mod use psi_serial_mod
implicit none implicit none
@ -410,6 +671,9 @@ contains
end subroutine z_base_scal end subroutine z_base_scal
!
! Norms 1, 2 and infinity
!
function z_base_nrm2(n,x) result(res) function z_base_nrm2(n,x) result(res)
implicit none implicit none
@ -442,53 +706,11 @@ contains
end function z_base_asum 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 ! Gather: Y = beta * Y + alpha * X(IDX(:))
! a placeholder.
! !
end subroutine z_base_sync
subroutine z_base_gthab(n,idx,alpha,x,beta,y) subroutine z_base_gthab(n,idx,alpha,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -499,7 +721,9 @@ contains
call psi_gth(n,idx,alpha,x%v,beta,y) call psi_gth(n,idx,alpha,x%v,beta,y)
end subroutine z_base_gthab end subroutine z_base_gthab
!
! shortcut alpha=1 beta=0
!
subroutine z_base_gthzv(n,idx,x,y) subroutine z_base_gthzv(n,idx,x,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -511,6 +735,10 @@ contains
end subroutine z_base_gthzv end subroutine z_base_gthzv
!
! Scatter:
! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
subroutine z_base_sctb(n,idx,x,beta,y) subroutine z_base_sctb(n,idx,x,beta,y)
use psi_serial_mod use psi_serial_mod
integer :: n, idx(:) integer :: n, idx(:)
@ -522,76 +750,4 @@ contains
end subroutine z_base_sctb 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 end module psb_z_base_vect_mod

Write Preview
Loading…
Cancel
Save