New method W%XYZW()

base/modules/auxil/psi_c_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_c_serial_mod
     end subroutine psi_cabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_cxyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_), intent(in)      :: m
+      complex(psb_spk_), intent (in)       :: x(:)
+      complex(psb_spk_), intent (inout)    :: y(:)
+      complex(psb_spk_), intent (inout)    :: z(:)
+      complex(psb_spk_), intent (inout)    :: w(:)
+      complex(psb_spk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_cxyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_spk_

base/modules/auxil/psi_d_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_d_serial_mod
     end subroutine psi_dabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_dxyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_), intent(in)      :: m
+      real(psb_dpk_), intent (in)       :: x(:)
+      real(psb_dpk_), intent (inout)    :: y(:)
+      real(psb_dpk_), intent (inout)    :: z(:)
+      real(psb_dpk_), intent (inout)    :: w(:)
+      real(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_dxyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_dpk_

base/modules/auxil/psi_e_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_e_serial_mod
     end subroutine psi_eabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_exyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+      integer(psb_epk_), intent (inout)    :: z(:)
+      integer(psb_epk_), intent (inout)    :: w(:)
+      integer(psb_epk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_exyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_egthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_

base/modules/auxil/psi_i2_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_i2_serial_mod
     end subroutine psi_i2abgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_i2xyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+      import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_
+      implicit none
+      integer(psb_ipk_), intent(in)      :: m
+      integer(psb_i2pk_), intent (in)       :: x(:)
+      integer(psb_i2pk_), intent (inout)    :: y(:)
+      integer(psb_i2pk_), intent (inout)    :: z(:)
+      integer(psb_i2pk_), intent (inout)    :: w(:)
+      integer(psb_i2pk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_i2xyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_i2gthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_

base/modules/auxil/psi_m_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_m_serial_mod
     end subroutine psi_mabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_mxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+      integer(psb_mpk_), intent (inout)    :: z(:)
+      integer(psb_mpk_), intent (inout)    :: w(:)
+      integer(psb_mpk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_mxyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_mgthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_lpk_,psb_mpk_, psb_epk_

base/modules/auxil/psi_s_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_s_serial_mod
     end subroutine psi_sabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_sxyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+      import :: psb_ipk_, psb_spk_
+      implicit none
+      integer(psb_ipk_), intent(in)      :: m
+      real(psb_spk_), intent (in)       :: x(:)
+      real(psb_spk_), intent (inout)    :: y(:)
+      real(psb_spk_), intent (inout)    :: z(:)
+      real(psb_spk_), intent (inout)    :: w(:)
+      real(psb_spk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_sxyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_spk_

base/modules/auxil/psi_z_serial_mod.f90

@@ -112,6 +112,20 @@ module psi_z_serial_mod
     end subroutine psi_zabgdxyz
   end interface psi_abgdxyz
   
+  interface psi_xyzw
+    subroutine psi_zxyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+      import :: psb_ipk_, psb_dpk_
+      implicit none
+      integer(psb_ipk_), intent(in)      :: m
+      complex(psb_dpk_), intent (in)       :: x(:)
+      complex(psb_dpk_), intent (inout)    :: y(:)
+      complex(psb_dpk_), intent (inout)    :: z(:)
+      complex(psb_dpk_), intent (inout)    :: w(:)
+      complex(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+      integer(psb_ipk_), intent(out)     :: info
+    end subroutine psi_zxyzw
+  end interface psi_xyzw
+  
   interface psi_gth
     subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y)
       import :: psb_ipk_, psb_dpk_

base/modules/serial/psb_c_base_vect_mod.F90

@@ -156,6 +156,7 @@ module psb_c_base_vect_mod
     procedure, pass(z) :: axpby_a2  => c_base_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => c_base_abgdxyz
+    procedure, pass(w) :: xyzw     => c_base_xyzw
     
     !
     ! Vector by vector multiplication. Need all variants
@@ -1155,22 +1156,37 @@ contains
     complex(psb_spk_), intent (in)       :: alpha, beta, gamma, delta
     integer(psb_ipk_), intent(out)              :: info
 
-    if (.false.) then 
-      if (x%is_dev()) call x%sync()
-      
-      call y%axpby(m,alpha,x,beta,info)
-      call z%axpby(m,gamma,y,delta,info)
-    else
     if (x%is_dev().and.(alpha/=czero)) call x%sync()
     if (y%is_dev().and.(beta/=czero))   call y%sync()
     if (z%is_dev().and.(delta/=czero))  call z%sync()
-      call psi_cabgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
+    call psi_abgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
     call y%set_host()
     call z%set_host()
-    end if
         
   end subroutine c_base_abgdxyz
 
+  subroutine c_base_xyzw(m,a,b,c,d,e,f,x, y, z, w,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
+    class(psb_c_base_vect_type), intent(inout)  :: w
+    complex(psb_spk_), intent (in)                :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+
+    if (x%is_dev().and.(a/=czero)) call x%sync()
+    if (y%is_dev().and.(b/=czero)) call y%sync()
+    if (z%is_dev().and.(d/=czero)) call z%sync()
+    if (w%is_dev().and.(f/=czero)) call w%sync()
+    call psi_xyzw(m,a,b,c,d,e,f,x%v, y%v, z%v, w%v, info)
+    call y%set_host()
+    call z%set_host()
+    call w%set_host()
+    
+  end subroutine c_base_xyzw
+
 
   !
   !  Multiple variants of two operations:

base/modules/serial/psb_c_vect_mod.F90

@@ -103,6 +103,7 @@ module psb_c_vect_mod
     procedure, pass(z) :: axpby_a2  => c_vect_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => c_vect_abgdxyz
+    procedure, pass(z) :: xyzw     => c_vect_xyzw
     
     procedure, pass(y) :: mlt_v    => c_vect_mlt_v
     procedure, pass(y) :: mlt_a    => c_vect_mlt_a
@@ -788,6 +789,22 @@ contains
     
   end subroutine c_vect_abgdxyz
 
+  subroutine c_vect_xyzw(m,a,b,c,d,e,f,x, y, z, w, 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
+    class(psb_c_vect_type), intent(inout)  :: w
+    complex(psb_spk_), intent (in)     :: a, b, c, d, e, f
+    integer(psb_ipk_), intent(out)   :: info
+
+    if (allocated(w%v)) &
+         call w%v%xyzw(m,a,b,c,d,e,f,x%v,y%v,z%v,info)
+    
+  end subroutine c_vect_xyzw
+
 
   subroutine c_vect_mlt_v(x, y, info)
     use psi_serial_mod

base/modules/serial/psb_d_base_vect_mod.F90

@@ -156,6 +156,7 @@ module psb_d_base_vect_mod
     procedure, pass(z) :: axpby_a2  => d_base_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => d_base_abgdxyz
+    procedure, pass(w) :: xyzw     => d_base_xyzw
     
     !
     ! Vector by vector multiplication. Need all variants
@@ -1162,22 +1163,37 @@ contains
     real(psb_dpk_), intent (in)       :: alpha, beta, gamma, delta
     integer(psb_ipk_), intent(out)              :: info
 
-    if (.false.) then 
-      if (x%is_dev()) call x%sync()
-      
-      call y%axpby(m,alpha,x,beta,info)
-      call z%axpby(m,gamma,y,delta,info)
-    else
     if (x%is_dev().and.(alpha/=dzero)) call x%sync()
     if (y%is_dev().and.(beta/=dzero))   call y%sync()
     if (z%is_dev().and.(delta/=dzero))  call z%sync()
-      call psi_dabgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
+    call psi_abgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
     call y%set_host()
     call z%set_host()
-    end if
         
   end subroutine d_base_abgdxyz
 
+  subroutine d_base_xyzw(m,a,b,c,d,e,f,x, y, z, w,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
+    class(psb_d_base_vect_type), intent(inout)  :: w
+    real(psb_dpk_), intent (in)                :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+
+    if (x%is_dev().and.(a/=dzero)) call x%sync()
+    if (y%is_dev().and.(b/=dzero)) call y%sync()
+    if (z%is_dev().and.(d/=dzero)) call z%sync()
+    if (w%is_dev().and.(f/=dzero)) call w%sync()
+    call psi_xyzw(m,a,b,c,d,e,f,x%v, y%v, z%v, w%v, info)
+    call y%set_host()
+    call z%set_host()
+    call w%set_host()
+    
+  end subroutine d_base_xyzw
+
 
   !
   !  Multiple variants of two operations:

base/modules/serial/psb_d_vect_mod.F90

@@ -103,6 +103,7 @@ module psb_d_vect_mod
     procedure, pass(z) :: axpby_a2  => d_vect_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => d_vect_abgdxyz
+    procedure, pass(z) :: xyzw     => d_vect_xyzw
     
     procedure, pass(y) :: mlt_v    => d_vect_mlt_v
     procedure, pass(y) :: mlt_a    => d_vect_mlt_a
@@ -795,6 +796,22 @@ contains
     
   end subroutine d_vect_abgdxyz
 
+  subroutine d_vect_xyzw(m,a,b,c,d,e,f,x, y, z, w, 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
+    class(psb_d_vect_type), intent(inout)  :: w
+    real(psb_dpk_), intent (in)     :: a, b, c, d, e, f
+    integer(psb_ipk_), intent(out)   :: info
+
+    if (allocated(w%v)) &
+         call w%v%xyzw(m,a,b,c,d,e,f,x%v,y%v,z%v,info)
+    
+  end subroutine d_vect_xyzw
+
 
   subroutine d_vect_mlt_v(x, y, info)
     use psi_serial_mod

base/modules/serial/psb_s_base_vect_mod.F90

@@ -156,6 +156,7 @@ module psb_s_base_vect_mod
     procedure, pass(z) :: axpby_a2  => s_base_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => s_base_abgdxyz
+    procedure, pass(w) :: xyzw     => s_base_xyzw
     
     !
     ! Vector by vector multiplication. Need all variants
@@ -1162,22 +1163,37 @@ contains
     real(psb_spk_), intent (in)       :: alpha, beta, gamma, delta
     integer(psb_ipk_), intent(out)              :: info
 
-    if (.false.) then 
-      if (x%is_dev()) call x%sync()
-      
-      call y%axpby(m,alpha,x,beta,info)
-      call z%axpby(m,gamma,y,delta,info)
-    else
     if (x%is_dev().and.(alpha/=szero)) call x%sync()
     if (y%is_dev().and.(beta/=szero))   call y%sync()
     if (z%is_dev().and.(delta/=szero))  call z%sync()
-      call psi_sabgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
+    call psi_abgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
     call y%set_host()
     call z%set_host()
-    end if
         
   end subroutine s_base_abgdxyz
 
+  subroutine s_base_xyzw(m,a,b,c,d,e,f,x, y, z, w,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
+    class(psb_s_base_vect_type), intent(inout)  :: w
+    real(psb_spk_), intent (in)                :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+
+    if (x%is_dev().and.(a/=szero)) call x%sync()
+    if (y%is_dev().and.(b/=szero)) call y%sync()
+    if (z%is_dev().and.(d/=szero)) call z%sync()
+    if (w%is_dev().and.(f/=szero)) call w%sync()
+    call psi_xyzw(m,a,b,c,d,e,f,x%v, y%v, z%v, w%v, info)
+    call y%set_host()
+    call z%set_host()
+    call w%set_host()
+    
+  end subroutine s_base_xyzw
+
 
   !
   !  Multiple variants of two operations:

base/modules/serial/psb_s_vect_mod.F90

@@ -103,6 +103,7 @@ module psb_s_vect_mod
     procedure, pass(z) :: axpby_a2  => s_vect_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => s_vect_abgdxyz
+    procedure, pass(z) :: xyzw     => s_vect_xyzw
     
     procedure, pass(y) :: mlt_v    => s_vect_mlt_v
     procedure, pass(y) :: mlt_a    => s_vect_mlt_a
@@ -795,6 +796,22 @@ contains
     
   end subroutine s_vect_abgdxyz
 
+  subroutine s_vect_xyzw(m,a,b,c,d,e,f,x, y, z, w, 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
+    class(psb_s_vect_type), intent(inout)  :: w
+    real(psb_spk_), intent (in)     :: a, b, c, d, e, f
+    integer(psb_ipk_), intent(out)   :: info
+
+    if (allocated(w%v)) &
+         call w%v%xyzw(m,a,b,c,d,e,f,x%v,y%v,z%v,info)
+    
+  end subroutine s_vect_xyzw
+
 
   subroutine s_vect_mlt_v(x, y, info)
     use psi_serial_mod

base/modules/serial/psb_z_base_vect_mod.F90

@@ -156,6 +156,7 @@ module psb_z_base_vect_mod
     procedure, pass(z) :: axpby_a2  => z_base_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => z_base_abgdxyz
+    procedure, pass(w) :: xyzw     => z_base_xyzw
     
     !
     ! Vector by vector multiplication. Need all variants
@@ -1155,22 +1156,37 @@ contains
     complex(psb_dpk_), intent (in)       :: alpha, beta, gamma, delta
     integer(psb_ipk_), intent(out)              :: info
 
-    if (.false.) then 
-      if (x%is_dev()) call x%sync()
-      
-      call y%axpby(m,alpha,x,beta,info)
-      call z%axpby(m,gamma,y,delta,info)
-    else
     if (x%is_dev().and.(alpha/=zzero)) call x%sync()
     if (y%is_dev().and.(beta/=zzero))   call y%sync()
     if (z%is_dev().and.(delta/=zzero))  call z%sync()
-      call psi_zabgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
+    call psi_abgdxyz(m,alpha, beta, gamma,delta,x%v, y%v, z%v, info)
     call y%set_host()
     call z%set_host()
-    end if
         
   end subroutine z_base_abgdxyz
 
+  subroutine z_base_xyzw(m,a,b,c,d,e,f,x, y, z, w,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
+    class(psb_z_base_vect_type), intent(inout)  :: w
+    complex(psb_dpk_), intent (in)                :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+
+    if (x%is_dev().and.(a/=zzero)) call x%sync()
+    if (y%is_dev().and.(b/=zzero)) call y%sync()
+    if (z%is_dev().and.(d/=zzero)) call z%sync()
+    if (w%is_dev().and.(f/=zzero)) call w%sync()
+    call psi_xyzw(m,a,b,c,d,e,f,x%v, y%v, z%v, w%v, info)
+    call y%set_host()
+    call z%set_host()
+    call w%set_host()
+    
+  end subroutine z_base_xyzw
+
 
   !
   !  Multiple variants of two operations:

base/modules/serial/psb_z_vect_mod.F90

@@ -103,6 +103,7 @@ module psb_z_vect_mod
     procedure, pass(z) :: axpby_a2  => z_vect_axpby_a2
     generic, public    :: axpby    => axpby_v, axpby_a, axpby_v2, axpby_a2
     procedure, pass(z) :: abgdxyz  => z_vect_abgdxyz
+    procedure, pass(z) :: xyzw     => z_vect_xyzw
     
     procedure, pass(y) :: mlt_v    => z_vect_mlt_v
     procedure, pass(y) :: mlt_a    => z_vect_mlt_a
@@ -788,6 +789,22 @@ contains
     
   end subroutine z_vect_abgdxyz
 
+  subroutine z_vect_xyzw(m,a,b,c,d,e,f,x, y, z, w, 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
+    class(psb_z_vect_type), intent(inout)  :: w
+    complex(psb_dpk_), intent (in)     :: a, b, c, d, e, f
+    integer(psb_ipk_), intent(out)   :: info
+
+    if (allocated(w%v)) &
+         call w%v%xyzw(m,a,b,c,d,e,f,x%v,y%v,z%v,info)
+    
+  end subroutine z_vect_xyzw
+
 
   subroutine z_vect_mlt_v(x, y, info)
     use psi_serial_mod

base/serial/psi_c_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_cabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_cabgdxyz
+
+subroutine psi_cxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  complex(psb_spk_), intent (inout)    :: z(:)
+  complex(psb_spk_), intent (inout)    :: w(:)
+  complex(psb_spk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='cabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==czero).or.(b==czero).or. &
+       & (c==czero).or.(d==czero).or.&
+       & (e==czero).or.(f==czero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_cxyzw

base/serial/psi_d_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_dabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_dabgdxyz
+
+subroutine psi_dxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  real(psb_dpk_), intent (inout)    :: z(:)
+  real(psb_dpk_), intent (inout)    :: w(:)
+  real(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='dabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==dzero).or.(b==dzero).or. &
+       & (c==dzero).or.(d==dzero).or.&
+       & (e==dzero).or.(f==dzero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_dxyzw

base/serial/psi_e_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_eabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_eabgdxyz
+
+subroutine psi_exyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  integer(psb_epk_), intent (inout)    :: z(:)
+  integer(psb_epk_), intent (inout)    :: w(:)
+  integer(psb_epk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='eabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==ezero).or.(b==ezero).or. &
+       & (c==ezero).or.(d==ezero).or.&
+       & (e==ezero).or.(f==ezero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_exyzw

base/serial/psi_i2_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_i2abgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_i2abgdxyz
+
+subroutine psi_i2xyzw(m,a,b,c,d,e,f,x, y, z,w, info)
+  use psb_const_mod
+  use psb_error_mod
+  implicit none
+  integer(psb_ipk_), intent(in)      :: m
+  integer(psb_i2pk_), intent (in)       :: x(:)
+  integer(psb_i2pk_), intent (inout)    :: y(:)
+  integer(psb_i2pk_), intent (inout)    :: z(:)
+  integer(psb_i2pk_), intent (inout)    :: w(:)
+  integer(psb_i2pk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='i2abgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==i2zero).or.(b==i2zero).or. &
+       & (c==i2zero).or.(d==i2zero).or.&
+       & (e==i2zero).or.(f==i2zero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_i2xyzw

base/serial/psi_m_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_mabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_mabgdxyz
+
+subroutine psi_mxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  integer(psb_mpk_), intent (inout)    :: z(:)
+  integer(psb_mpk_), intent (inout)    :: w(:)
+  integer(psb_mpk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='mabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==mzero).or.(b==mzero).or. &
+       & (c==mzero).or.(d==mzero).or.&
+       & (e==mzero).or.(f==mzero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_mxyzw

base/serial/psi_s_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_sabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_sabgdxyz
+
+subroutine psi_sxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  real(psb_spk_), intent (inout)    :: z(:)
+  real(psb_spk_), intent (inout)    :: w(:)
+  real(psb_spk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='sabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==szero).or.(b==szero).or. &
+       & (c==szero).or.(d==szero).or.&
+       & (e==szero).or.(f==szero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_sxyzw

base/serial/psi_z_serial_impl.F90

@@ -1792,3 +1792,75 @@ subroutine psi_zabgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
   return
 
 end subroutine psi_zabgdxyz
+
+subroutine psi_zxyzw(m,a,b,c,d,e,f,x, y, z,w, 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 (inout)    :: y(:)
+  complex(psb_dpk_), intent (inout)    :: z(:)
+  complex(psb_dpk_), intent (inout)    :: w(:)
+  complex(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+  integer(psb_ipk_), intent(out)     :: info
+
+  integer(psb_ipk_) :: i
+  integer(psb_ipk_) :: int_err(5)
+  character  name*20
+  name='zabgdxyz'
+
+  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 (size(x).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=6
+    int_err(2)=1
+    int_err(3)=size(x)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(y).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=7
+    int_err(2)=1
+    int_err(3)=size(y)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  else if (size(z).lt.max(1,m)) then
+    info=psb_err_iarg_not_gtia_ii_
+    int_err(1)=8
+    int_err(2)=1
+    int_err(3)=size(z)
+    int_err(4)=m
+    call fcpsb_errpush(info,name,int_err)
+    goto 9999
+  endif
+
+  if ((a==zzero).or.(b==zzero).or. &
+       & (c==zzero).or.(d==zzero).or.&
+       & (e==zzero).or.(f==zzero)) then
+    write(0,*) 'XYZW assumes  a,b,c,d,e,f are all nonzero'
+  else
+    !$omp parallel do private(i)
+    do i=1,m
+      y(i) = a*x(i)+b*y(i)
+      z(i) = c*y(i)+d*z(i)
+      w(i) = e*z(i)+f*w(i)
+    end do
+    
+  end if
+
+  return
+
+9999 continue
+  call fcpsb_serror()
+  return
+
+end subroutine psi_zxyzw