diff --git a/plugins/local/gpu_nvidia/gpu.c b/plugins/local/gpu_nvidia/gpu.c
index e77847a6..a775ab95 100644
--- a/plugins/local/gpu_nvidia/gpu.c
+++ b/plugins/local/gpu_nvidia/gpu.c
@@ -116,11 +116,6 @@ void gpu_ddot(cublasHandle_t handle, const int64_t n, const double* x, const int
   assert ( (int64_t) incy_ == incy);
 
   cublasStatus_t rc = cublasDdot(handle, n_, x, incx_, y, incy_, result);
-/*
-  double alpha = 1.0;
-  double beta = 0.0;
-  cublasStatus_t rc = cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 1, 1, n_, &alpha, x, 1, y, n_, &beta, &result_, 1);
-*/
   assert (rc == CUBLAS_STATUS_SUCCESS);
 }
 
@@ -149,8 +144,8 @@ void gpu_sdot(cublasHandle_t handle, const int64_t n, const float* x, const int6
 
 
 
-void gpu_dgemv(cublasHandle_t handle, const char* transa, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double* x, const int64_t incx, const double beta, double* y, const int64_t incy) {
+void gpu_dgemv(cublasHandle_t handle, const char* transa, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* x, const int64_t incx, const double* beta, double* y, const int64_t incy) {
 
   assert (handle != NULL);
 
@@ -173,13 +168,13 @@ void gpu_dgemv(cublasHandle_t handle, const char* transa, const int64_t m, const
   cublasOperation_t transa_ = CUBLAS_OP_N;
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
 
-  cublasDgemv(handle, transa_, m_, n_, &alpha, a, lda_, x, incx_, &beta, y, incy_);
+  cublasDgemv(handle, transa_, m_, n_, alpha, a, lda_, x, incx_, beta, y, incy_);
 }
 
 
 
-void gpu_sgemv(cublasHandle_t handle, const char* transa, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float* x, const int64_t incx, const float beta, float* y, const int64_t incy) {
+void gpu_sgemv(cublasHandle_t handle, const char* transa, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* x, const int64_t incx, const float* beta, float* y, const int64_t incy) {
 
   assert (handle != NULL);
 
@@ -202,12 +197,12 @@ void gpu_sgemv(cublasHandle_t handle, const char* transa, const int64_t m, const
   cublasOperation_t transa_ = CUBLAS_OP_N;
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
 
-  cublasSgemv(handle, transa_, m_, n_, &alpha, a, lda_, x, incx_, &beta, y, incy_);
+  cublasSgemv(handle, transa_, m_, n_, alpha, a, lda_, x, incx_, beta, y, incy_);
 }
 
 
-void gpu_dgemm(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const double alpha,
-               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double beta, double* c, const int64_t ldc) {
+void gpu_dgemm(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const double* alpha,
+               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double* beta, double* c, const int64_t ldc) {
 
   assert (handle != NULL);
 
@@ -234,13 +229,13 @@ void gpu_dgemm(cublasHandle_t handle, const char* transa, const char* transb, co
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
   if (*transb == 'T' || *transb == 't') transb_ = CUBLAS_OP_T;
 
-  cublasDgemm(handle, transa_, transb_, m_, n_, k_, &alpha, a, lda_, b, ldb_, &beta, c, ldc_);
+  cublasDgemm(handle, transa_, transb_, m_, n_, k_, alpha, a, lda_, b, ldb_, beta, c, ldc_);
 }
 
 
 
-void gpu_sgemm(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const float alpha,
-               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float beta, float* c, const int64_t ldc) {
+void gpu_sgemm(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const float* alpha,
+               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float* beta, float* c, const int64_t ldc) {
 
   assert (handle != NULL);
 
@@ -267,12 +262,12 @@ void gpu_sgemm(cublasHandle_t handle, const char* transa, const char* transb, co
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
   if (*transb == 'T' || *transb == 't') transb_ = CUBLAS_OP_T;
 
-  cublasSgemm(handle, transa_, transb_, m_, n_, k_, &alpha, a, lda_, b, ldb_, &beta, c, ldc_);
+  cublasSgemm(handle, transa_, transb_, m_, n_, k_, alpha, a, lda_, b, ldb_, beta, c, ldc_);
 }
 
 
-void gpu_dgeam(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double beta, const double* b, const int64_t ldb, double* c, const int64_t ldc) {
+void gpu_dgeam(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* beta, const double* b, const int64_t ldb, double* c, const int64_t ldc) {
   assert (handle != NULL);
 
   /* Convert to int */
@@ -296,13 +291,13 @@ void gpu_dgeam(cublasHandle_t handle, const char* transa, const char* transb, co
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
   if (*transb == 'T' || *transb == 't') transb_ = CUBLAS_OP_T;
 
-  cublasDgeam(handle, transa_, transb_, m_, n_, &alpha, a, lda_, &beta, b, ldb_, c, ldc_);
+  cublasDgeam(handle, transa_, transb_, m_, n_, alpha, a, lda_, beta, b, ldb_, c, ldc_);
 
 }
 
 
-void gpu_sgeam(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float beta, const float* b, const int64_t ldb, float* c, const int64_t ldc) {
+void gpu_sgeam(cublasHandle_t handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* beta, const float* b, const int64_t ldb, float* c, const int64_t ldc) {
   assert (handle != NULL);
 
   /* Convert to int */
@@ -326,6 +321,6 @@ void gpu_sgeam(cublasHandle_t handle, const char* transa, const char* transb, co
   if (*transa == 'T' || *transa == 't') transa_ = CUBLAS_OP_T;
   if (*transb == 'T' || *transb == 't') transb_ = CUBLAS_OP_T;
 
-  cublasSgeam(handle, transa_, transb_, m_, n_, &alpha, a, lda_, &beta, b, ldb_, c, ldc_);
+  cublasSgeam(handle, transa_, transb_, m_, n_, alpha, a, lda_, beta, b, ldb_, c, ldc_);
 
 }
diff --git a/plugins/local/gpu_x86/gpu.c b/plugins/local/gpu_x86/gpu.c
index fe3cadc5..49aec9d3 100644
--- a/plugins/local/gpu_x86/gpu.c
+++ b/plugins/local/gpu_x86/gpu.c
@@ -124,8 +124,8 @@ void gpu_sdot(void* handle, const int64_t n, const float* x, const int64_t incx,
 void dgemv_(const char* transa, const int32_t* m, const int32_t* n, const double* alpha,
             const double* a, const int32_t* lda, const double* x, const int32_t* incx, const double* beta, double* y, const int32_t* incy);
 
-void gpu_dgemv(void* handle, const char* transa, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double* x, const int64_t incx, const double beta, double* y, const int64_t incy) {
+void gpu_dgemv(void* handle, const char* transa, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* x, const int64_t incx, const double* beta, double* y, const int64_t incy) {
 
   assert (handle != NULL);
 
@@ -145,15 +145,15 @@ void gpu_dgemv(void* handle, const char* transa, const int64_t m, const int64_t
   assert ( (int64_t) incx_ == incx);
   assert ( (int64_t) incy_ == incy);
 
-  dgemv_(transa, &m_, &n_, &alpha, a, &lda_, x, &incx_, &beta, y, &incy_);
+  dgemv_(transa, &m_, &n_, alpha, a, &lda_, x, &incx_, beta, y, &incy_);
 }
 
 
 void sgemv_(const char* transa, const int32_t* m, const int32_t* n, const float* alpha,
                const float* a, const int32_t* lda, const float* x, const int32_t* incx, const float* beta, float* y, const int32_t* incy);
 
-void gpu_sgemv(void* handle, const char* transa, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float* x, const int64_t incx, const float beta, float* y, const int64_t incy) {
+void gpu_sgemv(void* handle, const char* transa, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* x, const int64_t incx, const float* beta, float* y, const int64_t incy) {
 
   assert (handle != NULL);
 
@@ -173,15 +173,15 @@ void gpu_sgemv(void* handle, const char* transa, const int64_t m, const int64_t
   assert ( (int64_t) incx_ == incx);
   assert ( (int64_t) incy_ == incy);
 
-  sgemv_(transa, &m_, &n_, &alpha, a, &lda_, x, &incx_, &beta, y, &incy_);
+  sgemv_(transa, &m_, &n_, alpha, a, &lda_, x, &incx_, beta, y, &incy_);
 }
 
 
 void dgemm_(const char* transa, const char* transb, const int32_t* m, const int32_t* n, const int32_t* k, const double* alpha,
             const double* a, const int32_t* lda, const double* b, const int32_t* ldb, const double* beta, double* c, const int32_t* ldc);
 
-void gpu_dgemm(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const double alpha,
-               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double beta, double* c, const int64_t ldc) {
+void gpu_dgemm(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const double* alpha,
+               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double* beta, double* c, const int64_t ldc) {
 
   assert (handle != NULL);
 
@@ -203,7 +203,7 @@ void gpu_dgemm(void* handle, const char* transa, const char* transb, const int64
   assert ( (int64_t) ldb_ == ldb);
   assert ( (int64_t) ldc_ == ldc);
 
-  dgemm_(transa, transb, &m_, &n_, &k_, &alpha, a, &lda_, b, &ldb_, &beta, c, &ldc_);
+  dgemm_(transa, transb, &m_, &n_, &k_, alpha, a, &lda_, b, &ldb_, beta, c, &ldc_);
 }
 
 
@@ -211,8 +211,8 @@ void gpu_dgemm(void* handle, const char* transa, const char* transb, const int64
 void sgemm_(const char* transa, const char* transb, const int32_t* m, const int32_t* n, const int32_t* k, const float* alpha,
             const float* a, const int32_t* lda, const float* b, const int32_t* ldb, const float* beta, float* c, const int32_t* ldc);
 
-void gpu_sgemm(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const float alpha,
-               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float beta, float* c, const int64_t ldc) {
+void gpu_sgemm(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const int64_t k, const float* alpha,
+               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float* beta, float* c, const int64_t ldc) {
 
   assert (handle != NULL);
 
@@ -234,12 +234,12 @@ void gpu_sgemm(void* handle, const char* transa, const char* transb, const int64
   assert ( (int64_t) ldb_ == ldb);
   assert ( (int64_t) ldc_ == ldc);
 
-  sgemm_(transa, transb, &m_, &n_, &k_, &alpha, a, &lda_, b, &ldb_, &beta, c, &ldc_);
+  sgemm_(transa, transb, &m_, &n_, &k_, alpha, a, &lda_, b, &ldb_, beta, c, &ldc_);
 }
 
 
-void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double beta, const double* b, const int64_t ldb, double* c, const int64_t ldc) {
+void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* beta, const double* b, const int64_t ldb, double* c, const int64_t ldc) {
   assert (handle != NULL);
 
   if ( (*transa == 'N' && *transb == 'N') ||
@@ -247,19 +247,19 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
        (*transa == 'N' && *transb == 'n') ||
        (*transa == 'n' && *transb == 'n') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[j*ldb+i];
+           c[j*ldc+i] = *beta * b[j*ldb+i];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i];
          }
        }
 
@@ -267,7 +267,7 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i] + beta * b[j*ldb+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i] + *beta * b[j*ldb+i];
          }
        }
 
@@ -278,19 +278,19 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'N' && *transb == 't') ||
               (*transa == 'n' && *transb == 't') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[i*ldb+j];
+           c[j*ldc+i] = *beta * b[i*ldb+j];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i];
          }
        }
 
@@ -298,7 +298,7 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i] + beta * b[i*ldb+j];
+           c[j*ldc+i] = *alpha * a[j*lda+i] + *beta * b[i*ldb+j];
          }
        }
 
@@ -309,19 +309,19 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'T' && *transb == 'n') ||
               (*transa == 't' && *transb == 'n') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[j*ldb+i];
+           c[j*ldc+i] = *beta * b[j*ldb+i];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j];
          }
        }
 
@@ -329,7 +329,7 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j] + beta * b[j*ldb+i];
+           c[j*ldc+i] = *alpha * a[i*lda+j] + *beta * b[j*ldb+i];
          }
        }
 
@@ -340,19 +340,19 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'T' && *transb == 't') ||
               (*transa == 't' && *transb == 't') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[i*ldb+j];
+           c[j*ldc+i] = *beta * b[i*ldb+j];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j];
          }
        }
 
@@ -360,7 +360,7 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j] + beta * b[i*ldb+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j] + *beta * b[i*ldb+j];
          }
        }
 
@@ -370,8 +370,8 @@ void gpu_dgeam(void* handle, const char* transa, const char* transb, const int64
 }
 
 
-void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float beta, const float* b, const int64_t ldb, float* c, const int64_t ldc) {
+void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* beta, const float* b, const int64_t ldb, float* c, const int64_t ldc) {
   assert (handle != NULL);
 
   if ( (*transa == 'N' && *transb == 'N') ||
@@ -379,19 +379,19 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
        (*transa == 'N' && *transb == 'n') ||
        (*transa == 'n' && *transb == 'n') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[j*ldb+i];
+           c[j*ldc+i] = *beta * b[j*ldb+i];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i];
          }
        }
 
@@ -399,7 +399,7 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i] + beta * b[j*ldb+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i] + *beta * b[j*ldb+i];
          }
        }
 
@@ -410,19 +410,19 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'N' && *transb == 't') ||
               (*transa == 'n' && *transb == 't') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[i*ldb+j];
+           c[j*ldc+i] = *beta * b[i*ldb+j];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i];
+           c[j*ldc+i] = *alpha * a[j*lda+i];
          }
        }
 
@@ -430,7 +430,7 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[j*lda+i] + beta * b[i*ldb+j];
+           c[j*ldc+i] = *alpha * a[j*lda+i] + *beta * b[i*ldb+j];
          }
        }
 
@@ -441,19 +441,19 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'T' && *transb == 'n') ||
               (*transa == 't' && *transb == 'n') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[j*ldb+i];
+           c[j*ldc+i] = *beta * b[j*ldb+i];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j];
          }
        }
 
@@ -461,7 +461,7 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j] + beta * b[j*ldb+i];
+           c[j*ldc+i] = *alpha * a[i*lda+j] + *beta * b[j*ldb+i];
          }
        }
 
@@ -472,19 +472,19 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
               (*transa == 'T' && *transb == 't') ||
               (*transa == 't' && *transb == 't') ) {
 
-     if (alpha == 0.) {
+     if (*alpha == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = beta * b[i*ldb+j];
+           c[j*ldc+i] = *beta * b[i*ldb+j];
          }
        }
 
-     } else if (beta == 0.) {
+     } else if (*beta == 0.) {
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j];
          }
        }
 
@@ -492,7 +492,7 @@ void gpu_sgeam(void* handle, const char* transa, const char* transb, const int64
 
        for (int64_t j=0 ; j<n ; ++j) {
          for (int64_t i=0 ; i<m ; ++i) {
-           c[j*ldc+i] = alpha * a[i*lda+j] + beta * b[i*ldb+j];
+           c[j*ldc+i] = *alpha * a[i*lda+j] + *beta * b[i*ldb+j];
          }
        }
 
diff --git a/src/ccsd/ccsd_space_orb_sub_chol.irp.f b/src/ccsd/ccsd_space_orb_sub_chol.irp.f
index e0048637..c34b390b 100644
--- a/src/ccsd/ccsd_space_orb_sub_chol.irp.f
+++ b/src/ccsd/ccsd_space_orb_sub_chol.irp.f
@@ -217,21 +217,18 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
   ! internal
   integer                       :: u,i,j,beta,a,b
 
-  call gpu_copy(d_cc_space_f_ov, r1)
+  type(gpu_stream) :: stream(nV)
+
+  do a=1,nV
+    call gpu_stream_create(stream(a))
+  enddo
 
   type(gpu_double2) :: X_oo
   call gpu_allocate(X_oo,nO,nO)
 
-  call gpu_dgemm(blas_handle, 'N','N', nO, nO, nV, &
-             -2d0, t1%f(1,1), size(t1%f,1), &
-                   d_cc_space_f_vo%f(1,1), size(d_cc_space_f_vo%f,1), &
-              0d0, X_oo%f(1,1), size(X_oo%f,1))
-
-  call gpu_dgemm(blas_handle, 'T','N', nO, nV, nO, &
-             1d0, X_oo%f(1,1), size(X_oo%f,2), &
-                  t1%f(1,1)  , size(t1%f,1), &
-             1d0, r1%f(1,1)  , size(r1%f,1))
+  call gpu_copy(d_cc_space_f_ov, r1)
 
+  call gpu_set_stream(blas_handle, stream(1))
   call gpu_dgemm(blas_handle, 'N','N', nO, nV, nV, &
              1d0, t1%f(1,1)  , size(t1%f,1), &
                   H_vv%f(1,1), size(H_vv%f,1), &
@@ -242,35 +239,32 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
                    t1%f(1,1)  , size(t1%f,1), &
               1d0, r1%f(1,1), size(r1%f,1))
 
+  call gpu_set_stream(blas_handle, stream(nV))
+  call gpu_dgemm(blas_handle, 'N','N', nO, nO, nV, &
+             -2d0, t1%f(1,1), size(t1%f,1), &
+                   d_cc_space_f_vo%f(1,1), size(d_cc_space_f_vo%f,1), &
+              0d0, X_oo%f(1,1), size(X_oo%f,1))
+
+  call gpu_synchronize()
+  call gpu_set_stream(blas_handle, gpu_default_stream)
+
+  call gpu_dgemm(blas_handle, 'T','N', nO, nV, nO, &
+             1d0, X_oo%f(1,1), size(X_oo%f,2), &
+                  t1%f(1,1)  , size(t1%f,1), &
+             1d0, r1%f(1,1)  , size(r1%f,1))
+
+
+
   type(gpu_double4) :: X_voov
   call gpu_allocate(X_voov, nV, nO, nO, nV)
 
-  type(gpu_stream) :: stream(nV)
-
-  do a=1,nV
-    call gpu_stream_create(stream(a))
-  enddo
-
-  call gpu_synchronize()
-
-!  do i=1,nO
-!    do beta=1,nV
-!      call gpu_set_stream(blas_handle, stream(beta))
-!      call gpu_dgeam(blas_handle, 'T', 'T', nV, nO,  -1.d0, t2%f(1,i,1,beta), &
-!         nO*nO, t1%f(i,beta), t1%f(1,1), nO, X_voov%f(1,i,1,beta), nV)
-!    enddo
-!  enddo
-
-  do beta = 1, nV
-    do u = 1, nO
-      do i = 1, nO
-        do a = 1, nV
-          X_voov%f(a,i,u,beta) = - t2%f(u,i,a,beta) + t1%f(u,a) * t1%f(i,beta)
-        enddo
-      enddo
+  do i=1,nO
+    do beta=1,nV
+      call gpu_set_stream(blas_handle, stream(beta))
+      call gpu_dgeam(blas_handle, 'T', 'T', nV, nO,  -1.d0, t2%f(1,i,1,beta), &
+         nO*nO, t1%f(i,beta), t1%f(1,1), nO, X_voov%f(1,i,1,beta), nV*nO)
     enddo
   enddo
-  call gpu_synchronize()
 
   do beta=1,nV
     call gpu_set_stream(blas_handle, stream(beta))
@@ -279,9 +273,8 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
   enddo
 
   call gpu_synchronize()
-  do a=1,nV
-    call gpu_stream_destroy(stream(a))
-  enddo
+  call gpu_deallocate(X_oo)
+
   call gpu_set_stream(blas_handle, gpu_default_stream)
 
   call gpu_dgemv(blas_handle, 'T', nV*nO, nO*nV, &
@@ -289,65 +282,38 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
                   H_vo%f(1,1)  , 1, &
              1d0, r1%f(1,1)    , 1)
 
-  call gpu_synchronize()
-  call gpu_deallocate(X_oo)
-  call gpu_deallocate(X_voov)
-
   type(gpu_double4) :: X_ovov
   call gpu_allocate(X_ovov, nO, nV, nO, nV)
 
-  !$omp parallel &
-  !$omp shared(nO,nV,d_cc_space_v_ovov,d_cc_space_v_voov,X_ovov) &
-  !$omp private(u,beta,i,a) &
-  !$omp default(none)
-  !$omp do
   do beta = 1, nV
-    do u = 1, nO
-      do a = 1, nv
-        do i = 1, nO
-          X_ovov%f(i,a,u,beta) = 2d0 * d_cc_space_v_voov%f(a,u,i,beta) - d_cc_space_v_ovov%f(u,a,i,beta)
-        enddo
-      enddo
+    call gpu_set_stream(blas_handle, stream(beta))
+    do u=1,nO
+      call gpu_dgeam(blas_handle, 'N', 'T', nO, nV, -1.d0, d_cc_space_v_ovov%f(1,1,u,beta), &
+      nO, 2.d0, d_cc_space_v_voov%f(1,u,1,beta), nV*nO, X_ovov%f(1,1,u,beta), nO)
     enddo
   enddo
-  !$omp end do
-  !$omp end parallel
 
-!  call dgemv('T', nO*nV, nO*nV, &
-!             1d0, X_ovov%f, size(X_ovov%f,1) * size(X_ovov%f,2), &
-!                  t1%f, 1, &
-!             1d0, r1%f, 1)
+  call gpu_set_stream(blas_handle, gpu_default_stream)
+  call gpu_synchronize()
+  call gpu_deallocate(X_voov)
+
   call gpu_dgemv(blas_handle, 'T', nO*nV, nO*nV, &
              1d0, X_ovov%f(1,1,1,1), size(X_ovov%f,1) * size(X_ovov%f,2), &
                   t1%f(1,1), 1, &
              1d0, r1%f(1,1), 1)
 
+
   integer :: iblock, block_size, nVmax
   type(gpu_double4) :: W_vvov, W_vvov_tmp, T_vvoo
 
   block_size = 16
-  call gpu_allocate(W_vvov,nV, nV,nO,block_size)
-  call gpu_allocate(W_vvov_tmp, nV,nO,nV,block_size)
   call gpu_allocate(T_vvoo, nV,nV,nO,nO)
 
-  call gpu_synchronize()
-  call gpu_deallocate(X_ovov)
+  call gpu_dgeam(blas_handle, 'T', 'N', nV*nV, nO*nO, 1.d0, tau%f(1,1,1,1), &
+    nO*nO, 0.d0, T_vvoo%f(1,1,1,1), nV*nV, T_vvoo%f(1,1,1,1), nV*nV)
 
-  !$omp parallel &
-  !$omp private(u,i,b,a) &
-  !$omp default(shared)
-  !$omp do
-  do u = 1, nO
-    do i = 1, nO
-      do b = 1, nV
-        do a = 1, nV
-          T_vvoo%f(a,b,i,u) = tau%f(i,u,a,b)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
+  call gpu_allocate(W_vvov,nV, nV,nO,block_size)
+  call gpu_allocate(W_vvov_tmp, nV,nO,nV,block_size)
 
   do iblock = 1, nV, block_size
     nVmax = min(block_size,nV-iblock+1)
@@ -358,22 +324,22 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
       0.d0, W_vvov_tmp%f(1,1,1,1), nV*nO)
 
     call gpu_synchronize()
-    !$omp parallel &
-    !$omp private(b,i,a,beta) &
-    !$omp default(shared)
-    do beta = 1,  nVmax
-      do i = 1, nO
-        !$omp do
-        do b = 1, nV
-          do a = 1, nV
-            W_vvov%f(a,b,i,beta) = 2d0 * W_vvov_tmp%f(a,i,b,beta) - W_vvov_tmp%f(b,i,a,beta)
-          enddo
-        enddo
-        !$omp end do nowait
+    do b=1,nV
+      call gpu_set_stream(blas_handle, stream(b))
+      do i=1,nO
+        call gpu_dgeam(blas_handle, 'N', 'N', nV, nVmax,  2.d0, W_vvov_tmp%f(1,i,b,1), &
+         nV*nO*nV, 0.d0, W_vvov_tmp%f(1,i,b,1), nV*nO*nV, W_vvov%f(1,b,i,1), nV*nV*nO)
       enddo
     enddo
-    !$omp barrier
-    !$omp end parallel
+
+    call gpu_synchronize()
+
+    do beta = 1,  nVmax
+      call gpu_set_stream(blas_handle, stream(beta))
+      call gpu_dgeam(blas_handle, 'N', 'T', nV, nV*nO,  1.d0, W_vvov%f(1,1,1,beta), &
+         nV, -1.d0, W_vvov_tmp%f(1,1,1,beta), nV*nO, W_vvov%f(1,1,1,beta), nV)
+    enddo
+    call gpu_synchronize()
 
     call gpu_dgemm(blas_handle, 'T','N',nO,nVmax,nO*nV*nV, &
              1d0, T_vvoo%f(1,1,1,1), nV*nV*nO, &
@@ -381,30 +347,24 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
              1d0, r1%f(1,iblock), nO)
   enddo
 
-  call gpu_synchronize()
-  call gpu_deallocate(W_vvov)
-  call gpu_deallocate(T_vvoo)
-
+  call gpu_deallocate(X_ovov)
 
   type(gpu_double4) :: W_oovo
   call gpu_allocate(W_oovo, nO,nO,nV,nO)
 
-  !$omp parallel &
-  !$omp shared(nO,nV,d_cc_space_v_oovo,W_oovo) &
-  !$omp private(u,a,i,j) &
-  !$omp default(none)
   do u = 1, nO
-    !$omp do
     do a = 1, nV
-      do j = 1, nO
-        do i = 1, nO
-          W_oovo%f(i,j,a,u) = 2d0 * d_cc_space_v_oovo%f(i,j,a,u) - d_cc_space_v_oovo%f(j,i,a,u)
-        enddo
-      enddo
+      call gpu_set_stream(blas_handle, stream(a))
+      call gpu_dgeam(blas_handle, 'N', 'T', nO, nO,  2.d0, d_cc_space_v_oovo%f(1,1,a,u), &
+         nO, -1.d0, d_cc_space_v_oovo%f(1,1,a,u), nO, W_oovo%f(1,1,a,u), nO)
     enddo
-    !$omp end do nowait
   enddo
-  !$omp end parallel
+
+  call gpu_set_stream(blas_handle, gpu_default_stream)
+  call gpu_synchronize()
+
+  call gpu_deallocate(W_vvov)
+  call gpu_deallocate(T_vvoo)
 
   ! Change the sign for consistency with the code in spin orbitals
   call gpu_dgemm(blas_handle, 'T','N', nO, nV, nO*nO*nV, &
@@ -422,6 +382,9 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1,d_cc_s
     enddo
   enddo
 
+  do a=1,nV
+    call gpu_stream_destroy(stream(a))
+  enddo
 
 end
 
diff --git a/src/gpu/gpu.h b/src/gpu/gpu.h
index 03aeeab7..ac70e21a 100644
--- a/src/gpu/gpu.h
+++ b/src/gpu/gpu.h
@@ -22,20 +22,20 @@ void gpu_ddot(const void* handle, const int64_t n, const double* x, const int64_
 
 void gpu_sdot(const void* handle, const int64_t n, const float* x, const int64_t incx, const float* y, const int64_t incy, float* result);
 
-void gpu_dgemv(const void* handle, const char transa, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double* x, const int64_t incx, const double beta, double* y, const int64_t incy);
+void gpu_dgemv(const void* handle, const char transa, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* x, const int64_t incx, const double* beta, double* y, const int64_t incy);
 
-void gpu_sgemv(const void* handle, const char transa, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float* x, const int64_t incx, const float beta, float* y, const int64_t incy);
+void gpu_sgemv(const void* handle, const char transa, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* x, const int64_t incx, const float* beta, float* y, const int64_t incy);
 
-void gpu_dgemm(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const int64_t k, const double alpha,
-               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double beta, double* c, const int64_t ldc);
+void gpu_dgemm(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const int64_t k, const double* alpha,
+               const double* a, const int64_t lda, const double* b, const int64_t ldb, const double* beta, double* c, const int64_t ldc);
 
-void gpu_sgemm(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const int64_t k, const float alpha,
-               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float beta, float* c, const int64_t ldc);
+void gpu_sgemm(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const int64_t k, const float* alpha,
+               const float* a, const int64_t lda, const float* b, const int64_t ldb, const float* beta, float* c, const int64_t ldc);
 
-void gpu_dgeam(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const double alpha,
-               const double* a, const int64_t lda, const double beta, const double* b, const int64_t ldb, double* c, const int64_t ldc);
+void gpu_dgeam(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const double* alpha,
+               const double* a, const int64_t lda, const double* beta, const double* b, const int64_t ldb, double* c, const int64_t ldc);
 
-void gpu_sgeam(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const float alpha,
-               const float* a, const int64_t lda, const float beta, const float* b, const int64_t ldb, float* c, const int64_t ldc);
+void gpu_sgeam(const void* handle, const char transa, const char transb, const int64_t m, const int64_t n, const float* alpha,
+               const float* a, const int64_t lda, const float* beta, const float* b, const int64_t ldb, float* c, const int64_t ldc);
diff --git a/src/gpu/gpu_module.F90 b/src/gpu/gpu_module.F90
index 949ae4fc..6050075f 100644
--- a/src/gpu/gpu_module.F90
+++ b/src/gpu/gpu_module.F90
@@ -145,7 +145,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in), value  :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc
-      real(c_double), intent(in), value      :: alpha, beta
+      real(c_double), intent(in)            :: alpha, beta
       type(c_ptr), value :: a, b, c
     end subroutine
 
@@ -155,7 +155,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in), value  :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc
-      real(c_float), intent(in), value      :: alpha, beta
+      real(c_float), intent(in)             :: alpha, beta
       real(c_float) :: a, b, c
     end subroutine
 
@@ -165,7 +165,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in)         :: transa
       integer(c_int64_t), intent(in), value :: m, n, lda, incx, incy
-      real(c_double), intent(in), value     :: alpha, beta
+      real(c_double), intent(in)            :: alpha, beta
       real(c_double)                        :: a, x, y
     end subroutine
 
@@ -175,7 +175,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in)         :: transa
       integer(c_int64_t), intent(in), value :: m, n, lda, incx, incy
-      real(c_float), intent(in), value      :: alpha, beta
+      real(c_float), intent(in)             :: alpha, beta
       real(c_float)                         :: a, x, y
     end subroutine
 
@@ -186,7 +186,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in)         :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, k, lda, ldb, ldc
-      real(c_double), intent(in), value :: alpha, beta
+      real(c_double), intent(in)            :: alpha, beta
       real(c_double) :: a, b, c
     end subroutine
 
@@ -196,7 +196,7 @@ module gpu
       type(c_ptr), value, intent(in)        :: handle
       character(c_char), intent(in), value  :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, k, lda, ldb, ldc
-      real(c_float), intent(in), value :: alpha, beta
+      real(c_float), intent(in)             :: alpha, beta
       real(c_float) :: a, b, c
     end subroutine