2 files changed, 155 insertions, 6 deletions
diff --git a/src/kernels/level2/level2.opencl b/src/kernels/level2/level2.opencl
index ad92595a..1b0efeab 100644
--- a/src/kernels/level2/level2.opencl
+++ b/src/kernels/level2/level2.opencl
@@ -39,10 +39,7 @@ inline real LoadVector(const int id, const int max,
   if (id < max) {
     real result = gm[id*inc + offset];
     if (do_conjugate) {
-      #if defined(ROUTINE_GERC)
-        COMPLEX_CONJUGATE(result);
-      #endif
-      #if defined(ROUTINE_HER) || defined(ROUTINE_HPR)
+      #if defined(ROUTINE_GERC) || defined(ROUTINE_HER) || defined(ROUTINE_HPR) || defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
         COMPLEX_CONJUGATE(result);
       #endif
     }
@@ -81,8 +78,16 @@ inline void MatrixUpdate(const int id1, const int id2, const int max1, const int
     const real avalue = agm[a_index];
 
     // Computes result = alpha * x[i] * y[j] + a[i][j]
-    real result;
-    GER(result, alpha, xvalue, yvalue, avalue);
+    #if PRECISION == 3232 || PRECISION == 6464
+      real ax;
+      ax.x = MulReal(alpha, xvalue);
+      ax.y = MulImag(alpha, xvalue);
+      real result;
+      result.x = MulReal(ax, yvalue) + avalue.x;
+      result.y = MulImag(ax, yvalue) + avalue.y;
+    #else
+      real result = alpha * xvalue * yvalue + avalue;
+    #endif
 
     // For hermetian matrices
     #if defined(ROUTINE_HER) || defined(ROUTINE_HPR)
@@ -94,6 +99,46 @@ inline void MatrixUpdate(const int id1, const int id2, const int max1, const int
   }
 }
 
+// Performs the rank-2 matrix update
+inline void MatrixUpdate2(const int id1, const int id2, const int max1, const int max2,
+                          __global real* agm, const int a_offset, const int a_ld,
+                          const real alpha1, const real xvalue, const real yvalue,
+                          const real alpha2, const real xtvalue, const real ytvalue,
+                          const int is_upper) {
+
+  // Bounds of a regular matrix
+  if (id1 < max1 && id2 < max2) {
+
+    const int a_index = id2*a_ld + id1 + a_offset;
+
+    // Loads the current value of the A matrix
+    const real avalue = agm[a_index];
+
+    // Computes result = alpha * x[i] * y[j] + alpha * x[j] * y[i] + a[i][j]
+    #if PRECISION == 3232 || PRECISION == 6464
+      real ax;
+      ax.x = MulReal(alpha2, xvalue);
+      ax.y = MulImag(alpha2, xvalue);
+      real atx;
+      atx.x = MulReal(alpha1, xtvalue);
+      atx.y = MulImag(alpha1, xtvalue);
+      real result;
+      result.x = MulReal(ax, yvalue) + MulReal(atx, ytvalue) + avalue.x;
+      result.y = MulImag(ax, yvalue) + MulImag(atx, ytvalue) + avalue.y;
+    #else
+      real result = alpha1 * xvalue * yvalue + alpha2 * xtvalue * ytvalue + avalue;
+    #endif
+
+    // For hermetian matrices
+    #if defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
+      if (id1 == id2) { result.y = ZERO; }
+    #endif
+
+    // Stores the final result
+    agm[a_index] = result;
+  }
+}
+
 // =================================================================================================
 
 // End of the C++11 raw string literal
diff --git a/src/kernels/level2/xher2.opencl b/src/kernels/level2/xher2.opencl
new file mode 100644
index 00000000..4a2edce8
--- /dev/null
+++ b/src/kernels/level2/xher2.opencl
@@ -0,0 +1,104 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file contains the Xher2 kernels for rank-2 matrix update.
+//
+// =================================================================================================
+
+// Enables loading of this file using the C++ pre-processor's #include (C++11 standard raw string
+// literal). Comment-out this line for syntax-highlighting when developing.
+R"(
+
+// =================================================================================================
+
+// Symmetric version of the rank-2 matrix update kernel (HER2, HPR2, SYR2, SPR2)
+__attribute__((reqd_work_group_size(WGS1, WGS2, 1)))
+__kernel void Xher2(const int n, const real alpha,
+                    const __global real* restrict xgm, const int x_offset, const int x_inc,
+                    const __global real* restrict ygm, const int y_offset, const int y_inc,
+                    __global real* restrict agm, const int a_offset, const int a_ld,
+                    const int is_upper, const int is_rowmajor) {
+
+  // Register storage for X and Y
+  real xvalues[WPT];
+  real yvalues[WPT];
+  real xtvalues[WPT];
+  real ytvalues[WPT];
+
+  // Loads the X-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id2 = w*get_global_size(1) + get_global_id(1);
+    xvalues[w] = LoadVector(id2, n, xgm, x_offset, x_inc, !is_rowmajor);
+  }
+
+  // Loads the X-transposed-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id1 = w*get_global_size(0) + get_global_id(0);
+    xtvalues[w] = LoadVector(id1, n, xgm, x_offset, x_inc, is_rowmajor);
+  }
+
+  // Loads the Y-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id1 = w*get_global_size(0) + get_global_id(0);
+    yvalues[w] = LoadVector(id1, n, ygm, y_offset, y_inc, is_rowmajor);
+  }
+
+  // Loads the Y-transposed-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id2 = w*get_global_size(1) + get_global_id(1);
+    ytvalues[w] = LoadVector(id2, n, ygm, y_offset, y_inc, !is_rowmajor);
+  }
+
+  // Sets the proper value of alpha in case conjugation is needed
+  real alpha1 = alpha;
+  real alpha2 = alpha;
+  #if defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
+    if (is_rowmajor) {
+      COMPLEX_CONJUGATE(alpha1);
+    }
+    else {
+      COMPLEX_CONJUGATE(alpha2);
+    }
+  #endif
+
+  // Loops over the work per thread twice
+  #pragma unroll
+  for (int w1=0; w1<WPT; ++w1) {
+    #pragma unroll
+    for (int w2=0; w2<WPT; ++w2) {
+
+      // Global thread IDs
+      const int id1 = w1*get_global_size(0) + get_global_id(0);
+      const int id2 = w2*get_global_size(1) + get_global_id(1);
+
+      // Skip these threads if they do not contain threads contributing to the matrix-triangle
+      if ((is_upper && (id1 > id2)) || (!is_upper && (id2 > id1))) {
+        // Do nothing
+      }
+
+      // Loads A, performs the operation, and stores the result into A
+      else {
+        MatrixUpdate2(id1, id2, n, n, agm, a_offset, a_ld,
+                      alpha1, xvalues[w2], yvalues[w1],
+                      alpha2, xtvalues[w1], ytvalues[w2], is_upper);
+      }
+    }
+  }
+}
+
+// =================================================================================================
+
+// End of the C++11 raw string literal
+)"
+
+// =================================================================================================