Added out-of-place axpby and relative c interface

base/modules/auxil/psi_c_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

base/modules/auxil/psi_d_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

base/modules/auxil/psi_e_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

base/modules/auxil/psi_m_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

base/modules/auxil/psi_s_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

base/modules/auxil/psi_z_serial_mod.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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

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_, &

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_, &

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_, &

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_, &

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:

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

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:

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

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:

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

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:

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

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,

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,

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,

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,

base/serial/psi_c_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == czero) then
-    if (alpha == 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 
+  else
-    if (beta == cone) then 
+    if (beta == cone) then
       ! Do nothing
-    else if (beta == -cone) then 
+    else if (beta == -cone) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == czero) then
-    if (alpha == 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 
+  else
-    if (beta == cone) then 
+    if (beta == cone) then
       ! Do nothing
-    else if (beta == -cone) then 
+    else if (beta == -cone) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.czero) then
-    if (beta.eq.czero) then 
+    if (beta.eq.czero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.czero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-cone) then 
+    else if (beta.eq.-cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.czero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-cone) then 
+    else if (beta.eq.-cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.czero) then 
+    if (beta.eq.czero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-cone) then 
+    else if (beta.eq.-cone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

base/serial/psi_d_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == dzero) then
-    if (alpha == 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 
+  else
-    if (beta == done) then 
+    if (beta == done) then
       ! Do nothing
-    else if (beta == -done) then 
+    else if (beta == -done) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == dzero) then
-    if (alpha == 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 
+  else
-    if (beta == done) then 
+    if (beta == done) then
       ! Do nothing
-    else if (beta == -done) then 
+    else if (beta == -done) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.dzero) then
-    if (beta.eq.dzero) then 
+    if (beta.eq.dzero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.dzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-done) then 
+    else if (beta.eq.-done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.dzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-done) then 
+    else if (beta.eq.-done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.dzero) then 
+    if (beta.eq.dzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-done) then 
+    else if (beta.eq.-done) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

base/serial/psi_e_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == ezero) then
-    if (alpha == 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 
+  else
-    if (beta == eone) then 
+    if (beta == eone) then
       ! Do nothing
-    else if (beta == -eone) then 
+    else if (beta == -eone) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == ezero) then
-    if (alpha == 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 
+  else
-    if (beta == eone) then 
+    if (beta == eone) then
       ! Do nothing
-    else if (beta == -eone) then 
+    else if (beta == -eone) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.ezero) then
-    if (beta.eq.ezero) then 
+    if (beta.eq.ezero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.ezero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-eone) then 
+    else if (beta.eq.-eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.ezero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-eone) then 
+    else if (beta.eq.-eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.ezero) then 
+    if (beta.eq.ezero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-eone) then 
+    else if (beta.eq.-eone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

base/serial/psi_m_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == mzero) then
-    if (alpha == 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 
+  else
-    if (beta == mone) then 
+    if (beta == mone) then
       ! Do nothing
-    else if (beta == -mone) then 
+    else if (beta == -mone) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == mzero) then
-    if (alpha == 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 
+  else
-    if (beta == mone) then 
+    if (beta == mone) then
       ! Do nothing
-    else if (beta == -mone) then 
+    else if (beta == -mone) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.mzero) then
-    if (beta.eq.mzero) then 
+    if (beta.eq.mzero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.mzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-mone) then 
+    else if (beta.eq.-mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.mzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-mone) then 
+    else if (beta.eq.-mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.mzero) then 
+    if (beta.eq.mzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-mone) then 
+    else if (beta.eq.-mone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

base/serial/psi_s_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == szero) then
-    if (alpha == 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 
+  else
-    if (beta == sone) then 
+    if (beta == sone) then
       ! Do nothing
-    else if (beta == -sone) then 
+    else if (beta == -sone) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == szero) then
-    if (alpha == 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 
+  else
-    if (beta == sone) then 
+    if (beta == sone) then
       ! Do nothing
-    else if (beta == -sone) then 
+    else if (beta == -sone) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.szero) then
-    if (beta.eq.szero) then 
+    if (beta.eq.szero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.szero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-sone) then 
+    else if (beta.eq.-sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.szero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-sone) then 
+    else if (beta.eq.-sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.szero) then 
+    if (beta.eq.szero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-sone) then 
+    else if (beta.eq.-sone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

base/serial/psi_z_serial_impl.f90

@@ -1,9 +1,9 @@
-!   
+!
 !                Parallel Sparse BLAS  version 3.5
 !      (C) Copyright 2006-2018
-!        Salvatore Filippone    
+!        Salvatore Filippone
-!        Alfredo Buttari      
+!        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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (lx < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 
+  if (ly < m) then
-    info = psb_err_input_asize_small_i_ 
+    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 (beta == zzero) then
-    if (alpha == 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 
+  else
-    if (beta == zone) then 
+    if (beta == zone) then
       ! Do nothing
-    else if (beta == -zone) then 
+    else if (beta == -zone) then
-      y(1:n*k) = -y(1:n*k) 
+      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 (beta == zzero) then
-    if (alpha == 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 
+  else
-    if (beta == zone) then 
+    if (beta == zone) then
       ! Do nothing
-    else if (beta == -zone) then 
+    else if (beta == -zone) then
-      y(1:n) = -y(1:n) 
+      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 (alpha.eq.zzero) then
-    if (beta.eq.zzero) then 
+    if (beta.eq.zzero) then
-      do j=1, n 
+      do j=1, n
-        do i=1,m 
+        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 
+    else if (beta.eq.-zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.zzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j)
         enddo
       enddo
     else if (beta.eq.zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-zone) then 
+    else if (beta.eq.-zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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 
+    if (beta.eq.zzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j)
         enddo
       enddo
     else if (beta.eq.zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-zone) then 
+    else if (beta.eq.-zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = -x(i,j) + beta*y(i,j)
         enddo
       enddo
     endif
 
-  else  
+  else
 
-    if (beta.eq.zzero) then 
+    if (beta.eq.zzero) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j)
         enddo
       enddo
     else if (beta.eq.zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) + y(i,j)
         enddo
       enddo
 
-    else if (beta.eq.-zone) then 
+    else if (beta.eq.-zone) then
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        do i=1,m
           y(i,j) = alpha*x(i,j) - y(i,j)
         enddo
       enddo
-    else  
+    else
-      do j=1,n 
+      do j=1,n
-        do i=1,m 
+        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

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,

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,

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

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,

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,

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

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

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