From 44f05510bc510a260151ab29e04a7c6d42f61ac6 Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Tue, 10 Mar 2020 17:34:08 +0100 Subject: [PATCH] Added out-of-place axpby and relative c interface --- base/modules/auxil/psi_c_serial_mod.f90 | 44 +- base/modules/auxil/psi_d_serial_mod.f90 | 44 +- base/modules/auxil/psi_e_serial_mod.f90 | 44 +- base/modules/auxil/psi_m_serial_mod.f90 | 44 +- base/modules/auxil/psi_s_serial_mod.f90 | 44 +- base/modules/auxil/psi_z_serial_mod.f90 | 44 +- base/modules/psblas/psb_c_psblas_mod.F90 | 22 + base/modules/psblas/psb_d_psblas_mod.F90 | 22 + base/modules/psblas/psb_s_psblas_mod.F90 | 22 + base/modules/psblas/psb_z_psblas_mod.F90 | 22 + base/modules/serial/psb_c_base_vect_mod.f90 | 66 ++- base/modules/serial/psb_c_vect_mod.F90 | 36 +- base/modules/serial/psb_d_base_vect_mod.f90 | 66 ++- base/modules/serial/psb_d_vect_mod.F90 | 36 +- base/modules/serial/psb_s_base_vect_mod.f90 | 66 ++- base/modules/serial/psb_s_vect_mod.F90 | 36 +- base/modules/serial/psb_z_base_vect_mod.f90 | 66 ++- base/modules/serial/psb_z_vect_mod.F90 | 36 +- base/psblas/psb_caxpby.f90 | 278 ++++++++++++ base/psblas/psb_daxpby.f90 | 278 ++++++++++++ base/psblas/psb_saxpby.f90 | 278 ++++++++++++ base/psblas/psb_zaxpby.f90 | 278 ++++++++++++ base/serial/psi_c_serial_impl.f90 | 453 +++++++++++++++----- base/serial/psi_d_serial_impl.f90 | 453 +++++++++++++++----- base/serial/psi_e_serial_impl.f90 | 453 +++++++++++++++----- base/serial/psi_m_serial_impl.f90 | 453 +++++++++++++++----- base/serial/psi_s_serial_impl.f90 | 453 +++++++++++++++----- base/serial/psi_z_serial_impl.f90 | 453 +++++++++++++++----- cbind/base/psb_c_cbase.h | 2 + cbind/base/psb_c_dbase.h | 2 + cbind/base/psb_c_psblas_cbind_mod.f90 | 43 ++ cbind/base/psb_c_sbase.h | 2 + cbind/base/psb_c_zbase.h | 2 + cbind/base/psb_d_psblas_cbind_mod.f90 | 43 ++ cbind/base/psb_s_psblas_cbind_mod.f90 | 43 ++ cbind/base/psb_z_psblas_cbind_mod.f90 | 43 ++ 36 files changed, 4018 insertions(+), 752 deletions(-) diff --git a/base/modules/auxil/psi_c_serial_mod.f90 b/base/modules/auxil/psi_c_serial_mod.f90 index 5dc25dc4..05145c1c 100644 --- a/base/modules/auxil/psi_c_serial_mod.f90 +++ b/base/modules/auxil/psi_c_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_c_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_, psb_spk_ @@ -36,7 +36,7 @@ module psi_c_serial_mod ! 2-D version subroutine psb_cgelp(trans,iperm,x,info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none complex(psb_spk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_c_serial_mod end subroutine psb_cgelp subroutine psb_cgelpv(trans,iperm,x,info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none complex(psb_spk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_cgelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_caxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n complex(psb_spk_), intent (in) :: x(:,:) complex(psb_spk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_c_serial_mod end subroutine psi_caxpby subroutine psi_caxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m complex(psb_spk_), intent (in) :: x(:) complex(psb_spk_), intent (inout) :: y(:) complex(psb_spk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_caxpbyv + subroutine psi_caxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_spk_ + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_spk_), intent (in) :: x(:) + complex(psb_spk_), intent (in) :: y(:) + complex(psb_spk_), intent (in) :: z(:) + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_caxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_c_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) complex(psb_spk_) :: x(:,:), y(:) - + end subroutine psi_cgthzmv subroutine psi_cgthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_spk_ implicit none integer(psb_ipk_) :: n, k, idx(:) complex(psb_spk_) :: x(:,:), y(:,:) - + end subroutine psi_cgthzmm subroutine psi_cgthzv(n,idx,x,y) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) complex(psb_spk_) :: x(:), y(:) end subroutine psi_cgthzv @@ -124,7 +134,7 @@ module psi_c_serial_mod subroutine psi_csctv(n,idx,x,beta,y) import :: psb_ipk_, psb_spk_ implicit none - + integer(psb_ipk_) :: n, idx(:) complex(psb_spk_) :: beta, x(:), y(:) end subroutine psi_csctv diff --git a/base/modules/auxil/psi_d_serial_mod.f90 b/base/modules/auxil/psi_d_serial_mod.f90 index 541446c7..0bea1bce 100644 --- a/base/modules/auxil/psi_d_serial_mod.f90 +++ b/base/modules/auxil/psi_d_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_d_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_, psb_dpk_ @@ -36,7 +36,7 @@ module psi_d_serial_mod ! 2-D version subroutine psb_dgelp(trans,iperm,x,info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none real(psb_dpk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_d_serial_mod end subroutine psb_dgelp subroutine psb_dgelpv(trans,iperm,x,info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none real(psb_dpk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_dgelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_daxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n real(psb_dpk_), intent (in) :: x(:,:) real(psb_dpk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_d_serial_mod end subroutine psi_daxpby subroutine psi_daxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m real(psb_dpk_), intent (in) :: x(:) real(psb_dpk_), intent (inout) :: y(:) real(psb_dpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_daxpbyv + subroutine psi_daxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_dpk_ + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_dpk_), intent (in) :: x(:) + real(psb_dpk_), intent (in) :: y(:) + real(psb_dpk_), intent (in) :: z(:) + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_daxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_d_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) real(psb_dpk_) :: x(:,:), y(:) - + end subroutine psi_dgthzmv subroutine psi_dgthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_dpk_ implicit none integer(psb_ipk_) :: n, k, idx(:) real(psb_dpk_) :: x(:,:), y(:,:) - + end subroutine psi_dgthzmm subroutine psi_dgthzv(n,idx,x,y) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) real(psb_dpk_) :: x(:), y(:) end subroutine psi_dgthzv @@ -124,7 +134,7 @@ module psi_d_serial_mod subroutine psi_dsctv(n,idx,x,beta,y) import :: psb_ipk_, psb_dpk_ implicit none - + integer(psb_ipk_) :: n, idx(:) real(psb_dpk_) :: beta, x(:), y(:) end subroutine psi_dsctv diff --git a/base/modules/auxil/psi_e_serial_mod.f90 b/base/modules/auxil/psi_e_serial_mod.f90 index f95d812b..f8b9694d 100644 --- a/base/modules/auxil/psi_e_serial_mod.f90 +++ b/base/modules/auxil/psi_e_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_e_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_ @@ -36,7 +36,7 @@ module psi_e_serial_mod ! 2-D version subroutine psb_egelp(trans,iperm,x,info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_epk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_e_serial_mod end subroutine psb_egelp subroutine psb_egelpv(trans,iperm,x,info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_epk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_egelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_eaxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n integer(psb_epk_), intent (in) :: x(:,:) integer(psb_epk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_e_serial_mod end subroutine psi_eaxpby subroutine psi_eaxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m integer(psb_epk_), intent (in) :: x(:) integer(psb_epk_), intent (inout) :: y(:) integer(psb_epk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_eaxpbyv + subroutine psi_eaxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ + implicit none + integer(psb_ipk_), intent(in) :: m + integer(psb_epk_), intent (in) :: x(:) + integer(psb_epk_), intent (in) :: y(:) + integer(psb_epk_), intent (in) :: z(:) + integer(psb_epk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_eaxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_e_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) integer(psb_epk_) :: x(:,:), y(:) - + end subroutine psi_egthzmv subroutine psi_egthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ implicit none integer(psb_ipk_) :: n, k, idx(:) integer(psb_epk_) :: x(:,:), y(:,:) - + end subroutine psi_egthzmm subroutine psi_egthzv(n,idx,x,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) integer(psb_epk_) :: x(:), y(:) end subroutine psi_egthzv @@ -124,7 +134,7 @@ module psi_e_serial_mod subroutine psi_esctv(n,idx,x,beta,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ implicit none - + integer(psb_ipk_) :: n, idx(:) integer(psb_epk_) :: beta, x(:), y(:) end subroutine psi_esctv diff --git a/base/modules/auxil/psi_m_serial_mod.f90 b/base/modules/auxil/psi_m_serial_mod.f90 index e9575fc7..2acb7482 100644 --- a/base/modules/auxil/psi_m_serial_mod.f90 +++ b/base/modules/auxil/psi_m_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_m_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_ @@ -36,7 +36,7 @@ module psi_m_serial_mod ! 2-D version subroutine psb_mgelp(trans,iperm,x,info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_mpk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_m_serial_mod end subroutine psb_mgelp subroutine psb_mgelpv(trans,iperm,x,info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_mpk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_mgelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_maxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n integer(psb_mpk_), intent (in) :: x(:,:) integer(psb_mpk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_m_serial_mod end subroutine psi_maxpby subroutine psi_maxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m integer(psb_mpk_), intent (in) :: x(:) integer(psb_mpk_), intent (inout) :: y(:) integer(psb_mpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_maxpbyv + subroutine psi_maxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ + implicit none + integer(psb_ipk_), intent(in) :: m + integer(psb_mpk_), intent (in) :: x(:) + integer(psb_mpk_), intent (in) :: y(:) + integer(psb_mpk_), intent (in) :: z(:) + integer(psb_mpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_maxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_m_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) integer(psb_mpk_) :: x(:,:), y(:) - + end subroutine psi_mgthzmv subroutine psi_mgthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ implicit none integer(psb_ipk_) :: n, k, idx(:) integer(psb_mpk_) :: x(:,:), y(:,:) - + end subroutine psi_mgthzmm subroutine psi_mgthzv(n,idx,x,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) integer(psb_mpk_) :: x(:), y(:) end subroutine psi_mgthzv @@ -124,7 +134,7 @@ module psi_m_serial_mod subroutine psi_msctv(n,idx,x,beta,y) import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_ implicit none - + integer(psb_ipk_) :: n, idx(:) integer(psb_mpk_) :: beta, x(:), y(:) end subroutine psi_msctv diff --git a/base/modules/auxil/psi_s_serial_mod.f90 b/base/modules/auxil/psi_s_serial_mod.f90 index 443b16fe..ac3dbb62 100644 --- a/base/modules/auxil/psi_s_serial_mod.f90 +++ b/base/modules/auxil/psi_s_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_s_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_, psb_spk_ @@ -36,7 +36,7 @@ module psi_s_serial_mod ! 2-D version subroutine psb_sgelp(trans,iperm,x,info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none real(psb_spk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_s_serial_mod end subroutine psb_sgelp subroutine psb_sgelpv(trans,iperm,x,info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none real(psb_spk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_sgelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_saxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n real(psb_spk_), intent (in) :: x(:,:) real(psb_spk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_s_serial_mod end subroutine psi_saxpby subroutine psi_saxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m real(psb_spk_), intent (in) :: x(:) real(psb_spk_), intent (inout) :: y(:) real(psb_spk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_saxpbyv + subroutine psi_saxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_spk_ + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_spk_), intent (in) :: x(:) + real(psb_spk_), intent (in) :: y(:) + real(psb_spk_), intent (in) :: z(:) + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_saxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_s_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) real(psb_spk_) :: x(:,:), y(:) - + end subroutine psi_sgthzmv subroutine psi_sgthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_spk_ implicit none integer(psb_ipk_) :: n, k, idx(:) real(psb_spk_) :: x(:,:), y(:,:) - + end subroutine psi_sgthzmm subroutine psi_sgthzv(n,idx,x,y) import :: psb_ipk_, psb_spk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) real(psb_spk_) :: x(:), y(:) end subroutine psi_sgthzv @@ -124,7 +134,7 @@ module psi_s_serial_mod subroutine psi_ssctv(n,idx,x,beta,y) import :: psb_ipk_, psb_spk_ implicit none - + integer(psb_ipk_) :: n, idx(:) real(psb_spk_) :: beta, x(:), y(:) end subroutine psi_ssctv diff --git a/base/modules/auxil/psi_z_serial_mod.f90 b/base/modules/auxil/psi_z_serial_mod.f90 index f0d7dd11..ee148ef2 100644 --- a/base/modules/auxil/psi_z_serial_mod.f90 +++ b/base/modules/auxil/psi_z_serial_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,8 +27,8 @@ ! 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. -! -! +! +! module psi_z_serial_mod use psb_const_mod, only : psb_ipk_, psb_lpk_, psb_mpk_, psb_epk_, psb_dpk_ @@ -36,7 +36,7 @@ module psi_z_serial_mod ! 2-D version subroutine psb_zgelp(trans,iperm,x,info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none complex(psb_dpk_), intent(inout) :: x(:,:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info @@ -44,18 +44,18 @@ module psi_z_serial_mod end subroutine psb_zgelp subroutine psb_zgelpv(trans,iperm,x,info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none complex(psb_dpk_), intent(inout) :: x(:) integer(psb_ipk_), intent(in) :: iperm(:) integer(psb_ipk_), intent(out) :: info character, intent(in) :: trans end subroutine psb_zgelpv end interface psb_gelp - - interface psb_geaxpby + + interface psb_geaxpby subroutine psi_zaxpby(m,n,alpha, x, beta, y, info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n complex(psb_dpk_), intent (in) :: x(:,:) complex(psb_dpk_), intent (inout) :: y(:,:) @@ -64,13 +64,23 @@ module psi_z_serial_mod end subroutine psi_zaxpby subroutine psi_zaxpbyv(m,alpha, x, beta, y, info) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_), intent(in) :: m complex(psb_dpk_), intent (in) :: x(:) complex(psb_dpk_), intent (inout) :: y(:) complex(psb_dpk_), intent (in) :: alpha, beta integer(psb_ipk_), intent(out) :: info end subroutine psi_zaxpbyv + subroutine psi_zaxpbyv2(m,alpha, x, beta, y, z, info) + import :: psb_ipk_, psb_dpk_ + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_dpk_), intent (in) :: x(:) + complex(psb_dpk_), intent (in) :: y(:) + complex(psb_dpk_), intent (in) :: z(:) + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + end subroutine psi_zaxpbyv2 end interface psb_geaxpby interface psi_gth @@ -91,18 +101,18 @@ module psi_z_serial_mod implicit none integer(psb_ipk_) :: n, k, idx(:) complex(psb_dpk_) :: x(:,:), y(:) - + end subroutine psi_zgthzmv subroutine psi_zgthzmm(n,k,idx,x,y) import :: psb_ipk_, psb_dpk_ implicit none integer(psb_ipk_) :: n, k, idx(:) complex(psb_dpk_) :: x(:,:), y(:,:) - + end subroutine psi_zgthzmm subroutine psi_zgthzv(n,idx,x,y) import :: psb_ipk_, psb_dpk_ - implicit none + implicit none integer(psb_ipk_) :: n, idx(:) complex(psb_dpk_) :: x(:), y(:) end subroutine psi_zgthzv @@ -124,7 +134,7 @@ module psi_z_serial_mod subroutine psi_zsctv(n,idx,x,beta,y) import :: psb_ipk_, psb_dpk_ implicit none - + integer(psb_ipk_) :: n, idx(:) complex(psb_dpk_) :: beta, x(:), y(:) end subroutine psi_zsctv diff --git a/base/modules/psblas/psb_c_psblas_mod.F90 b/base/modules/psblas/psb_c_psblas_mod.F90 index d20e1a48..632bc043 100644 --- a/base/modules/psblas/psb_c_psblas_mod.F90 +++ b/base/modules/psblas/psb_c_psblas_mod.F90 @@ -98,6 +98,17 @@ module psb_c_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_caxpby_vect + subroutine psb_caxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_c_vect_type, psb_cspmat_type + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + type(psb_c_vect_type), intent (inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_caxpby_vect_out subroutine psb_caxpbyv(alpha, x, beta, y,& & desc_a, info) import :: psb_desc_type, psb_spk_, psb_ipk_, & @@ -108,6 +119,17 @@ module psb_c_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_caxpbyv + subroutine psb_caxpbyvout(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_c_vect_type, psb_cspmat_type + complex(psb_spk_), intent (in) :: x(:) + complex(psb_spk_), intent (in) :: y(:) + complex(psb_spk_), intent (inout) :: z(:) + complex(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_caxpbyvout subroutine psb_caxpby(alpha, x, beta, y,& & desc_a, info, n, jx, jy) import :: psb_desc_type, psb_spk_, psb_ipk_, & diff --git a/base/modules/psblas/psb_d_psblas_mod.F90 b/base/modules/psblas/psb_d_psblas_mod.F90 index fe0d7658..aafe6931 100644 --- a/base/modules/psblas/psb_d_psblas_mod.F90 +++ b/base/modules/psblas/psb_d_psblas_mod.F90 @@ -98,6 +98,17 @@ module psb_d_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_daxpby_vect + subroutine psb_daxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_vect_type, psb_dspmat_type + type(psb_d_vect_type), intent (inout) :: x + type(psb_d_vect_type), intent (inout) :: y + type(psb_d_vect_type), intent (inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_daxpby_vect_out subroutine psb_daxpbyv(alpha, x, beta, y,& & desc_a, info) import :: psb_desc_type, psb_dpk_, psb_ipk_, & @@ -108,6 +119,17 @@ module psb_d_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_daxpbyv + subroutine psb_daxpbyvout(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_d_vect_type, psb_dspmat_type + real(psb_dpk_), intent (in) :: x(:) + real(psb_dpk_), intent (in) :: y(:) + real(psb_dpk_), intent (inout) :: z(:) + real(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_daxpbyvout subroutine psb_daxpby(alpha, x, beta, y,& & desc_a, info, n, jx, jy) import :: psb_desc_type, psb_dpk_, psb_ipk_, & diff --git a/base/modules/psblas/psb_s_psblas_mod.F90 b/base/modules/psblas/psb_s_psblas_mod.F90 index 4544e011..48ccc212 100644 --- a/base/modules/psblas/psb_s_psblas_mod.F90 +++ b/base/modules/psblas/psb_s_psblas_mod.F90 @@ -98,6 +98,17 @@ module psb_s_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_saxpby_vect + subroutine psb_saxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_s_vect_type, psb_sspmat_type + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + type(psb_s_vect_type), intent (inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_saxpby_vect_out subroutine psb_saxpbyv(alpha, x, beta, y,& & desc_a, info) import :: psb_desc_type, psb_spk_, psb_ipk_, & @@ -108,6 +119,17 @@ module psb_s_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_saxpbyv + subroutine psb_saxpbyvout(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_spk_, psb_ipk_, & + & psb_s_vect_type, psb_sspmat_type + real(psb_spk_), intent (in) :: x(:) + real(psb_spk_), intent (in) :: y(:) + real(psb_spk_), intent (inout) :: z(:) + real(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_saxpbyvout subroutine psb_saxpby(alpha, x, beta, y,& & desc_a, info, n, jx, jy) import :: psb_desc_type, psb_spk_, psb_ipk_, & diff --git a/base/modules/psblas/psb_z_psblas_mod.F90 b/base/modules/psblas/psb_z_psblas_mod.F90 index 94d31e00..3e920f83 100644 --- a/base/modules/psblas/psb_z_psblas_mod.F90 +++ b/base/modules/psblas/psb_z_psblas_mod.F90 @@ -98,6 +98,17 @@ module psb_z_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_zaxpby_vect + subroutine psb_zaxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_z_vect_type, psb_zspmat_type + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + type(psb_z_vect_type), intent (inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_zaxpby_vect_out subroutine psb_zaxpbyv(alpha, x, beta, y,& & desc_a, info) import :: psb_desc_type, psb_dpk_, psb_ipk_, & @@ -108,6 +119,17 @@ module psb_z_psblas_mod type(psb_desc_type), intent (in) :: desc_a integer(psb_ipk_), intent(out) :: info end subroutine psb_zaxpbyv + subroutine psb_zaxpbyvout(alpha, x, beta, y,& + & z, desc_a, info) + import :: psb_desc_type, psb_dpk_, psb_ipk_, & + & psb_z_vect_type, psb_zspmat_type + complex(psb_dpk_), intent (in) :: x(:) + complex(psb_dpk_), intent (in) :: y(:) + complex(psb_dpk_), intent (inout) :: z(:) + complex(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + end subroutine psb_zaxpbyvout subroutine psb_zaxpby(alpha, x, beta, y,& & desc_a, info, n, jx, jy) import :: psb_desc_type, psb_dpk_, psb_ipk_, & diff --git a/base/modules/serial/psb_c_base_vect_mod.f90 b/base/modules/serial/psb_c_base_vect_mod.f90 index 32f4e22e..887b9e3b 100644 --- a/base/modules/serial/psb_c_base_vect_mod.f90 +++ b/base/modules/serial/psb_c_base_vect_mod.f90 @@ -151,7 +151,9 @@ module psb_c_base_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => c_base_axpby_v procedure, pass(y) :: axpby_a => c_base_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => c_base_axpby_v2 + procedure, pass(z) :: axpby_a2 => c_base_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 ! ! Vector by vector multiplication. Need all variants ! to handle multiple requirements from preconditioners @@ -974,6 +976,38 @@ contains end subroutine c_base_axpby_v + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + ! + !> Function base_axpby_v2 + !! \memberof psb_c_base_vect_type + !! \brief AXPBY by a (base_vect) z=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x The class(base_vect) to be added + !! \param beta scalar alpha + !! \param y The class(base_vect) to be added + !! \param z The class(base_vect) to be returned + !! \param info return code + !! + subroutine c_base_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_c_base_vect_type), intent(inout) :: x + class(psb_c_base_vect_type), intent(inout) :: y + class(psb_c_base_vect_type), intent(inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (x%is_dev()) call x%sync() + + call z%axpby(m,alpha,x%v,beta,y%v,info) + + end subroutine c_base_axpby_v2 + ! ! AXPBY is invoked via Y, hence the structure below. ! @@ -1002,6 +1036,36 @@ contains end subroutine c_base_axpby_a + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + !> Function base_axpby_a2 + !! \memberof psb_c_base_vect_type + !! \brief AXPBY by a normal array y=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x(:) The array to be added + !! \param beta scalar beta + !! \param y(:) The array to be added + !! \param info return code + !! + subroutine c_base_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_spk_), intent(in) :: x(:) + complex(psb_spk_), intent(in) :: y(:) + class(psb_c_base_vect_type), intent(inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (z%is_dev()) call z%sync() + call psb_geaxpby(m,alpha,x,beta,y,z%v,info) + call z%set_host() + + end subroutine c_base_axpby_a2 + ! ! Multiple variants of two operations: diff --git a/base/modules/serial/psb_c_vect_mod.F90 b/base/modules/serial/psb_c_vect_mod.F90 index afa2f0a9..cbfdc116 100644 --- a/base/modules/serial/psb_c_vect_mod.F90 +++ b/base/modules/serial/psb_c_vect_mod.F90 @@ -85,7 +85,9 @@ module psb_c_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => c_vect_axpby_v procedure, pass(y) :: axpby_a => c_vect_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => c_vect_axpby_v2 + procedure, pass(z) :: axpby_a2 => c_vect_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 procedure, pass(y) :: mlt_v => c_vect_mlt_v procedure, pass(y) :: mlt_a => c_vect_mlt_a procedure, pass(z) :: mlt_a_2 => c_vect_mlt_a_2 @@ -640,6 +642,24 @@ contains end subroutine c_vect_axpby_v + subroutine c_vect_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_c_vect_type), intent(inout) :: x + class(psb_c_vect_type), intent(inout) :: y + class(psb_c_vect_type), intent(inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v).and.allocated(y%v)) then + call z%v%axpby(m,alpha,x%v,beta,y%v,info) + else + info = psb_err_invalid_vect_state_ + end if + + end subroutine c_vect_axpby_v2 + subroutine c_vect_axpby_a(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -654,6 +674,20 @@ contains end subroutine c_vect_axpby_a + subroutine c_vect_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_spk_), intent(in) :: x(:) + complex(psb_spk_), intent(in) :: y(:) + class(psb_c_vect_type), intent(inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(z%v)) & + & call z%v%axpby(m,alpha,x,beta,y,info) + + end subroutine c_vect_axpby_a2 subroutine c_vect_mlt_v(x, y, info) use psi_serial_mod diff --git a/base/modules/serial/psb_d_base_vect_mod.f90 b/base/modules/serial/psb_d_base_vect_mod.f90 index fbdcc0cc..0b6ea129 100644 --- a/base/modules/serial/psb_d_base_vect_mod.f90 +++ b/base/modules/serial/psb_d_base_vect_mod.f90 @@ -151,7 +151,9 @@ module psb_d_base_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => d_base_axpby_v procedure, pass(y) :: axpby_a => d_base_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => d_base_axpby_v2 + procedure, pass(z) :: axpby_a2 => d_base_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 ! ! Vector by vector multiplication. Need all variants ! to handle multiple requirements from preconditioners @@ -981,6 +983,38 @@ contains end subroutine d_base_axpby_v + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + ! + !> Function base_axpby_v2 + !! \memberof psb_d_base_vect_type + !! \brief AXPBY by a (base_vect) z=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x The class(base_vect) to be added + !! \param beta scalar alpha + !! \param y The class(base_vect) to be added + !! \param z The class(base_vect) to be returned + !! \param info return code + !! + subroutine d_base_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_d_base_vect_type), intent(inout) :: x + class(psb_d_base_vect_type), intent(inout) :: y + class(psb_d_base_vect_type), intent(inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (x%is_dev()) call x%sync() + + call z%axpby(m,alpha,x%v,beta,y%v,info) + + end subroutine d_base_axpby_v2 + ! ! AXPBY is invoked via Y, hence the structure below. ! @@ -1009,6 +1043,36 @@ contains end subroutine d_base_axpby_a + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + !> Function base_axpby_a2 + !! \memberof psb_d_base_vect_type + !! \brief AXPBY by a normal array y=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x(:) The array to be added + !! \param beta scalar beta + !! \param y(:) The array to be added + !! \param info return code + !! + subroutine d_base_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(in) :: y(:) + class(psb_d_base_vect_type), intent(inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (z%is_dev()) call z%sync() + call psb_geaxpby(m,alpha,x,beta,y,z%v,info) + call z%set_host() + + end subroutine d_base_axpby_a2 + ! ! Multiple variants of two operations: diff --git a/base/modules/serial/psb_d_vect_mod.F90 b/base/modules/serial/psb_d_vect_mod.F90 index f49db541..edd99bf5 100644 --- a/base/modules/serial/psb_d_vect_mod.F90 +++ b/base/modules/serial/psb_d_vect_mod.F90 @@ -85,7 +85,9 @@ module psb_d_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => d_vect_axpby_v procedure, pass(y) :: axpby_a => d_vect_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => d_vect_axpby_v2 + procedure, pass(z) :: axpby_a2 => d_vect_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 procedure, pass(y) :: mlt_v => d_vect_mlt_v procedure, pass(y) :: mlt_a => d_vect_mlt_a procedure, pass(z) :: mlt_a_2 => d_vect_mlt_a_2 @@ -647,6 +649,24 @@ contains end subroutine d_vect_axpby_v + subroutine d_vect_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_d_vect_type), intent(inout) :: x + class(psb_d_vect_type), intent(inout) :: y + class(psb_d_vect_type), intent(inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v).and.allocated(y%v)) then + call z%v%axpby(m,alpha,x%v,beta,y%v,info) + else + info = psb_err_invalid_vect_state_ + end if + + end subroutine d_vect_axpby_v2 + subroutine d_vect_axpby_a(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -661,6 +681,20 @@ contains end subroutine d_vect_axpby_a + subroutine d_vect_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(in) :: y(:) + class(psb_d_vect_type), intent(inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(z%v)) & + & call z%v%axpby(m,alpha,x,beta,y,info) + + end subroutine d_vect_axpby_a2 subroutine d_vect_mlt_v(x, y, info) use psi_serial_mod diff --git a/base/modules/serial/psb_s_base_vect_mod.f90 b/base/modules/serial/psb_s_base_vect_mod.f90 index 451baef5..1a951177 100644 --- a/base/modules/serial/psb_s_base_vect_mod.f90 +++ b/base/modules/serial/psb_s_base_vect_mod.f90 @@ -151,7 +151,9 @@ module psb_s_base_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => s_base_axpby_v procedure, pass(y) :: axpby_a => s_base_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => s_base_axpby_v2 + procedure, pass(z) :: axpby_a2 => s_base_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 ! ! Vector by vector multiplication. Need all variants ! to handle multiple requirements from preconditioners @@ -981,6 +983,38 @@ contains end subroutine s_base_axpby_v + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + ! + !> Function base_axpby_v2 + !! \memberof psb_s_base_vect_type + !! \brief AXPBY by a (base_vect) z=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x The class(base_vect) to be added + !! \param beta scalar alpha + !! \param y The class(base_vect) to be added + !! \param z The class(base_vect) to be returned + !! \param info return code + !! + subroutine s_base_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_s_base_vect_type), intent(inout) :: x + class(psb_s_base_vect_type), intent(inout) :: y + class(psb_s_base_vect_type), intent(inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (x%is_dev()) call x%sync() + + call z%axpby(m,alpha,x%v,beta,y%v,info) + + end subroutine s_base_axpby_v2 + ! ! AXPBY is invoked via Y, hence the structure below. ! @@ -1009,6 +1043,36 @@ contains end subroutine s_base_axpby_a + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + !> Function base_axpby_a2 + !! \memberof psb_s_base_vect_type + !! \brief AXPBY by a normal array y=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x(:) The array to be added + !! \param beta scalar beta + !! \param y(:) The array to be added + !! \param info return code + !! + subroutine s_base_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(in) :: y(:) + class(psb_s_base_vect_type), intent(inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (z%is_dev()) call z%sync() + call psb_geaxpby(m,alpha,x,beta,y,z%v,info) + call z%set_host() + + end subroutine s_base_axpby_a2 + ! ! Multiple variants of two operations: diff --git a/base/modules/serial/psb_s_vect_mod.F90 b/base/modules/serial/psb_s_vect_mod.F90 index aad99c9c..edeedacb 100644 --- a/base/modules/serial/psb_s_vect_mod.F90 +++ b/base/modules/serial/psb_s_vect_mod.F90 @@ -85,7 +85,9 @@ module psb_s_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => s_vect_axpby_v procedure, pass(y) :: axpby_a => s_vect_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => s_vect_axpby_v2 + procedure, pass(z) :: axpby_a2 => s_vect_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 procedure, pass(y) :: mlt_v => s_vect_mlt_v procedure, pass(y) :: mlt_a => s_vect_mlt_a procedure, pass(z) :: mlt_a_2 => s_vect_mlt_a_2 @@ -647,6 +649,24 @@ contains end subroutine s_vect_axpby_v + subroutine s_vect_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_s_vect_type), intent(inout) :: x + class(psb_s_vect_type), intent(inout) :: y + class(psb_s_vect_type), intent(inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v).and.allocated(y%v)) then + call z%v%axpby(m,alpha,x%v,beta,y%v,info) + else + info = psb_err_invalid_vect_state_ + end if + + end subroutine s_vect_axpby_v2 + subroutine s_vect_axpby_a(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -661,6 +681,20 @@ contains end subroutine s_vect_axpby_a + subroutine s_vect_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(in) :: y(:) + class(psb_s_vect_type), intent(inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(z%v)) & + & call z%v%axpby(m,alpha,x,beta,y,info) + + end subroutine s_vect_axpby_a2 subroutine s_vect_mlt_v(x, y, info) use psi_serial_mod diff --git a/base/modules/serial/psb_z_base_vect_mod.f90 b/base/modules/serial/psb_z_base_vect_mod.f90 index 502b7795..a7e6978a 100644 --- a/base/modules/serial/psb_z_base_vect_mod.f90 +++ b/base/modules/serial/psb_z_base_vect_mod.f90 @@ -151,7 +151,9 @@ module psb_z_base_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => z_base_axpby_v procedure, pass(y) :: axpby_a => z_base_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => z_base_axpby_v2 + procedure, pass(z) :: axpby_a2 => z_base_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 ! ! Vector by vector multiplication. Need all variants ! to handle multiple requirements from preconditioners @@ -974,6 +976,38 @@ contains end subroutine z_base_axpby_v + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + ! + !> Function base_axpby_v2 + !! \memberof psb_z_base_vect_type + !! \brief AXPBY by a (base_vect) z=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x The class(base_vect) to be added + !! \param beta scalar alpha + !! \param y The class(base_vect) to be added + !! \param z The class(base_vect) to be returned + !! \param info return code + !! + subroutine z_base_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_z_base_vect_type), intent(inout) :: x + class(psb_z_base_vect_type), intent(inout) :: y + class(psb_z_base_vect_type), intent(inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (x%is_dev()) call x%sync() + + call z%axpby(m,alpha,x%v,beta,y%v,info) + + end subroutine z_base_axpby_v2 + ! ! AXPBY is invoked via Y, hence the structure below. ! @@ -1002,6 +1036,36 @@ contains end subroutine z_base_axpby_a + ! + ! AXPBY is invoked via Z, hence the structure below. + ! + ! + !> Function base_axpby_a2 + !! \memberof psb_z_base_vect_type + !! \brief AXPBY by a normal array y=alpha*x+beta*y + !! \param m Number of entries to be considered + !! \param alpha scalar alpha + !! \param x(:) The array to be added + !! \param beta scalar beta + !! \param y(:) The array to be added + !! \param info return code + !! + subroutine z_base_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_dpk_), intent(in) :: x(:) + complex(psb_dpk_), intent(in) :: y(:) + class(psb_z_base_vect_type), intent(inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (z%is_dev()) call z%sync() + call psb_geaxpby(m,alpha,x,beta,y,z%v,info) + call z%set_host() + + end subroutine z_base_axpby_a2 + ! ! Multiple variants of two operations: diff --git a/base/modules/serial/psb_z_vect_mod.F90 b/base/modules/serial/psb_z_vect_mod.F90 index 7b01aba0..77aba066 100644 --- a/base/modules/serial/psb_z_vect_mod.F90 +++ b/base/modules/serial/psb_z_vect_mod.F90 @@ -85,7 +85,9 @@ module psb_z_vect_mod generic, public :: dot => dot_v, dot_a procedure, pass(y) :: axpby_v => z_vect_axpby_v procedure, pass(y) :: axpby_a => z_vect_axpby_a - generic, public :: axpby => axpby_v, axpby_a + procedure, pass(z) :: axpby_v2 => z_vect_axpby_v2 + procedure, pass(z) :: axpby_a2 => z_vect_axpby_a2 + generic, public :: axpby => axpby_v, axpby_a, axpby_v2, axpby_a2 procedure, pass(y) :: mlt_v => z_vect_mlt_v procedure, pass(y) :: mlt_a => z_vect_mlt_a procedure, pass(z) :: mlt_a_2 => z_vect_mlt_a_2 @@ -640,6 +642,24 @@ contains end subroutine z_vect_axpby_v + subroutine z_vect_axpby_v2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + class(psb_z_vect_type), intent(inout) :: x + class(psb_z_vect_type), intent(inout) :: y + class(psb_z_vect_type), intent(inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(x%v).and.allocated(y%v)) then + call z%v%axpby(m,alpha,x%v,beta,y%v,info) + else + info = psb_err_invalid_vect_state_ + end if + + end subroutine z_vect_axpby_v2 + subroutine z_vect_axpby_a(m,alpha, x, beta, y, info) use psi_serial_mod implicit none @@ -654,6 +674,20 @@ contains end subroutine z_vect_axpby_a + subroutine z_vect_axpby_a2(m,alpha, x, beta, y, z, info) + use psi_serial_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_dpk_), intent(in) :: x(:) + complex(psb_dpk_), intent(in) :: y(:) + class(psb_z_vect_type), intent(inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + + if (allocated(z%v)) & + & call z%v%axpby(m,alpha,x,beta,y,info) + + end subroutine z_vect_axpby_a2 subroutine z_vect_mlt_v(x, y, info) use psi_serial_mod diff --git a/base/psblas/psb_caxpby.f90 b/base/psblas/psb_caxpby.f90 index 280022b1..b92072ab 100644 --- a/base/psblas/psb_caxpby.f90 +++ b/base/psblas/psb_caxpby.f90 @@ -129,6 +129,152 @@ subroutine psb_caxpby_vect(alpha, x, beta, y,& end subroutine psb_caxpby_vect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! +! +! File: psb_caxpby.f90 + +! +! Subroutine: psb_caxpby_vect_out +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - complex,input The scalar used to multiply each component of X +! x - type(psb_c_vect_type) The input vector containing the entries of X +! beta - complex,input The scalar used to multiply each component of Y +! y - type(psb_c_vect_type) The input vector Y +! z - type(psb_c_vect_type) The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! Note: from a functional point of view, X is input, but here +! it's declared INOUT because of the sync() methods. +! +subroutine psb_caxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + use psb_base_mod, psb_protect_name => psb_caxpby_vect_out + implicit none + type(psb_c_vect_type), intent (inout) :: x + type(psb_c_vect_type), intent (inout) :: y + type(psb_c_vect_type), intent (inout) :: z + complex(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + + name='psb_cgeaxpby' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,y%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 3' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call z%axpby(desc_a%get_local_rows(),& + & alpha,x,beta,z,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_caxpby_vect_out + ! ! Subroutine: psb_caxpby ! Adds one distributed matrix to another, @@ -372,6 +518,138 @@ subroutine psb_caxpbyv(alpha, x, beta,y,desc_a,info) return end subroutine psb_caxpbyv + +!!$ +!!$ Parallel Sparse BLAS version 3.5 +!!$ (C) Copyright 2006-2018 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari +!!$ +!!$ 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. +!!$ +!!$ +! +! Subroutine: psb_caxpbyvout +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - complex,input The scalar used to multiply each component of X +! x(:) - complex,input The input vector containing the entries of X +! beta - complex,input The scalar used to multiply each component of Y +! y(:) - complex,input The input vector Y containing the entries of Y +! Z(:) - complex,inout The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! +subroutine psb_caxpbyvout(alpha, x, beta,y, z, desc_a,info) + use psb_base_mod, psb_protect_name => psb_caxpbyvout + implicit none + + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(in) :: desc_a + complex(psb_spk_), intent(in) :: alpha, beta + complex(psb_spk_), intent(in) :: x(:) + complex(psb_spk_), intent(in) :: y(:) + complex(psb_spk_), intent(inout) :: z(:) + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz, lldx, lldy, lldz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + lldx = size(x,1) + lldy = size(y,1) + lldz = size(z,1) + ! check vector correctness + call psb_chkvect(m,lone,lldx,ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldy,iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldz,iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call caxpby(desc_a%get_local_cols(),ione,& + & alpha,x,lldx,beta,& + & y,lldy,z,lldz,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end subroutine psb_caxpbyvout + ! ! Subroutine: psb_caddconst_vect ! Adds one distributed vector to another, diff --git a/base/psblas/psb_daxpby.f90 b/base/psblas/psb_daxpby.f90 index 4bdb24a8..4ee12d46 100644 --- a/base/psblas/psb_daxpby.f90 +++ b/base/psblas/psb_daxpby.f90 @@ -129,6 +129,152 @@ subroutine psb_daxpby_vect(alpha, x, beta, y,& end subroutine psb_daxpby_vect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! +! +! File: psb_daxpby.f90 + +! +! Subroutine: psb_daxpby_vect_out +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - real,input The scalar used to multiply each component of X +! x - type(psb_d_vect_type) The input vector containing the entries of X +! beta - real,input The scalar used to multiply each component of Y +! y - type(psb_d_vect_type) The input vector Y +! z - type(psb_d_vect_type) The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! Note: from a functional point of view, X is input, but here +! it's declared INOUT because of the sync() methods. +! +subroutine psb_daxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + use psb_base_mod, psb_protect_name => psb_daxpby_vect_out + implicit none + type(psb_d_vect_type), intent (inout) :: x + type(psb_d_vect_type), intent (inout) :: y + type(psb_d_vect_type), intent (inout) :: z + real(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + + name='psb_dgeaxpby' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,y%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 3' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call z%axpby(desc_a%get_local_rows(),& + & alpha,x,beta,z,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_daxpby_vect_out + ! ! Subroutine: psb_daxpby ! Adds one distributed matrix to another, @@ -372,6 +518,138 @@ subroutine psb_daxpbyv(alpha, x, beta,y,desc_a,info) return end subroutine psb_daxpbyv + +!!$ +!!$ Parallel Sparse BLAS version 3.5 +!!$ (C) Copyright 2006-2018 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari +!!$ +!!$ 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. +!!$ +!!$ +! +! Subroutine: psb_daxpbyvout +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - real,input The scalar used to multiply each component of X +! x(:) - real,input The input vector containing the entries of X +! beta - real,input The scalar used to multiply each component of Y +! y(:) - real,input The input vector Y containing the entries of Y +! Z(:) - real,inout The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! +subroutine psb_daxpbyvout(alpha, x, beta,y, z, desc_a,info) + use psb_base_mod, psb_protect_name => psb_daxpbyvout + implicit none + + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: alpha, beta + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(in) :: y(:) + real(psb_dpk_), intent(inout) :: z(:) + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz, lldx, lldy, lldz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + lldx = size(x,1) + lldy = size(y,1) + lldz = size(z,1) + ! check vector correctness + call psb_chkvect(m,lone,lldx,ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldy,iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldz,iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call daxpby(desc_a%get_local_cols(),ione,& + & alpha,x,lldx,beta,& + & y,lldy,z,lldz,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end subroutine psb_daxpbyvout + ! ! Subroutine: psb_daddconst_vect ! Adds one distributed vector to another, diff --git a/base/psblas/psb_saxpby.f90 b/base/psblas/psb_saxpby.f90 index 2d3be43e..0c619984 100644 --- a/base/psblas/psb_saxpby.f90 +++ b/base/psblas/psb_saxpby.f90 @@ -129,6 +129,152 @@ subroutine psb_saxpby_vect(alpha, x, beta, y,& end subroutine psb_saxpby_vect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! +! +! File: psb_saxpby.f90 + +! +! Subroutine: psb_saxpby_vect_out +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - real,input The scalar used to multiply each component of X +! x - type(psb_s_vect_type) The input vector containing the entries of X +! beta - real,input The scalar used to multiply each component of Y +! y - type(psb_s_vect_type) The input vector Y +! z - type(psb_s_vect_type) The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! Note: from a functional point of view, X is input, but here +! it's declared INOUT because of the sync() methods. +! +subroutine psb_saxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + use psb_base_mod, psb_protect_name => psb_saxpby_vect_out + implicit none + type(psb_s_vect_type), intent (inout) :: x + type(psb_s_vect_type), intent (inout) :: y + type(psb_s_vect_type), intent (inout) :: z + real(psb_spk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + + name='psb_sgeaxpby' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,y%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 3' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call z%axpby(desc_a%get_local_rows(),& + & alpha,x,beta,z,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_saxpby_vect_out + ! ! Subroutine: psb_saxpby ! Adds one distributed matrix to another, @@ -372,6 +518,138 @@ subroutine psb_saxpbyv(alpha, x, beta,y,desc_a,info) return end subroutine psb_saxpbyv + +!!$ +!!$ Parallel Sparse BLAS version 3.5 +!!$ (C) Copyright 2006-2018 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari +!!$ +!!$ 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. +!!$ +!!$ +! +! Subroutine: psb_saxpbyvout +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - real,input The scalar used to multiply each component of X +! x(:) - real,input The input vector containing the entries of X +! beta - real,input The scalar used to multiply each component of Y +! y(:) - real,input The input vector Y containing the entries of Y +! Z(:) - real,inout The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! +subroutine psb_saxpbyvout(alpha, x, beta,y, z, desc_a,info) + use psb_base_mod, psb_protect_name => psb_saxpbyvout + implicit none + + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: alpha, beta + real(psb_spk_), intent(in) :: x(:) + real(psb_spk_), intent(in) :: y(:) + real(psb_spk_), intent(inout) :: z(:) + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz, lldx, lldy, lldz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + lldx = size(x,1) + lldy = size(y,1) + lldz = size(z,1) + ! check vector correctness + call psb_chkvect(m,lone,lldx,ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldy,iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldz,iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call saxpby(desc_a%get_local_cols(),ione,& + & alpha,x,lldx,beta,& + & y,lldy,z,lldz,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end subroutine psb_saxpbyvout + ! ! Subroutine: psb_saddconst_vect ! Adds one distributed vector to another, diff --git a/base/psblas/psb_zaxpby.f90 b/base/psblas/psb_zaxpby.f90 index b2cd7a0e..f60e8739 100644 --- a/base/psblas/psb_zaxpby.f90 +++ b/base/psblas/psb_zaxpby.f90 @@ -129,6 +129,152 @@ subroutine psb_zaxpby_vect(alpha, x, beta, y,& end subroutine psb_zaxpby_vect +! +! Parallel Sparse BLAS version 3.5 +! (C) Copyright 2006-2018 +! Salvatore Filippone +! Alfredo Buttari +! +! 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. +! +! +! File: psb_zaxpby.f90 + +! +! Subroutine: psb_zaxpby_vect_out +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - complex,input The scalar used to multiply each component of X +! x - type(psb_z_vect_type) The input vector containing the entries of X +! beta - complex,input The scalar used to multiply each component of Y +! y - type(psb_z_vect_type) The input vector Y +! z - type(psb_z_vect_type) The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! Note: from a functional point of view, X is input, but here +! it's declared INOUT because of the sync() methods. +! +subroutine psb_zaxpby_vect_out(alpha, x, beta, y,& + & z, desc_a, info) + use psb_base_mod, psb_protect_name => psb_zaxpby_vect_out + implicit none + type(psb_z_vect_type), intent (inout) :: x + type(psb_z_vect_type), intent (inout) :: y + type(psb_z_vect_type), intent (inout) :: z + complex(psb_dpk_), intent (in) :: alpha, beta + type(psb_desc_type), intent (in) :: desc_a + integer(psb_ipk_), intent(out) :: info + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + + name='psb_zgeaxpby' + if (psb_errstatus_fatal()) return + info=psb_success_ + call psb_erractionsave(err_act) + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(x%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(y%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + if (.not.allocated(z%v)) then + info = psb_err_invalid_vect_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + + ! check vector correctness + call psb_chkvect(m,lone,x%get_nrows(),ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,y%get_nrows(),iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,z%get_nrows(),iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 3' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call z%axpby(desc_a%get_local_rows(),& + & alpha,x,beta,z,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return + +end subroutine psb_zaxpby_vect_out + ! ! Subroutine: psb_zaxpby ! Adds one distributed matrix to another, @@ -372,6 +518,138 @@ subroutine psb_zaxpbyv(alpha, x, beta,y,desc_a,info) return end subroutine psb_zaxpbyv + +!!$ +!!$ Parallel Sparse BLAS version 3.5 +!!$ (C) Copyright 2006-2018 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari +!!$ +!!$ 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. +!!$ +!!$ +! +! Subroutine: psb_zaxpbyvout +! Adds one distributed vector to another, +! +! Z := beta * Y + alpha * X +! +! Arguments: +! alpha - complex,input The scalar used to multiply each component of X +! x(:) - complex,input The input vector containing the entries of X +! beta - complex,input The scalar used to multiply each component of Y +! y(:) - complex,input The input vector Y containing the entries of Y +! Z(:) - complex,inout The output vector Z +! desc_a - type(psb_desc_type) The communication descriptor. +! info - integer Return code +! +! +subroutine psb_zaxpbyvout(alpha, x, beta,y, z, desc_a,info) + use psb_base_mod, psb_protect_name => psb_zaxpbyvout + implicit none + + integer(psb_ipk_), intent(out) :: info + type(psb_desc_type), intent(in) :: desc_a + complex(psb_dpk_), intent(in) :: alpha, beta + complex(psb_dpk_), intent(in) :: x(:) + complex(psb_dpk_), intent(in) :: y(:) + complex(psb_dpk_), intent(inout) :: z(:) + + ! locals + integer(psb_ipk_) :: ictxt, np, me,& + & err_act, iix, jjx, iiy, jjy, iiz, jjz, lldx, lldy, lldz + integer(psb_lpk_) :: ix, ijx, iy, ijy, iz, ijz, m + character(len=20) :: name, ch_err + logical, parameter :: debug=.false. + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + ictxt=desc_a%get_context() + + call psb_info(ictxt, me, np) + if (np == -ione) then + info = psb_err_context_error_ + call psb_errpush(info,name) + goto 9999 + endif + + ix = ione + iy = ione + iz = ione + + m = desc_a%get_global_rows() + lldx = size(x,1) + lldy = size(y,1) + lldz = size(z,1) + ! check vector correctness + call psb_chkvect(m,lone,lldx,ix,lone,desc_a,info,iix,jjx) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 1' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldy,iy,lone,desc_a,info,iiy,jjy) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + call psb_chkvect(m,lone,lldz,iz,lone,desc_a,info,iiz,jjz) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_chkvect 2' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + if ((iix /= ione).or.(iiy /= ione).or.(iiz /= ione)) then + info=psb_err_ix_n1_iy_n1_unsupported_ + call psb_errpush(info,name) + end if + + if(desc_a%get_local_rows() > 0) then + call zaxpby(desc_a%get_local_cols(),ione,& + & alpha,x,lldx,beta,& + & y,lldy,z,lldz,info) + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(ictxt,err_act) + + return +end subroutine psb_zaxpbyvout + ! ! Subroutine: psb_zaddconst_vect ! Adds one distributed vector to another, diff --git a/base/serial/psi_c_serial_impl.f90 b/base/serial/psi_c_serial_impl.f90 index 77841b76..2120683d 100644 --- a/base/serial/psi_c_serial_impl.f90 +++ b/base/serial/psi_c_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_caxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n complex(psb_spk_), intent (in) :: x(:,:) complex(psb_spk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_caxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call caxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_caxpby(m,n,alpha, x, beta, y, info) end subroutine psi_caxpby subroutine psi_caxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m complex(psb_spk_), intent (in) :: x(:) complex(psb_spk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_caxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call caxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_caxpbyv(m,alpha, x, beta, y, info) end subroutine psi_caxpbyv +subroutine psi_caxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_spk_), intent (in) :: x(:) + complex(psb_spk_), intent (in) :: y(:) + complex(psb_spk_), intent (inout) :: z(:) + complex(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call caxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_caxpbyv2 subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == czero) then - if (alpha == czero) then + if (beta == czero) then + if (alpha == czero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -cone) then + else if (alpha == -cone) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == cone) then + else + if (beta == cone) then ! Do nothing - else if (beta == -cone) then - y(1:n*k) = -y(1:n*k) + else if (beta == -cone) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == czero) then + if (alpha == czero) then ! do nothing - else if (alpha == cone) then + else if (alpha == cone) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_cgthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == czero) then - if (alpha == czero) then + if (beta == czero) then + if (alpha == czero) then do i=1,n y(i) = czero end do - else if (alpha == cone) then + else if (alpha == cone) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -cone) then + else if (alpha == -cone) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == cone) then + else + if (beta == cone) then ! Do nothing - else if (beta == -cone) then - y(1:n) = -y(1:n) + else if (beta == -cone) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == czero) then + if (alpha == czero) then ! do nothing - else if (alpha == cone) then + else if (alpha == cone) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -cone) then + else if (alpha == -cone) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_cgthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_cgthzv(n,idx,x,y) end subroutine psi_cgthzv subroutine psi_csctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_csctmm(n,k,idx,x,beta,y) end subroutine psi_csctmm subroutine psi_csctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_csctv subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info complex(psb_spk_) X(lldx,*), Y(lldy,*) complex(psb_spk_) alpha, beta @@ -447,19 +508,19 @@ subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.czero) then - if (beta.eq.czero) then - do j=1, n - do i=1,m + if (alpha.eq.czero) then + if (beta.eq.czero) then + do j=1, n + do i=1,m y(i,j) = czero enddo enddo else if (beta.eq.cone) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-cone) then - do j=1,n - do i=1,m + else if (beta.eq.-cone) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.cone) then - if (beta.eq.czero) then - do j=1,n - do i=1,m + if (beta.eq.czero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.cone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-cone) then - do j=1,n - do i=1,m + else if (beta.eq.-cone) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-cone) then + else if (alpha.eq.-cone) then - if (beta.eq.czero) then - do j=1,n - do i=1,m + if (beta.eq.czero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.cone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-cone) then - do j=1,n - do i=1,m + else if (beta.eq.-cone) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.czero) then - do j=1,n - do i=1,m + if (beta.eq.czero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.cone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-cone) then - do j=1,n - do i=1,m + else if (beta.eq.-cone) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine caxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine caxpby + +subroutine caxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + complex(psb_spk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + complex(psb_spk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='caxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.czero) then + if (beta.eq.czero) then + do j=1, n + do i=1,m + Z(i,j) = czero + enddo + enddo + else if (beta.eq.cone) then + ! + ! Do nothing! + ! + + else if (beta.eq.-cone) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.cone) then + + if (beta.eq.czero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.cone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-cone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-cone) then + + if (beta.eq.czero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.cone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-cone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.czero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.cone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-cone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine caxpbyv2 diff --git a/base/serial/psi_d_serial_impl.f90 b/base/serial/psi_d_serial_impl.f90 index 0e1904f1..0d80f459 100644 --- a/base/serial/psi_d_serial_impl.f90 +++ b/base/serial/psi_d_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_daxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n real(psb_dpk_), intent (in) :: x(:,:) real(psb_dpk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_daxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call daxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_daxpby(m,n,alpha, x, beta, y, info) end subroutine psi_daxpby subroutine psi_daxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m real(psb_dpk_), intent (in) :: x(:) real(psb_dpk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_daxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call daxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_daxpbyv(m,alpha, x, beta, y, info) end subroutine psi_daxpbyv +subroutine psi_daxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_dpk_), intent (in) :: x(:) + real(psb_dpk_), intent (in) :: y(:) + real(psb_dpk_), intent (inout) :: z(:) + real(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call daxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_daxpbyv2 subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == dzero) then - if (alpha == dzero) then + if (beta == dzero) then + if (alpha == dzero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -done) then + else if (alpha == -done) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == done) then + else + if (beta == done) then ! Do nothing - else if (beta == -done) then - y(1:n*k) = -y(1:n*k) + else if (beta == -done) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == dzero) then + if (alpha == dzero) then ! do nothing - else if (alpha == done) then + else if (alpha == done) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_dgthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == dzero) then - if (alpha == dzero) then + if (beta == dzero) then + if (alpha == dzero) then do i=1,n y(i) = dzero end do - else if (alpha == done) then + else if (alpha == done) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -done) then + else if (alpha == -done) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == done) then + else + if (beta == done) then ! Do nothing - else if (beta == -done) then - y(1:n) = -y(1:n) + else if (beta == -done) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == dzero) then + if (alpha == dzero) then ! do nothing - else if (alpha == done) then + else if (alpha == done) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -done) then + else if (alpha == -done) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_dgthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_dgthzv(n,idx,x,y) end subroutine psi_dgthzv subroutine psi_dsctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_dsctmm(n,k,idx,x,beta,y) end subroutine psi_dsctmm subroutine psi_dsctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_dsctv subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info real(psb_dpk_) X(lldx,*), Y(lldy,*) real(psb_dpk_) alpha, beta @@ -447,19 +508,19 @@ subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.dzero) then - if (beta.eq.dzero) then - do j=1, n - do i=1,m + if (alpha.eq.dzero) then + if (beta.eq.dzero) then + do j=1, n + do i=1,m y(i,j) = dzero enddo enddo else if (beta.eq.done) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-done) then - do j=1,n - do i=1,m + else if (beta.eq.-done) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.done) then - if (beta.eq.dzero) then - do j=1,n - do i=1,m + if (beta.eq.dzero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.done) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-done) then - do j=1,n - do i=1,m + else if (beta.eq.-done) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-done) then + else if (alpha.eq.-done) then - if (beta.eq.dzero) then - do j=1,n - do i=1,m + if (beta.eq.dzero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.done) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-done) then - do j=1,n - do i=1,m + else if (beta.eq.-done) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.dzero) then - do j=1,n - do i=1,m + if (beta.eq.dzero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.done) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-done) then - do j=1,n - do i=1,m + else if (beta.eq.-done) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine daxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine daxpby + +subroutine daxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + real(psb_dpk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + real(psb_dpk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='daxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.dzero) then + if (beta.eq.dzero) then + do j=1, n + do i=1,m + Z(i,j) = dzero + enddo + enddo + else if (beta.eq.done) then + ! + ! Do nothing! + ! + + else if (beta.eq.-done) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.done) then + + if (beta.eq.dzero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.done) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-done) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-done) then + + if (beta.eq.dzero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.done) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-done) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.dzero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.done) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-done) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine daxpbyv2 diff --git a/base/serial/psi_e_serial_impl.f90 b/base/serial/psi_e_serial_impl.f90 index 9f198651..0595d87e 100644 --- a/base/serial/psi_e_serial_impl.f90 +++ b/base/serial/psi_e_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_eaxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n integer(psb_epk_), intent (in) :: x(:,:) integer(psb_epk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_eaxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call eaxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_eaxpby(m,n,alpha, x, beta, y, info) end subroutine psi_eaxpby subroutine psi_eaxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m integer(psb_epk_), intent (in) :: x(:) integer(psb_epk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_eaxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call eaxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_eaxpbyv(m,alpha, x, beta, y, info) end subroutine psi_eaxpbyv +subroutine psi_eaxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + integer(psb_epk_), intent (in) :: x(:) + integer(psb_epk_), intent (in) :: y(:) + integer(psb_epk_), intent (inout) :: z(:) + integer(psb_epk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call eaxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_eaxpbyv2 subroutine psi_egthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_egthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == ezero) then - if (alpha == ezero) then + if (beta == ezero) then + if (alpha == ezero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_egthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -eone) then + else if (alpha == -eone) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_egthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == eone) then + else + if (beta == eone) then ! Do nothing - else if (beta == -eone) then - y(1:n*k) = -y(1:n*k) + else if (beta == -eone) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == ezero) then + if (alpha == ezero) then ! do nothing - else if (alpha == eone) then + else if (alpha == eone) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_egthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_egthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == ezero) then - if (alpha == ezero) then + if (beta == ezero) then + if (alpha == ezero) then do i=1,n y(i) = ezero end do - else if (alpha == eone) then + else if (alpha == eone) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -eone) then + else if (alpha == -eone) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == eone) then + else + if (beta == eone) then ! Do nothing - else if (beta == -eone) then - y(1:n) = -y(1:n) + else if (beta == -eone) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == ezero) then + if (alpha == ezero) then ! do nothing - else if (alpha == eone) then + else if (alpha == eone) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -eone) then + else if (alpha == -eone) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_egthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_egthzv(n,idx,x,y) end subroutine psi_egthzv subroutine psi_esctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_esctmm(n,k,idx,x,beta,y) end subroutine psi_esctmm subroutine psi_esctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_esctv subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info integer(psb_epk_) X(lldx,*), Y(lldy,*) integer(psb_epk_) alpha, beta @@ -447,19 +508,19 @@ subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.ezero) then - if (beta.eq.ezero) then - do j=1, n - do i=1,m + if (alpha.eq.ezero) then + if (beta.eq.ezero) then + do j=1, n + do i=1,m y(i,j) = ezero enddo enddo else if (beta.eq.eone) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-eone) then - do j=1,n - do i=1,m + else if (beta.eq.-eone) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.eone) then - if (beta.eq.ezero) then - do j=1,n - do i=1,m + if (beta.eq.ezero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.eone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-eone) then - do j=1,n - do i=1,m + else if (beta.eq.-eone) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-eone) then + else if (alpha.eq.-eone) then - if (beta.eq.ezero) then - do j=1,n - do i=1,m + if (beta.eq.ezero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.eone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-eone) then - do j=1,n - do i=1,m + else if (beta.eq.-eone) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.ezero) then - do j=1,n - do i=1,m + if (beta.eq.ezero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.eone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-eone) then - do j=1,n - do i=1,m + else if (beta.eq.-eone) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine eaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine eaxpby + +subroutine eaxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + integer(psb_epk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + integer(psb_epk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='eaxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.ezero) then + if (beta.eq.ezero) then + do j=1, n + do i=1,m + Z(i,j) = ezero + enddo + enddo + else if (beta.eq.eone) then + ! + ! Do nothing! + ! + + else if (beta.eq.-eone) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.eone) then + + if (beta.eq.ezero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.eone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-eone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-eone) then + + if (beta.eq.ezero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.eone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-eone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.ezero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.eone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-eone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine eaxpbyv2 diff --git a/base/serial/psi_m_serial_impl.f90 b/base/serial/psi_m_serial_impl.f90 index a885f2bd..cc8b9f4f 100644 --- a/base/serial/psi_m_serial_impl.f90 +++ b/base/serial/psi_m_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_maxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n integer(psb_mpk_), intent (in) :: x(:,:) integer(psb_mpk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_maxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call maxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_maxpby(m,n,alpha, x, beta, y, info) end subroutine psi_maxpby subroutine psi_maxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m integer(psb_mpk_), intent (in) :: x(:) integer(psb_mpk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_maxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call maxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_maxpbyv(m,alpha, x, beta, y, info) end subroutine psi_maxpbyv +subroutine psi_maxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + integer(psb_mpk_), intent (in) :: x(:) + integer(psb_mpk_), intent (in) :: y(:) + integer(psb_mpk_), intent (inout) :: z(:) + integer(psb_mpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call maxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_maxpbyv2 subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == mzero) then - if (alpha == mzero) then + if (beta == mzero) then + if (alpha == mzero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -mone) then + else if (alpha == -mone) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == mone) then + else + if (beta == mone) then ! Do nothing - else if (beta == -mone) then - y(1:n*k) = -y(1:n*k) + else if (beta == -mone) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == mzero) then + if (alpha == mzero) then ! do nothing - else if (alpha == mone) then + else if (alpha == mone) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_mgthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == mzero) then - if (alpha == mzero) then + if (beta == mzero) then + if (alpha == mzero) then do i=1,n y(i) = mzero end do - else if (alpha == mone) then + else if (alpha == mone) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -mone) then + else if (alpha == -mone) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == mone) then + else + if (beta == mone) then ! Do nothing - else if (beta == -mone) then - y(1:n) = -y(1:n) + else if (beta == -mone) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == mzero) then + if (alpha == mzero) then ! do nothing - else if (alpha == mone) then + else if (alpha == mone) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -mone) then + else if (alpha == -mone) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_mgthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_mgthzv(n,idx,x,y) end subroutine psi_mgthzv subroutine psi_msctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_msctmm(n,k,idx,x,beta,y) end subroutine psi_msctmm subroutine psi_msctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_msctv subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info integer(psb_mpk_) X(lldx,*), Y(lldy,*) integer(psb_mpk_) alpha, beta @@ -447,19 +508,19 @@ subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.mzero) then - if (beta.eq.mzero) then - do j=1, n - do i=1,m + if (alpha.eq.mzero) then + if (beta.eq.mzero) then + do j=1, n + do i=1,m y(i,j) = mzero enddo enddo else if (beta.eq.mone) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-mone) then - do j=1,n - do i=1,m + else if (beta.eq.-mone) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.mone) then - if (beta.eq.mzero) then - do j=1,n - do i=1,m + if (beta.eq.mzero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.mone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-mone) then - do j=1,n - do i=1,m + else if (beta.eq.-mone) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-mone) then + else if (alpha.eq.-mone) then - if (beta.eq.mzero) then - do j=1,n - do i=1,m + if (beta.eq.mzero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.mone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-mone) then - do j=1,n - do i=1,m + else if (beta.eq.-mone) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.mzero) then - do j=1,n - do i=1,m + if (beta.eq.mzero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.mone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-mone) then - do j=1,n - do i=1,m + else if (beta.eq.-mone) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine maxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine maxpby + +subroutine maxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + integer(psb_mpk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + integer(psb_mpk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='maxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.mzero) then + if (beta.eq.mzero) then + do j=1, n + do i=1,m + Z(i,j) = mzero + enddo + enddo + else if (beta.eq.mone) then + ! + ! Do nothing! + ! + + else if (beta.eq.-mone) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.mone) then + + if (beta.eq.mzero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.mone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-mone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-mone) then + + if (beta.eq.mzero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.mone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-mone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.mzero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.mone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-mone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine maxpbyv2 diff --git a/base/serial/psi_s_serial_impl.f90 b/base/serial/psi_s_serial_impl.f90 index f9b83727..dfe2559b 100644 --- a/base/serial/psi_s_serial_impl.f90 +++ b/base/serial/psi_s_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_saxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n real(psb_spk_), intent (in) :: x(:,:) real(psb_spk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_saxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call saxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_saxpby(m,n,alpha, x, beta, y, info) end subroutine psi_saxpby subroutine psi_saxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m real(psb_spk_), intent (in) :: x(:) real(psb_spk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_saxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call saxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_saxpbyv(m,alpha, x, beta, y, info) end subroutine psi_saxpbyv +subroutine psi_saxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + real(psb_spk_), intent (in) :: x(:) + real(psb_spk_), intent (in) :: y(:) + real(psb_spk_), intent (inout) :: z(:) + real(psb_spk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call saxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_saxpbyv2 subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == szero) then - if (alpha == szero) then + if (beta == szero) then + if (alpha == szero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -sone) then + else if (alpha == -sone) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == sone) then + else + if (beta == sone) then ! Do nothing - else if (beta == -sone) then - y(1:n*k) = -y(1:n*k) + else if (beta == -sone) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == szero) then + if (alpha == szero) then ! do nothing - else if (alpha == sone) then + else if (alpha == sone) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_sgthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == szero) then - if (alpha == szero) then + if (beta == szero) then + if (alpha == szero) then do i=1,n y(i) = szero end do - else if (alpha == sone) then + else if (alpha == sone) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -sone) then + else if (alpha == -sone) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == sone) then + else + if (beta == sone) then ! Do nothing - else if (beta == -sone) then - y(1:n) = -y(1:n) + else if (beta == -sone) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == szero) then + if (alpha == szero) then ! do nothing - else if (alpha == sone) then + else if (alpha == sone) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -sone) then + else if (alpha == -sone) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_sgthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_sgthzv(n,idx,x,y) end subroutine psi_sgthzv subroutine psi_ssctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_ssctmm(n,k,idx,x,beta,y) end subroutine psi_ssctmm subroutine psi_ssctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_ssctv subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info real(psb_spk_) X(lldx,*), Y(lldy,*) real(psb_spk_) alpha, beta @@ -447,19 +508,19 @@ subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.szero) then - if (beta.eq.szero) then - do j=1, n - do i=1,m + if (alpha.eq.szero) then + if (beta.eq.szero) then + do j=1, n + do i=1,m y(i,j) = szero enddo enddo else if (beta.eq.sone) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-sone) then - do j=1,n - do i=1,m + else if (beta.eq.-sone) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.sone) then - if (beta.eq.szero) then - do j=1,n - do i=1,m + if (beta.eq.szero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.sone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-sone) then - do j=1,n - do i=1,m + else if (beta.eq.-sone) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-sone) then + else if (alpha.eq.-sone) then - if (beta.eq.szero) then - do j=1,n - do i=1,m + if (beta.eq.szero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.sone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-sone) then - do j=1,n - do i=1,m + else if (beta.eq.-sone) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.szero) then - do j=1,n - do i=1,m + if (beta.eq.szero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.sone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-sone) then - do j=1,n - do i=1,m + else if (beta.eq.-sone) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine saxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine saxpby + +subroutine saxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + real(psb_spk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + real(psb_spk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='saxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.szero) then + if (beta.eq.szero) then + do j=1, n + do i=1,m + Z(i,j) = szero + enddo + enddo + else if (beta.eq.sone) then + ! + ! Do nothing! + ! + + else if (beta.eq.-sone) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.sone) then + + if (beta.eq.szero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.sone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-sone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-sone) then + + if (beta.eq.szero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.sone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-sone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.szero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.sone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-sone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine saxpbyv2 diff --git a/base/serial/psi_z_serial_impl.f90 b/base/serial/psi_z_serial_impl.f90 index 8d945430..5b7036e6 100644 --- a/base/serial/psi_z_serial_impl.f90 +++ b/base/serial/psi_z_serial_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 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 @@ -27,13 +27,13 @@ ! 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. -! -! +! +! subroutine psi_zaxpby(m,n,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m, n complex(psb_dpk_), intent (in) :: x(:,:) complex(psb_dpk_), intent (inout) :: y(:,:) @@ -55,27 +55,27 @@ subroutine psi_zaxpby(m,n,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (n < 0) then info = psb_err_iarg_neg_ ierr(1) = 2; ierr(2) = n call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then + if (lx < m) then info = psb_err_input_asize_small_i_ ierr(1) = 4; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 6; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if ((m>0).and.(n>0)) call zaxpby(m,n,alpha,x,lx,beta,y,ly,info) @@ -89,10 +89,10 @@ subroutine psi_zaxpby(m,n,alpha, x, beta, y, info) end subroutine psi_zaxpby subroutine psi_zaxpbyv(m,alpha, x, beta, y, info) - + use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_), intent(in) :: m complex(psb_dpk_), intent (in) :: x(:) complex(psb_dpk_), intent (inout) :: y(:) @@ -114,21 +114,21 @@ subroutine psi_zaxpbyv(m,alpha, x, beta, y, info) info = psb_err_iarg_neg_ ierr(1) = 1; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if lx = size(x,1) ly = size(y,1) - if (lx < m) then - info = psb_err_input_asize_small_i_ + if (lx < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 3; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if - if (ly < m) then - info = psb_err_input_asize_small_i_ + if (ly < m) then + info = psb_err_input_asize_small_i_ ierr(1) = 5; ierr(2) = m call psb_errpush(info,name,i_err=ierr) - goto 9999 + goto 9999 end if if (m>0) call zaxpby(m,ione,alpha,x,lx,beta,y,ly,info) @@ -142,6 +142,67 @@ subroutine psi_zaxpbyv(m,alpha, x, beta, y, info) end subroutine psi_zaxpbyv +subroutine psi_zaxpbyv2(m,alpha, x, beta, y, z, info) + + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_), intent(in) :: m + complex(psb_dpk_), intent (in) :: x(:) + complex(psb_dpk_), intent (in) :: y(:) + complex(psb_dpk_), intent (inout) :: z(:) + complex(psb_dpk_), intent (in) :: alpha, beta + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: lx, ly, lz + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name, ch_err + + name='psb_geaxpby' + info=psb_success_ + call psb_erractionsave(err_act) + if (psb_errstatus_fatal()) then + info = psb_err_internal_error_ ; goto 9999 + end if + + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 1; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + lx = size(x,1) + ly = size(y,1) + lz = size(z,1) + if (lx < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 3; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (ly < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (lz < m) then + info = psb_err_input_asize_small_i_ + ierr(1) = 5; ierr(2) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (m>0) call zaxpbyv2(m,ione,alpha,x,lx,beta,y,ly,z,lz,info) + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine psi_zaxpbyv2 subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y) @@ -154,8 +215,8 @@ subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i, j, pt - if (beta == zzero) then - if (alpha == zzero) then + if (beta == zzero) then + if (alpha == zzero) then pt=0 do j=1,k do i=1,n @@ -171,11 +232,11 @@ subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y) y(pt) = x(idx(i),j) end do end do - else if (alpha == -zone) then + else if (alpha == -zone) then pt=0 do j=1,k do i=1,n - pt=pt+1 + pt=pt+1 y(pt) = -x(idx(i),j) end do end do @@ -188,18 +249,18 @@ subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y) end do end do end if - else - if (beta == zone) then + else + if (beta == zone) then ! Do nothing - else if (beta == -zone) then - y(1:n*k) = -y(1:n*k) + else if (beta == -zone) then + y(1:n*k) = -y(1:n*k) else - y(1:n*k) = beta*y(1:n*k) + y(1:n*k) = beta*y(1:n*k) end if - if (alpha == zzero) then + if (alpha == zzero) then ! do nothing - else if (alpha == zone) then + else if (alpha == zone) then pt=0 do j=1,k do i=1,n @@ -215,7 +276,7 @@ subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y) y(pt) = y(pt) - x(idx(i),j) end do end do - else + else pt=0 do j=1,k do i=1,n @@ -238,44 +299,44 @@ subroutine psi_zgthv(n,idx,alpha,x,beta,y) ! Locals integer(psb_ipk_) :: i - if (beta == zzero) then - if (alpha == zzero) then + if (beta == zzero) then + if (alpha == zzero) then do i=1,n y(i) = zzero end do - else if (alpha == zone) then + else if (alpha == zone) then do i=1,n y(i) = x(idx(i)) end do - else if (alpha == -zone) then + else if (alpha == -zone) then do i=1,n y(i) = -x(idx(i)) end do - else + else do i=1,n y(i) = alpha*x(idx(i)) end do end if - else - if (beta == zone) then + else + if (beta == zone) then ! Do nothing - else if (beta == -zone) then - y(1:n) = -y(1:n) + else if (beta == -zone) then + y(1:n) = -y(1:n) else - y(1:n) = beta*y(1:n) + y(1:n) = beta*y(1:n) end if - if (alpha == zzero) then + if (alpha == zzero) then ! do nothing - else if (alpha == zone) then + else if (alpha == zone) then do i=1,n y(i) = y(i) + x(idx(i)) end do - else if (alpha == -zone) then + else if (alpha == -zone) then do i=1,n y(i) = y(i) - x(idx(i)) end do - else + else do i=1,n y(i) = y(i) + alpha*x(idx(i)) end do @@ -295,7 +356,7 @@ subroutine psi_zgthzmm(n,k,idx,x,y) ! Locals integer(psb_ipk_) :: i - + do i=1,n y(i,1:k)=x(idx(i),1:k) end do @@ -341,7 +402,7 @@ subroutine psi_zgthzv(n,idx,x,y) end subroutine psi_zgthzv subroutine psi_zsctmm(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -367,7 +428,7 @@ subroutine psi_zsctmm(n,k,idx,x,beta,y) end subroutine psi_zsctmm subroutine psi_zsctmv(n,k,idx,x,beta,y) - + use psb_const_mod implicit none @@ -433,7 +494,7 @@ end subroutine psi_zsctv subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) use psb_const_mod use psb_error_mod - implicit none + implicit none integer(psb_ipk_) :: n, m, lldx, lldy, info complex(psb_dpk_) X(lldx,*), Y(lldy,*) complex(psb_dpk_) alpha, beta @@ -447,19 +508,19 @@ subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) ! Error handling ! info = psb_success_ - if (m.lt.0) then + if (m.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (n.lt.0) then + else if (n.lt.0) then info=psb_err_iarg_neg_ int_err(1)=1 int_err(2)=n call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldx.lt.max(1,m)) then + else if (lldx.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=5 int_err(2)=1 @@ -467,7 +528,7 @@ subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) int_err(4)=m call fcpsb_errpush(info,name,int_err) goto 9999 - else if (lldy.lt.max(1,m)) then + else if (lldy.lt.max(1,m)) then info=psb_err_iarg_not_gtia_ii_ int_err(1)=8 int_err(2)=1 @@ -477,27 +538,27 @@ subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) goto 9999 endif - if (alpha.eq.zzero) then - if (beta.eq.zzero) then - do j=1, n - do i=1,m + if (alpha.eq.zzero) then + if (beta.eq.zzero) then + do j=1, n + do i=1,m y(i,j) = zzero enddo enddo else if (beta.eq.zone) then - ! - ! Do nothing! - ! + ! + ! Do nothing! + ! - else if (beta.eq.-zone) then - do j=1,n - do i=1,m + else if (beta.eq.-zone) then + do j=1,n + do i=1,m y(i,j) = - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = beta*y(i,j) enddo enddo @@ -505,86 +566,86 @@ subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) else if (alpha.eq.zone) then - if (beta.eq.zzero) then - do j=1,n - do i=1,m + if (beta.eq.zzero) then + do j=1,n + do i=1,m y(i,j) = x(i,j) enddo enddo else if (beta.eq.zone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-zone) then - do j=1,n - do i=1,m + else if (beta.eq.-zone) then + do j=1,n + do i=1,m y(i,j) = x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = x(i,j) + beta*y(i,j) enddo enddo endif - else if (alpha.eq.-zone) then + else if (alpha.eq.-zone) then - if (beta.eq.zzero) then - do j=1,n - do i=1,m + if (beta.eq.zzero) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) enddo enddo else if (beta.eq.zone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = -x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-zone) then - do j=1,n - do i=1,m + else if (beta.eq.-zone) then + do j=1,n + do i=1,m y(i,j) = -x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = -x(i,j) + beta*y(i,j) enddo enddo endif - else + else - if (beta.eq.zzero) then - do j=1,n - do i=1,m + if (beta.eq.zzero) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) enddo enddo else if (beta.eq.zone) then - do j=1,n - do i=1,m + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + y(i,j) enddo enddo - else if (beta.eq.-zone) then - do j=1,n - do i=1,m + else if (beta.eq.-zone) then + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) - y(i,j) enddo enddo - else - do j=1,n - do i=1,m + else + do j=1,n + do i=1,m y(i,j) = alpha*x(i,j) + beta*y(i,j) enddo enddo @@ -599,3 +660,181 @@ subroutine zaxpby(m, n, alpha, X, lldx, beta, Y, lldy, info) return end subroutine zaxpby + +subroutine zaxpbyv2(m, n, alpha, X, lldx, beta, Y, lldy, Z, lldz, info) + use psb_const_mod + use psb_error_mod + implicit none + integer(psb_ipk_) :: n, m, lldx, lldy, lldz, info + complex(psb_dpk_) X(lldx,*), Y(lldy,*), Z(lldy,*) + complex(psb_dpk_) alpha, beta + integer(psb_ipk_) :: i, j + integer(psb_ipk_) :: int_err(5) + character name*20 + name='zaxpby' + + + ! + ! Error handling + ! + info = psb_success_ + if (m.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (n.lt.0) then + info=psb_err_iarg_neg_ + int_err(1)=1 + int_err(2)=n + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldx.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=5 + int_err(2)=1 + int_err(3)=lldx + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldy.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldy + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + else if (lldz.lt.max(1,m)) then + info=psb_err_iarg_not_gtia_ii_ + int_err(1)=8 + int_err(2)=1 + int_err(3)=lldz + int_err(4)=m + call fcpsb_errpush(info,name,int_err) + goto 9999 + endif + + if (alpha.eq.zzero) then + if (beta.eq.zzero) then + do j=1, n + do i=1,m + Z(i,j) = zzero + enddo + enddo + else if (beta.eq.zone) then + ! + ! Do nothing! + ! + + else if (beta.eq.-zone) then + do j=1,n + do i=1,m + Z(i,j) = - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.zone) then + + if (beta.eq.zzero) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + enddo + enddo + else if (beta.eq.zone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-zone) then + do j=1,n + do i=1,m + Z(i,j) = x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else if (alpha.eq.-zone) then + + if (beta.eq.zzero) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + enddo + enddo + else if (beta.eq.zone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-zone) then + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = -x(i,j) + beta*y(i,j) + enddo + enddo + endif + + else + + if (beta.eq.zzero) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + enddo + enddo + else if (beta.eq.zone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + y(i,j) + enddo + enddo + + else if (beta.eq.-zone) then + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) - y(i,j) + enddo + enddo + else + do j=1,n + do i=1,m + Z(i,j) = alpha*x(i,j) + beta*y(i,j) + enddo + enddo + endif + + endif + + return + +9999 continue + call fcpsb_serror() + return + +end subroutine zaxpbyv2 diff --git a/cbind/base/psb_c_cbase.h b/cbind/base/psb_c_cbase.h index 9ae6ef1f..f0258419 100644 --- a/cbind/base/psb_c_cbase.h +++ b/cbind/base/psb_c_cbase.h @@ -53,6 +53,8 @@ psb_s_t psb_c_cgeasum(psb_c_cvector *xh, psb_c_descriptor *cdh); psb_s_t psb_c_cspnrmi(psb_c_cspmat *ah, psb_c_descriptor *cdh); psb_i_t psb_c_cgeaxpby(psb_c_t alpha, psb_c_cvector *xh, psb_c_t beta, psb_c_cvector *yh, psb_c_descriptor *cdh); +psb_i_t psb_c_cgeaxpbyz(psb_c_t alpha, psb_c_cvector *xh, + psb_c_t beta, psb_c_cvector *yh, psb_c_cvector *zh, psb_c_descriptor *cdh); psb_i_t psb_c_cspmm(psb_c_t alpha, psb_c_cspmat *ah, psb_c_cvector *xh, psb_c_t beta, psb_c_cvector *yh, psb_c_descriptor *cdh); psb_i_t psb_c_cspmm_opt(psb_c_t alpha, psb_c_cspmat *ah, psb_c_cvector *xh, diff --git a/cbind/base/psb_c_dbase.h b/cbind/base/psb_c_dbase.h index 04f5d232..a3104107 100644 --- a/cbind/base/psb_c_dbase.h +++ b/cbind/base/psb_c_dbase.h @@ -53,6 +53,8 @@ psb_d_t psb_c_dgeasum(psb_c_dvector *xh, psb_c_descriptor *cdh); psb_d_t psb_c_dspnrmi(psb_c_dvector *xh, psb_c_descriptor *cdh); psb_i_t psb_c_dgeaxpby(psb_d_t alpha, psb_c_dvector *xh, psb_d_t beta, psb_c_dvector *yh, psb_c_descriptor *cdh); +psb_i_t psb_c_dgeaxpbyz(psb_d_t alpha, psb_c_dvector *xh, + psb_d_t beta, psb_c_dvector *yh, psb_c_dvector *zh, psb_c_descriptor *cdh); psb_i_t psb_c_dspmm(psb_d_t alpha, psb_c_dspmat *ah, psb_c_dvector *xh, psb_d_t beta, psb_c_dvector *yh, psb_c_descriptor *cdh); psb_i_t psb_c_dspmm_opt(psb_d_t alpha, psb_c_dspmat *ah, psb_c_dvector *xh, diff --git a/cbind/base/psb_c_psblas_cbind_mod.f90 b/cbind/base/psb_c_psblas_cbind_mod.f90 index 608db58b..6bbee920 100644 --- a/cbind/base/psb_c_psblas_cbind_mod.f90 +++ b/cbind/base/psb_c_psblas_cbind_mod.f90 @@ -43,6 +43,49 @@ contains end function psb_c_cgeaxpby + function psb_c_cgeaxpbyz(alpha,xh,beta,yh,zh,cdh) bind(c) result(res) + implicit none + integer(psb_c_ipk_) :: res + + type(psb_c_cvector) :: xh,yh,zh + type(psb_c_descriptor) :: cdh + complex(c_float_complex), value :: alpha,beta + + type(psb_desc_type), pointer :: descp + type(psb_c_vect_type), pointer :: xp,yp,zp + integer(psb_c_ipk_) :: info + + + res = -1 + + if (c_associated(cdh%item)) then + call c_f_pointer(cdh%item,descp) + else + return + end if + if (c_associated(xh%item)) then + call c_f_pointer(xh%item,xp) + else + return + end if + if (c_associated(yh%item)) then + call c_f_pointer(yh%item,yp) + else + return + end if + + if (c_associated(zh%item)) then + call c_f_pointer(zh%item,zp) + else + return + end if + + call psb_geaxpby(alpha,xp,beta,yp,zp,descp,info) + + res = info + + end function psb_c_cgeaxpbyz + function psb_c_cgemlt(xh,yh,cdh) bind(c) result(res) implicit none integer(psb_c_ipk_) :: res diff --git a/cbind/base/psb_c_sbase.h b/cbind/base/psb_c_sbase.h index 6aa4db18..55f6b0e8 100644 --- a/cbind/base/psb_c_sbase.h +++ b/cbind/base/psb_c_sbase.h @@ -53,6 +53,8 @@ psb_s_t psb_c_sgeasum(psb_c_svector *xh, psb_c_descriptor *cdh); psb_s_t psb_c_sspnrmi(psb_c_sspmat *ah, psb_c_descriptor *cdh); psb_i_t psb_c_sgeaxpby(psb_s_t alpha, psb_c_svector *xh, psb_s_t beta, psb_c_svector *yh, psb_c_descriptor *cdh); +psb_i_t psb_c_sgeaxpbyz(psb_s_t alpha, psb_c_svector *xh, + psb_s_t beta, psb_c_svector *yh, psb_c_svector *zh, psb_c_descriptor *cdh); psb_i_t psb_c_sspmm(psb_s_t alpha, psb_c_sspmat *ah, psb_c_svector *xh, psb_s_t beta, psb_c_svector *yh, psb_c_descriptor *cdh); psb_i_t psb_c_sspmm_opt(psb_s_t alpha, psb_c_sspmat *ah, psb_c_svector *xh, diff --git a/cbind/base/psb_c_zbase.h b/cbind/base/psb_c_zbase.h index c3c9e342..fadd23a9 100644 --- a/cbind/base/psb_c_zbase.h +++ b/cbind/base/psb_c_zbase.h @@ -53,6 +53,8 @@ psb_d_t psb_c_zgeasum(psb_c_zvector *xh, psb_c_descriptor *cdh); psb_d_t psb_c_zspnrmi(psb_c_zspmat *ah, psb_c_descriptor *cdh); psb_i_t psb_c_zgeaxpby(psb_z_t alpha, psb_c_zvector *xh, psb_z_t beta, psb_c_zvector *yh, psb_c_descriptor *cdh); +psb_i_t psb_c_zgeaxpbyz(psb_z_t alpha, psb_c_zvector *xh, + psb_z_t beta, psb_c_zvector *yh, psb_c_zvector *zh, psb_c_descriptor *cdh); psb_i_t psb_c_zspmm(psb_z_t alpha, psb_c_zspmat *ah, psb_c_zvector *xh, psb_z_t beta, psb_c_zvector *yh, psb_c_descriptor *cdh); psb_i_t psb_c_zspmm_opt(psb_z_t alpha, psb_c_zspmat *ah, psb_c_zvector *xh, diff --git a/cbind/base/psb_d_psblas_cbind_mod.f90 b/cbind/base/psb_d_psblas_cbind_mod.f90 index 37d315c0..fd2d670a 100644 --- a/cbind/base/psb_d_psblas_cbind_mod.f90 +++ b/cbind/base/psb_d_psblas_cbind_mod.f90 @@ -43,6 +43,49 @@ contains end function psb_c_dgeaxpby + function psb_c_dgeaxpbyz(alpha,xh,beta,yh,zh,cdh) bind(c) result(res) + implicit none + integer(psb_c_ipk_) :: res + + type(psb_c_dvector) :: xh,yh,zh + type(psb_c_descriptor) :: cdh + real(c_double), value :: alpha,beta + + type(psb_desc_type), pointer :: descp + type(psb_d_vect_type), pointer :: xp,yp,zp + integer(psb_c_ipk_) :: info + + + res = -1 + + if (c_associated(cdh%item)) then + call c_f_pointer(cdh%item,descp) + else + return + end if + if (c_associated(xh%item)) then + call c_f_pointer(xh%item,xp) + else + return + end if + if (c_associated(yh%item)) then + call c_f_pointer(yh%item,yp) + else + return + end if + + if (c_associated(zh%item)) then + call c_f_pointer(zh%item,zp) + else + return + end if + + call psb_geaxpby(alpha,xp,beta,yp,zp,descp,info) + + res = info + + end function psb_c_dgeaxpbyz + function psb_c_dgemlt(xh,yh,cdh) bind(c) result(res) implicit none integer(psb_c_ipk_) :: res diff --git a/cbind/base/psb_s_psblas_cbind_mod.f90 b/cbind/base/psb_s_psblas_cbind_mod.f90 index dd42f6f5..9a7a3571 100644 --- a/cbind/base/psb_s_psblas_cbind_mod.f90 +++ b/cbind/base/psb_s_psblas_cbind_mod.f90 @@ -43,6 +43,49 @@ contains end function psb_c_sgeaxpby + function psb_c_sgeaxpbyz(alpha,xh,beta,yh,zh,cdh) bind(c) result(res) + implicit none + integer(psb_c_ipk_) :: res + + type(psb_c_svector) :: xh,yh,zh + type(psb_c_descriptor) :: cdh + real(c_float), value :: alpha,beta + + type(psb_desc_type), pointer :: descp + type(psb_s_vect_type), pointer :: xp,yp,zp + integer(psb_c_ipk_) :: info + + + res = -1 + + if (c_associated(cdh%item)) then + call c_f_pointer(cdh%item,descp) + else + return + end if + if (c_associated(xh%item)) then + call c_f_pointer(xh%item,xp) + else + return + end if + if (c_associated(yh%item)) then + call c_f_pointer(yh%item,yp) + else + return + end if + + if (c_associated(zh%item)) then + call c_f_pointer(zh%item,zp) + else + return + end if + + call psb_geaxpby(alpha,xp,beta,yp,zp,descp,info) + + res = info + + end function psb_c_sgeaxpbyz + function psb_c_sgemlt(xh,yh,cdh) bind(c) result(res) implicit none integer(psb_c_ipk_) :: res diff --git a/cbind/base/psb_z_psblas_cbind_mod.f90 b/cbind/base/psb_z_psblas_cbind_mod.f90 index b2506853..1b7f83a6 100644 --- a/cbind/base/psb_z_psblas_cbind_mod.f90 +++ b/cbind/base/psb_z_psblas_cbind_mod.f90 @@ -43,6 +43,49 @@ contains end function psb_c_zgeaxpby + function psb_c_zgeaxpbyz(alpha,xh,beta,yh,zh,cdh) bind(c) result(res) + implicit none + integer(psb_c_ipk_) :: res + + type(psb_c_zvector) :: xh,yh,zh + type(psb_c_descriptor) :: cdh + complex(c_double_complex), value :: alpha,beta + + type(psb_desc_type), pointer :: descp + type(psb_z_vect_type), pointer :: xp,yp,zp + integer(psb_c_ipk_) :: info + + + res = -1 + + if (c_associated(cdh%item)) then + call c_f_pointer(cdh%item,descp) + else + return + end if + if (c_associated(xh%item)) then + call c_f_pointer(xh%item,xp) + else + return + end if + if (c_associated(yh%item)) then + call c_f_pointer(yh%item,yp) + else + return + end if + + if (c_associated(zh%item)) then + call c_f_pointer(zh%item,zp) + else + return + end if + + call psb_geaxpby(alpha,xp,beta,yp,zp,descp,info) + + res = info + + end function psb_c_zgeaxpbyz + function psb_c_zgemlt(xh,yh,cdh) bind(c) result(res) implicit none integer(psb_c_ipk_) :: res