From 6eb6fa305ccd5aff3f0b3df0ecc5294aa8e75152 Mon Sep 17 00:00:00 2001
From: Salvatore Filippone <salvatore.filippone@cranfield.ac.uk>
Date: Thu, 2 Jul 2015 18:18:14 +0000
Subject: [PATCH] psblas:  base/serial/psi_c_serial_impl.f90 
 base/serial/psi_d_serial_impl.f90  base/serial/psi_i_serial_impl.f90 
 base/serial/psi_s_serial_impl.f90  base/serial/psi_z_serial_impl.f90

New gth for multivectors.
---
 base/serial/psi_c_serial_impl.f90 | 97 +++++++++++++++++++++++++++++--
 base/serial/psi_d_serial_impl.f90 | 97 +++++++++++++++++++++++++++++--
 base/serial/psi_i_serial_impl.f90 | 97 +++++++++++++++++++++++++++++--
 base/serial/psi_s_serial_impl.f90 | 97 +++++++++++++++++++++++++++++--
 base/serial/psi_z_serial_impl.f90 | 97 +++++++++++++++++++++++++++++--
 5 files changed, 460 insertions(+), 25 deletions(-)

diff --git a/base/serial/psi_c_serial_impl.f90 b/base/serial/psi_c_serial_impl.f90
index bacc90f2..13795f4a 100644
--- a/base/serial/psi_c_serial_impl.f90
+++ b/base/serial/psi_c_serial_impl.f90
@@ -54,6 +54,7 @@ subroutine psi_caxpby(m,n,alpha, x, beta, y, info)
 
   return
 end subroutine psi_caxpby
+
 subroutine psi_caxpbyv(m,alpha, x, beta, y, info)
   
   use psb_const_mod
@@ -105,8 +106,94 @@ subroutine psi_caxpbyv(m,alpha, x, beta, y, info)
   return
 
 end subroutine psi_caxpbyv
+
+
+subroutine psi_cgthmv(n,k,idx,alpha,x,beta,y)
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_spk_) :: x(:,:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == czero) then 
+    if (alpha == czero) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = czero
+        end do
+      end do
+    else if (alpha == cone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -cone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1                
+          y(pt) = -x(idx(i),j)
+        end do
+      end do
+    else
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  else 
+    if (beta == cone) then 
+      ! Do nothing
+    else if (beta == -cone) then 
+      y(1:n*k) = -y(1:n*k) 
+    else
+      y(1:n*k) = beta*y(1:n*k) 
+    end if
+
+    if (alpha == czero) then 
+      ! do nothing
+    else if (alpha == cone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -cone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) - x(idx(i),j)
+        end do
+      end do
+    else  
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  end if
+
+end subroutine psi_cgthmv
+
 subroutine psi_cgthv(n,idx,alpha,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
@@ -161,8 +248,8 @@ subroutine psi_cgthv(n,idx,alpha,x,beta,y)
   end if
 
 end subroutine psi_cgthv
+
 subroutine psi_cgthzmv(n,k,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -182,8 +269,8 @@ subroutine psi_cgthzmv(n,k,idx,x,y)
   end do
 
 end subroutine psi_cgthzmv
+
 subroutine psi_cgthzv(n,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -199,9 +286,9 @@ subroutine psi_cgthzv(n,idx,x,y)
   end do
 
 end subroutine psi_cgthzv
+
 subroutine psi_csctmv(n,k,idx,x,beta,y)
   
-
   use psb_const_mod
   implicit none
 
@@ -237,8 +324,8 @@ subroutine psi_csctmv(n,k,idx,x,beta,y)
     end do
   end if
 end subroutine psi_csctmv
+
 subroutine psi_csctv(n,idx,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
diff --git a/base/serial/psi_d_serial_impl.f90 b/base/serial/psi_d_serial_impl.f90
index 38031e3a..490ad385 100644
--- a/base/serial/psi_d_serial_impl.f90
+++ b/base/serial/psi_d_serial_impl.f90
@@ -54,6 +54,7 @@ subroutine psi_daxpby(m,n,alpha, x, beta, y, info)
 
   return
 end subroutine psi_daxpby
+
 subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
   
   use psb_const_mod
@@ -105,8 +106,94 @@ subroutine psi_daxpbyv(m,alpha, x, beta, y, info)
   return
 
 end subroutine psi_daxpbyv
+
+
+subroutine psi_dgthmv(n,k,idx,alpha,x,beta,y)
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_dpk_) :: x(:,:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == dzero) then 
+    if (alpha == dzero) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = dzero
+        end do
+      end do
+    else if (alpha == done) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -done) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1                
+          y(pt) = -x(idx(i),j)
+        end do
+      end do
+    else
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  else 
+    if (beta == done) then 
+      ! Do nothing
+    else if (beta == -done) then 
+      y(1:n*k) = -y(1:n*k) 
+    else
+      y(1:n*k) = beta*y(1:n*k) 
+    end if
+
+    if (alpha == dzero) then 
+      ! do nothing
+    else if (alpha == done) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -done) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) - x(idx(i),j)
+        end do
+      end do
+    else  
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  end if
+
+end subroutine psi_dgthmv
+
 subroutine psi_dgthv(n,idx,alpha,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
@@ -161,8 +248,8 @@ subroutine psi_dgthv(n,idx,alpha,x,beta,y)
   end if
 
 end subroutine psi_dgthv
+
 subroutine psi_dgthzmv(n,k,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -182,8 +269,8 @@ subroutine psi_dgthzmv(n,k,idx,x,y)
   end do
 
 end subroutine psi_dgthzmv
+
 subroutine psi_dgthzv(n,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -199,9 +286,9 @@ subroutine psi_dgthzv(n,idx,x,y)
   end do
 
 end subroutine psi_dgthzv
+
 subroutine psi_dsctmv(n,k,idx,x,beta,y)
   
-
   use psb_const_mod
   implicit none
 
@@ -237,8 +324,8 @@ subroutine psi_dsctmv(n,k,idx,x,beta,y)
     end do
   end if
 end subroutine psi_dsctmv
+
 subroutine psi_dsctv(n,idx,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
diff --git a/base/serial/psi_i_serial_impl.f90 b/base/serial/psi_i_serial_impl.f90
index 6b227983..3c3bc1c6 100644
--- a/base/serial/psi_i_serial_impl.f90
+++ b/base/serial/psi_i_serial_impl.f90
@@ -54,6 +54,7 @@ subroutine psi_iaxpby(m,n,alpha, x, beta, y, info)
 
   return
 end subroutine psi_iaxpby
+
 subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)
   
   use psb_const_mod
@@ -105,8 +106,94 @@ subroutine psi_iaxpbyv(m,alpha, x, beta, y, info)
   return
 
 end subroutine psi_iaxpbyv
+
+
+subroutine psi_igthmv(n,k,idx,alpha,x,beta,y)
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  integer(psb_ipk_) :: x(:,:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == izero) then 
+    if (alpha == izero) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = izero
+        end do
+      end do
+    else if (alpha == ione) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -ione) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1                
+          y(pt) = -x(idx(i),j)
+        end do
+      end do
+    else
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  else 
+    if (beta == ione) then 
+      ! Do nothing
+    else if (beta == -ione) then 
+      y(1:n*k) = -y(1:n*k) 
+    else
+      y(1:n*k) = beta*y(1:n*k) 
+    end if
+
+    if (alpha == izero) then 
+      ! do nothing
+    else if (alpha == ione) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -ione) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) - x(idx(i),j)
+        end do
+      end do
+    else  
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  end if
+
+end subroutine psi_igthmv
+
 subroutine psi_igthv(n,idx,alpha,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
@@ -161,8 +248,8 @@ subroutine psi_igthv(n,idx,alpha,x,beta,y)
   end if
 
 end subroutine psi_igthv
+
 subroutine psi_igthzmv(n,k,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -182,8 +269,8 @@ subroutine psi_igthzmv(n,k,idx,x,y)
   end do
 
 end subroutine psi_igthzmv
+
 subroutine psi_igthzv(n,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -199,9 +286,9 @@ subroutine psi_igthzv(n,idx,x,y)
   end do
 
 end subroutine psi_igthzv
+
 subroutine psi_isctmv(n,k,idx,x,beta,y)
   
-
   use psb_const_mod
   implicit none
 
@@ -237,8 +324,8 @@ subroutine psi_isctmv(n,k,idx,x,beta,y)
     end do
   end if
 end subroutine psi_isctmv
+
 subroutine psi_isctv(n,idx,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
diff --git a/base/serial/psi_s_serial_impl.f90 b/base/serial/psi_s_serial_impl.f90
index aaae65f6..287520a7 100644
--- a/base/serial/psi_s_serial_impl.f90
+++ b/base/serial/psi_s_serial_impl.f90
@@ -54,6 +54,7 @@ subroutine psi_saxpby(m,n,alpha, x, beta, y, info)
 
   return
 end subroutine psi_saxpby
+
 subroutine psi_saxpbyv(m,alpha, x, beta, y, info)
   
   use psb_const_mod
@@ -105,8 +106,94 @@ subroutine psi_saxpbyv(m,alpha, x, beta, y, info)
   return
 
 end subroutine psi_saxpbyv
+
+
+subroutine psi_sgthmv(n,k,idx,alpha,x,beta,y)
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  real(psb_spk_) :: x(:,:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == szero) then 
+    if (alpha == szero) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = szero
+        end do
+      end do
+    else if (alpha == sone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -sone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1                
+          y(pt) = -x(idx(i),j)
+        end do
+      end do
+    else
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  else 
+    if (beta == sone) then 
+      ! Do nothing
+    else if (beta == -sone) then 
+      y(1:n*k) = -y(1:n*k) 
+    else
+      y(1:n*k) = beta*y(1:n*k) 
+    end if
+
+    if (alpha == szero) then 
+      ! do nothing
+    else if (alpha == sone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -sone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) - x(idx(i),j)
+        end do
+      end do
+    else  
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  end if
+
+end subroutine psi_sgthmv
+
 subroutine psi_sgthv(n,idx,alpha,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
@@ -161,8 +248,8 @@ subroutine psi_sgthv(n,idx,alpha,x,beta,y)
   end if
 
 end subroutine psi_sgthv
+
 subroutine psi_sgthzmv(n,k,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -182,8 +269,8 @@ subroutine psi_sgthzmv(n,k,idx,x,y)
   end do
 
 end subroutine psi_sgthzmv
+
 subroutine psi_sgthzv(n,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -199,9 +286,9 @@ subroutine psi_sgthzv(n,idx,x,y)
   end do
 
 end subroutine psi_sgthzv
+
 subroutine psi_ssctmv(n,k,idx,x,beta,y)
   
-
   use psb_const_mod
   implicit none
 
@@ -237,8 +324,8 @@ subroutine psi_ssctmv(n,k,idx,x,beta,y)
     end do
   end if
 end subroutine psi_ssctmv
+
 subroutine psi_ssctv(n,idx,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
diff --git a/base/serial/psi_z_serial_impl.f90 b/base/serial/psi_z_serial_impl.f90
index a3842126..259c04aa 100644
--- a/base/serial/psi_z_serial_impl.f90
+++ b/base/serial/psi_z_serial_impl.f90
@@ -54,6 +54,7 @@ subroutine psi_zaxpby(m,n,alpha, x, beta, y, info)
 
   return
 end subroutine psi_zaxpby
+
 subroutine psi_zaxpbyv(m,alpha, x, beta, y, info)
   
   use psb_const_mod
@@ -105,8 +106,94 @@ subroutine psi_zaxpbyv(m,alpha, x, beta, y, info)
   return
 
 end subroutine psi_zaxpbyv
+
+
+subroutine psi_zgthmv(n,k,idx,alpha,x,beta,y)
+
+  use psb_const_mod
+  implicit none
+
+  integer(psb_ipk_) :: n, idx(:)
+  complex(psb_dpk_) :: x(:,:), y(:),alpha,beta
+
+  ! Locals
+  integer(psb_ipk_) :: i, j, pt
+
+  if (beta == zzero) then 
+    if (alpha == zzero) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = zzero
+        end do
+      end do
+    else if (alpha == zone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -zone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1                
+          y(pt) = -x(idx(i),j)
+        end do
+      end do
+    else
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  else 
+    if (beta == zone) then 
+      ! Do nothing
+    else if (beta == -zone) then 
+      y(1:n*k) = -y(1:n*k) 
+    else
+      y(1:n*k) = beta*y(1:n*k) 
+    end if
+
+    if (alpha == zzero) then 
+      ! do nothing
+    else if (alpha == zone) then 
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + x(idx(i),j)
+        end do
+      end do
+    else if (alpha == -zone) then
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) - x(idx(i),j)
+        end do
+      end do
+    else  
+      pt=0
+      do j=1,k
+        do i=1,n
+          pt=pt+1
+          y(pt) = y(pt) + alpha*x(idx(i),j)
+        end do
+      end do
+    end if
+  end if
+
+end subroutine psi_zgthmv
+
 subroutine psi_zgthv(n,idx,alpha,x,beta,y)
-  
 
   use psb_const_mod
   implicit none
@@ -161,8 +248,8 @@ subroutine psi_zgthv(n,idx,alpha,x,beta,y)
   end if
 
 end subroutine psi_zgthv
+
 subroutine psi_zgthzmv(n,k,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -182,8 +269,8 @@ subroutine psi_zgthzmv(n,k,idx,x,y)
   end do
 
 end subroutine psi_zgthzmv
+
 subroutine psi_zgthzv(n,idx,x,y)
-  
 
   use psb_const_mod
   implicit none
@@ -199,9 +286,9 @@ subroutine psi_zgthzv(n,idx,x,y)
   end do
 
 end subroutine psi_zgthzv
+
 subroutine psi_zsctmv(n,k,idx,x,beta,y)
   
-
   use psb_const_mod
   implicit none
 
@@ -237,8 +324,8 @@ subroutine psi_zsctmv(n,k,idx,x,beta,y)
     end do
   end if
 end subroutine psi_zsctmv
+
 subroutine psi_zsctv(n,idx,x,beta,y)
-  
 
   use psb_const_mod
   implicit none