diff --git a/cuda/cvectordev.c b/cuda/cvectordev.c
index 9db5202e..cdfda481 100644
--- a/cuda/cvectordev.c
+++ b/cuda/cvectordev.c
@@ -273,6 +273,28 @@ int abgdxyzMultiVecDeviceFloatComplex(int n,cuFloatComplex  alpha,cuFloatComplex
   return(i);
 }
 
+int xyzwMultiVecDeviceFloatComplex(int n,cuFloatComplex  a,cuFloatComplex  b,
+				   cuFloatComplex c, cuFloatComplex  d,
+				   cuFloatComplex e, cuFloatComplex  f, 
+				   void* devMultiVecX, void* devMultiVecY,
+				   void* devMultiVecZ, void* devMultiVecW)
+{ int j=0, i=0;
+  int pitch = 0;
+  struct MultiVectDevice *devVecX = (struct MultiVectDevice *) devMultiVecX;
+  struct MultiVectDevice *devVecY = (struct MultiVectDevice *) devMultiVecY;
+  struct MultiVectDevice *devVecZ = (struct MultiVectDevice *) devMultiVecZ;
+  struct MultiVectDevice *devVecW = (struct MultiVectDevice *) devMultiVecW;
+  spgpuHandle_t handle=psb_cudaGetHandle();
+  pitch = devVecY->pitch_;
+  if ((n > devVecY->size_) || (n>devVecX->size_ )) 
+    return SPGPU_UNSUPPORTED;
+
+  spgpuCxyzw(handle,n, a,b,c,d,e,f,
+	     (cuFloatComplex *)devVecX->v_,(cuFloatComplex *) devVecY->v_,
+	     (cuFloatComplex *) devVecZ->v_,(cuFloatComplex *) devVecW->v_);
+  return(i);
+}
+
 int axyMultiVecDeviceFloatComplex(int n, cuFloatComplex alpha, 
 				  void *deviceVecA, void *deviceVecB)
 { int i = 0; 
diff --git a/cuda/cvectordev.h b/cuda/cvectordev.h
index fc18e328..62693e27 100644
--- a/cuda/cvectordev.h
+++ b/cuda/cvectordev.h
@@ -72,6 +72,11 @@ int axpbyMultiVecDeviceFloatComplex(int n, cuFloatComplex alpha, void* devVecX,
 int abgdxyzMultiVecDeviceFloatComplex(int n,cuFloatComplex  alpha,cuFloatComplex  beta,
 				      cuFloatComplex  gamma, cuFloatComplex  delta, 
 				      void* devMultiVecX, void* devMultiVecY, void* devMultiVecZ);
+int xyzwMultiVecDeviceFloatComplex(int n,cuFloatComplex  a,cuFloatComplex  b,
+				   cuFloatComplex c, cuFloatComplex  d,
+				   cuFloatComplex e, cuFloatComplex  f, 
+				   void* devMultiVecX, void* devMultiVecY,
+				   void* devMultiVecZ, void* devMultiVecW);
 int axyMultiVecDeviceFloatComplex(int n, cuFloatComplex alpha, void *deviceVecA, void *deviceVecB);
 int axybzMultiVecDeviceFloatComplex(int n, cuFloatComplex alpha, void *deviceVecA,
 			      void *deviceVecB, cuFloatComplex beta, void *deviceVecZ);
diff --git a/cuda/dvectordev.c b/cuda/dvectordev.c
index b4ca95f4..723f48d8 100644
--- a/cuda/dvectordev.c
+++ b/cuda/dvectordev.c
@@ -258,6 +258,25 @@ int abgdxyzMultiVecDeviceDouble(int n,double alpha,double beta, double gamma, do
   return(i);
 }
  
+int xyzwMultiVecDeviceDouble(int n,double a, double b, double c, double d, double e, double f, 
+			     void* devMultiVecX, void* devMultiVecY,
+			     void* devMultiVecZ, void* devMultiVecW)
+{ int j=0, i=0;
+  int pitch = 0;
+  struct MultiVectDevice *devVecX = (struct MultiVectDevice *) devMultiVecX;
+  struct MultiVectDevice *devVecY = (struct MultiVectDevice *) devMultiVecY;
+  struct MultiVectDevice *devVecZ = (struct MultiVectDevice *) devMultiVecZ;
+  struct MultiVectDevice *devVecW = (struct MultiVectDevice *) devMultiVecW;
+  spgpuHandle_t handle=psb_cudaGetHandle();
+  pitch = devVecY->pitch_;
+  if ((n > devVecY->size_) || (n>devVecX->size_ )) 
+    return SPGPU_UNSUPPORTED;
+
+  spgpuDxyzw(handle,n, a,b,c,d,e,f,
+	      (double*)devVecX->v_,(double*) devVecY->v_,(double*) devVecZ->v_,(double*) devVecW->v_);
+  return(i);
+}
+ 
 int axyMultiVecDeviceDouble(int n, double alpha, void *deviceVecA, void *deviceVecB)
 { int i = 0; 
   struct MultiVectDevice *devVecA = (struct MultiVectDevice *) deviceVecA;
diff --git a/cuda/dvectordev.h b/cuda/dvectordev.h
index 81a2e8f6..c2bfa1b5 100644
--- a/cuda/dvectordev.h
+++ b/cuda/dvectordev.h
@@ -69,6 +69,9 @@ int dotMultiVecDeviceDouble(double* y_res, int n, void* devVecA, void* devVecB);
 int axpbyMultiVecDeviceDouble(int n, double alpha, void* devVecX, double beta, void* devVecY);
 int abgdxyzMultiVecDeviceDouble(int n,double alpha,double beta, double gamma, double delta, 
 				void* devMultiVecX, void* devMultiVecY, void* devMultiVecZ);
+int xyzwMultiVecDeviceDouble(int n,double a, double b, double c, double d, double e, double f, 
+			     void* devMultiVecX, void* devMultiVecY,
+			     void* devMultiVecZ, void* devMultiVecW);
 int axyMultiVecDeviceDouble(int n, double alpha, void *deviceVecA, void *deviceVecB);
 int axybzMultiVecDeviceDouble(int n, double alpha, void *deviceVecA,
 			      void *deviceVecB, double beta, void *deviceVecZ);
diff --git a/cuda/psb_c_cuda_vect_mod.F90 b/cuda/psb_c_cuda_vect_mod.F90
index 7eee128f..727249df 100644
--- a/cuda/psb_c_cuda_vect_mod.F90
+++ b/cuda/psb_c_cuda_vect_mod.F90
@@ -914,7 +914,6 @@ contains
 
   end subroutine c_cuda_axpby_v
 
-  
   subroutine c_cuda_abgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
     use psi_serial_mod
     implicit none
@@ -975,9 +974,70 @@ contains
       call z%axpby(m,gamma,y,delta,info)
     end if
 
-    
   end subroutine c_cuda_abgdxyz
 
+  subroutine c_cuda_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_vect_cuda), intent(inout)  :: w
+    complex(psb_spk_), intent (in)       :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_) :: nx, ny, nz, nw
+    logical :: gpu_done
+
+    info = psb_success_
+
+    gpu_done = .false.
+    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
+      select type(xx => x)
+      class is (psb_c_vect_cuda)
+        select type(yy => y)
+        class is (psb_c_vect_cuda)
+          select type(zz => z)
+          class is (psb_c_vect_cuda)
+            ! Do something different here 
+            if (xx%is_host()) call xx%sync()
+            if (yy%is_host()) call yy%sync()
+            if (zz%is_host()) call zz%sync()
+            if (w%is_host())  call w%sync()
+            nx = getMultiVecDeviceSize(xx%deviceVect)
+            ny = getMultiVecDeviceSize(yy%deviceVect)
+            nz = getMultiVecDeviceSize(zz%deviceVect)
+            nw = getMultiVecDeviceSize(w%deviceVect)
+            if ((nx<m).or.(ny<m).or.(nz<m).or.(nw<m)) then
+              info = psb_err_internal_error_
+            else
+              info = xyzwMultiVecDevice(m,a,b,c,d,e,f,&
+                   & xx%deviceVect,yy%deviceVect,zz%deviceVect,w%deviceVect)
+            end if
+            call yy%set_dev()
+            call zz%set_dev()
+            call w%set_dev()
+            gpu_done = .true.
+          end select
+        end select
+      end select
+      
+      if (.not.gpu_done) then 
+        if (x%is_host()) call x%sync()
+        if (y%is_host()) call y%sync()
+        if (z%is_host()) call z%sync()
+        if (w%is_host()) call w%sync()
+        call y%axpby(m,a,x,b,info)
+        call z%axpby(m,c,y,d,info)
+        call w%axpby(m,e,z,f,info)
+      end if
+    end if
+
+  end subroutine c_cuda_xyzw
   
   subroutine c_cuda_axpby_a(m,alpha, x, beta, y, info)
     use psi_serial_mod
diff --git a/cuda/psb_c_vectordev_mod.F90 b/cuda/psb_c_vectordev_mod.F90
index 20a4ac3f..a9dd1ee0 100644
--- a/cuda/psb_c_vectordev_mod.F90
+++ b/cuda/psb_c_vectordev_mod.F90
@@ -325,6 +325,18 @@ module psb_c_vectordev_mod
     end function abgdxyzMultiVecDeviceFloatComplex
   end interface
 
+  interface xyzwMultiVecDevice
+    function xyzwMultiVecDeviceFloatComplex(n,a,b,c,d,e,f,deviceVecX,&
+         & deviceVecY,deviceVecZ,deviceVecW) &
+         & result(res) bind(c,name='xyzwMultiVecDeviceFloatComplex')
+      use iso_c_binding
+      integer(c_int)      :: res
+      integer(c_int), value :: n
+      complex(c_float_complex), value :: a,b,c,d,e,f
+      type(c_ptr), value  :: deviceVecX, deviceVecY, deviceVecZ, deviceVecW
+    end function xyzwMultiVecDeviceFloatComplex
+  end interface
+
   interface axyMultiVecDevice
     function axyMultiVecDeviceFloatComplex(n,alpha,deviceVecA,deviceVecB) &
          & result(res) bind(c,name='axyMultiVecDeviceFloatComplex')
diff --git a/cuda/psb_d_cuda_vect_mod.F90 b/cuda/psb_d_cuda_vect_mod.F90
index 1e6e9f2a..b480fdf7 100644
--- a/cuda/psb_d_cuda_vect_mod.F90
+++ b/cuda/psb_d_cuda_vect_mod.F90
@@ -914,7 +914,6 @@ contains
 
   end subroutine d_cuda_axpby_v
 
-  
   subroutine d_cuda_abgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
     use psi_serial_mod
     implicit none
@@ -975,9 +974,70 @@ contains
       call z%axpby(m,gamma,y,delta,info)
     end if
 
-    
   end subroutine d_cuda_abgdxyz
 
+  subroutine d_cuda_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_vect_cuda), intent(inout)  :: w
+    real(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_) :: nx, ny, nz, nw
+    logical :: gpu_done
+
+    info = psb_success_
+
+    gpu_done = .false.
+    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
+      select type(xx => x)
+      class is (psb_d_vect_cuda)
+        select type(yy => y)
+        class is (psb_d_vect_cuda)
+          select type(zz => z)
+          class is (psb_d_vect_cuda)
+            ! Do something different here 
+            if (xx%is_host()) call xx%sync()
+            if (yy%is_host()) call yy%sync()
+            if (zz%is_host()) call zz%sync()
+            if (w%is_host())  call w%sync()
+            nx = getMultiVecDeviceSize(xx%deviceVect)
+            ny = getMultiVecDeviceSize(yy%deviceVect)
+            nz = getMultiVecDeviceSize(zz%deviceVect)
+            nw = getMultiVecDeviceSize(w%deviceVect)
+            if ((nx<m).or.(ny<m).or.(nz<m).or.(nw<m)) then
+              info = psb_err_internal_error_
+            else
+              info = xyzwMultiVecDevice(m,a,b,c,d,e,f,&
+                   & xx%deviceVect,yy%deviceVect,zz%deviceVect,w%deviceVect)
+            end if
+            call yy%set_dev()
+            call zz%set_dev()
+            call w%set_dev()
+            gpu_done = .true.
+          end select
+        end select
+      end select
+      
+      if (.not.gpu_done) then 
+        if (x%is_host()) call x%sync()
+        if (y%is_host()) call y%sync()
+        if (z%is_host()) call z%sync()
+        if (w%is_host()) call w%sync()
+        call y%axpby(m,a,x,b,info)
+        call z%axpby(m,c,y,d,info)
+        call w%axpby(m,e,z,f,info)
+      end if
+    end if
+
+  end subroutine d_cuda_xyzw
   
   subroutine d_cuda_axpby_a(m,alpha, x, beta, y, info)
     use psi_serial_mod
diff --git a/cuda/psb_d_vectordev_mod.F90 b/cuda/psb_d_vectordev_mod.F90
index 080b27fe..b6c851dd 100644
--- a/cuda/psb_d_vectordev_mod.F90
+++ b/cuda/psb_d_vectordev_mod.F90
@@ -325,6 +325,18 @@ module psb_d_vectordev_mod
     end function abgdxyzMultiVecDeviceDouble
   end interface
 
+  interface xyzwMultiVecDevice
+    function xyzwMultiVecDeviceDouble(n,a,b,c,d,e,f,deviceVecX,&
+         & deviceVecY,deviceVecZ,deviceVecW) &
+         & result(res) bind(c,name='xyzwMultiVecDeviceDouble')
+      use iso_c_binding
+      integer(c_int)      :: res
+      integer(c_int), value :: n
+      real(c_double), value :: a,b,c,d,e,f
+      type(c_ptr), value  :: deviceVecX, deviceVecY, deviceVecZ, deviceVecW
+    end function xyzwMultiVecDeviceDouble
+  end interface
+
   interface axyMultiVecDevice
     function axyMultiVecDeviceDouble(n,alpha,deviceVecA,deviceVecB) &
          & result(res) bind(c,name='axyMultiVecDeviceDouble')
diff --git a/cuda/psb_s_cuda_vect_mod.F90 b/cuda/psb_s_cuda_vect_mod.F90
index 3497c33e..244c3b72 100644
--- a/cuda/psb_s_cuda_vect_mod.F90
+++ b/cuda/psb_s_cuda_vect_mod.F90
@@ -914,7 +914,6 @@ contains
 
   end subroutine s_cuda_axpby_v
 
-  
   subroutine s_cuda_abgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
     use psi_serial_mod
     implicit none
@@ -975,9 +974,70 @@ contains
       call z%axpby(m,gamma,y,delta,info)
     end if
 
-    
   end subroutine s_cuda_abgdxyz
 
+  subroutine s_cuda_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_vect_cuda), intent(inout)  :: w
+    real(psb_spk_), intent (in)       :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_) :: nx, ny, nz, nw
+    logical :: gpu_done
+
+    info = psb_success_
+
+    gpu_done = .false.
+    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
+      select type(xx => x)
+      class is (psb_s_vect_cuda)
+        select type(yy => y)
+        class is (psb_s_vect_cuda)
+          select type(zz => z)
+          class is (psb_s_vect_cuda)
+            ! Do something different here 
+            if (xx%is_host()) call xx%sync()
+            if (yy%is_host()) call yy%sync()
+            if (zz%is_host()) call zz%sync()
+            if (w%is_host())  call w%sync()
+            nx = getMultiVecDeviceSize(xx%deviceVect)
+            ny = getMultiVecDeviceSize(yy%deviceVect)
+            nz = getMultiVecDeviceSize(zz%deviceVect)
+            nw = getMultiVecDeviceSize(w%deviceVect)
+            if ((nx<m).or.(ny<m).or.(nz<m).or.(nw<m)) then
+              info = psb_err_internal_error_
+            else
+              info = xyzwMultiVecDevice(m,a,b,c,d,e,f,&
+                   & xx%deviceVect,yy%deviceVect,zz%deviceVect,w%deviceVect)
+            end if
+            call yy%set_dev()
+            call zz%set_dev()
+            call w%set_dev()
+            gpu_done = .true.
+          end select
+        end select
+      end select
+      
+      if (.not.gpu_done) then 
+        if (x%is_host()) call x%sync()
+        if (y%is_host()) call y%sync()
+        if (z%is_host()) call z%sync()
+        if (w%is_host()) call w%sync()
+        call y%axpby(m,a,x,b,info)
+        call z%axpby(m,c,y,d,info)
+        call w%axpby(m,e,z,f,info)
+      end if
+    end if
+
+  end subroutine s_cuda_xyzw
   
   subroutine s_cuda_axpby_a(m,alpha, x, beta, y, info)
     use psi_serial_mod
diff --git a/cuda/psb_s_vectordev_mod.F90 b/cuda/psb_s_vectordev_mod.F90
index 19776cbc..a43fe0b6 100644
--- a/cuda/psb_s_vectordev_mod.F90
+++ b/cuda/psb_s_vectordev_mod.F90
@@ -325,6 +325,18 @@ module psb_s_vectordev_mod
     end function abgdxyzMultiVecDeviceFloat
   end interface
 
+  interface xyzwMultiVecDevice
+    function xyzwMultiVecDeviceFloat(n,a,b,c,d,e,f,deviceVecX,&
+         & deviceVecY,deviceVecZ,deviceVecW) &
+         & result(res) bind(c,name='xyzwMultiVecDeviceFloat')
+      use iso_c_binding
+      integer(c_int)      :: res
+      integer(c_int), value :: n
+      real(c_float), value :: a,b,c,d,e,f
+      type(c_ptr), value  :: deviceVecX, deviceVecY, deviceVecZ, deviceVecW
+    end function xyzwMultiVecDeviceFloat
+  end interface
+
   interface axyMultiVecDevice
     function axyMultiVecDeviceFloat(n,alpha,deviceVecA,deviceVecB) &
          & result(res) bind(c,name='axyMultiVecDeviceFloat')
diff --git a/cuda/psb_z_cuda_vect_mod.F90 b/cuda/psb_z_cuda_vect_mod.F90
index 8483544c..20768dd2 100644
--- a/cuda/psb_z_cuda_vect_mod.F90
+++ b/cuda/psb_z_cuda_vect_mod.F90
@@ -914,7 +914,6 @@ contains
 
   end subroutine z_cuda_axpby_v
 
-  
   subroutine z_cuda_abgdxyz(m,alpha, beta, gamma,delta,x, y, z, info)
     use psi_serial_mod
     implicit none
@@ -975,9 +974,70 @@ contains
       call z%axpby(m,gamma,y,delta,info)
     end if
 
-    
   end subroutine z_cuda_abgdxyz
 
+  subroutine z_cuda_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_vect_cuda), intent(inout)  :: w
+    complex(psb_dpk_), intent (in)       :: a,b,c,d,e,f
+    integer(psb_ipk_), intent(out)              :: info
+    integer(psb_ipk_) :: nx, ny, nz, nw
+    logical :: gpu_done
+
+    info = psb_success_
+
+    gpu_done = .false.
+    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
+      select type(xx => x)
+      class is (psb_z_vect_cuda)
+        select type(yy => y)
+        class is (psb_z_vect_cuda)
+          select type(zz => z)
+          class is (psb_z_vect_cuda)
+            ! Do something different here 
+            if (xx%is_host()) call xx%sync()
+            if (yy%is_host()) call yy%sync()
+            if (zz%is_host()) call zz%sync()
+            if (w%is_host())  call w%sync()
+            nx = getMultiVecDeviceSize(xx%deviceVect)
+            ny = getMultiVecDeviceSize(yy%deviceVect)
+            nz = getMultiVecDeviceSize(zz%deviceVect)
+            nw = getMultiVecDeviceSize(w%deviceVect)
+            if ((nx<m).or.(ny<m).or.(nz<m).or.(nw<m)) then
+              info = psb_err_internal_error_
+            else
+              info = xyzwMultiVecDevice(m,a,b,c,d,e,f,&
+                   & xx%deviceVect,yy%deviceVect,zz%deviceVect,w%deviceVect)
+            end if
+            call yy%set_dev()
+            call zz%set_dev()
+            call w%set_dev()
+            gpu_done = .true.
+          end select
+        end select
+      end select
+      
+      if (.not.gpu_done) then 
+        if (x%is_host()) call x%sync()
+        if (y%is_host()) call y%sync()
+        if (z%is_host()) call z%sync()
+        if (w%is_host()) call w%sync()
+        call y%axpby(m,a,x,b,info)
+        call z%axpby(m,c,y,d,info)
+        call w%axpby(m,e,z,f,info)
+      end if
+    end if
+
+  end subroutine z_cuda_xyzw
   
   subroutine z_cuda_axpby_a(m,alpha, x, beta, y, info)
     use psi_serial_mod
diff --git a/cuda/psb_z_vectordev_mod.F90 b/cuda/psb_z_vectordev_mod.F90
index 07e4ba37..320f362a 100644
--- a/cuda/psb_z_vectordev_mod.F90
+++ b/cuda/psb_z_vectordev_mod.F90
@@ -325,6 +325,18 @@ module psb_z_vectordev_mod
     end function abgdxyzMultiVecDeviceDoubleComplex
   end interface
 
+  interface xyzwMultiVecDevice
+    function xyzwMultiVecDeviceDoubleComplex(n,a,b,c,d,e,f,deviceVecX,&
+         & deviceVecY,deviceVecZ,deviceVecW) &
+         & result(res) bind(c,name='xyzwMultiVecDeviceDoubleComplex')
+      use iso_c_binding
+      integer(c_int)      :: res
+      integer(c_int), value :: n
+      complex(c_double_complex), value :: a,b,c,d,e,f
+      type(c_ptr), value  :: deviceVecX, deviceVecY, deviceVecZ, deviceVecW
+    end function xyzwMultiVecDeviceDoubleComplex
+  end interface
+
   interface axyMultiVecDevice
     function axyMultiVecDeviceDoubleComplex(n,alpha,deviceVecA,deviceVecB) &
          & result(res) bind(c,name='axyMultiVecDeviceDoubleComplex')
diff --git a/cuda/spgpu/kernels/Makefile b/cuda/spgpu/kernels/Makefile
index 3e668b4e..c4ae38a0 100644
--- a/cuda/spgpu/kernels/Makefile
+++ b/cuda/spgpu/kernels/Makefile
@@ -19,7 +19,8 @@ OBJS=cabs.o  camax.o  casum.o  caxpby.o  caxy.o  cdot.o  cgath.o \
   hdia_cspmv.o hdia_dspmv.o hdia_sspmv.o hdia_zspmv.o hell_cspmv.o hell_dspmv.o \
   hell_sspmv.o hell_zspmv.o igath.o iscat.o isetscal.o sabs.o samax.o sasum.o \
   saxpby.o saxy.o sdot.o sgath.o snrm2.o sscal.o sscat.o ssetscal.o zabs.o zamax.o sabgdxyz.o\
-  zasum.o zaxpby.o zaxy.o zdot.o zgath.o znrm2.o zscal.o zscat.o zsetscal.o zabgdxyz.o
+  zasum.o zaxpby.o zaxy.o zdot.o zgath.o znrm2.o zscal.o zscat.o zsetscal.o zabgdxyz.o \
+  sxyzw.o cxyzw.o dxyzw.o zxyzw.o
 
 objs: $(OBJS)
 lib: objs
diff --git a/cuda/spgpu/vector.h b/cuda/spgpu/vector.h
index 9fc3e658..d08e6edd 100644
--- a/cuda/spgpu/vector.h
+++ b/cuda/spgpu/vector.h
@@ -192,6 +192,17 @@ void spgpuSabgdxyz(spgpuHandle_t handle,
 		   __device float *y,
 		   __device float *z)
 ;
+
+  void spgpuSxyzw(spgpuHandle_t handle,
+		  int n,
+		  float a, float b,
+		  float c, float d,
+		  float e, float f,
+		   __device float* x,
+		   __device float *y,
+		  __device float *z,
+		  __device float *w)
+;
 	
 /** 
 * \fn void spgpuSmaxpby(spgpuHandle_t handle, __device float *z, int n, float beta, __device float *y, float alpha, __device float* x, int count, int pitch)
@@ -487,6 +498,16 @@ void spgpuDabgdxyz(spgpuHandle_t handle,
 		   __device double *z)
 ;
 	
+  void spgpuDxyzw(spgpuHandle_t handle,
+		  int n,
+		  double a, double b,
+		  double c, double d,
+		  double e, double f,
+		  __device double* x,
+		  __device double *y,
+		  __device double *z,
+		  __device double *w)
+;
 
   /** 
 * \fn void spgpuDmaxpby(spgpuHandle_t handle, __device double *z, int n, double beta, __device double *y, double alpha, __device double* x, int count, int pitch)
@@ -778,6 +799,18 @@ void spgpuCabgdxyz(spgpuHandle_t handle,
 		   __device cuFloatComplex *y,
 		   __device cuFloatComplex *z)
 ;
+  
+  void spgpuCxyzw(spgpuHandle_t handle,
+		  int n,
+		  cuFloatComplex a, cuFloatComplex b,
+		  cuFloatComplex c, cuFloatComplex d,
+		  cuFloatComplex e, cuFloatComplex f,
+		   __device cuFloatComplex* x,
+		   __device cuFloatComplex *y,
+		  __device cuFloatComplex *z,
+		  __device cuFloatComplex *w)
+;
+
 	
 /** 
 * \fn void spgpuCmaxpby(spgpuHandle_t handle, __device cuFloatComplex *z, int n, cuFloatComplex beta, __device cuFloatComplex *y, cuFloatComplex alpha, __device cuFloatComplex* x, int count, int pitch)
@@ -1069,6 +1102,19 @@ void spgpuZabgdxyz(spgpuHandle_t handle,
 		   __device cuDoubleComplex *y,
 		   __device cuDoubleComplex *z)
 ;
+
+  
+  void spgpuZxyzw(spgpuHandle_t handle,
+		  int n,
+		  cuDoubleComplex a, cuDoubleComplex b,
+		  cuDoubleComplex c, cuDoubleComplex d,
+		  cuDoubleComplex e, cuDoubleComplex f,
+		   __device cuDoubleComplex* x,
+		   __device cuDoubleComplex *y,
+		  __device cuDoubleComplex *z,
+		  __device cuDoubleComplex *w)
+;
+  
 	
 /** 
 * \fn void spgpuZmaxpby(spgpuHandle_t handle, __device cuDoubleComplex *z, int n, cuDoubleComplex beta, __device cuDoubleComplex *y, cuDoubleComplex alpha, __device cuDoubleComplex* x, int count, int pitch)
diff --git a/cuda/svectordev.c b/cuda/svectordev.c
index b84718f5..bf7545bb 100644
--- a/cuda/svectordev.c
+++ b/cuda/svectordev.c
@@ -258,6 +258,27 @@ int abgdxyzMultiVecDeviceFloat(int n,float alpha,float beta, float gamma, float
   return(i);
 }
  
+
+int xyzwMultiVecDeviceFloat(int n,float a,float b, float c, float d, float e, float f,
+			    void* devMultiVecX, void* devMultiVecY,
+			    void* devMultiVecZ, void* devMultiVecW)
+{ int j=0, i=0;
+  int pitch = 0;
+  struct MultiVectDevice *devVecX = (struct MultiVectDevice *) devMultiVecX;
+  struct MultiVectDevice *devVecY = (struct MultiVectDevice *) devMultiVecY;
+  struct MultiVectDevice *devVecZ = (struct MultiVectDevice *) devMultiVecZ;
+  struct MultiVectDevice *devVecW = (struct MultiVectDevice *) devMultiVecW;
+  spgpuHandle_t handle=psb_cudaGetHandle();
+  pitch = devVecY->pitch_;
+  if ((n > devVecY->size_) || (n>devVecX->size_ )) 
+    return SPGPU_UNSUPPORTED;
+
+  spgpuSxyzw(handle,n, a,b,c,d,e,f,
+	      (float*)devVecX->v_,(float*) devVecY->v_,
+	      (float*) devVecZ->v_,(float*) devVecW->v_);
+  return(i);
+}
+ 
 int axyMultiVecDeviceFloat(int n, float alpha, void *deviceVecA, void *deviceVecB)
 { int i = 0; 
   struct MultiVectDevice *devVecA = (struct MultiVectDevice *) deviceVecA;
diff --git a/cuda/svectordev.h b/cuda/svectordev.h
index 730f929a..363c0108 100644
--- a/cuda/svectordev.h
+++ b/cuda/svectordev.h
@@ -69,6 +69,9 @@ int dotMultiVecDeviceFloat(float* y_res, int n, void* devVecA, void* devVecB);
 int axpbyMultiVecDeviceFloat(int n, float alpha, void* devVecX, float beta, void* devVecY);
 int abgdxyzMultiVecDeviceFloat(int n,float alpha,float beta, float gamma, float delta, 
 			       void* devMultiVecX, void* devMultiVecY, void* devMultiVecZ);
+int xyzwMultiVecDeviceFloat(int n,float a,float b, float c, float d, float e, float f, 
+			    void* devMultiVecX, void* devMultiVecY,
+			    void* devMultiVecZ, void* devMultiVecW);
 int axyMultiVecDeviceFloat(int n, float alpha, void *deviceVecA, void *deviceVecB);
 int axybzMultiVecDeviceFloat(int n, float alpha, void *deviceVecA,
 			      void *deviceVecB, float beta, void *deviceVecZ);
diff --git a/cuda/zvectordev.c b/cuda/zvectordev.c
index d1f23f2a..e9f0cec7 100644
--- a/cuda/zvectordev.c
+++ b/cuda/zvectordev.c
@@ -251,7 +251,29 @@ int abgdxyzMultiVecDeviceDoubleComplex(int n,cuDoubleComplex  alpha,
 	      (cuDoubleComplex *)devVecX->v_,(cuDoubleComplex *) devVecY->v_,(cuDoubleComplex *) devVecZ->v_);
   return(i);
 }
- 
+
+int xyzwMultiVecDeviceDoubleComplex(int n,cuDoubleComplex  a,   cuDoubleComplex  b,
+				    cuDoubleComplex  c, cuDoubleComplex  d,
+				    cuDoubleComplex  e, cuDoubleComplex  f, 
+				    void* devMultiVecX, void* devMultiVecY,
+				    void* devMultiVecZ, void* devMultiVecW)  
+{ int j=0, i=0;
+  int pitch = 0;
+  struct MultiVectDevice *devVecX = (struct MultiVectDevice *) devMultiVecX;
+  struct MultiVectDevice *devVecY = (struct MultiVectDevice *) devMultiVecY;
+  struct MultiVectDevice *devVecZ = (struct MultiVectDevice *) devMultiVecZ;
+  struct MultiVectDevice *devVecW = (struct MultiVectDevice *) devMultiVecW;
+  spgpuHandle_t handle=psb_cudaGetHandle();
+  pitch = devVecY->pitch_;
+  if ((n > devVecY->size_) || (n>devVecX->size_ )) 
+    return SPGPU_UNSUPPORTED;
+
+  spgpuZxyzw(handle,n, a,b,c,d,e,f,
+	      (cuDoubleComplex *)devVecX->v_,(cuDoubleComplex *) devVecY->v_,
+	      (cuDoubleComplex *) devVecZ->v_,(cuDoubleComplex *) devVecW->v_);
+  return(i);
+}
+
 int axpbyMultiVecDeviceDoubleComplex(int n,cuDoubleComplex alpha, void* devMultiVecX, 
 			      cuDoubleComplex beta, void* devMultiVecY)
 { int j=0, i=0;
diff --git a/cuda/zvectordev.h b/cuda/zvectordev.h
index 4c32f11c..ae623bdb 100644
--- a/cuda/zvectordev.h
+++ b/cuda/zvectordev.h
@@ -80,6 +80,11 @@ int axpbyMultiVecDeviceDoubleComplex(int n, cuDoubleComplex alpha, void* devVecX
 int abgdxyzMultiVecDeviceDoubleComplex(int n,cuDoubleComplex  alpha,
 				       cuDoubleComplex  beta, cuDoubleComplex  gamma, cuDoubleComplex  delta, 
 				       void* devMultiVecX, void* devMultiVecY, void* devMultiVecZ);
+int xyzwMultiVecDeviceDoubleComplex(int n,cuDoubleComplex  a, cuDoubleComplex  b,
+				    cuDoubleComplex  c, cuDoubleComplex  d,
+				    cuDoubleComplex  e, cuDoubleComplex  f, 
+				    void* devMultiVecX, void* devMultiVecY,
+				    void* devMultiVecZ, void* devMultiVecW);
 int axyMultiVecDeviceDoubleComplex(int n, cuDoubleComplex alpha, 
 				   void *deviceVecA, void *deviceVecB);
 int axybzMultiVecDeviceDoubleComplex(int n, cuDoubleComplex alpha, void *deviceVecA,