From 095053333a85f9ed7bde695f1b19844809b86c51 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 21 Aug 2023 13:17:39 +0200
Subject: [PATCH] Added SP files

---
 devel/ccsd_gpu/gpu_init_sp.c |  344 ++++++
 devel/ccsd_gpu/gpu_sp.c      | 2134 ++++++++++++++++++++++++++++++++++
 2 files changed, 2478 insertions(+)
 create mode 100644 devel/ccsd_gpu/gpu_init_sp.c
 create mode 100644 devel/ccsd_gpu/gpu_sp.c

diff --git a/devel/ccsd_gpu/gpu_init_sp.c b/devel/ccsd_gpu/gpu_init_sp.c
new file mode 100644
index 0000000..709f817
--- /dev/null
+++ b/devel/ccsd_gpu/gpu_init_sp.c
@@ -0,0 +1,344 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <omp.h>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+#include <assert.h>
+#include "gpu.h"
+
+gpu_data_sp* gpu_init_sp(
+   int nO, int nV, int cholesky_mo_num,
+   double* cc_space_v_oo_chol_, double* cc_space_v_ov_chol_,
+   double* cc_space_v_vo_chol_, double* cc_space_v_vv_chol_,
+   double* cc_space_v_oooo_, double* cc_space_v_vooo_, double* cc_space_v_voov_,
+   double* cc_space_v_oovv_, double* cc_space_v_vvoo_,
+   double* cc_space_v_oovo_, double* cc_space_v_ovvo_,
+   double* cc_space_v_ovov_, double* cc_space_v_ovoo_,
+   double* cc_space_f_oo_, double* cc_space_f_ov_,
+   double* cc_space_f_vo_, double* cc_space_f_vv_)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    gpu_data_sp* data = (gpu_data_sp*) malloc ((size_t) ngpus*sizeof(gpu_data_sp));
+    assert (data != NULL);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      cudaError_t cudaStat = cudaSuccess;
+      size_t lda;
+
+      int igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+      cublasHandle_t handle;
+      cublasCreate(&handle);
+
+/*
+      size_t limit;
+      size_t mtotal, mfree;
+      cudaStat = cudaMemGetInfo(&mfree, &mtotal);
+      fprintf(stderr, "%s\nTotal:%lu\nFree :%lu\n", cudaGetErrorString(cudaStat), mtotal, mfree);
+*/
+
+      float* d_cc_space_v_oo_chol = NULL;
+      lda = cholesky_mo_num * nO;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_oo_chol, lda * nO * sizeof(float));
+
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_oo_chol = malloc((size_t) lda * nO * sizeof(float));
+      assert (cc_space_v_oo_chol != NULL);
+      for (size_t i=0 ; i<(size_t)  lda * nO ; ++i) {
+        cc_space_v_oo_chol[i] = cc_space_v_oo_chol_[i];
+      }
+      cublasSetMatrix(cholesky_mo_num*nO, nO, sizeof(float), cc_space_v_oo_chol, lda, d_cc_space_v_oo_chol, lda);
+      free(cc_space_v_oo_chol);
+
+      float* d_cc_space_v_ov_chol = NULL;
+      lda = cholesky_mo_num * nO;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_ov_chol = malloc((size_t) lda * nV * sizeof(float));
+      assert (cc_space_v_ov_chol != NULL);
+      for (size_t i=0 ; i<(size_t)  lda * nV ; ++i) {
+        cc_space_v_ov_chol[i] = cc_space_v_ov_chol_[i];
+      }
+      cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(float), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
+      free(cc_space_v_ov_chol);
+
+      float* d_cc_space_v_vo_chol = NULL;
+      lda = cholesky_mo_num * nV;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_vo_chol = malloc((size_t) lda * nO * sizeof(float));
+      assert (cc_space_v_vo_chol != NULL);
+      for (size_t i=0 ; i<(size_t)  lda * nO ; ++i) {
+        cc_space_v_vo_chol[i] = cc_space_v_vo_chol_[i];
+      }
+      cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(float), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
+      free(cc_space_v_vo_chol);
+
+      float* d_cc_space_v_vv_chol = NULL;
+      lda = cholesky_mo_num * nV;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_vv_chol = malloc((size_t) lda * nV * sizeof(float));
+      assert (cc_space_v_vv_chol != NULL);
+      for (size_t i=0 ; i< (size_t) lda * nV ; ++i) {
+        cc_space_v_vv_chol[i] = cc_space_v_vv_chol_[i];
+      }
+      cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(float), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
+      free(cc_space_v_vv_chol);
+
+      float* d_cc_space_v_oooo = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_oooo, nO*nO*nO*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_oooo = malloc((size_t) nO*nO*nO*nO*sizeof(float));
+      assert (cc_space_v_oooo != NULL);
+      for (size_t i=0 ; i<(size_t) nO*nO*nO*nO ; ++i) {
+        cc_space_v_oooo[i] = cc_space_v_oooo_[i];
+      }
+      cublasSetMatrix(nO*nO, nO*nO, sizeof(float), cc_space_v_oooo, nO*nO, d_cc_space_v_oooo, nO*nO);
+      free(cc_space_v_oooo);
+
+      float* d_cc_space_v_vooo = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_vooo, nV*nO*nO*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_vooo = malloc((size_t) nV*nO*nO*nO*sizeof(float));
+      assert (cc_space_v_vooo != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nO*nO*nO; ++i) {
+        cc_space_v_vooo[i] = cc_space_v_vooo_[i];
+      }
+      cublasSetMatrix(nV*nO, nO*nO, sizeof(float), cc_space_v_vooo, nV*nO, d_cc_space_v_vooo, nV*nO);
+      free(cc_space_v_vooo);
+
+      float* d_cc_space_v_voov = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_voov, nV*nO*nO*nV*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_voov = malloc((size_t) nV*nO*nO*nV*sizeof(float));
+      assert (cc_space_v_voov != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nO*nO*nV; ++i) {
+        cc_space_v_voov[i] = cc_space_v_voov_[i];
+      }
+      cublasSetMatrix(nV*nO, nO*nV, sizeof(float), cc_space_v_voov, nV*nO, d_cc_space_v_voov, nV*nO);
+      free(cc_space_v_voov);
+
+      float* d_cc_space_v_oovv = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_oovv, nO*nO*nV*nV*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_oovv = malloc((size_t) nO*nO*nV*nV*sizeof(float));
+      assert (cc_space_v_oovv != NULL);
+      for (size_t i=0 ; i<(size_t) nO*nO*nV*nV; ++i) {
+        cc_space_v_oovv[i] = cc_space_v_oovv_[i];
+      }
+      cublasSetMatrix(nO*nO, nV*nV, sizeof(float), cc_space_v_oovv, nO*nO, d_cc_space_v_oovv, nO*nO);
+      free(cc_space_v_oovv);
+
+      float* d_cc_space_v_vvoo = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_vvoo, nV*nV*nO*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_vvoo = malloc((size_t) nV*nV*nO*nO*sizeof(float));
+      assert (cc_space_v_vvoo != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nV*nO*nO; ++i) {
+        cc_space_v_vvoo[i] = cc_space_v_vvoo_[i];
+      }
+      cublasSetMatrix(nV*nV, nO*nO, sizeof(float), cc_space_v_vvoo, nV*nV, d_cc_space_v_vvoo, nV*nV);
+      free(cc_space_v_vvoo);
+
+      float* d_cc_space_v_oovo = NULL;
+      lda = nO*nO;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_oovo, nO*nO*nV*nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_oovo = malloc((size_t) nO*nO*nV*nO * sizeof(float));
+      assert (cc_space_v_oovo != NULL);
+      for (size_t i=0 ; i<(size_t)  nO*nO*nV*nO ; ++i) {
+        cc_space_v_oovo[i] = cc_space_v_oovo_[i];
+      }
+      cublasSetMatrix(lda, nV*nO, sizeof(float), cc_space_v_oovo, lda, d_cc_space_v_oovo, lda);
+      free(cc_space_v_oovo);
+
+      float* d_cc_space_v_ovvo = NULL;
+      lda = nO*nV;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovvo, nO*nV*nV*nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_ovvo = malloc((size_t) nO*nV*nV*nO * sizeof(float));
+      assert (cc_space_v_ovvo != NULL);
+      for (size_t i=0 ; i<(size_t)  nO*nV*nV*nO ; ++i) {
+        cc_space_v_ovvo[i] = cc_space_v_ovvo_[i];
+      }
+      cublasSetMatrix(lda, nV*nO, sizeof(float), cc_space_v_ovvo, lda, d_cc_space_v_ovvo, lda);
+      free(cc_space_v_ovvo);
+
+      float* d_cc_space_v_ovov = NULL;
+      lda = nO*nV;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovov, nO*nV*nV*nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_ovov = malloc((size_t) nO*nV*nV*nO * sizeof(float));
+      assert (cc_space_v_ovov != NULL);
+      for (size_t i=0 ; i<(size_t)  nO*nV*nV*nO ; ++i) {
+        cc_space_v_ovov[i] = cc_space_v_ovov_[i];
+      }
+      cublasSetMatrix(lda, nV*nO, sizeof(float), cc_space_v_ovov, lda, d_cc_space_v_ovov, lda);
+      free(cc_space_v_ovov);
+
+      float* d_cc_space_v_ovoo = NULL;
+      lda = nO*nV;
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovoo, nO*nV*nO*nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_v_ovoo = malloc((size_t) nO*nV*nO*nO * sizeof(float));
+      assert (cc_space_v_ovoo != NULL);
+      for (size_t i=0 ; i<(size_t)  nO*nV*nO*nO ; ++i) {
+        cc_space_v_ovoo[i] = cc_space_v_ovoo_[i];
+      }
+      cublasSetMatrix(lda, nO*nO, sizeof(float), cc_space_v_ovoo, lda, d_cc_space_v_ovoo, lda);
+      free(cc_space_v_ovoo);
+
+      float* d_cc_space_f_oo = NULL;
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_oo, nO*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_f_oo = malloc((size_t) nO*nO*sizeof(float));
+      assert (cc_space_f_oo != NULL);
+      for (size_t i=0 ; i<(size_t) nO*nO; ++i) {
+        cc_space_f_oo[i] = cc_space_f_oo_[i];
+      }
+      cublasSetMatrix(nO, nO, sizeof(float), cc_space_f_oo, nO, d_cc_space_f_oo, nO);
+      free(cc_space_f_oo);
+
+      float* d_cc_space_f_vo = NULL;
+
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_f_vo = malloc((size_t) nV*nO*sizeof(float));
+      assert (cc_space_f_vo != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nO; ++i) {
+        cc_space_f_vo[i] = cc_space_f_vo_[i];
+      }
+      cublasSetMatrix(nV, nO, sizeof(float), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
+      free(cc_space_f_vo);
+
+      float* d_cc_space_f_ov = NULL;
+
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_ov, nV*nO*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_f_ov = malloc((size_t) nV*nO*sizeof(float));
+      assert (cc_space_f_ov != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nO; ++i) {
+        cc_space_f_ov[i] = cc_space_f_ov_[i];
+      }
+      cublasSetMatrix(nO, nV, sizeof(float), cc_space_f_ov, nO, d_cc_space_f_ov, nO);
+      free(cc_space_f_ov);
+
+      float* d_cc_space_f_vv = NULL;
+
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_vv, nV*nV*sizeof(float));
+      assert (cudaStat == cudaSuccess);
+      float* cc_space_f_vv = malloc((size_t) nV*nV*sizeof(float));
+      assert (cc_space_f_vv != NULL);
+      for (size_t i=0 ; i<(size_t) nV*nV; ++i) {
+        cc_space_f_vv[i] = cc_space_f_vv_[i];
+      }
+      cublasSetMatrix(nV, nV, sizeof(float), cc_space_f_vv, nV, d_cc_space_f_vv, nV);
+      free(cc_space_f_vv);
+
+      float* d_tau = NULL;
+      lda = nO * nO;
+      cudaStat = cudaMalloc((void **)&d_tau, lda * nV * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_tau_x = NULL;
+      lda = nO * nO;
+      cudaStat = cudaMalloc((void **)&d_tau_x, lda * nV * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_t1 = NULL;
+      cudaStat = cudaMalloc((void **)&d_t1, nO * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_t2 = NULL;
+      cudaStat = cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_H_oo = NULL;
+      cudaStat = cudaMalloc((void **)&d_H_oo, nO * nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_H_vo = NULL;
+      cudaStat = cudaMalloc((void **)&d_H_vo, nV * nO * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      float* d_H_vv = NULL;
+      cudaStat = cudaMalloc((void **)&d_H_vv, nV * nV * sizeof(float));
+      assert (cudaStat == cudaSuccess);
+
+      data[igpu].cc_space_v_oo_chol = d_cc_space_v_oo_chol;
+      data[igpu].cc_space_v_ov_chol = d_cc_space_v_ov_chol;
+      data[igpu].cc_space_v_vo_chol = d_cc_space_v_vo_chol;
+      data[igpu].cc_space_v_vv_chol = d_cc_space_v_vv_chol;
+      data[igpu].cc_space_v_oooo = d_cc_space_v_oooo;
+      data[igpu].cc_space_v_vooo = d_cc_space_v_vooo;
+      data[igpu].cc_space_v_voov = d_cc_space_v_voov;
+      data[igpu].cc_space_v_oovv = d_cc_space_v_oovv;
+      data[igpu].cc_space_v_vvoo = d_cc_space_v_vvoo;
+      data[igpu].cc_space_v_oovo = d_cc_space_v_oovo;
+      data[igpu].cc_space_v_ovvo = d_cc_space_v_ovvo;
+      data[igpu].cc_space_v_ovov = d_cc_space_v_ovov;
+      data[igpu].cc_space_v_ovoo = d_cc_space_v_ovoo;
+      data[igpu].cc_space_f_oo   = d_cc_space_f_oo;
+      data[igpu].cc_space_f_ov   = d_cc_space_f_ov;
+      data[igpu].cc_space_f_vo   = d_cc_space_f_vo;
+      data[igpu].cc_space_f_vv   = d_cc_space_f_vv;
+      data[igpu].tau   = d_tau;
+      data[igpu].tau_x = d_tau_x;
+      data[igpu].t1    = d_t1;
+      data[igpu].t2    = d_t2;
+      data[igpu].H_oo  = d_H_oo;
+      data[igpu].H_vo  = d_H_vo;
+      data[igpu].H_vv  = d_H_vv;
+
+      data[igpu].nO = nO;
+      data[igpu].nV = nV;
+      data[igpu].cholesky_mo_num = cholesky_mo_num;
+
+    }
+    return data;
+}
+
+
+void gpu_deinit_sp(gpu_data_sp* data)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      size_t lda;
+      int igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+
+      free(data[igpu].cc_space_v_oo_chol);
+      free(data[igpu].cc_space_v_ov_chol);
+      free(data[igpu].cc_space_v_vo_chol);
+      free(data[igpu].cc_space_v_vv_chol);
+      free(data[igpu].cc_space_v_oooo);
+      free(data[igpu].cc_space_v_vooo);
+      free(data[igpu].cc_space_v_voov);
+      free(data[igpu].cc_space_v_oovv);
+      free(data[igpu].cc_space_v_vvoo);
+      free(data[igpu].cc_space_v_oovo);
+      free(data[igpu].cc_space_v_ovvo);
+      free(data[igpu].cc_space_v_ovov);
+      free(data[igpu].cc_space_v_ovoo);
+      free(data[igpu].cc_space_f_oo);
+      free(data[igpu].cc_space_f_ov);
+      free(data[igpu].cc_space_f_vo);
+      free(data[igpu].cc_space_f_vv);
+      free(data[igpu].tau);
+      free(data[igpu].tau_x);
+      free(data[igpu].t1);
+      free(data[igpu].t2);
+      free(data[igpu].H_oo);
+      free(data[igpu].H_vo);
+      free(data[igpu].H_vv);
+    }
+}
+
diff --git a/devel/ccsd_gpu/gpu_sp.c b/devel/ccsd_gpu/gpu_sp.c
new file mode 100644
index 0000000..2dc84f3
--- /dev/null
+++ b/devel/ccsd_gpu/gpu_sp.c
@@ -0,0 +1,2134 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <omp.h>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+#include <assert.h>
+#include "gpu.h"
+
+void gpu_upload_sp(gpu_data_sp* data,
+                int nO, int nV,
+                double* t1,
+                double* t2)
+{
+    size_t lda;
+
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    float* tau = malloc((size_t) nO*nO*nV*nV * sizeof(float));
+    assert (tau != NULL);
+
+    float* tau_x = malloc((size_t) nO*nO*nV*nV * sizeof(float));
+    assert (tau_x != NULL);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+        cudaError_t cudaStat = cudaSuccess;
+        float* AA;
+        size_t igpu = omp_get_thread_num();
+        cudaSetDevice(igpu);
+
+        float* d_t1 = data[igpu].t1;
+        lda = nO;
+        AA = malloc((size_t) nO*nV*sizeof(float));
+        assert (AA != NULL);
+        for (size_t i=0 ; i<(size_t) nO*nV ; ++i) {
+          AA[i] = t1[i];
+        }
+        cudaStat = cublasSetMatrix(nO, nV, sizeof(float), AA, lda, d_t1, lda);
+        assert (cudaStat == cudaSuccess);
+        free(AA);
+
+        float* d_t2 = data[igpu].t2;
+        lda = nO*nO;
+        AA = malloc((size_t) nO*nO*nV*nV*sizeof(float));
+        assert (AA != NULL);
+        for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
+          AA[i] = t2[i];
+        }
+        cudaStat = cublasSetMatrix(nO*nO, nV*nV, sizeof(float), AA, lda, d_t2, lda);
+        assert (cudaStat == cudaSuccess);
+        free(AA);
+
+        size_t lda, ldb, ldc;
+        float alpha, beta;
+        float* A;
+        float* B;
+        float* C;
+
+        cublasHandle_t handle;
+        cublasCreate(&handle);
+        cudaStream_t stream[nV];
+
+        float* d_tau = data[igpu].tau;
+
+        float* d_tau_x = data[igpu].tau_x;
+        lda = nO * nO;
+        cublasSetMatrix(nO*nO, nV*nV, sizeof(float), tau_x, lda, d_tau_x, lda);
+
+        if (igpu == 0) {
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              beta = t1[j+b*nO];
+              A = &(d_t2[nO*(j + nO*nV*b)]); lda = nO*nO;
+              B = d_t1; ldb = nO;
+              C = &(d_tau[nO*(j + nO*nV*b)]); ldc = nO*nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          cudaDeviceSynchronize();
+          alpha = 2.0;
+          beta = -1.0;
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t a=0 ; a<nV ; ++a) {
+              cublasSetStream(handle, stream[a]);
+              A = &(d_tau[nO*nO*(a + nV*b)]); lda = nO;
+              B = &(d_tau[nO*nO*(b + nV*a)]); ldb = nO;
+              C = &(d_tau_x[nO*nO*(a + nV*b)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          lda = nO*nO;
+          cublasGetMatrix(nO*nO, nV*nV, sizeof(float), d_tau, lda, tau, lda);
+          cublasGetMatrix(nO*nO, nV*nV, sizeof(float), d_tau_x, lda, tau_x, lda);
+        }
+
+        #pragma omp barrier
+
+        if (igpu > 0) {
+          lda = nO * nO;
+          cublasSetMatrix(nO*nO, nV*nV, sizeof(float), tau, lda, d_tau, lda);
+          cublasSetMatrix(nO*nO, nV*nV, sizeof(float), tau_x, lda, d_tau_x, lda);
+        }
+        cublasDestroy(handle);
+    }
+    free(tau);
+    free(tau_x);
+}
+
+
+
+
+void compute_h_oo_chol_gpu_sp(gpu_data_sp* data, int igpu)
+{
+    cudaError_t cudaStat;
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
+    cudaSetDevice(igpu);
+
+    size_t m,n,k, lda, ldb, ldc;
+    float alpha, beta;
+    float* A;
+    float* B;
+    float* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    float* d_H_oo               = data[igpu].H_oo;
+    float* d_tau_x              = data[igpu].tau_x;
+    float* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
+    float* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+    float* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+
+    float* d_tau_kau;
+    cudaStat = cudaMalloc((void **)&d_tau_kau, cholesky_mo_num*nV*nO * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    float* d_tmp_ovv;
+    cudaStat = cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    float* d_tmp_vov;
+    cudaStat = cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    alpha = 1.0;
+    beta = 0.0;
+    for (size_t u=0 ; u<nO ; ++u) {
+      cublasScopy(handle, nO*nV*nV, &(d_tau_x[u]), nO, d_tmp_ovv, 1);
+      for (size_t b=0 ; b<nV ; ++b) {
+            cublasSetStream(handle, stream[b]);
+            A = &(d_tmp_ovv[nO*nV*b]); lda = nO;
+            B = &(d_tmp_ovv[nO*nV*b]); ldb = nO;
+            C = &(d_tmp_vov[nV*nO*b]); ldc = nV;
+            cublasSgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+      }
+      cudaDeviceSynchronize();
+      cublasSetStream(handle, NULL);
+      alpha = 1.0;
+      beta = 0.0;
+      m=cholesky_mo_num;  n=nV;  k=nO*nV;
+      A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
+      B=d_tmp_vov;  ldb=nV;
+      C=&(d_tau_kau[cholesky_mo_num*nV*u]);  ldc=cholesky_mo_num;
+      cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+    }
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamDestroy(stream[i]);
+    }
+
+    cudaFree(d_tmp_vov);
+    cudaFree(d_tmp_ovv);
+
+    cublasScopy(handle, nO*nO, d_cc_space_f_oo, 1, d_H_oo, 1);
+    alpha = 1.0;
+    beta  = 1.0;
+    m=nO;  n=nO;  k=cholesky_mo_num*nV;
+    A=d_tau_kau;  lda=cholesky_mo_num*nV;
+    B=d_cc_space_v_vo_chol;  ldb=cholesky_mo_num*nV;
+    C=d_H_oo;  ldc=nO;
+    cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    cudaFree(d_tau_kau);
+
+    float* H_oo = malloc((size_t) nO*nO*sizeof(float));
+    assert (H_oo != NULL);
+    cublasGetMatrix(nO, nO, sizeof(float), d_H_oo, nO, H_oo, nO);
+    for (size_t i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         float* d_H_oo = data[i].H_oo;
+         cudaSetDevice(i);
+         cublasSetMatrix(nO, nO, sizeof(float), H_oo, nO, d_H_oo, nO);
+       }
+    }
+    free(H_oo);
+
+    cublasDestroy(handle);
+
+}
+
+
+void compute_h_vo_chol_gpu_sp(gpu_data_sp* data, int igpu)
+{
+    cudaError_t cudaStat;
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
+    cudaSetDevice(igpu);
+
+    size_t m,n,k, lda, ldb, ldc;
+    float alpha, beta;
+    float* A;
+    float* B;
+    float* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    float* d_t1                 = data[igpu].t1;
+    float* d_H_vo               = data[igpu].H_vo;
+    float* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+    float* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+    float* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+
+    cublasScopy(handle, nV*nO, d_cc_space_f_vo, 1, d_H_vo, 1);
+
+    float* d_tmp_k;
+    cudaStat = cudaMalloc((void **)&d_tmp_k, cholesky_mo_num * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    alpha = 2.0;
+    beta  = 0.0;
+    m=cholesky_mo_num;  n=1;  k=nO*nV;
+    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
+    B=d_t1;  ldb=nO*nV;
+    C=d_tmp_k;  ldc=cholesky_mo_num;
+    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    alpha = 1.0;
+    beta  = 1.0;
+    m=nV*nO;  n=1;  k=cholesky_mo_num;
+    A=d_cc_space_v_vo_chol;  lda=cholesky_mo_num;
+    B=d_tmp_k;  ldb=cholesky_mo_num;
+    C=d_H_vo;  ldc=nV*nO;
+    cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    cudaFree(d_tmp_k);
+
+    float* d_tmp;
+    cudaStat = cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    alpha = 1.0;
+    beta  = 0.0;
+    m=cholesky_mo_num*nO;  n=nO;  k=nV;
+    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
+    B=d_t1;  ldb=nO;
+    C=d_tmp;  ldc=cholesky_mo_num*nO;
+    cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    float* d_tmp2;
+    cudaStat = cudaMalloc((void **)&d_tmp2, cholesky_mo_num*nO*nO * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    for (size_t i=0 ; i<nO ; ++i) {
+      for (size_t j=0 ; j<nO ; ++j) {
+         cublasSetStream(handle, stream[j]);
+         cublasScopy(handle, cholesky_mo_num, &(d_tmp [cholesky_mo_num*(i+nO*j)]), 1,
+                                              &(d_tmp2[cholesky_mo_num*(j+nO*i)]), 1);
+      }
+    }
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamDestroy(stream[i]);
+    }
+    cublasSetStream(handle, NULL);
+
+    alpha = -1.0;
+    beta  = 1.0;
+    m=nV;  n=nO;  k=cholesky_mo_num*nO;
+    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
+    B=d_tmp2;  ldb=cholesky_mo_num*nO;
+    C=d_H_vo;  ldc=nV;
+    cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    float* H_vo = malloc((size_t) nV*nO*sizeof(float));
+    assert (H_vo != NULL);
+    cublasGetMatrix(nV, nO, sizeof(float), d_H_vo, nV, H_vo, nV);
+    for (size_t i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         float* d_H_vo = data[i].H_vo;
+         cudaSetDevice(i);
+         cublasSetMatrix(nV, nO, sizeof(float), H_vo, nV, d_H_vo, nV);
+       }
+    }
+    free(H_vo);
+
+    cublasDestroy(handle);
+
+}
+
+
+
+void compute_h_vv_chol_gpu_sp(gpu_data_sp* data, int igpu)
+{
+    cudaError_t cudaStat;
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
+    cudaSetDevice(igpu);
+
+    size_t m,n,k, lda, ldb, ldc;
+    float alpha, beta;
+    float* A;
+    float* B;
+    float* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    float* d_H_vv               = data[igpu].H_vv;
+    float* d_tau_x              = data[igpu].tau_x;
+    float* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
+    float* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+
+    float* d_tau_kia;
+    cudaStat = cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    float* d_tmp_oov;
+    cudaStat = cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(float));
+    assert(cudaStat == cudaSuccess);
+
+    alpha = 1.0;
+    beta = 0.0;
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    for (size_t a=0 ; a<nV ; ++a) {
+      for (size_t b=0 ; b<nV ; ++b) {
+         cublasSetStream(handle, stream[b]);
+         cublasScopy(handle, nO*nO, &(d_tau_x[nO*nO*(a+nV*b)]), 1, &(d_tmp_oov[nO*nO*b]), 1);
+      }
+      cudaDeviceSynchronize();
+      cublasSetStream(handle, NULL);
+      alpha = 1.0;
+      beta = 0.0;
+      m=cholesky_mo_num;  n=nO;  k=nO*nV;
+      A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
+      B=d_tmp_oov;  ldb=nO;
+      C=&(d_tau_kia[cholesky_mo_num*nO*a]);  ldc=cholesky_mo_num;
+      cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+    }
+    for (size_t i=0 ; i<nV ; ++i) {
+      cudaStreamDestroy(stream[i]);
+    }
+    cudaFree(d_tmp_oov);
+
+    cublasScopy(handle, nV*nV, d_cc_space_f_vv, 1, d_H_vv, 1);
+    alpha = -1.0;
+    beta  = 1.0;
+    m=nV;  n=nV;  k=cholesky_mo_num*nO;
+    A=d_tau_kia;  lda=cholesky_mo_num*nO;
+    B=d_cc_space_v_ov_chol;  ldb=cholesky_mo_num*nO;
+    C=d_H_vv;  ldc=nV;
+    cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    cudaFree(d_tau_kia);
+
+    float* H_vv = malloc((size_t) nV*nV*sizeof(float));
+    assert (H_vv != NULL);
+    cublasGetMatrix(nV, nV, sizeof(float), d_H_vv, nV, H_vv, nV);
+    for (size_t i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         float* d_H_vv = data[i].H_vv;
+         cudaSetDevice(i);
+         cublasSetMatrix(nV, nV, sizeof(float), H_vv, nV, d_H_vv, nV);
+       }
+    }
+    free(H_vv);
+
+    cublasDestroy(handle);
+
+}
+
+
+
+
+void compute_r2_space_chol_gpu_sp(gpu_data_sp* data, int nO, int nV, double* t1, double* r2, double* max_r2)
+{
+    const size_t cholesky_mo_num = data->cholesky_mo_num;
+
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    float* J1 = malloc((size_t) nO*nV*nV*nO*sizeof(float));
+    assert (J1 != NULL);
+
+    float* K1 = malloc((size_t) nO*nV*nV*nO*sizeof(float));
+    assert (K1 != NULL);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      cudaError_t cudaStat;
+      size_t m,n,k, lda, ldb, ldc;
+      float alpha, beta;
+      float* A;
+      float* B;
+      float* C;
+      cudaStream_t stream[nV];
+
+      size_t igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+
+      cublasHandle_t handle;
+      cublasCreate(&handle);
+
+      float* d_r2;
+      lda = nO * nO;
+      cudaStat = cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(float));
+      assert(cudaStat == cudaSuccess);
+      cudaMemset(d_r2, 0, nO*nO*nV*nV*sizeof(float));
+      memset(r2, 0, nO*nO*nV*nV*sizeof(double));
+
+      float* d_cc_space_v_oo_chol = data[igpu].cc_space_v_oo_chol;
+      float* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+      float* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+      float* d_cc_space_v_vv_chol = data[igpu].cc_space_v_vv_chol;
+      float* d_cc_space_v_oooo    = data[igpu].cc_space_v_oooo;
+      float* d_cc_space_v_vooo    = data[igpu].cc_space_v_vooo;
+      float* d_cc_space_v_oovv    = data[igpu].cc_space_v_oovv;
+      float* d_cc_space_v_vvoo    = data[igpu].cc_space_v_vvoo;
+      float* d_cc_space_v_oovo    = data[igpu].cc_space_v_oovo;
+      float* d_cc_space_v_ovvo    = data[igpu].cc_space_v_ovvo;
+      float* d_cc_space_v_ovov    = data[igpu].cc_space_v_ovov;
+      float* d_cc_space_v_ovoo    = data[igpu].cc_space_v_ovoo;
+      float* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+      float* d_tau                = data[igpu].tau;
+      float* d_t1                 = data[igpu].t1;
+      float* d_t2                 = data[igpu].t2;
+      float* d_H_oo               = data[igpu].H_oo;
+      float* d_H_vv               = data[igpu].H_vv;
+
+      float* d_K1;
+      cudaStat = cudaMalloc((void **)&d_K1, nO*nV*nO*nV * sizeof(float));
+      assert(cudaStat == cudaSuccess);
+
+      #pragma omp sections
+      {
+
+        #pragma omp section
+        {
+          float* d_J1;
+          cudaStat = cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta =  0.0;
+          A = d_cc_space_v_ovvo; lda = nO*nV;
+          B = d_cc_space_v_ovvo; ldb = nO*nV;
+          C = d_J1; ldc = nO*nV;
+          cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nV*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+
+          float* d_X_ovoo;
+          cudaStat = cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+          alpha = 0.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
+              cublasSetStream(handle, stream[i]);
+              A = &(d_X_ovoo[nO*nV*(i+nO*j)]); lda = nO;
+              B = &(d_cc_space_v_ovoo[nO*nV*(j+nO*i)]); ldb = nO;
+              C = &(d_X_ovoo[nO*nV*(i+nO*j)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nO ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+
+          float* d_Y_ovov;
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV*nO;  n=nV;  k=nO;
+          A=d_X_ovoo;  lda=nO*nV*nO;
+          B=d_t1;  ldb=nO;
+          C=d_Y_ovov;  ldc=nO*nV*nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_X_ovoo);
+
+          alpha = 1.0;
+          beta =  -1.0;
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
+              cublasSetStream(handle, stream[i]);
+              A = &(d_J1[nO*nV*(i+nV*j)]); lda = nO;
+              B = &(d_Y_ovov[nO*nV*(j+nO*i)]); ldb = nO;
+              C = &(d_J1[nO*nV*(i+nV*j)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_tmp_cc;
+          cudaStat = cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=cholesky_mo_num*nV;  n=nO;  k=nV;
+          A=d_cc_space_v_vv_chol;  lda=cholesky_mo_num*nV;
+          B=d_t1;  ldb=nO;
+          C=d_tmp_cc;  ldc=cholesky_mo_num*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_J1_tmp;
+          cudaStat = cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nV*nO;  n=nV*nO;  k=cholesky_mo_num;
+          A=d_tmp_cc;  lda=cholesky_mo_num;
+          B=d_cc_space_v_vo_chol;  ldb=cholesky_mo_num;
+          C=d_J1_tmp;  ldc=nV*nO;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_tmp_cc);
+
+          alpha = 1.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t i=0 ; i<nO ; ++i) {
+              cublasSetStream(handle, stream[i]);
+              A = &(d_J1[nO*nV*nV*i]); lda = nO*nV;
+              B = &(d_J1_tmp[nV*nO*nV*i]); ldb = nV;
+              C = &(d_J1[nO*nV*nV*i]); ldc = nO*nV;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO*nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          }
+          for (size_t i=0 ; i<nO ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_J1_tmp);
+
+          float* d_X_voov;
+          cudaStat = cudaMalloc((void **)&d_X_voov, nV*nO*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 0.5;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              beta = t1[j+b*nO];
+              A = &(d_t2[nO*(j + nO*nV*b)]); lda = nO*nO;
+              B = d_t1; ldb = nO;
+              C = &(d_Y_ovov[nO*(b+nV*j)]); ldc = nO*nV*nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          alpha = 1.0;
+          beta  = 0.0;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); lda = nV*nV;
+              B = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); ldb = nV*nV;
+              C = &(d_X_voov[nV*nO*(j+nO*b)]); ldc = nV;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Z_ovvo;
+          cudaStat = cudaMalloc((void **)&d_Z_ovvo, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = -1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nV*nO;  k=nO*nV;
+          A=d_Y_ovov;  lda=nO*nV;
+          B=d_X_voov;  ldb=nV*nO;
+          C=d_Z_ovvo;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_X_voov);
+          cudaFree(d_Y_ovov);
+
+          alpha = 1.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
+              B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb=nO*nV;
+              C = &(d_J1[nO*nV*(b+nV*i)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          float* d_Y_vovo;
+          cudaStat = cudaMalloc((void **)&d_Y_vovo, nV*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta = -0.5;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
+              cublasSetStream(handle, stream[i]);
+              A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
+              B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
+              C = &(d_Y_vovo[nV*(i+nO*nV*j)]); ldc = nV*nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          float* d_X_ovvo;
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta =  0.0;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
+              B = &(d_t2[nO*(j+nO*nV*b)]); ldb = nO*nO;
+              C = &(d_X_ovvo[nO*nV*(b+nV*j)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nV*nO;  k=nV*nO;
+          A=d_X_ovvo;  lda=nO*nV;
+          B=d_Y_vovo;  ldb=nV*nO;
+          C=d_Z_ovvo;  ldc=nO*nV;
+         cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_X_ovvo);
+          cudaFree(d_Y_vovo);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
+              B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb = nO*nV;
+              C = &(d_J1[nO*nV*(b+nV*i)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          cudaFree(d_Z_ovvo);
+
+          lda = nO*nV;
+          cublasGetMatrix(nO*nV, nV*nO, sizeof(float), d_J1, lda, J1, lda);
+          cudaFree(d_J1);
+
+        }
+
+        #pragma omp section
+        {
+
+          alpha = 1.0;
+          beta =  0.0;
+          A = d_cc_space_v_ovov; lda = nO*nV;
+          B = d_cc_space_v_ovov; ldb = nO*nV;
+          C = d_K1; ldc = nO*nV;
+          cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nO*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+          alpha = -1.0;
+          beta  = 1.0;
+          m=nO*nV*nO;  n=nV;  k=nO;
+          A=d_cc_space_v_ovoo;  lda=nO*nV*nO;
+          B=d_t1;  ldb=nO;
+          C=d_K1;  ldc=nO*nV*nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_X;
+          cudaStat = cudaMalloc((void **)&d_X, nV*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          float* d_Y;
+          cudaStat = cudaMalloc((void **)&d_Y, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha =-1.0;
+          beta =  0.0;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
+              cublasSetStream(handle, stream[i]);
+              A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
+              B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
+              C = &(d_X[nV*(j+nO*nV*i)]); ldc = nV*nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          alpha = 0.5;
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              beta = t1[j+b*nO];
+              A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
+              B = d_t1; ldb = nO;
+              C = &(d_Y[nO*(b+nV*nV*j)]); ldc = nO*nV;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Z;
+          cudaStat = cudaMalloc((void **)&d_Z, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nV*nO;  n=nO*nV;  k=nV*nO;
+          A=d_Y;  lda=nO*nV;
+          B=d_X;  ldb=nV*nO;
+          C=d_Z;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X);
+          cudaFree(d_Y);
+
+
+          float* d_t1v;
+          cudaStat = cudaMalloc((void **)&d_t1v, cholesky_mo_num*nO*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=cholesky_mo_num*nO;  n=nO;  k=nV;
+          A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
+          B=d_t1;  ldb=nO;
+          C=d_t1v;  ldc=cholesky_mo_num*nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_K1tmp;
+          cudaStat = cudaMalloc((void **)&d_K1tmp, nO*nO*nV*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nO;  n=nV*nV;  k=cholesky_mo_num;
+          A=d_t1v;  lda=cholesky_mo_num;
+          B=d_cc_space_v_vv_chol;  ldb=cholesky_mo_num;
+          C=d_K1tmp;  ldc=nO*nO;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_t1v);
+
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
+              B = &(d_K1tmp[nO*(i+nO*nV*b)]); ldb = nO*nO;
+              C = &(d_K1[nO*nV*(i+nO*b)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
+              B = &(d_Z[nO*(b+nV*nV*i)]); ldb = nO*nV;
+              C = &(d_K1[nO*nV*(i+nO*b)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_K1tmp);
+          cudaFree(d_Z);
+
+          lda = nO*nV;
+          cublasGetMatrix(nO*nV, nO*nV, sizeof(float), d_K1, lda, K1, lda);
+
+        }
+
+        #pragma omp section
+        {
+           float* d_Y_oooo;
+           cudaStat = cudaMalloc((void**)&d_Y_oooo, nO*nO*nO*nO*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+
+           alpha = 1.0;
+           beta =  0.0;
+           m=nO ; n=nO*nO*nO; k=nV;
+           A = d_t1 ; lda = nO;
+           B = d_cc_space_v_vooo ; ldb = nV;
+           C = d_Y_oooo; ldc = nO;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+           float* d_A1;
+           cudaStat = cudaMalloc((void**)&d_A1, nO*nO*nO*nO*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+
+           alpha = 1.0;
+           beta =  1.0;
+           A = d_cc_space_v_oooo; lda = nO*nO;
+           B = d_Y_oooo; ldb = nO*nO;
+           C = d_A1; ldc = nO*nO;
+           cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+           for (size_t i=0 ; i<nO ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (size_t j=0 ; j<nO ; ++j) {
+             for (size_t i=0 ; i<nO ; ++i) {
+                cublasSetStream(handle, stream[i]);
+                alpha = 1.0;
+                beta =  1.0;
+                A = &(d_A1[nO*nO*(i+nO*j)]); lda = nO;
+                B = &(d_Y_oooo[nO*nO*(j+nO*i)]); ldb = nO;
+                C = &(d_A1[nO*nO*(i+nO*j)]); ldc = nO;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+             }
+           }
+           for (size_t i=0 ; i<nO ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
+           cudaFree(d_Y_oooo);
+
+           alpha = 1.0;
+           beta =  1.0;
+           m=nO*nO ; n=nO*nO; k=nV*nV;
+           A = d_tau ; lda = nO*nO;
+           B = d_cc_space_v_vvoo ; ldb = nV*nV;
+           C = d_A1; ldc = nO*nO;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+           alpha = 1.0;
+           beta =  0.0;
+           m=nO*nO ; n=nV*nV; k=nO*nO;
+           A = d_A1  ; lda = nO*nO;
+           B = d_tau ; ldb = nO*nO;
+           C = d_r2; ldc = nO*nO;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           cudaFree(d_A1);
+        }
+
+        // g_vir
+        #pragma omp section
+        {
+           float* d_g_vir;
+           cudaStat = cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+           cublasScopy(handle, nV*nV, d_H_vv, 1, d_g_vir, 1);
+
+           alpha = -1.0;
+           beta =  1.0;
+           m=nV ; n=nV; k=nO;
+           A = d_cc_space_f_vo  ; lda = nV;
+           B = d_t1 ; ldb = nO;
+           C = d_g_vir; ldc = nV;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+           float* d_tmp_k;
+           cudaStat = cudaMalloc((void**)&d_tmp_k, cholesky_mo_num*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+           alpha = 1.0;
+           beta =  0.0;
+           m=cholesky_mo_num ; n=1; k=nO*nV;
+           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num;
+           B = d_t1 ; ldb = nO*nV;
+           C = d_tmp_k; ldc = cholesky_mo_num;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+           alpha = 2.0;
+           beta =  1.0;
+           m=nV*nV; n=1; k=cholesky_mo_num;
+           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num;
+           B = d_tmp_k ; ldb = cholesky_mo_num;
+           C = d_g_vir; ldc = nV*nV;
+           cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           cudaFree(d_tmp_k);
+
+           float* d_tmp_vo;
+           cudaStat = cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+           alpha = 1.0;
+           beta =  0.0;
+           m=cholesky_mo_num*nV ; n=nO; k=nV;
+           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num*nV;
+           B = d_t1 ; ldb = nO;
+           C = d_tmp_vo; ldc = cholesky_mo_num*nV;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+           float* d_tmp_vo2;
+           cudaStat = cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+           for (size_t i=0 ; i<nO ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (size_t i=0 ; i<nO ; ++i) {
+             cublasSetStream(handle, stream[i]);
+             alpha = -1.0;
+             beta =  0.0;
+             A = &(d_tmp_vo[cholesky_mo_num*nV*i]); lda = cholesky_mo_num;
+             B = &(d_tmp_vo[cholesky_mo_num*nV*i]); ldb = cholesky_mo_num;
+             C = &(d_tmp_vo2[cholesky_mo_num*i]); ldc = cholesky_mo_num*nO;
+             cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+           }
+           for (size_t i=0 ; i<nO ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
+           cudaFree(d_tmp_vo);
+
+           alpha = 1.0;
+           beta =  1.0;
+           m=nV ; n=nV; k=nO*cholesky_mo_num;
+           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num*nO;
+           B = d_tmp_vo2 ; ldb = cholesky_mo_num*nO;
+           C = d_g_vir; ldc = nV;
+           cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           cudaFree(d_tmp_vo2);
+
+           float* d_Y_oovv;
+           cudaStat = cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(float));
+           assert(cudaStat == cudaSuccess);
+           alpha = 1.0;
+           beta =  0.0;
+           m=nO*nO*nV ; n=nV; k=nV;
+           A = d_t2; lda = nO*nO*nV;
+           B = d_g_vir; ldb = nV;
+           C = d_Y_oovv; ldc = nO*nO*nV;
+           cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           cudaFree(d_g_vir);
+
+           alpha = 1.0;
+           beta =  1.0;
+           A = d_r2; lda = nO*nO;
+           B = d_Y_oovv; ldb = nO*nO;
+           C = d_r2; ldc = nO*nO;
+           cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
+                cublasSetStream(handle, stream[i]);
+                alpha = 1.0;
+                beta =  1.0;
+                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+                B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
+                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+             }
+           }
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
+           cudaFree(d_Y_oovv);
+        }
+
+        // g_occ
+        #pragma omp section
+        {
+          float* d_g_occ;
+          lda = nO;
+          cudaStat = cudaMalloc((void **)&d_g_occ, nO*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+          cublasScopy(handle, nO*nO, d_H_oo, 1, d_g_occ, 1);
+
+          float* d_X;
+          cudaStat = cudaMalloc((void **)&d_X, cholesky_mo_num*sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 2.0;
+          beta  = 0.0;
+          m=cholesky_mo_num;  n=nO*nV;
+          A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
+          B=d_t1;  ldb=1;
+          C=d_X;  ldc=1;
+          cublasSgemv(handle, CUBLAS_OP_N, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          alpha = 1.0;
+          beta  = 1.0;
+          m=cholesky_mo_num;  n=nO*nO;
+          A=d_cc_space_v_oo_chol;  lda=cholesky_mo_num;
+          B=d_X;  ldb=1;
+          C=d_g_occ;  ldc=1;
+          cublasSgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X);
+          alpha = -1.0;
+          beta  = 1.0;
+          m=nO*nV;  n=nO*nO;
+          A=d_cc_space_v_ovoo;  lda=nO*nV;
+          B=d_t1;  ldb=1;
+          C=d_g_occ;  ldc=1;
+          cublasSgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          alpha = 1.0;
+          beta  = 1.0;
+          m=nO;  n=nO;  k=nV;
+          A=d_t1;  lda=nO;
+          B=d_cc_space_f_vo;  ldb=nV;
+          C=d_g_occ;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_X_oovv;
+          cudaStat = cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO;  n=nO*nV*nV;  k=nO;
+          A=d_g_occ;  lda=nO;
+          B=d_t2;  ldb=nO;
+          C=d_X_oovv;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_g_occ);
+
+           alpha = 1.0;
+           beta = -1.0;
+           A = d_r2; lda = nO*nO;
+           B = d_X_oovv; ldb = nO*nO;
+           C = d_r2; ldc = nO*nO;
+           cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
+                cublasSetStream(handle, stream[i]);
+                alpha = 1.0;
+                beta = -1.0;
+                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+                B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
+                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+             }
+           }
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
+           cudaFree(d_X_oovv);
+        }
+
+      } // end sections
+
+      lda = nO*nV;
+      cublasSetMatrix(lda, nO*nV, sizeof(float), K1, lda, d_K1, lda);
+
+#define BLOCK_SIZE 16
+
+      #pragma omp sections
+      {
+
+        #pragma omp section
+        {
+           float* d_X_vovv;
+           cudaStat = cudaMalloc((void **)&d_X_vovv, nV*nO*nV*BLOCK_SIZE * sizeof(float));
+           assert(cudaStat == cudaSuccess);
+
+           float* d_Y_oovv;
+           cudaStat = cudaMalloc((void **)&d_Y_oovv, nO*nO*nV*nV * sizeof(float));
+           assert(cudaStat == cudaSuccess);
+
+           for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
+             size_t mbs = nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE;
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
+               cudaStreamCreate(&(stream[gam]));
+             }
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
+               cublasSetStream(handle, stream[gam]);
+               alpha = 1.0;
+               beta  = 0.0;
+               m=nV;  n=nO*nV;  k=cholesky_mo_num;
+               A=&(d_cc_space_v_vv_chol[cholesky_mo_num*nV*gam]);  lda=cholesky_mo_num;
+               B=d_cc_space_v_ov_chol;  ldb=cholesky_mo_num;
+               C=&(d_X_vovv[nV*nO*nV*(gam-iblock)]);  ldc=nV;
+               cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+             }
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
+               cudaStreamDestroy(stream[gam]);
+             }
+             cublasSetStream(handle, NULL);
+             mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+             alpha = 1.0;
+             beta  = 0.0;
+             m=nO;  n=nO*nV*mbs;  k=nV;
+             A=d_t1;  lda=nO;
+             B=d_X_vovv;  ldb=nV;
+             C=&(d_Y_oovv[nO*nO*nV*iblock]);  ldc=nO;
+             cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           }
+           cudaFree(d_X_vovv);
+           alpha = 1.0;
+           beta =  1.0;
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
+                cublasSetStream(handle, stream[i]);
+                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+                B = &(d_Y_oovv[nO*nO*(i+nV*j)]); ldb = nO;
+                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+             }
+             for (size_t i=0 ; i<nV ; ++i) {
+                cublasSetStream(handle, stream[i]);
+                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+                B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
+                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+             }
+           }
+           for (size_t i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
+           cudaFree(d_Y_oovv);
+        }
+
+        #pragma omp section
+        {
+          float* d_tcc2;
+          cudaStat = cudaMalloc((void **)&d_tcc2, cholesky_mo_num*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=cholesky_mo_num*nV;  n=nO;  k=nV;
+          A=d_cc_space_v_vv_chol;  lda=cholesky_mo_num*nV;
+          B=d_t1;  ldb=nO;
+          C=d_tcc2;  ldc=cholesky_mo_num*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_tcc;
+          cudaStat = cudaMalloc((void **)&d_tcc, cholesky_mo_num*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=cholesky_mo_num*nO;  n=nV;  k=nO;
+          A=d_cc_space_v_oo_chol;  lda=cholesky_mo_num*nO;
+          B=d_t1;  ldb=nO;
+          C=d_tcc;  ldc=cholesky_mo_num*nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          float* d_X_ovvo;
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nV*nO;  k=cholesky_mo_num;
+          A=d_tcc;  lda=cholesky_mo_num;
+          B=d_tcc2;  ldb=cholesky_mo_num;
+          C=d_X_ovvo;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_tcc);
+          cudaFree(d_tcc2);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta = -1.0;
+          for(size_t gam = 0; gam < nV; gam++){
+            for(size_t bet = 0; bet < nV; bet++){
+               cublasSetStream(handle, stream[bet]);
+               A = &(d_r2[nO*nO*(bet+nV*gam)]); lda = nO;
+               B = &(d_X_ovvo[nO*(bet+nV*gam)]); ldb = nO*nV*nV;
+               C = &(d_r2[nO*nO*(bet+nV*gam)]); ldc = nO;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for(size_t bet = 0; bet < nV; bet++){
+               cublasSetStream(handle, stream[bet]);
+               A = &(d_r2[nO*nO*(bet+nV*gam)]); lda = nO;
+               B = &(d_X_ovvo[nO*(gam+nV*bet)]); ldb = nO*nV*nV;
+               C = &(d_r2[nO*nO*(bet+nV*gam)]); ldc = nO;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_X_ovvo);
+        }
+
+        #pragma omp section
+        {
+          float* d_X_oovv;
+          cudaStat = cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nO*nV;  n=nV;  k=nO;
+          A=d_cc_space_v_oovo;  lda=nO*nO*nV;
+          B=d_t1;  ldb=nO;
+          C=d_X_oovv;  ldc=nO*nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          alpha = 1.0;
+          beta = -1.0;
+          A = d_r2; lda = nO*nO;
+          B = d_X_oovv; ldb = nO*nO;
+          C = d_r2; ldc = nO*nO;
+          cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t j=0 ; j<nV ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
+               cublasSetStream(handle, stream[i]);
+               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+               B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
+               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+
+          float* d_X_vovo;
+          cudaStat = cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 0.0;
+          beta =  1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t gam=0 ; gam<nV ; ++gam) {
+               cublasSetStream(handle, stream[gam]);
+               A = &(d_X_vovo[nV*nO*(gam+nV*i)]); lda = nV;
+               B = &(d_cc_space_v_ovvo[nO*nV*(gam+nV*i)]); ldb = nO;
+               C = &(d_X_vovo[nV*nO*(gam+nV*i)]); ldc = nV;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Y_oovo;
+          cudaStat = cudaMalloc((void **)&d_Y_oovo, nO*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO;  n=nO*nV*nO;  k=nV;
+          A=d_t1;  lda=nO;
+          B=d_X_vovo;  ldb=nV;
+          C=d_Y_oovo;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_vovo);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nO*nV;  n=nV;  k=nO;
+          A=d_Y_oovo;  lda=nO*nO*nV;
+          B=d_t1;  ldb=nO;
+          C=d_X_oovv;  ldc=nO*nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_Y_oovo);
+
+          alpha = 1.0;
+          beta = -1.0;
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t j=0 ; j<nV ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
+               cublasSetStream(handle, stream[i]);
+               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+               B = &(d_X_oovv[nO*nO*(i+nV*j)]); ldb = nO;
+               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for (size_t i=0 ; i<nV ; ++i) {
+               cublasSetStream(handle, stream[i]);
+               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
+               B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
+               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
+               cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_X_oovv);
+        }
+
+        #pragma omp section
+        {
+          alpha = 1.0;
+          beta =  1.0;
+          A = d_r2; lda = nO*nO;
+          B = d_cc_space_v_oovv; ldb = nO*nO;
+          C = d_r2; ldc = nO*nO;
+          cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+        }
+
+        #pragma omp section
+        {
+          float* d_J1;
+          lda = nO*nV;
+          cudaStat = cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+          cublasSetMatrix(lda, nV*nO, sizeof(float), J1, lda, d_J1, lda);
+
+          float* d_X_ovvo;
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta = -0.5;
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
+              B = &(d_K1[nO*nV*(i+nO*b)]); ldb = nO;
+              C = &(d_X_ovvo[nO*(b+nV*nV*i)]); ldc = nO*nV;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_J1);
+
+          float* d_Y_voov;
+          cudaStat = cudaMalloc((void **)&d_Y_voov, nV*nO*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 2.0;
+          beta = -1.0;
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_t2[nO*(v+nO*nV*g)]); lda = nO*nO;
+              B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
+              C = &(d_Y_voov[nV*nO*(v+nO*g)]); ldc = nV;
+              cublasSgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Z_ovov;
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nO*nV;  k=nV*nO;
+          A=d_X_ovvo;  lda=nO*nV;
+          B=d_Y_voov;  ldb=nV*nO;
+          C=d_Z_ovov;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+          cudaFree(d_X_ovvo);
+          cudaFree(d_Y_voov);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta =  1.0;
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_Z_ovov);
+
+        }
+
+        #pragma omp section
+        {
+          float* d_X_ovov;
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          float* d_Y_ovov;
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 0.5;
+          beta =  0.0;
+          for (size_t a=0 ; a<nV ; ++a) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_K1[nO*(a+nV*nO*b)]); lda = nO*nV;
+              B = &(d_K1[nO*(a+nV*nO*b)]); ldb = nO*nV;
+              C = &(d_X_ovov[nO*(a+nV*nO*b)]); ldc = nO*nV;
+              cublasSgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          alpha = 1.0;
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_t2[nO*(v+nO*g)]); lda = nO*nO*nV;
+              B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
+              C = &(d_Y_ovov[nO*nV*(v+nO*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Z_ovov;
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nO*nV;  k=nO*nV;
+          A=d_X_ovov;  lda=nO*nV;
+          B=d_Y_ovov;  ldb=nO*nV;
+          C=d_Z_ovov;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_ovov);
+          cudaFree(d_Y_ovov);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta = -1.0;
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_Z_ovov);
+
+        }
+
+        #pragma omp section
+        {
+          float* d_X_ovov;
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          float* d_Y_ovov;
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta =  0.0;
+          for (size_t a=0 ; a<nV ; ++a) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_K1[nO*(a+nV*nO*g)]); lda = nO*nV;
+              B = &(d_K1[nO*(a+nV*nO*g)]); ldb = nO*nV;
+              C = &(d_X_ovov[nO*(g+nV*nO*a)]); ldc = nO*nV;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          alpha = 1.0;
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t b=0 ; b<nV ; ++b) {
+              cublasSetStream(handle, stream[b]);
+              A = &(d_t2[nO*(v+nO*b)]); lda = nO*nO*nV;
+              B = &(d_t2[nO*(v+nO*b)]); ldb = nO*nO*nV;
+              C = &(d_Y_ovov[nO*nV*(v+nO*b)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          float* d_Z_ovov;
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 1.0;
+          beta  = 0.0;
+          m=nO*nV;  n=nO*nV;  k=nO*nV;
+          A=d_X_ovov;  lda=nO*nV;
+          B=d_Y_ovov;  ldb=nO*nV;
+          C=d_Z_ovov;  ldc=nO*nV;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_ovov);
+          cudaFree(d_Y_ovov);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta = -1.0;
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for (size_t g=0 ; g<nV ; ++g) {
+              cublasSetStream(handle, stream[g]);
+              A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
+              B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
+              C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_Z_ovov);
+
+        }
+      } // end sections
+
+      cudaFree(d_K1);
+
+      float* d_tmp_cc;
+      lda = cholesky_mo_num * nV;
+      cudaStat = cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(float));
+      assert(cudaStat == cudaSuccess);
+
+      alpha=1.0; beta=0.0;
+      m=cholesky_mo_num*nV;  n=nV;  k=nO;
+      A = d_cc_space_v_vo_chol; B = d_t1; C = d_tmp_cc;
+      cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, m, B, k, &beta, C, m);
+
+      float* d_tmp_cc2;
+      cudaStat = cudaMalloc((void **)&d_tmp_cc2, cholesky_mo_num*nV*sizeof(float));
+      assert(cudaStat == cudaSuccess);
+
+      float* d_B1;
+      cudaStat = cudaMalloc((void**)&d_B1, nV*nV*BLOCK_SIZE*sizeof(float));
+      assert(cudaStat == cudaSuccess);
+
+      float* d_tmpB1;
+      cudaStat = cudaMalloc((void**)&d_tmpB1, nV*BLOCK_SIZE*nV*sizeof(float));
+      assert(cudaStat == cudaSuccess);
+
+      #pragma omp for
+      for (size_t gam=0 ; gam<nV ; ++gam)
+      {
+         float* d_tmp_cc_ = &(d_tmp_cc[gam*nV*cholesky_mo_num]);
+         float* d_cc_space_v_vv_chol_ = &(d_cc_space_v_vv_chol[gam*nV*cholesky_mo_num]);
+
+         alpha = 1.0;
+         beta = -1.0;
+         A = d_cc_space_v_vv_chol_; lda = cholesky_mo_num;
+         B = d_tmp_cc_; ldb = cholesky_mo_num;
+         C = d_tmp_cc2 ; ldc = cholesky_mo_num;
+         cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+         for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE)
+         {
+              const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+
+              alpha=-1.0; beta=0.0;
+              m=nV*mbs;  n=nV;  k=cholesky_mo_num;
+
+              A=&(d_tmp_cc[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
+              B=d_cc_space_v_vv_chol_; ldb=cholesky_mo_num;
+              C=d_tmpB1 ; ldc=nV*BLOCK_SIZE;
+              cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, lda, &beta, C, ldc);
+
+              alpha=1.0; beta=1.0;
+              m=nV*mbs;  n=nV;  k=cholesky_mo_num;
+
+              A=&(d_cc_space_v_vv_chol[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
+              B=d_tmp_cc2; ldb=cholesky_mo_num;
+              C=d_tmpB1 ; ldc=nV*BLOCK_SIZE;
+              cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, lda, &beta, C, ldc);
+
+              for (size_t bet=iblock ; bet<(nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE) ; ++bet)
+              {
+                alpha = 1.0;
+                beta = 0.0;
+                A = &(d_tmpB1[nV*(bet-iblock)]); lda = nV*BLOCK_SIZE;
+                B = d_tmpB1; ldb = nV;
+                C = &(d_B1[nV*nV*(bet-iblock)]) ; ldc = nV;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+             }
+
+             alpha=1.0; beta=1.0;
+             m=nO*nO;  n=mbs;  k=nV*nV;
+             A=d_tau; lda=nO*nO;
+             B=d_B1 ; ldb=nV*nV;
+             C=&(d_r2[nO*nO*(iblock + nV*gam)]); ldc=nO*nO;
+             cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+         }
+      }
+      cudaFree(d_tmpB1);
+      cudaFree(d_B1);
+      cudaFree(d_tmp_cc2);
+      cudaFree(d_tmp_cc);
+
+      float* r2_tmp = malloc((size_t) nO*nO*nV*nV*sizeof(float));
+      assert (r2_tmp != NULL);
+      lda=nO*nO;
+      cublasGetMatrix(nO*nO, nV*nV, sizeof(float), d_r2, lda, r2_tmp, lda);
+      #pragma omp critical
+      {
+         for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
+           r2[i] -= r2_tmp[i];
+         }
+      }
+      free(r2_tmp);
+
+      cudaFree(d_r2);
+
+      cublasDestroy(handle);
+  }
+  free(K1);
+  free(J1);
+
+  *max_r2 = 0.;
+  for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
+    const double x = r2[i] > 0. ? r2[i] : -r2[i];
+    *max_r2 = *max_r2 > x ? *max_r2 : x;
+  }
+
+}
+
+
+
+
+
+
+void compute_r1_space_chol_gpu_sp(gpu_data_sp* data, int nO, int nV, double* t1, double* r1, double* max_r1)
+{
+    const size_t cholesky_mo_num = data->cholesky_mo_num;
+
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      cudaError_t cudaStat;
+      size_t m,n,k, lda, ldb, ldc;
+      float alpha, beta;
+      float* A;
+      float* B;
+      float* C;
+      cudaStream_t stream[nV];
+
+      size_t igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+
+      cublasHandle_t handle;
+      cublasCreate(&handle);
+
+      float* d_r1;
+      lda = nO ;
+      cudaStat = cudaMalloc((void **)&d_r1, lda * nV * sizeof(float));
+      assert(cudaStat == cudaSuccess);
+      cudaMemset(d_r1, 0, nO*nV*sizeof(float));
+      memset(r1, 0, nO*nV*sizeof(double));
+
+      float* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+      float* d_cc_space_v_vv_chol = data[igpu].cc_space_v_vv_chol;
+      float* d_cc_space_v_oovo    = data[igpu].cc_space_v_oovo;
+      float* d_cc_space_v_ovov    = data[igpu].cc_space_v_ovov;
+      float* d_cc_space_v_voov    = data[igpu].cc_space_v_voov;
+      float* d_cc_space_f_ov      = data[igpu].cc_space_f_ov;
+      float* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+      float* d_tau                = data[igpu].tau;
+      float* d_t1                 = data[igpu].t1;
+      float* d_t2                 = data[igpu].t2;
+      float* d_H_oo               = data[igpu].H_oo;
+      float* d_H_vo               = data[igpu].H_vo;
+      float* d_H_vv               = data[igpu].H_vv;
+
+      #pragma omp sections
+      {
+
+        #pragma omp section
+        {
+          cublasScopy(handle, nO*nV, d_cc_space_f_ov, 1, d_r1, 1);
+
+          float* d_X_oo;
+          cudaStat = cudaMalloc((void **)&d_X_oo, nO*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = -2.0;
+          beta  = 0.0;
+          m=nO;  n=nO;  k=nV;
+          A=d_t1;  lda=nO;
+          B=d_cc_space_f_vo;  ldb=nV;
+          C=d_X_oo;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          alpha = 1.0;
+          beta  = 1.0;
+          m=nO;  n=nV;  k=nO;
+          A=d_X_oo;  lda=nO;
+          B=d_t1;  ldb=nO;
+          C=d_r1;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_oo);
+        }
+
+        #pragma omp section
+        {
+          alpha = 1.0;
+          beta  = 1.0;
+          m=nO;  n=nV;  k=nV;
+          A=d_t1;  lda=nO;
+          B=d_H_vv;  ldb=nV;
+          C=d_r1;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+        }
+
+        #pragma omp section
+        {
+          alpha = -1.0;
+          beta  = 1.0;
+          m=nO;  n=nV;  k=nO;
+          A=d_H_oo;  lda=nO;
+          B=d_t1;  ldb=nO;
+          C=d_r1;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+        }
+
+        #pragma omp section
+        {
+          float* d_X_voov;
+          cudaStat = cudaMalloc((void **)&d_X_voov, nV* nO* nO* nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = -1.0;
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t bet=0 ; bet<nV ; ++bet) {
+              cublasSetStream(handle, stream[bet]);
+              beta = t1[i+bet*nO];
+              A = &(d_t2[nO*(i+nO*nV*bet)]); lda = nO*nO;
+              B = &(d_t1[0]); ldb = nO;
+              C = &(d_X_voov[nV*(i+nO*nO*bet)]); ldc = nV*nO;
+              cublasSgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          cudaDeviceSynchronize();
+          alpha = 1.0;
+          beta  = 2.0;
+          for (size_t bet=0 ; bet<nV ; ++bet) {
+            cublasSetStream(handle, stream[bet]);
+            A = &(d_X_voov[nV*nO*nO*bet]); lda = nV;
+            B = &(d_t2[nO*nO*nV*bet]); ldb = nO*nO;
+            C = A ; ldc = lda;
+            cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          alpha = 1.0;
+          beta  = 1.0;
+          m=nV*nO;  n=nO*nV;
+          A=d_X_voov;  lda=nV * nO;
+          B=d_H_vo;  ldb=1;
+          C=d_r1;  ldc=1;
+          cublasSgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_voov);
+        }
+
+        #pragma omp section
+        {
+          float* d_X_ovov;
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO* nV* nO* nV * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+          cublasScopy(handle, nO*nV*nO*nV, d_cc_space_v_ovov, 1, d_X_ovov, 1);
+
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = -1.0;
+          beta  = 2.0;
+          for (size_t u=0 ; u<nO ; ++u) {
+            for (size_t bet=0 ; bet<nV ; ++bet) {
+              cublasSetStream(handle, stream[bet]);
+              A = &(d_X_ovov[nO*nV*(u+nO*bet)]); lda = nO;
+              B = &(d_cc_space_v_voov[(nV*(u+nO*nO*bet))]); ldb = nV*nO;
+              C = A ; ldc = lda;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          alpha = 1.0;
+          beta  = 1.0;
+          m=nO*nV;  n=nO*nV;
+          A=d_X_ovov;  lda=nO * nV;
+          B=d_t1;  ldb=1;
+          C=d_r1;  ldc=1;
+          cublasSgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_ovov);
+        }
+
+
+        #pragma omp section
+        {
+          float* d_T_vvoo;
+          cudaStat = cudaMalloc((void **)&d_T_vvoo, nV*nV*nO*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 0.0;
+          beta  = 1.0;
+          A = d_T_vvoo; lda = nV*nV;
+          B = d_tau; ldb = nO*nO;
+          C = A ; ldc = lda;
+          cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV*nV, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+          float* d_W_vvov;
+          cudaStat = cudaMalloc((void **)&d_W_vvov, nV*nV*nO*BLOCK_SIZE * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          float* d_W_vvov_tmp;
+          cudaStat = cudaMalloc((void **)&d_W_vvov_tmp, nV*nO*nV*BLOCK_SIZE * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+
+          for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
+            const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+
+            alpha = 1.0;
+            beta  = 0.0;
+            m=nV*nO;  n=nV*mbs;  k=cholesky_mo_num;
+            A=d_cc_space_v_vo_chol;  lda=cholesky_mo_num;
+            B=&(d_cc_space_v_vv_chol[cholesky_mo_num*nV*iblock]);  ldb=cholesky_mo_num;
+            C=d_W_vvov_tmp;  ldc=nV*nO;
+            cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+            alpha = 2.0;
+            beta = -1.0;
+            size_t kk=0;
+            for (size_t i=0 ; i<nV ; ++i) {
+              cudaStreamCreate(&(stream[i]));
+            }
+            for (size_t i=0 ; i<nO ; ++i) {
+              for (size_t bet=0 ; bet<mbs ; ++bet) {
+                cublasSetStream(handle, stream[kk]);
+                ++kk;
+                if (kk >= nV) kk = 0;
+                A = &(d_W_vvov_tmp[nV*(i+nO*nV*bet)]); lda = nV*nO;
+                B = &(d_W_vvov_tmp[nV*(i+nO*nV*bet)]); ldb = nV*nO;
+                C = &(d_W_vvov[nV*nV*(i+nO*bet)]); ldc = nV;
+                cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+              }
+            }
+            for (size_t i=0 ; i<nV ; ++i) {
+              cudaStreamDestroy(stream[i]);
+            }
+            cublasSetStream(handle, NULL);
+
+            alpha = 1.0;
+            beta  = 1.0;
+            m=nO;  n=mbs;  k=nO*nV*nV;
+            A=d_T_vvoo;  lda=nV*nV*nO;
+            B=d_W_vvov;  ldb=nO*nV*nV;
+            C=&(d_r1[nO*iblock]);  ldc=nO;
+            cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          }
+          cudaFree(d_W_vvov);
+          cudaFree(d_W_vvov_tmp);
+          cudaFree(d_T_vvoo);
+        }
+
+        #pragma omp section
+        {
+          float* d_W_oovo;
+          cudaStat = cudaMalloc((void **)&d_W_oovo, nO*nO*nV*nO * sizeof(float));
+          assert(cudaStat == cudaSuccess);
+
+          alpha = 2.0;
+          beta  = -1.0;
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (size_t u=0 ; u<nO ; ++u) {
+            for (size_t a=0 ; a<nV ; ++a) {
+              cublasSetStream(handle, stream[a]);
+              A = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); lda = nO;
+              B = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); ldb = nO;
+              C = &(d_W_oovo[nO*nO*(a+nV*u)]); ldc = nO;
+              cublasSgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (size_t i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          alpha = -1.0;
+          beta  = 1.0;
+          m=nO;  n=nV;  k=nO*nO*nV;
+          A=d_W_oovo;  lda=nO * nO * nV;
+          B=d_tau;  ldb=nO * nO * nV;
+          C=d_r1;  ldc=nO;
+          cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+        }
+      }
+
+      float* r1_tmp = malloc((size_t) nO*nV*sizeof(float));
+      assert (r1_tmp != NULL);
+      lda=nO;
+      cublasGetMatrix(nO, nV, sizeof(float), d_r1, lda, r1_tmp, lda);
+      #pragma omp critical
+      {
+         for (size_t i=0 ; i<(size_t) nO*nV ; ++i) {
+           r1[i] -= r1_tmp[i];
+         }
+      }
+      free(r1_tmp);
+
+      cudaFree(d_r1);
+
+      cublasDestroy(handle);
+  }
+
+  *max_r1 = 0.;
+  for (size_t i=0 ; i<(size_t) nO*nV ; ++i) {
+    const double x = r1[i] > 0. ? r1[i] : -r1[i];
+    *max_r1 = *max_r1 > x ? *max_r1 : x;
+  }
+}
+
+
+
+double ccsd_energy_space_gpu_sp(gpu_data_sp* data)
+{
+    double result = 0.0;
+
+    const size_t nO = data->nO;
+    const size_t nV = data->nV;
+
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      cudaError_t cudaStat;
+      size_t igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+
+      cublasHandle_t handle;
+      cublasCreate(&handle);
+
+      double result_local = 0.0;
+      #pragma omp sections
+      {
+        #pragma omp section
+        {
+          float* d_cc_space_f_ov = data[igpu].cc_space_f_ov;
+          float* d_t1 = data[igpu].t1;
+          float x;
+          cublasSdot(handle, nO*nV, d_cc_space_f_ov, 1, d_t1, 1, &x);
+          result_local += 2.0*x;
+        }
+        #pragma omp section
+        {
+          float* d_tau_x = data[igpu].tau_x;
+          float* d_cc_space_v_oovv = data[igpu].cc_space_v_oovv;
+          float x;
+          cublasSdot(handle, nO*nO*nV*nV, d_tau_x, 1, d_cc_space_v_oovv, 1, &x);
+          result_local += x;
+        }
+      }
+      cublasDestroy(handle);
+      #pragma omp critical
+      {
+        result += result_local;
+      }
+  }
+  return result;
+}