fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-17 03:12:15 +00:00
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_i_base_vect_mod.f90
 base/modules/psb_s_base_vect_mod.f90
 base/modules/psb_z_base_vect_mod.f90

Additional methods for multivectors.
This commit is contained in:
Salvatore Filippone
2015-07-07 17:59:18 +00:00
parent 585cf6f735
commit 1917e2212a
5 changed files with 1092 additions and 1032 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/modules/psb_c_base_vect_mod.f90
+268 -253
View File
@@ -925,7 +925,6 @@ contains
if (z%is_dev()) call z%sync()
n = min(size(z%v), size(x), size(y))
!!$ write(0,*) 'Mlt_a_2: ',n
if (alpha == czero) then
if (beta == cone) then
return
@@ -1192,26 +1191,6 @@ contains
call x%gth(n,idx%v(i:),x%combuf(i:))
end subroutine c_base_gthzbuf
subroutine c_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex(psb_spk_) :: beta
class(psb_c_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine c_base_sctb_buf
!
!> Function base_device_wait:
!! \memberof psb_c_base_vect_type
@@ -1329,6 +1308,26 @@ contains
end subroutine c_base_sctb_x
subroutine c_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex(psb_spk_) :: beta
class(psb_c_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine c_base_sctb_buf
end module psb_c_base_vect_mod
@@ -1415,20 +1414,23 @@ module psb_c_base_multivect_mod
procedure, pass(y) :: axpby_v => c_base_mlv_axpby_v
procedure, pass(y) :: axpby_a => c_base_mlv_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_mlv_mlt_v
!!$ procedure, pass(y) :: mlt_a => c_base_mlv_mlt_a
!!$ procedure, pass(z) :: mlt_a_2 => c_base_mlv_mlt_a_2
!!$ procedure, pass(z) :: mlt_v_2 => c_base_mlv_mlt_v_2
!
! MultiVector by vector/multivector multiplication. Need all variants
! to handle multiple requirements from preconditioners
!
procedure, pass(y) :: mlt_mv => c_base_mlv_mlt_mv
procedure, pass(y) :: mlt_mv_v => c_base_mlv_mlt_mv_v
procedure, pass(y) :: mlt_ar1 => c_base_mlv_mlt_ar1
procedure, pass(y) :: mlt_ar2 => c_base_mlv_mlt_ar2
procedure, pass(z) :: mlt_a_2 => c_base_mlv_mlt_a_2
procedure, pass(z) :: mlt_v_2 => c_base_mlv_mlt_v_2
!!$ procedure, pass(z) :: mlt_va => c_base_mlv_mlt_va
!!$ procedure, pass(z) :: mlt_av => c_base_mlv_mlt_av
!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
!!$ !
!!$ ! Scaling and norms
!!$ !
generic, public :: mlt => mlt_mv, mlt_mv_v, mlt_ar1, mlt_ar2, &
& mlt_a_2, mlt_v_2 !, mlt_av, mlt_va
!
! Scaling and norms
!
procedure, pass(x) :: scal => c_base_mlv_scal
procedure, pass(x) :: nrm2 => c_base_mlv_nrm2
procedure, pass(x) :: amax => c_base_mlv_amax
@@ -1444,9 +1446,9 @@ module psb_c_base_multivect_mod
procedure, pass(x) :: gthzv => c_base_mlv_gthzv
procedure, pass(x) :: gthzv_x => c_base_mlv_gthzv_x
generic, public :: gth => gthab, gthzv, gthzv_x
!!$ procedure, pass(y) :: sctb => c_base_mlv_sctb
!!$ procedure, pass(y) :: sctb_x => c_base_mlv_sctb_x
!!$ generic, public :: sct => sctb, sctb_x
procedure, pass(y) :: sctb => c_base_mlv_sctb
procedure, pass(y) :: sctb_x => c_base_mlv_sctb_x
generic, public :: sct => sctb, sctb_x
end type psb_c_base_multivect_type
interface psb_c_base_multivect
@@ -2082,183 +2084,194 @@ contains
end subroutine c_base_mlv_axpby_a
!!$ !
!!$ ! Multiple variants of two operations:
!!$ ! Simple multiplication Y(:) = X(:)*Y(:)
!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!!$ !
!!$ ! Variants expanded according to the dynamic type
!!$ ! of the involved entities
!!$ !
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_c_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a base_mlv_vect array y=x*y
!!$ !! \param x The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine c_base_mlv_mlt_v(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ class(psb_c_base_multivect_type), intent(inout) :: x
!!$ class(psb_c_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ select type(xx => x)
!!$ type is (psb_c_base_multivect_type)
!!$ n = min(size(y%v), size(xx%v))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*xx%v(i)
!!$ end do
!!$ class default
!!$ call y%mlt(x%v,info)
!!$ end select
!!$
!!$ end subroutine c_base_mlv_mlt_v
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_c_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y
!!$ !! \param x(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine c_base_mlv_mlt_a(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ complex(psb_spk_), intent(in) :: x(:)
!!$ class(psb_c_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(y%v), size(x))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*x(i)
!!$ end do
!!$
!!$ end subroutine c_base_mlv_mlt_a
!!$
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a_2
!!$ !! \memberof psb_c_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x(:) The array to be multiplied b
!!$ !! \param y(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine c_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ complex(psb_spk_), intent(in) :: alpha,beta
!!$ complex(psb_spk_), intent(in) :: y(:)
!!$ complex(psb_spk_), intent(in) :: x(:)
!!$ class(psb_c_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(z%v), size(x), size(y))
!!$ if (alpha == izero) then
!!$ if (beta == ione) then
!!$ return
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (alpha == ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + y(i)*x(i)
!!$ end do
!!$ end if
!!$ else if (alpha == -ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = -y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) - y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) - y(i)*x(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = alpha*y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ end if
!!$ end if
!!$ end if
!!$ end subroutine c_base_mlv_mlt_a_2
!!$
!!$ !
!!$ !> Function base_mlv_mlt_v_2
!!$ !! \memberof psb_c_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x The class(base_mlv_vect) to be multiplied b
!!$ !! \param y The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine c_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
!!$ use psi_serial_mod
!!$ use psb_string_mod
!!$ implicit none
!!$ complex(psb_spk_), intent(in) :: alpha,beta
!!$ class(psb_c_base_multivect_type), intent(inout) :: x
!!$ class(psb_c_base_multivect_type), intent(inout) :: y
!!$ class(psb_c_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ character(len=1), intent(in), optional :: conjgx, conjgy
!!$ integer(psb_ipk_) :: i, n
!!$ logical :: conjgx_, conjgy_
!!$
!!$ info = 0
!!$ if (.not.psb_i_is_complex_) then
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ else
!!$ conjgx_=.false.
!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
!!$ conjgy_=.false.
!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ end if
!!$ end subroutine c_base_mlv_mlt_v_2
!
! Multiple variants of two operations:
! Simple multiplication Y(:.:) = X(:,:)*Y(:,:)
! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!
! Variants expanded according to the dynamic type
! of the involved entities
!
!
!> Function base_mlv_mlt_mv
!! \memberof psb_c_base_multivect_type
!! \brief Multivector entry-by-entry multiply by a base_mlv_multivect y=x*y
!! \param x The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine c_base_mlv_mlt_mv(x, y, info)
use psi_serial_mod
implicit none
class(psb_c_base_multivect_type), intent(inout) :: x
class(psb_c_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine c_base_mlv_mlt_mv
subroutine c_base_mlv_mlt_mv_v(x, y, info)
use psi_serial_mod
implicit none
class(psb_c_base_vect_type), intent(inout) :: x
class(psb_c_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine c_base_mlv_mlt_mv_v
!
!> Function base_mlv_mlt_ar1
!! \memberof psb_c_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 1 array y=x*y
!! \param x(:) The array to be multiplied by
!! \param info return code
!!
subroutine c_base_mlv_mlt_ar1(x, y, info)
use psi_serial_mod
implicit none
complex(psb_spk_), intent(in) :: x(:)
class(psb_c_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, n
info = 0
n = min(psb_size(y%v,1), size(x))
do i=1, n
y%v(i,:) = y%v(i,:)*x(i)
end do
end subroutine c_base_mlv_mlt_ar1
!
!> Function base_mlv_mlt_ar2
!! \memberof psb_c_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 2 array y=x*y
!! \param x(:,:) The array to be multiplied by
!! \param info return code
!!
subroutine c_base_mlv_mlt_ar2(x, y, info)
use psi_serial_mod
implicit none
complex(psb_spk_), intent(in) :: x(:,:)
class(psb_c_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr,nc
info = 0
nr = min(psb_size(y%v,1), size(x,1))
nc = min(psb_size(y%v,2), size(x,2))
y%v(1:nr,1:nc) = y%v(1:nr,1:nc)*x(1:nr,1:nc)
end subroutine c_base_mlv_mlt_ar2
!
!> Function base_mlv_mlt_a_2
!! \memberof psb_c_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x(:) The array to be multiplied b
!! \param y(:) The array to be multiplied by
!! \param info return code
!!
subroutine c_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
use psi_serial_mod
implicit none
complex(psb_spk_), intent(in) :: alpha,beta
complex(psb_spk_), intent(in) :: y(:,:)
complex(psb_spk_), intent(in) :: x(:,:)
class(psb_c_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr, nc
info = 0
nr = min(psb_size(z%v,1), size(x,1), size(y,1))
nc = min(psb_size(z%v,2), size(x,2), size(y,2))
if (alpha == czero) then
if (beta == cone) then
return
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc)
end if
else
if (alpha == cone) then
if (beta == czero) then
z%v(1:nr,1:nc) = y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == cone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else if (alpha == -cone) then
if (beta == czero) then
z%v(1:nr,1:nc) = -y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == cone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else
if (beta == czero) then
z%v(1:nr,1:nc) = alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == cone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
end if
end if
end if
end subroutine c_base_mlv_mlt_a_2
!
!> Function base_mlv_mlt_v_2
!! \memberof psb_c_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x The class(base_mlv_vect) to be multiplied b
!! \param y The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine c_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
use psi_serial_mod
use psb_string_mod
implicit none
complex(psb_spk_), intent(in) :: alpha,beta
class(psb_c_base_multivect_type), intent(inout) :: x
class(psb_c_base_multivect_type), intent(inout) :: y
class(psb_c_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
character(len=1), intent(in), optional :: conjgx, conjgy
integer(psb_ipk_) :: i, n
logical :: conjgx_, conjgy_
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
if (z%is_dev()) call z%sync()
if (.not.psb_c_is_complex_) then
call z%mlt(alpha,x%v,y%v,beta,info)
else
conjgx_=.false.
if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
conjgy_=.false.
if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
if (conjgx_) x%v=conjg(x%v)
if (conjgy_) y%v=conjg(y%v)
call z%mlt(alpha,x%v,y%v,beta,info)
if (conjgx_) x%v=conjg(x%v)
if (conjgy_) y%v=conjg(y%v)
end if
end subroutine c_base_mlv_mlt_v_2
!!$
!!$ subroutine c_base_mlv_mlt_av(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
@@ -2400,7 +2413,7 @@ contains
class(psb_c_base_multivect_type), intent(inout) :: x
class(psb_c_base_multivect_type), intent(inout) :: y
if (.not.x%is_host()) call x%sync()
if (x%is_dev()) call x%sync()
if (allocated(x%v)) then
call y%axpby(min(x%get_nrows(),y%get_nrows()),cone,x,czero,info)
call y%absval()
@@ -2428,7 +2441,7 @@ contains
class(psb_c_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2452,7 +2465,7 @@ contains
complex(psb_spk_) :: y(:)
class(psb_c_base_multivect_type) :: x
call x%sync()
if (x%is_dev()) call x%sync()
call x%gth(n,idx%v(i:),y)
end subroutine c_base_mlv_gthzv_x
@@ -2472,8 +2485,8 @@ contains
complex(psb_spk_) :: y(:)
class(psb_c_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2483,40 +2496,42 @@ contains
end subroutine c_base_mlv_gthzv
!!$ !
!!$ ! Scatter:
!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
!!$ !
!!$ !
!!$ !> Function base_mlv_sctb
!!$ !! \memberof psb_c_base_multivect_type
!!$ !! \brief scatter into a class(base_mlv_vect)
!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!!$ !! \param n how many entries to consider
!!$ !! \param idx(:) indices
!!$ !! \param beta
!!$ !! \param x(:)
!!$ subroutine c_base_mlv_sctb(n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: n, idx(:)
!!$ complex(psb_spk_) :: beta, x(:)
!!$ class(psb_c_base_multivect_type) :: y
!!$
!!$ call y%sync()
!!$ call psi_sct(n,idx,x,beta,y%v)
!!$ call y%set_host()
!!$
!!$ end subroutine c_base_mlv_sctb
!!$
!!$ subroutine c_base_mlv_sctb_x(i,n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: i, n
!!$ class(psb_c_base_multivect_type) :: idx
!!$ complex( psb_spk_) :: beta, x(:)
!!$ class(psb_c_base_multivect_type) :: y
!!$
!!$ call y%sct(n,idx%v(i:),x,beta)
!!$
!!$ end subroutine c_base_mlv_sctb_x
!
! Scatter:
! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
!
!
!> Function base_mlv_sctb
!! \memberof psb_c_base_multivect_type
!! \brief scatter into a class(base_mlv_vect)
!! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!! \param n how many entries to consider
!! \param idx(:) indices
!! \param beta
!! \param x(:)
subroutine c_base_mlv_sctb(n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: n, idx(:)
complex(psb_spk_) :: beta, x(:)
class(psb_c_base_multivect_type) :: y
integer(psb_ipk_) :: nc
if (y%is_dev()) call y%sync()
nc = psb_size(y%v,2)
call psi_sct(n,nc,idx,x,beta,y%v)
call y%set_host()
end subroutine c_base_mlv_sctb
subroutine c_base_mlv_sctb_x(i,n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex( psb_spk_) :: beta, x(:)
class(psb_c_base_multivect_type) :: y
call y%sct(n,idx%v(i:),x,beta)
end subroutine c_base_mlv_sctb_x
end module psb_c_base_multivect_mod
base/modules/psb_d_base_vect_mod.f90
+268 -253
View File
@@ -925,7 +925,6 @@ contains
if (z%is_dev()) call z%sync()
n = min(size(z%v), size(x), size(y))
!!$ write(0,*) 'Mlt_a_2: ',n
if (alpha == dzero) then
if (beta == done) then
return
@@ -1192,26 +1191,6 @@ contains
call x%gth(n,idx%v(i:),x%combuf(i:))
end subroutine d_base_gthzbuf
subroutine d_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real(psb_dpk_) :: beta
class(psb_d_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine d_base_sctb_buf
!
!> Function base_device_wait:
!! \memberof psb_d_base_vect_type
@@ -1329,6 +1308,26 @@ contains
end subroutine d_base_sctb_x
subroutine d_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real(psb_dpk_) :: beta
class(psb_d_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine d_base_sctb_buf
end module psb_d_base_vect_mod
@@ -1415,20 +1414,23 @@ module psb_d_base_multivect_mod
procedure, pass(y) :: axpby_v => d_base_mlv_axpby_v
procedure, pass(y) :: axpby_a => d_base_mlv_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_mlv_mlt_v
!!$ procedure, pass(y) :: mlt_a => d_base_mlv_mlt_a
!!$ procedure, pass(z) :: mlt_a_2 => d_base_mlv_mlt_a_2
!!$ procedure, pass(z) :: mlt_v_2 => d_base_mlv_mlt_v_2
!
! MultiVector by vector/multivector multiplication. Need all variants
! to handle multiple requirements from preconditioners
!
procedure, pass(y) :: mlt_mv => d_base_mlv_mlt_mv
procedure, pass(y) :: mlt_mv_v => d_base_mlv_mlt_mv_v
procedure, pass(y) :: mlt_ar1 => d_base_mlv_mlt_ar1
procedure, pass(y) :: mlt_ar2 => d_base_mlv_mlt_ar2
procedure, pass(z) :: mlt_a_2 => d_base_mlv_mlt_a_2
procedure, pass(z) :: mlt_v_2 => d_base_mlv_mlt_v_2
!!$ procedure, pass(z) :: mlt_va => d_base_mlv_mlt_va
!!$ procedure, pass(z) :: mlt_av => d_base_mlv_mlt_av
!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
!!$ !
!!$ ! Scaling and norms
!!$ !
generic, public :: mlt => mlt_mv, mlt_mv_v, mlt_ar1, mlt_ar2, &
& mlt_a_2, mlt_v_2 !, mlt_av, mlt_va
!
! Scaling and norms
!
procedure, pass(x) :: scal => d_base_mlv_scal
procedure, pass(x) :: nrm2 => d_base_mlv_nrm2
procedure, pass(x) :: amax => d_base_mlv_amax
@@ -1444,9 +1446,9 @@ module psb_d_base_multivect_mod
procedure, pass(x) :: gthzv => d_base_mlv_gthzv
procedure, pass(x) :: gthzv_x => d_base_mlv_gthzv_x
generic, public :: gth => gthab, gthzv, gthzv_x
!!$ procedure, pass(y) :: sctb => d_base_mlv_sctb
!!$ procedure, pass(y) :: sctb_x => d_base_mlv_sctb_x
!!$ generic, public :: sct => sctb, sctb_x
procedure, pass(y) :: sctb => d_base_mlv_sctb
procedure, pass(y) :: sctb_x => d_base_mlv_sctb_x
generic, public :: sct => sctb, sctb_x
end type psb_d_base_multivect_type
interface psb_d_base_multivect
@@ -2082,183 +2084,194 @@ contains
end subroutine d_base_mlv_axpby_a
!!$ !
!!$ ! Multiple variants of two operations:
!!$ ! Simple multiplication Y(:) = X(:)*Y(:)
!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!!$ !
!!$ ! Variants expanded according to the dynamic type
!!$ ! of the involved entities
!!$ !
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_d_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a base_mlv_vect array y=x*y
!!$ !! \param x The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine d_base_mlv_mlt_v(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ class(psb_d_base_multivect_type), intent(inout) :: x
!!$ class(psb_d_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ select type(xx => x)
!!$ type is (psb_d_base_multivect_type)
!!$ n = min(size(y%v), size(xx%v))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*xx%v(i)
!!$ end do
!!$ class default
!!$ call y%mlt(x%v,info)
!!$ end select
!!$
!!$ end subroutine d_base_mlv_mlt_v
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_d_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y
!!$ !! \param x(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine d_base_mlv_mlt_a(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ real(psb_dpk_), intent(in) :: x(:)
!!$ class(psb_d_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(y%v), size(x))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*x(i)
!!$ end do
!!$
!!$ end subroutine d_base_mlv_mlt_a
!!$
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a_2
!!$ !! \memberof psb_d_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x(:) The array to be multiplied b
!!$ !! \param y(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine d_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ real(psb_dpk_), intent(in) :: alpha,beta
!!$ real(psb_dpk_), intent(in) :: y(:)
!!$ real(psb_dpk_), intent(in) :: x(:)
!!$ class(psb_d_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(z%v), size(x), size(y))
!!$ if (alpha == izero) then
!!$ if (beta == ione) then
!!$ return
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (alpha == ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + y(i)*x(i)
!!$ end do
!!$ end if
!!$ else if (alpha == -ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = -y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) - y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) - y(i)*x(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = alpha*y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ end if
!!$ end if
!!$ end if
!!$ end subroutine d_base_mlv_mlt_a_2
!!$
!!$ !
!!$ !> Function base_mlv_mlt_v_2
!!$ !! \memberof psb_d_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x The class(base_mlv_vect) to be multiplied b
!!$ !! \param y The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine d_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
!!$ use psi_serial_mod
!!$ use psb_string_mod
!!$ implicit none
!!$ real(psb_dpk_), intent(in) :: alpha,beta
!!$ class(psb_d_base_multivect_type), intent(inout) :: x
!!$ class(psb_d_base_multivect_type), intent(inout) :: y
!!$ class(psb_d_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ character(len=1), intent(in), optional :: conjgx, conjgy
!!$ integer(psb_ipk_) :: i, n
!!$ logical :: conjgx_, conjgy_
!!$
!!$ info = 0
!!$ if (.not.psb_i_is_complex_) then
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ else
!!$ conjgx_=.false.
!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
!!$ conjgy_=.false.
!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ end if
!!$ end subroutine d_base_mlv_mlt_v_2
!
! Multiple variants of two operations:
! Simple multiplication Y(:.:) = X(:,:)*Y(:,:)
! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!
! Variants expanded according to the dynamic type
! of the involved entities
!
!
!> Function base_mlv_mlt_mv
!! \memberof psb_d_base_multivect_type
!! \brief Multivector entry-by-entry multiply by a base_mlv_multivect y=x*y
!! \param x The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine d_base_mlv_mlt_mv(x, y, info)
use psi_serial_mod
implicit none
class(psb_d_base_multivect_type), intent(inout) :: x
class(psb_d_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine d_base_mlv_mlt_mv
subroutine d_base_mlv_mlt_mv_v(x, y, info)
use psi_serial_mod
implicit none
class(psb_d_base_vect_type), intent(inout) :: x
class(psb_d_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine d_base_mlv_mlt_mv_v
!
!> Function base_mlv_mlt_ar1
!! \memberof psb_d_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 1 array y=x*y
!! \param x(:) The array to be multiplied by
!! \param info return code
!!
subroutine d_base_mlv_mlt_ar1(x, y, info)
use psi_serial_mod
implicit none
real(psb_dpk_), intent(in) :: x(:)
class(psb_d_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, n
info = 0
n = min(psb_size(y%v,1), size(x))
do i=1, n
y%v(i,:) = y%v(i,:)*x(i)
end do
end subroutine d_base_mlv_mlt_ar1
!
!> Function base_mlv_mlt_ar2
!! \memberof psb_d_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 2 array y=x*y
!! \param x(:,:) The array to be multiplied by
!! \param info return code
!!
subroutine d_base_mlv_mlt_ar2(x, y, info)
use psi_serial_mod
implicit none
real(psb_dpk_), intent(in) :: x(:,:)
class(psb_d_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr,nc
info = 0
nr = min(psb_size(y%v,1), size(x,1))
nc = min(psb_size(y%v,2), size(x,2))
y%v(1:nr,1:nc) = y%v(1:nr,1:nc)*x(1:nr,1:nc)
end subroutine d_base_mlv_mlt_ar2
!
!> Function base_mlv_mlt_a_2
!! \memberof psb_d_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x(:) The array to be multiplied b
!! \param y(:) The array to be multiplied by
!! \param info return code
!!
subroutine d_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
use psi_serial_mod
implicit none
real(psb_dpk_), intent(in) :: alpha,beta
real(psb_dpk_), intent(in) :: y(:,:)
real(psb_dpk_), intent(in) :: x(:,:)
class(psb_d_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr, nc
info = 0
nr = min(psb_size(z%v,1), size(x,1), size(y,1))
nc = min(psb_size(z%v,2), size(x,2), size(y,2))
if (alpha == dzero) then
if (beta == done) then
return
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc)
end if
else
if (alpha == done) then
if (beta == dzero) then
z%v(1:nr,1:nc) = y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == done) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else if (alpha == -done) then
if (beta == dzero) then
z%v(1:nr,1:nc) = -y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == done) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else
if (beta == dzero) then
z%v(1:nr,1:nc) = alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == done) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
end if
end if
end if
end subroutine d_base_mlv_mlt_a_2
!
!> Function base_mlv_mlt_v_2
!! \memberof psb_d_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x The class(base_mlv_vect) to be multiplied b
!! \param y The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine d_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
use psi_serial_mod
use psb_string_mod
implicit none
real(psb_dpk_), intent(in) :: alpha,beta
class(psb_d_base_multivect_type), intent(inout) :: x
class(psb_d_base_multivect_type), intent(inout) :: y
class(psb_d_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
character(len=1), intent(in), optional :: conjgx, conjgy
integer(psb_ipk_) :: i, n
logical :: conjgx_, conjgy_
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
if (z%is_dev()) call z%sync()
if (.not.psb_d_is_complex_) then
call z%mlt(alpha,x%v,y%v,beta,info)
else
conjgx_=.false.
if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
conjgy_=.false.
if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
if (conjgx_) x%v=(x%v)
if (conjgy_) y%v=(y%v)
call z%mlt(alpha,x%v,y%v,beta,info)
if (conjgx_) x%v=(x%v)
if (conjgy_) y%v=(y%v)
end if
end subroutine d_base_mlv_mlt_v_2
!!$
!!$ subroutine d_base_mlv_mlt_av(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
@@ -2400,7 +2413,7 @@ contains
class(psb_d_base_multivect_type), intent(inout) :: x
class(psb_d_base_multivect_type), intent(inout) :: y
if (.not.x%is_host()) call x%sync()
if (x%is_dev()) call x%sync()
if (allocated(x%v)) then
call y%axpby(min(x%get_nrows(),y%get_nrows()),done,x,dzero,info)
call y%absval()
@@ -2428,7 +2441,7 @@ contains
class(psb_d_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2452,7 +2465,7 @@ contains
real(psb_dpk_) :: y(:)
class(psb_d_base_multivect_type) :: x
call x%sync()
if (x%is_dev()) call x%sync()
call x%gth(n,idx%v(i:),y)
end subroutine d_base_mlv_gthzv_x
@@ -2472,8 +2485,8 @@ contains
real(psb_dpk_) :: y(:)
class(psb_d_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2483,40 +2496,42 @@ contains
end subroutine d_base_mlv_gthzv
!!$ !
!!$ ! Scatter:
!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
!!$ !
!!$ !
!!$ !> Function base_mlv_sctb
!!$ !! \memberof psb_d_base_multivect_type
!!$ !! \brief scatter into a class(base_mlv_vect)
!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!!$ !! \param n how many entries to consider
!!$ !! \param idx(:) indices
!!$ !! \param beta
!!$ !! \param x(:)
!!$ subroutine d_base_mlv_sctb(n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: n, idx(:)
!!$ real(psb_dpk_) :: beta, x(:)
!!$ class(psb_d_base_multivect_type) :: y
!!$
!!$ call y%sync()
!!$ call psi_sct(n,idx,x,beta,y%v)
!!$ call y%set_host()
!!$
!!$ end subroutine d_base_mlv_sctb
!!$
!!$ subroutine d_base_mlv_sctb_x(i,n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: i, n
!!$ class(psb_d_base_multivect_type) :: idx
!!$ real( psb_dpk_) :: beta, x(:)
!!$ class(psb_d_base_multivect_type) :: y
!!$
!!$ call y%sct(n,idx%v(i:),x,beta)
!!$
!!$ end subroutine d_base_mlv_sctb_x
!
! Scatter:
! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
!
!
!> Function base_mlv_sctb
!! \memberof psb_d_base_multivect_type
!! \brief scatter into a class(base_mlv_vect)
!! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!! \param n how many entries to consider
!! \param idx(:) indices
!! \param beta
!! \param x(:)
subroutine d_base_mlv_sctb(n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: n, idx(:)
real(psb_dpk_) :: beta, x(:)
class(psb_d_base_multivect_type) :: y
integer(psb_ipk_) :: nc
if (y%is_dev()) call y%sync()
nc = psb_size(y%v,2)
call psi_sct(n,nc,idx,x,beta,y%v)
call y%set_host()
end subroutine d_base_mlv_sctb
subroutine d_base_mlv_sctb_x(i,n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real( psb_dpk_) :: beta, x(:)
class(psb_d_base_multivect_type) :: y
call y%sct(n,idx%v(i:),x,beta)
end subroutine d_base_mlv_sctb_x
end module psb_d_base_multivect_mod
base/modules/psb_i_base_vect_mod.f90
+20 -20
View File
@@ -732,26 +732,6 @@ contains
call x%gth(n,idx%v(i:),x%combuf(i:))
end subroutine i_base_gthzbuf
subroutine i_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
integer(psb_ipk_) :: beta
class(psb_i_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine i_base_sctb_buf
!
!> Function base_device_wait:
!! \memberof psb_i_base_vect_type
@@ -869,6 +849,26 @@ contains
end subroutine i_base_sctb_x
subroutine i_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
integer(psb_ipk_) :: beta
class(psb_i_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine i_base_sctb_buf
end module psb_i_base_vect_mod
base/modules/psb_s_base_vect_mod.f90
+268 -253
View File
@@ -925,7 +925,6 @@ contains
if (z%is_dev()) call z%sync()
n = min(size(z%v), size(x), size(y))
!!$ write(0,*) 'Mlt_a_2: ',n
if (alpha == szero) then
if (beta == sone) then
return
@@ -1192,26 +1191,6 @@ contains
call x%gth(n,idx%v(i:),x%combuf(i:))
end subroutine s_base_gthzbuf
subroutine s_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real(psb_spk_) :: beta
class(psb_s_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine s_base_sctb_buf
!
!> Function base_device_wait:
!! \memberof psb_s_base_vect_type
@@ -1329,6 +1308,26 @@ contains
end subroutine s_base_sctb_x
subroutine s_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real(psb_spk_) :: beta
class(psb_s_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine s_base_sctb_buf
end module psb_s_base_vect_mod
@@ -1415,20 +1414,23 @@ module psb_s_base_multivect_mod
procedure, pass(y) :: axpby_v => s_base_mlv_axpby_v
procedure, pass(y) :: axpby_a => s_base_mlv_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_mlv_mlt_v
!!$ procedure, pass(y) :: mlt_a => s_base_mlv_mlt_a
!!$ procedure, pass(z) :: mlt_a_2 => s_base_mlv_mlt_a_2
!!$ procedure, pass(z) :: mlt_v_2 => s_base_mlv_mlt_v_2
!
! MultiVector by vector/multivector multiplication. Need all variants
! to handle multiple requirements from preconditioners
!
procedure, pass(y) :: mlt_mv => s_base_mlv_mlt_mv
procedure, pass(y) :: mlt_mv_v => s_base_mlv_mlt_mv_v
procedure, pass(y) :: mlt_ar1 => s_base_mlv_mlt_ar1
procedure, pass(y) :: mlt_ar2 => s_base_mlv_mlt_ar2
procedure, pass(z) :: mlt_a_2 => s_base_mlv_mlt_a_2
procedure, pass(z) :: mlt_v_2 => s_base_mlv_mlt_v_2
!!$ procedure, pass(z) :: mlt_va => s_base_mlv_mlt_va
!!$ procedure, pass(z) :: mlt_av => s_base_mlv_mlt_av
!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
!!$ !
!!$ ! Scaling and norms
!!$ !
generic, public :: mlt => mlt_mv, mlt_mv_v, mlt_ar1, mlt_ar2, &
& mlt_a_2, mlt_v_2 !, mlt_av, mlt_va
!
! Scaling and norms
!
procedure, pass(x) :: scal => s_base_mlv_scal
procedure, pass(x) :: nrm2 => s_base_mlv_nrm2
procedure, pass(x) :: amax => s_base_mlv_amax
@@ -1444,9 +1446,9 @@ module psb_s_base_multivect_mod
procedure, pass(x) :: gthzv => s_base_mlv_gthzv
procedure, pass(x) :: gthzv_x => s_base_mlv_gthzv_x
generic, public :: gth => gthab, gthzv, gthzv_x
!!$ procedure, pass(y) :: sctb => s_base_mlv_sctb
!!$ procedure, pass(y) :: sctb_x => s_base_mlv_sctb_x
!!$ generic, public :: sct => sctb, sctb_x
procedure, pass(y) :: sctb => s_base_mlv_sctb
procedure, pass(y) :: sctb_x => s_base_mlv_sctb_x
generic, public :: sct => sctb, sctb_x
end type psb_s_base_multivect_type
interface psb_s_base_multivect
@@ -2082,183 +2084,194 @@ contains
end subroutine s_base_mlv_axpby_a
!!$ !
!!$ ! Multiple variants of two operations:
!!$ ! Simple multiplication Y(:) = X(:)*Y(:)
!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!!$ !
!!$ ! Variants expanded according to the dynamic type
!!$ ! of the involved entities
!!$ !
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_s_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a base_mlv_vect array y=x*y
!!$ !! \param x The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine s_base_mlv_mlt_v(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ class(psb_s_base_multivect_type), intent(inout) :: x
!!$ class(psb_s_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ select type(xx => x)
!!$ type is (psb_s_base_multivect_type)
!!$ n = min(size(y%v), size(xx%v))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*xx%v(i)
!!$ end do
!!$ class default
!!$ call y%mlt(x%v,info)
!!$ end select
!!$
!!$ end subroutine s_base_mlv_mlt_v
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_s_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y
!!$ !! \param x(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine s_base_mlv_mlt_a(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ real(psb_spk_), intent(in) :: x(:)
!!$ class(psb_s_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(y%v), size(x))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*x(i)
!!$ end do
!!$
!!$ end subroutine s_base_mlv_mlt_a
!!$
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a_2
!!$ !! \memberof psb_s_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x(:) The array to be multiplied b
!!$ !! \param y(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine s_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ real(psb_spk_), intent(in) :: alpha,beta
!!$ real(psb_spk_), intent(in) :: y(:)
!!$ real(psb_spk_), intent(in) :: x(:)
!!$ class(psb_s_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(z%v), size(x), size(y))
!!$ if (alpha == izero) then
!!$ if (beta == ione) then
!!$ return
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (alpha == ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + y(i)*x(i)
!!$ end do
!!$ end if
!!$ else if (alpha == -ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = -y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) - y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) - y(i)*x(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = alpha*y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ end if
!!$ end if
!!$ end if
!!$ end subroutine s_base_mlv_mlt_a_2
!!$
!!$ !
!!$ !> Function base_mlv_mlt_v_2
!!$ !! \memberof psb_s_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x The class(base_mlv_vect) to be multiplied b
!!$ !! \param y The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine s_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
!!$ use psi_serial_mod
!!$ use psb_string_mod
!!$ implicit none
!!$ real(psb_spk_), intent(in) :: alpha,beta
!!$ class(psb_s_base_multivect_type), intent(inout) :: x
!!$ class(psb_s_base_multivect_type), intent(inout) :: y
!!$ class(psb_s_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ character(len=1), intent(in), optional :: conjgx, conjgy
!!$ integer(psb_ipk_) :: i, n
!!$ logical :: conjgx_, conjgy_
!!$
!!$ info = 0
!!$ if (.not.psb_i_is_complex_) then
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ else
!!$ conjgx_=.false.
!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
!!$ conjgy_=.false.
!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ end if
!!$ end subroutine s_base_mlv_mlt_v_2
!
! Multiple variants of two operations:
! Simple multiplication Y(:.:) = X(:,:)*Y(:,:)
! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!
! Variants expanded according to the dynamic type
! of the involved entities
!
!
!> Function base_mlv_mlt_mv
!! \memberof psb_s_base_multivect_type
!! \brief Multivector entry-by-entry multiply by a base_mlv_multivect y=x*y
!! \param x The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine s_base_mlv_mlt_mv(x, y, info)
use psi_serial_mod
implicit none
class(psb_s_base_multivect_type), intent(inout) :: x
class(psb_s_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine s_base_mlv_mlt_mv
subroutine s_base_mlv_mlt_mv_v(x, y, info)
use psi_serial_mod
implicit none
class(psb_s_base_vect_type), intent(inout) :: x
class(psb_s_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine s_base_mlv_mlt_mv_v
!
!> Function base_mlv_mlt_ar1
!! \memberof psb_s_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 1 array y=x*y
!! \param x(:) The array to be multiplied by
!! \param info return code
!!
subroutine s_base_mlv_mlt_ar1(x, y, info)
use psi_serial_mod
implicit none
real(psb_spk_), intent(in) :: x(:)
class(psb_s_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, n
info = 0
n = min(psb_size(y%v,1), size(x))
do i=1, n
y%v(i,:) = y%v(i,:)*x(i)
end do
end subroutine s_base_mlv_mlt_ar1
!
!> Function base_mlv_mlt_ar2
!! \memberof psb_s_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 2 array y=x*y
!! \param x(:,:) The array to be multiplied by
!! \param info return code
!!
subroutine s_base_mlv_mlt_ar2(x, y, info)
use psi_serial_mod
implicit none
real(psb_spk_), intent(in) :: x(:,:)
class(psb_s_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr,nc
info = 0
nr = min(psb_size(y%v,1), size(x,1))
nc = min(psb_size(y%v,2), size(x,2))
y%v(1:nr,1:nc) = y%v(1:nr,1:nc)*x(1:nr,1:nc)
end subroutine s_base_mlv_mlt_ar2
!
!> Function base_mlv_mlt_a_2
!! \memberof psb_s_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x(:) The array to be multiplied b
!! \param y(:) The array to be multiplied by
!! \param info return code
!!
subroutine s_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
use psi_serial_mod
implicit none
real(psb_spk_), intent(in) :: alpha,beta
real(psb_spk_), intent(in) :: y(:,:)
real(psb_spk_), intent(in) :: x(:,:)
class(psb_s_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr, nc
info = 0
nr = min(psb_size(z%v,1), size(x,1), size(y,1))
nc = min(psb_size(z%v,2), size(x,2), size(y,2))
if (alpha == szero) then
if (beta == sone) then
return
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc)
end if
else
if (alpha == sone) then
if (beta == szero) then
z%v(1:nr,1:nc) = y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == sone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else if (alpha == -sone) then
if (beta == szero) then
z%v(1:nr,1:nc) = -y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == sone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else
if (beta == szero) then
z%v(1:nr,1:nc) = alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == sone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
end if
end if
end if
end subroutine s_base_mlv_mlt_a_2
!
!> Function base_mlv_mlt_v_2
!! \memberof psb_s_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x The class(base_mlv_vect) to be multiplied b
!! \param y The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine s_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
use psi_serial_mod
use psb_string_mod
implicit none
real(psb_spk_), intent(in) :: alpha,beta
class(psb_s_base_multivect_type), intent(inout) :: x
class(psb_s_base_multivect_type), intent(inout) :: y
class(psb_s_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
character(len=1), intent(in), optional :: conjgx, conjgy
integer(psb_ipk_) :: i, n
logical :: conjgx_, conjgy_
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
if (z%is_dev()) call z%sync()
if (.not.psb_s_is_complex_) then
call z%mlt(alpha,x%v,y%v,beta,info)
else
conjgx_=.false.
if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
conjgy_=.false.
if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
if (conjgx_) x%v=(x%v)
if (conjgy_) y%v=(y%v)
call z%mlt(alpha,x%v,y%v,beta,info)
if (conjgx_) x%v=(x%v)
if (conjgy_) y%v=(y%v)
end if
end subroutine s_base_mlv_mlt_v_2
!!$
!!$ subroutine s_base_mlv_mlt_av(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
@@ -2400,7 +2413,7 @@ contains
class(psb_s_base_multivect_type), intent(inout) :: x
class(psb_s_base_multivect_type), intent(inout) :: y
if (.not.x%is_host()) call x%sync()
if (x%is_dev()) call x%sync()
if (allocated(x%v)) then
call y%axpby(min(x%get_nrows(),y%get_nrows()),sone,x,szero,info)
call y%absval()
@@ -2428,7 +2441,7 @@ contains
class(psb_s_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2452,7 +2465,7 @@ contains
real(psb_spk_) :: y(:)
class(psb_s_base_multivect_type) :: x
call x%sync()
if (x%is_dev()) call x%sync()
call x%gth(n,idx%v(i:),y)
end subroutine s_base_mlv_gthzv_x
@@ -2472,8 +2485,8 @@ contains
real(psb_spk_) :: y(:)
class(psb_s_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2483,40 +2496,42 @@ contains
end subroutine s_base_mlv_gthzv
!!$ !
!!$ ! Scatter:
!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
!!$ !
!!$ !
!!$ !> Function base_mlv_sctb
!!$ !! \memberof psb_s_base_multivect_type
!!$ !! \brief scatter into a class(base_mlv_vect)
!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!!$ !! \param n how many entries to consider
!!$ !! \param idx(:) indices
!!$ !! \param beta
!!$ !! \param x(:)
!!$ subroutine s_base_mlv_sctb(n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: n, idx(:)
!!$ real(psb_spk_) :: beta, x(:)
!!$ class(psb_s_base_multivect_type) :: y
!!$
!!$ call y%sync()
!!$ call psi_sct(n,idx,x,beta,y%v)
!!$ call y%set_host()
!!$
!!$ end subroutine s_base_mlv_sctb
!!$
!!$ subroutine s_base_mlv_sctb_x(i,n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: i, n
!!$ class(psb_s_base_multivect_type) :: idx
!!$ real( psb_spk_) :: beta, x(:)
!!$ class(psb_s_base_multivect_type) :: y
!!$
!!$ call y%sct(n,idx%v(i:),x,beta)
!!$
!!$ end subroutine s_base_mlv_sctb_x
!
! Scatter:
! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
!
!
!> Function base_mlv_sctb
!! \memberof psb_s_base_multivect_type
!! \brief scatter into a class(base_mlv_vect)
!! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!! \param n how many entries to consider
!! \param idx(:) indices
!! \param beta
!! \param x(:)
subroutine s_base_mlv_sctb(n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: n, idx(:)
real(psb_spk_) :: beta, x(:)
class(psb_s_base_multivect_type) :: y
integer(psb_ipk_) :: nc
if (y%is_dev()) call y%sync()
nc = psb_size(y%v,2)
call psi_sct(n,nc,idx,x,beta,y%v)
call y%set_host()
end subroutine s_base_mlv_sctb
subroutine s_base_mlv_sctb_x(i,n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
real( psb_spk_) :: beta, x(:)
class(psb_s_base_multivect_type) :: y
call y%sct(n,idx%v(i:),x,beta)
end subroutine s_base_mlv_sctb_x
end module psb_s_base_multivect_mod
base/modules/psb_z_base_vect_mod.f90
+268 -253
View File
@@ -925,7 +925,6 @@ contains
if (z%is_dev()) call z%sync()
n = min(size(z%v), size(x), size(y))
!!$ write(0,*) 'Mlt_a_2: ',n
if (alpha == zzero) then
if (beta == zone) then
return
@@ -1192,26 +1191,6 @@ contains
call x%gth(n,idx%v(i:),x%combuf(i:))
end subroutine z_base_gthzbuf
subroutine z_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex(psb_dpk_) :: beta
class(psb_z_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine z_base_sctb_buf
!
!> Function base_device_wait:
!! \memberof psb_z_base_vect_type
@@ -1329,6 +1308,26 @@ contains
end subroutine z_base_sctb_x
subroutine z_base_sctb_buf(i,n,idx,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex(psb_dpk_) :: beta
class(psb_z_base_vect_type) :: y
if (.not.allocated(y%combuf)) then
call psb_errpush(psb_err_alloc_dealloc_,'sctb_buf')
return
end if
if (y%is_dev()) call y%sync()
if (idx%is_dev()) call idx%sync()
call y%sct(n,idx%v(i:),y%combuf(i:),beta)
call y%set_host()
end subroutine z_base_sctb_buf
end module psb_z_base_vect_mod
@@ -1415,20 +1414,23 @@ module psb_z_base_multivect_mod
procedure, pass(y) :: axpby_v => z_base_mlv_axpby_v
procedure, pass(y) :: axpby_a => z_base_mlv_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_mlv_mlt_v
!!$ procedure, pass(y) :: mlt_a => z_base_mlv_mlt_a
!!$ procedure, pass(z) :: mlt_a_2 => z_base_mlv_mlt_a_2
!!$ procedure, pass(z) :: mlt_v_2 => z_base_mlv_mlt_v_2
!
! MultiVector by vector/multivector multiplication. Need all variants
! to handle multiple requirements from preconditioners
!
procedure, pass(y) :: mlt_mv => z_base_mlv_mlt_mv
procedure, pass(y) :: mlt_mv_v => z_base_mlv_mlt_mv_v
procedure, pass(y) :: mlt_ar1 => z_base_mlv_mlt_ar1
procedure, pass(y) :: mlt_ar2 => z_base_mlv_mlt_ar2
procedure, pass(z) :: mlt_a_2 => z_base_mlv_mlt_a_2
procedure, pass(z) :: mlt_v_2 => z_base_mlv_mlt_v_2
!!$ procedure, pass(z) :: mlt_va => z_base_mlv_mlt_va
!!$ procedure, pass(z) :: mlt_av => z_base_mlv_mlt_av
!!$ generic, public :: mlt => mlt_v, mlt_a, mlt_a_2, mlt_v_2, mlt_av, mlt_va
!!$ !
!!$ ! Scaling and norms
!!$ !
generic, public :: mlt => mlt_mv, mlt_mv_v, mlt_ar1, mlt_ar2, &
& mlt_a_2, mlt_v_2 !, mlt_av, mlt_va
!
! Scaling and norms
!
procedure, pass(x) :: scal => z_base_mlv_scal
procedure, pass(x) :: nrm2 => z_base_mlv_nrm2
procedure, pass(x) :: amax => z_base_mlv_amax
@@ -1444,9 +1446,9 @@ module psb_z_base_multivect_mod
procedure, pass(x) :: gthzv => z_base_mlv_gthzv
procedure, pass(x) :: gthzv_x => z_base_mlv_gthzv_x
generic, public :: gth => gthab, gthzv, gthzv_x
!!$ procedure, pass(y) :: sctb => z_base_mlv_sctb
!!$ procedure, pass(y) :: sctb_x => z_base_mlv_sctb_x
!!$ generic, public :: sct => sctb, sctb_x
procedure, pass(y) :: sctb => z_base_mlv_sctb
procedure, pass(y) :: sctb_x => z_base_mlv_sctb_x
generic, public :: sct => sctb, sctb_x
end type psb_z_base_multivect_type
interface psb_z_base_multivect
@@ -2082,183 +2084,194 @@ contains
end subroutine z_base_mlv_axpby_a
!!$ !
!!$ ! Multiple variants of two operations:
!!$ ! Simple multiplication Y(:) = X(:)*Y(:)
!!$ ! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!!$ !
!!$ ! Variants expanded according to the dynamic type
!!$ ! of the involved entities
!!$ !
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_z_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a base_mlv_vect array y=x*y
!!$ !! \param x The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine z_base_mlv_mlt_v(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ class(psb_z_base_multivect_type), intent(inout) :: x
!!$ class(psb_z_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ select type(xx => x)
!!$ type is (psb_z_base_multivect_type)
!!$ n = min(size(y%v), size(xx%v))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*xx%v(i)
!!$ end do
!!$ class default
!!$ call y%mlt(x%v,info)
!!$ end select
!!$
!!$ end subroutine z_base_mlv_mlt_v
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a
!!$ !! \memberof psb_z_base_multivect_type
!!$ !! \brief Vector entry-by-entry multiply by a normal array y=x*y
!!$ !! \param x(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine z_base_mlv_mlt_a(x, y, info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ complex(psb_dpk_), intent(in) :: x(:)
!!$ class(psb_z_base_multivect_type), intent(inout) :: y
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(y%v), size(x))
!!$ do i=1, n
!!$ y%v(i) = y%v(i)*x(i)
!!$ end do
!!$
!!$ end subroutine z_base_mlv_mlt_a
!!$
!!$
!!$ !
!!$ !> Function base_mlv_mlt_a_2
!!$ !! \memberof psb_z_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x(:) The array to be multiplied b
!!$ !! \param y(:) The array to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine z_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
!!$ implicit none
!!$ complex(psb_dpk_), intent(in) :: alpha,beta
!!$ complex(psb_dpk_), intent(in) :: y(:)
!!$ complex(psb_dpk_), intent(in) :: x(:)
!!$ class(psb_z_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ integer(psb_ipk_) :: i, n
!!$
!!$ info = 0
!!$ n = min(size(z%v), size(x), size(y))
!!$ if (alpha == izero) then
!!$ if (beta == ione) then
!!$ return
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (alpha == ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + y(i)*x(i)
!!$ end do
!!$ end if
!!$ else if (alpha == -ione) then
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = -y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) - y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) - y(i)*x(i)
!!$ end do
!!$ end if
!!$ else
!!$ if (beta == izero) then
!!$ do i=1, n
!!$ z%v(i) = alpha*y(i)*x(i)
!!$ end do
!!$ else if (beta == ione) then
!!$ do i=1, n
!!$ z%v(i) = z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ else
!!$ do i=1, n
!!$ z%v(i) = beta*z%v(i) + alpha*y(i)*x(i)
!!$ end do
!!$ end if
!!$ end if
!!$ end if
!!$ end subroutine z_base_mlv_mlt_a_2
!!$
!!$ !
!!$ !> Function base_mlv_mlt_v_2
!!$ !! \memberof psb_z_base_multivect_type
!!$ !! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!!$ !! z=beta*z+alpha*x*y
!!$ !! \param alpha
!!$ !! \param beta
!!$ !! \param x The class(base_mlv_vect) to be multiplied b
!!$ !! \param y The class(base_mlv_vect) to be multiplied by
!!$ !! \param info return code
!!$ !!
!!$ subroutine z_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
!!$ use psi_serial_mod
!!$ use psb_string_mod
!!$ implicit none
!!$ complex(psb_dpk_), intent(in) :: alpha,beta
!!$ class(psb_z_base_multivect_type), intent(inout) :: x
!!$ class(psb_z_base_multivect_type), intent(inout) :: y
!!$ class(psb_z_base_multivect_type), intent(inout) :: z
!!$ integer(psb_ipk_), intent(out) :: info
!!$ character(len=1), intent(in), optional :: conjgx, conjgy
!!$ integer(psb_ipk_) :: i, n
!!$ logical :: conjgx_, conjgy_
!!$
!!$ info = 0
!!$ if (.not.psb_i_is_complex_) then
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ else
!!$ conjgx_=.false.
!!$ if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
!!$ conjgy_=.false.
!!$ if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ call z%mlt(alpha,x%v,y%v,beta,info)
!!$ if (conjgx_) x%v=(x%v)
!!$ if (conjgy_) y%v=(y%v)
!!$ end if
!!$ end subroutine z_base_mlv_mlt_v_2
!
! Multiple variants of two operations:
! Simple multiplication Y(:.:) = X(:,:)*Y(:,:)
! blas-like: Z(:) = alpha*X(:)*Y(:)+beta*Z(:)
!
! Variants expanded according to the dynamic type
! of the involved entities
!
!
!> Function base_mlv_mlt_mv
!! \memberof psb_z_base_multivect_type
!! \brief Multivector entry-by-entry multiply by a base_mlv_multivect y=x*y
!! \param x The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine z_base_mlv_mlt_mv(x, y, info)
use psi_serial_mod
implicit none
class(psb_z_base_multivect_type), intent(inout) :: x
class(psb_z_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine z_base_mlv_mlt_mv
subroutine z_base_mlv_mlt_mv_v(x, y, info)
use psi_serial_mod
implicit none
class(psb_z_base_vect_type), intent(inout) :: x
class(psb_z_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
info = 0
if (x%is_dev()) call x%sync()
call y%mlt(x%v,info)
end subroutine z_base_mlv_mlt_mv_v
!
!> Function base_mlv_mlt_ar1
!! \memberof psb_z_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 1 array y=x*y
!! \param x(:) The array to be multiplied by
!! \param info return code
!!
subroutine z_base_mlv_mlt_ar1(x, y, info)
use psi_serial_mod
implicit none
complex(psb_dpk_), intent(in) :: x(:)
class(psb_z_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, n
info = 0
n = min(psb_size(y%v,1), size(x))
do i=1, n
y%v(i,:) = y%v(i,:)*x(i)
end do
end subroutine z_base_mlv_mlt_ar1
!
!> Function base_mlv_mlt_ar2
!! \memberof psb_z_base_multivect_type
!! \brief MultiVector entry-by-entry multiply by a rank 2 array y=x*y
!! \param x(:,:) The array to be multiplied by
!! \param info return code
!!
subroutine z_base_mlv_mlt_ar2(x, y, info)
use psi_serial_mod
implicit none
complex(psb_dpk_), intent(in) :: x(:,:)
class(psb_z_base_multivect_type), intent(inout) :: y
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr,nc
info = 0
nr = min(psb_size(y%v,1), size(x,1))
nc = min(psb_size(y%v,2), size(x,2))
y%v(1:nr,1:nc) = y%v(1:nr,1:nc)*x(1:nr,1:nc)
end subroutine z_base_mlv_mlt_ar2
!
!> Function base_mlv_mlt_a_2
!! \memberof psb_z_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by normal arrays
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x(:) The array to be multiplied b
!! \param y(:) The array to be multiplied by
!! \param info return code
!!
subroutine z_base_mlv_mlt_a_2(alpha,x,y,beta,z,info)
use psi_serial_mod
implicit none
complex(psb_dpk_), intent(in) :: alpha,beta
complex(psb_dpk_), intent(in) :: y(:,:)
complex(psb_dpk_), intent(in) :: x(:,:)
class(psb_z_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
integer(psb_ipk_) :: i, nr, nc
info = 0
nr = min(psb_size(z%v,1), size(x,1), size(y,1))
nc = min(psb_size(z%v,2), size(x,2), size(y,2))
if (alpha == zzero) then
if (beta == zone) then
return
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc)
end if
else
if (alpha == zone) then
if (beta == zzero) then
z%v(1:nr,1:nc) = y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == zone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else if (alpha == -zone) then
if (beta == zzero) then
z%v(1:nr,1:nc) = -y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == zone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) - y(1:nr,1:nc)*x(1:nr,1:nc)
end if
else
if (beta == zzero) then
z%v(1:nr,1:nc) = alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else if (beta == zone) then
z%v(1:nr,1:nc) = z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
else
z%v(1:nr,1:nc) = beta*z%v(1:nr,1:nc) + alpha*y(1:nr,1:nc)*x(1:nr,1:nc)
end if
end if
end if
end subroutine z_base_mlv_mlt_a_2
!
!> Function base_mlv_mlt_v_2
!! \memberof psb_z_base_multivect_type
!! \brief AXPBY-like Vector entry-by-entry multiply by class(base_mlv_vect)
!! z=beta*z+alpha*x*y
!! \param alpha
!! \param beta
!! \param x The class(base_mlv_vect) to be multiplied b
!! \param y The class(base_mlv_vect) to be multiplied by
!! \param info return code
!!
subroutine z_base_mlv_mlt_v_2(alpha,x,y,beta,z,info,conjgx,conjgy)
use psi_serial_mod
use psb_string_mod
implicit none
complex(psb_dpk_), intent(in) :: alpha,beta
class(psb_z_base_multivect_type), intent(inout) :: x
class(psb_z_base_multivect_type), intent(inout) :: y
class(psb_z_base_multivect_type), intent(inout) :: z
integer(psb_ipk_), intent(out) :: info
character(len=1), intent(in), optional :: conjgx, conjgy
integer(psb_ipk_) :: i, n
logical :: conjgx_, conjgy_
info = 0
if (x%is_dev()) call x%sync()
if (y%is_dev()) call y%sync()
if (z%is_dev()) call z%sync()
if (.not.psb_z_is_complex_) then
call z%mlt(alpha,x%v,y%v,beta,info)
else
conjgx_=.false.
if (present(conjgx)) conjgx_ = (psb_toupper(conjgx)=='C')
conjgy_=.false.
if (present(conjgy)) conjgy_ = (psb_toupper(conjgy)=='C')
if (conjgx_) x%v=conjg(x%v)
if (conjgy_) y%v=conjg(y%v)
call z%mlt(alpha,x%v,y%v,beta,info)
if (conjgx_) x%v=conjg(x%v)
if (conjgy_) y%v=conjg(y%v)
end if
end subroutine z_base_mlv_mlt_v_2
!!$
!!$ subroutine z_base_mlv_mlt_av(alpha,x,y,beta,z,info)
!!$ use psi_serial_mod
@@ -2400,7 +2413,7 @@ contains
class(psb_z_base_multivect_type), intent(inout) :: x
class(psb_z_base_multivect_type), intent(inout) :: y
if (.not.x%is_host()) call x%sync()
if (x%is_dev()) call x%sync()
if (allocated(x%v)) then
call y%axpby(min(x%get_nrows(),y%get_nrows()),zone,x,zzero,info)
call y%absval()
@@ -2428,7 +2441,7 @@ contains
class(psb_z_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2452,7 +2465,7 @@ contains
complex(psb_dpk_) :: y(:)
class(psb_z_base_multivect_type) :: x
call x%sync()
if (x%is_dev()) call x%sync()
call x%gth(n,idx%v(i:),y)
end subroutine z_base_mlv_gthzv_x
@@ -2472,8 +2485,8 @@ contains
complex(psb_dpk_) :: y(:)
class(psb_z_base_multivect_type) :: x
integer(psb_ipk_) :: nc
call x%sync()
if (x%is_dev()) call x%sync()
if (.not.allocated(x%v)) then
return
end if
@@ -2483,40 +2496,42 @@ contains
end subroutine z_base_mlv_gthzv
!!$ !
!!$ ! Scatter:
!!$ ! Y(IDX(:)) = beta*Y(IDX(:)) + X(:)
!!$ !
!!$ !
!!$ !> Function base_mlv_sctb
!!$ !! \memberof psb_z_base_multivect_type
!!$ !! \brief scatter into a class(base_mlv_vect)
!!$ !! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!!$ !! \param n how many entries to consider
!!$ !! \param idx(:) indices
!!$ !! \param beta
!!$ !! \param x(:)
!!$ subroutine z_base_mlv_sctb(n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: n, idx(:)
!!$ complex(psb_dpk_) :: beta, x(:)
!!$ class(psb_z_base_multivect_type) :: y
!!$
!!$ call y%sync()
!!$ call psi_sct(n,idx,x,beta,y%v)
!!$ call y%set_host()
!!$
!!$ end subroutine z_base_mlv_sctb
!!$
!!$ subroutine z_base_mlv_sctb_x(i,n,idx,x,beta,y)
!!$ use psi_serial_mod
!!$ integer(psb_ipk_) :: i, n
!!$ class(psb_z_base_multivect_type) :: idx
!!$ complex( psb_dpk_) :: beta, x(:)
!!$ class(psb_z_base_multivect_type) :: y
!!$
!!$ call y%sct(n,idx%v(i:),x,beta)
!!$
!!$ end subroutine z_base_mlv_sctb_x
!
! Scatter:
! Y(IDX(:),:) = beta*Y(IDX(:),:) + X(:)
!
!
!> Function base_mlv_sctb
!! \memberof psb_z_base_multivect_type
!! \brief scatter into a class(base_mlv_vect)
!! Y(IDX(:)) = beta * Y(IDX(:)) + X(:)
!! \param n how many entries to consider
!! \param idx(:) indices
!! \param beta
!! \param x(:)
subroutine z_base_mlv_sctb(n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: n, idx(:)
complex(psb_dpk_) :: beta, x(:)
class(psb_z_base_multivect_type) :: y
integer(psb_ipk_) :: nc
if (y%is_dev()) call y%sync()
nc = psb_size(y%v,2)
call psi_sct(n,nc,idx,x,beta,y%v)
call y%set_host()
end subroutine z_base_mlv_sctb
subroutine z_base_mlv_sctb_x(i,n,idx,x,beta,y)
use psi_serial_mod
integer(psb_ipk_) :: i, n
class(psb_i_base_vect_type) :: idx
complex( psb_dpk_) :: beta, x(:)
class(psb_z_base_multivect_type) :: y
call y%sct(n,idx%v(i:),x,beta)
end subroutine z_base_mlv_sctb_x
end module psb_z_base_multivect_mod