More dgemm

devel/ccsd_gpu/ccsd_space_orb_sub.irp.f

@@ -1,5 +1,5 @@
 subroutine run_ccsd_space_orb
-
+  use gpu_module
   implicit none
 
   integer :: i,j,k,l,a,b,c,d,tmp_a,tmp_b,tmp_c,tmp_d
@@ -115,6 +115,13 @@ subroutine run_ccsd_space_orb
   write(*,'(A77)') ' -----------------------------------------------------------------------------'
   call wall_time(ta)
 
+  type(c_ptr) :: gpu_data
+
+  gpu_data = gpu_init(nO, nV, cholesky_mo_num, &
+      cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
+      cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
+      cc_space_v_oovo, cc_space_v_ovvo, cc_space_v_ovoo, cc_space_f_vo)
+
   do while (not_converged)
 
     ! Residue
@@ -125,7 +132,7 @@ subroutine run_ccsd_space_orb
       call compute_H_vo_chol(nO,nV,t1,H_vo)
 
       call compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
-      call compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
+      call compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2,gpu_data)
     else
       call compute_H_oo(nO,nV,t1,t2,tau,H_oo)
       call compute_H_vv(nO,nV,t1,t2,tau,H_vv)

devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f

@@ -441,7 +441,7 @@ end
 
 ! R2
 
-subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
+subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2,gpu_data)
   use gpu_module
   implicit none
 
@@ -449,6 +449,7 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   integer, intent(in)           :: nO, nV
   double precision, intent(in)  :: t1(nO,nV), t2(nO,nO,nV,nV), tau(nO,nO,nV,nV)
   double precision, intent(in)  :: H_oo(nO,nO), H_vv(nV,nV), H_vo(nV,nO)
+  type(c_ptr), intent(in) :: gpu_data
 
   ! out
   double precision, intent(out) :: r2(nO,nO,nV,nV), max_r2
@@ -489,10 +490,12 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   !$omp end do
   !$omp end parallel
 
-  call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,tau, &
+  double precision, allocatable :: J1(:,:,:,:)
-      cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
+  allocate(J1(nO,nV,nV,nO))
-      cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
+  J1 = 0.d0
-      cc_space_f_vo, H_vv, g_occ, r2)
+
+  call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,gpu_data,t1,t2,tau, &
+      H_vv, g_occ, J1, r2)
 
 
 !---
@@ -501,86 +504,8 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   double precision, allocatable :: X_ovvo(:,:,:,:)
   double precision, allocatable :: tcc(:,:,:), tcc2(:,:,:)
 
-  allocate(X_oovv(nO,nO,nV,nV))
-
-  call dgemm('N','N',nO*nO*nV,nV,nO, &
-             1d0, cc_space_v_oovo, size(cc_space_v_oovo,1) * size(cc_space_v_oovo,2) * size(cc_space_v_oovo,3), &
-                  t1 , size(t1,1), &
-             0d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,r2,X_oovv) &
-  !$omp private(u,v,gam,beta) &
-  !$omp default(none)
-  !$omp do
-  do gam = 1, nV
-    do beta = 1, nV
-      do v = 1, nO
-        do u = 1, nO
-           r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(u,v,beta,gam) - X_oovv(v,u,gam,beta)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-  deallocate(X_oovv)
-
   double precision, allocatable :: X_vovo(:,:,:,:), Y_oovo(:,:,:,:)
-  allocate(X_vovo(nV,nO,nV,nO))
 
-  !$omp parallel &
-  !$omp shared(nO,nV,X_vovo,cc_space_v_ovvo) &
-  !$omp private(a,v,gam,i) &
-  !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do gam = 1, nV
-      do v = 1, nO
-        do a = 1, nV
-          X_vovo(a,v,gam,i) = cc_space_v_ovvo(v,a,gam,i)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
-  !$omp end parallel
-
-  allocate(Y_oovo(nO,nO,nV,nO))
-  call dgemm('N','N',nO,nO*nV*nO,nV, &
-             1d0, t1, size(t1,1), &
-                  X_vovo, size(X_vovo,1), &
-             0d0, Y_oovo, size(Y_oovo,1))
-
-  deallocate(X_vovo)
-  allocate(X_oovv(nO,nO,nV,nV))
-  call dgemm('N','N',nO*nO*nV, nV, nO, &
-             1d0, Y_oovo, size(Y_oovo,1) * size(Y_oovo,2) * size(Y_oovo,3), &
-                  t1    , size(t1,1), &
-             0d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3))
-  deallocate(Y_oovo)
-
-  !$omp parallel &
-  !$omp shared(nO,nV,r2,X_oovv) &
-  !$omp private(u,v,gam,beta) &
-  !$omp default(none)
-  !$omp do
-  do gam = 1, nV
-    do beta = 1, nV
-      do v = 1, nO
-        do u = 1, nO
-           r2(u,v,beta,gam) = r2(u,v,beta,gam) - X_oovv(u,v,gam,beta) - X_oovv(v,u,beta,gam)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-  deallocate(X_oovv)
-
-
-  double precision, allocatable :: J1(:,:,:,:)
-  allocate(J1(nO,nV,nV,nO))
   call compute_J1_chol(nO,nV,t1,t2,cc_space_v_ovvo,cc_space_v_ovoo, &
        cc_space_v_vvoo,J1)
 
@@ -810,146 +735,28 @@ subroutine compute_J1_chol(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvoo,J1)
   double precision, intent(in)  :: t2(nO, nO, nV, nV)
   double precision, intent(in)  :: v_ovvo(nO,nV,nV,nO), v_ovoo(nO,nV,nO,nO)
   double precision, intent(in)  :: v_vvoo(nV,nV,nO,nO)
-  double precision, intent(out) :: J1(nO, nV, nV, nO)
+  double precision, intent(inout) :: J1(nO, nV, nV, nO)
 
   integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
 
   double precision, allocatable :: X_ovoo(:,:,:,:), Y_ovov(:,:,:,:)
   allocate(X_ovoo(nO,nV,nO,nO),Y_ovov(nO,nV,nO,nV))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,J1,v_ovvo,v_ovoo,X_ovoo) &
-  !$omp private(i,j,a,u,beta) &
-  !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do beta = 1, nV
-      do a = 1, nV
-        do u = 1, nO
-          J1(u,a,beta,i) = v_ovvo(u,a,beta,i)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
-
-  !$omp do collapse(2)
-  do j = 1, nO
-    do i = 1, nO
-      do a = 1, nV
-        do u = 1, nO
-          X_ovoo(u,a,i,j) = v_ovoo(u,a,j,i)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-  call dgemm('N','N',nO*nV*nO,nV,nO, &
-            -1d0, X_ovoo, size(X_ovoo,1) * size(X_ovoo,2) * size(X_ovoo,3), &
-                  t1    , size(t1,1), &
-             0d0, Y_ovov, size(Y_ovov,1) * size(Y_ovov,2) * size(Y_ovov,3))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,J1,Y_ovov) &
-  !$omp private(i,beta,a,u) &
-  !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do beta = 1, nV
-      do a = 1, nV
-        do u = 1, nO
-          J1(u,a,beta,i) = J1(u,a,beta,i) + Y_ovov(u,a,i,beta)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
-  !$omp end parallel
   deallocate(X_ovoo)
 
   double precision, allocatable :: tmp_cc(:,:,:), J1_tmp(:,:,:,:)
-  allocate(tmp_cc(cholesky_mo_num,nV,nO), J1_tmp(nV,nO,nV,nO))
 
-  call dgemm('N','T', cholesky_mo_num*nV, nO, nV, 1.d0, &
-      cc_space_v_vv_chol, cholesky_mo_num*nV, &
-      t1, nO, &
-      0.d0, tmp_cc, cholesky_mo_num*nV)
-
-  call dgemm('T','N', nV*nO, nV*nO, cholesky_mo_num, 1.d0, &
-      tmp_cc, cholesky_mo_num, cc_space_v_vo_chol, cholesky_mo_num, &
-      0.d0, J1_tmp, nV*nO)
-
-  deallocate(tmp_cc)
-
-  do i=1,nO
-    do b=1,nV
-      do a=1,nV
-        do u=1,nO
-          J1(u,a,b,i) = J1(u,a,b,i) + J1_tmp(b,u,a,i)
-        enddo
-      enddo
-    enddo
-  enddo
-
-  deallocate(J1_tmp)
 
+!TODO: I am here
   !- cc_space_v_vvoo(a,b,i,j) * (0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)) &
   double precision, allocatable :: X_voov(:,:,:,:), Z_ovvo(:,:,:,:)
   allocate(X_voov(nV,nO,nO,nV), Z_ovvo(nO,nV,nV,nO))
-  !$omp parallel &
-  !$omp shared(nO,nV,t2,t1,Y_ovov,X_voov,v_vvoo) &
-  !$omp private(i,beta,a,u,b,j) &
-  !$omp default(none)
-  !$omp do
-  do b = 1, nV
-    do j = 1, nO
-      do beta = 1, nV
-        do u = 1, nO
-          Y_ovov(u,beta,j,b) = 0.5d0 * t2(u,j,b,beta) + t1(u,b) * t1(j,beta)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do nowait
-
-  !$omp do
-  do b = 1, nV
-    do j = 1, nO
-      do i = 1, nO
-        do a = 1, nV
-          X_voov(a,i,j,b) = v_vvoo(a,b,i,j)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-  call dgemm('N','T',nO*nV,nV*nO,nO*nV, &
-             -1d0, Y_ovov, size(Y_ovov,1) * size(Y_ovov,2), &
-                   X_voov, size(X_voov,1) * size(X_voov,2), &
-              0d0, Z_ovvo, size(Z_ovvo,1) * size(Z_ovvo,2))
-  deallocate(X_voov)
-
   double precision, allocatable :: X_ovvo(:,:,:,:), Y_vovo(:,:,:,:)
-  allocate(X_ovvo(nO,nV,nV,nO),Y_vovo(nV,nO,nV,nO))
+  allocate(X_ovvo(nO,nV,nV,nO))
+  allocate(Y_vovo(nV,nO,nV,nO))
   !$omp parallel &
   !$omp shared(nO,nV,J1,Z_ovvo,t2,Y_vovo,v_vvoo,X_ovvo) &
   !$omp private(i,beta,a,u,j,b) &
   !$omp default(none)
-  do i = 1, nO
-    !$omp do
-    do beta = 1, nV
-      do a = 1, nV
-        do u = 1, nO
-          J1(u,a,beta,i) = J1(u,a,beta,i) + Z_ovvo(u,beta,a,i)
-        enddo
-      enddo
-    enddo
-    !$omp end do nowait
-  enddo
 
   !+ 0.5d0 * (2d0 * cc_space_v_vvoo(a,b,i,j) - cc_space_v_vvoo(b,a,i,j)) * t2(u,j,beta,b)
   do j = 1, nO
@@ -999,7 +806,7 @@ subroutine compute_J1_chol(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvoo,J1)
   enddo
   !$omp end parallel
 
-  deallocate(X_ovvo,Z_ovvo,Y_ovov)
+  deallocate(X_ovvo,Z_ovvo)
 
 end
 

devel/ccsd_gpu/gpu.c

@@ -4,81 +4,25 @@
 #include <omp.h>
 #include <cublas_v2.h>
 #include <cuda_runtime.h>
+#include "gpu.h"
 
 #define BLOCK_SIZE 16
 
-void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*,
-		double*, double*, int*);
-
-
-
-
-void gpu_dgemm(char transa, char transb, int m, int n, int k, double alpha,
-                  double* A, int lda, double* B, int ldb, double beta, double* C, int ldc)
-{
-    cublasHandle_t handle;
-    cublasCreate(&handle);
-
-    double * d_A;
-    double * d_B;
-    double * d_C;
-    cublasOperation_t ta, tb;
-
-    if (transa == 'N') {
-      cudaMalloc((void**)&d_A, lda*k*sizeof(double));
-      cublasSetMatrix(m, k, sizeof(double), A, lda, d_A, lda);
-      ta = CUBLAS_OP_N;
-    } else {
-      cudaMalloc((void**)&d_A, lda*m*sizeof(double));
-      cublasSetMatrix(k, m, sizeof(double), A, lda, d_A, lda);
-      ta = CUBLAS_OP_T;
-    }
-
-    if (transb == 'N') {
-      cudaMalloc((void**)&d_B, ldb*n*sizeof(double));
-      cublasSetMatrix(k, n, sizeof(double), B, ldb, d_B, ldb);
-      tb = CUBLAS_OP_N;
-    } else {
-      cudaMalloc((void**)&d_B, ldb*k*sizeof(double));
-      cublasSetMatrix(n, k, sizeof(double), B, ldb, d_B, ldb);
-      tb = CUBLAS_OP_T;
-    }
-
-    cudaMalloc((void**)&d_C, ldc*n*sizeof(double));
-    if (beta != 0.) {
-      cublasSetMatrix(m, n, sizeof(double), C, ldc, d_C, ldc);
-    }
-
-    cublasDgemm(handle, ta, tb, m, n, k, &alpha, d_A, lda, d_B, ldb, &beta, d_C, ldc);
-
-    cublasGetMatrix(m, n,  sizeof(double), d_C, ldc, C, ldc);
-
-    cudaFree(d_A);
-    cudaFree(d_B);
-    cudaFree(d_C);
-    cublasDestroy(handle);
-}
 
 void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
+                gpu_data* data,
                 double* t1,
                 double* t2,
                 double* tau,
-                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_oovv,
-                double* cc_space_v_vvoo,
-                double* cc_space_f_vo,
                 double* H_vv,
                 double* g_occ,
+                double* J1,
                 double* r2)
 {
 
-    int ngpus;
+    int ngpus = 1;
     cudaGetDeviceCount(&ngpus);
+
     #pragma omp parallel num_threads(ngpus)
     {
       int m,n,k, lda, ldb, ldc;
@@ -87,20 +31,30 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       double* B;
       double* C;
 
-      int ithread = omp_get_thread_num();
+      int igpu = omp_get_thread_num();
-      int igpu = ithread ;
+
-//igpu=1;
       cudaSetDevice(igpu);
       cublasHandle_t handle;
 
       cublasCreate(&handle);
 
+      cudaStream_t stream[nV];
+
+      double* d_cc_space_v_oo_chol = data[igpu].cc_space_v_oo_chol;
+      double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+      double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+      double* d_cc_space_v_vv_chol = data[igpu].cc_space_v_vv_chol;
+      double* d_cc_space_v_oooo    = data[igpu].cc_space_v_oooo;
+      double* d_cc_space_v_vooo    = data[igpu].cc_space_v_vooo;
+      double* d_cc_space_v_oovv = data[igpu].cc_space_v_oovv;
+      double* d_cc_space_v_vvoo    = data[igpu].cc_space_v_vvoo;
+      double* d_cc_space_v_oovo    = data[igpu].cc_space_v_oovo;
+      double* d_cc_space_v_ovvo    = data[igpu].cc_space_v_ovvo;
+      double* d_cc_space_v_ovoo    = data[igpu].cc_space_v_ovoo;
+      double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+
       double* d_tau;
       double* d_r2;
-      double* d_cc_space_v_oo_chol;
-      double* d_cc_space_v_ov_chol;
-      double* d_cc_space_v_vo_chol;
-      double* d_cc_space_v_vv_chol;
       double* d_t1;
       double* d_t2;
       double* d_tmp_cc;
@@ -114,22 +68,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       memset(r2, 0, nO*nO*nV*nV*sizeof(double));
       cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
 
-      lda = cholesky_mo_num * nV;
-      cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
-
-      lda = cholesky_mo_num * nO;
-      cudaMalloc((void **)&d_cc_space_v_oo_chol, lda * nO * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nO, nO, sizeof(double), cc_space_v_oo_chol, lda, d_cc_space_v_oo_chol, lda);
-
-      lda = cholesky_mo_num * nO;
-      cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
-
-      lda = cholesky_mo_num * nV;
-      cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(double), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
-
       lda = nO;
       cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
       cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
@@ -138,19 +76,12 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
       cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
 
-      double* d_cc_space_f_vo;
-      cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
-      cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
 
       #pragma omp sections
       {
 
         #pragma omp section
         {
-           double* d_cc_space_v_vooo;
-           cudaMalloc((void**)&d_cc_space_v_vooo, nV*nO*nO*nO*sizeof(double));
-           cublasSetMatrix(nV*nO, nO*nO, sizeof(double), cc_space_v_vooo, nV*nO, d_cc_space_v_vooo, nV*nO);
-
            double* d_Y_oooo;
            cudaMalloc((void**)&d_Y_oooo, nO*nO*nO*nO*sizeof(double));
 
@@ -161,23 +92,22 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            B = d_cc_space_v_vooo ; ldb = nV;
            C = d_Y_oooo; ldc = nO;
            cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_cc_space_v_vooo);
 
            double* d_A1;
            cudaMalloc((void**)&d_A1, nO*nO*nO*nO*sizeof(double));
 
-           double* d_cc_space_v_oooo;
-           cudaMalloc((void**)&d_cc_space_v_oooo, nO*nO*nO*nO*sizeof(double));
-           cublasSetMatrix(nO*nO, nO*nO, sizeof(double), cc_space_v_oooo, nO*nO, d_cc_space_v_oooo, nO*nO);
-
            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;
            cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+           for (int i=0 ; i<nO ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
            for (int j=0 ; j<nO ; ++j) {
              for (int 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;
@@ -187,12 +117,12 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
 
              }
            }
+           for (int i=0 ; i<nO ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
            cudaFree(d_Y_oooo);
 
-           double* d_cc_space_v_vvoo;
-           cudaMalloc((void**)&d_cc_space_v_vvoo, nV*nV*nO*nO*sizeof(double));
-           cublasSetMatrix(nV*nV, nO*nO, sizeof(double), cc_space_v_vvoo, nV*nV, d_cc_space_v_vvoo, nV*nV);
-
            alpha = 1.0;
            beta =  1.0;
            m=nO*nO ; n=nO*nO; k=nV*nV;
@@ -200,7 +130,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            B = d_cc_space_v_vvoo ; ldb = nV*nV;
            C = d_A1; ldc = nO*nO;
            cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_cc_space_v_vvoo);
 
            alpha = 1.0;
            beta =  0.0;
@@ -259,6 +188,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            double* d_tmp_vo2;
            cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(double));
            for (int i=0 ; i<nO ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (int 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;
@@ -266,6 +199,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
              C = &(d_tmp_vo2[cholesky_mo_num*i]); ldc = cholesky_mo_num*nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
+           for (int i=0 ; i<nO ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
            cudaFree(d_tmp_vo);
 
            alpha = 1.0;
@@ -294,8 +231,43 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            B = d_Y_oovv; ldb = nO*nO;
            C = d_r2; ldc = nO*nO;
            cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+/*
+           double * Y_oovv = malloc(nO*nO*nV*nV*sizeof(double));
+           lda=nO*nO;
+           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_Y_oovv, lda, Y_oovv, lda);
+           cudaFree(d_Y_oovv);
+
+           double * r2_tmp = malloc(nO*nO*nV*nV*sizeof(double));
+           lda=nO*nO;
+           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, lda, r2_tmp, lda);
+
            for (int j=0 ; j<nV ; ++j) {
              for (int i=0 ; i<nV ; ++i) {
+               double * xx = &(r2_tmp[nO*nO*(i+nV*j)]);
+               const double * yy = &(Y_oovv[nO*nO*(j+nV*i)]);
+               for (int k=0 ; k<nO ; ++k) {
+                 for (int l=0 ; l<nO ; ++l) {
+                   xx[l + k*nO] += yy[k + l*nO];
+                 }
+               }
+             }
+           }
+           free(Y_oovv);
+
+           lda=nO*nO;
+           cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2_tmp, lda, d_r2, lda);
+           free(r2_tmp);
+
+*/
+
+//--
+           for (int i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (int j=0 ; j<nV ; ++j) {
+             for (int 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;
@@ -305,7 +277,12 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
 
              }
            }
+           for (int i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
            cudaFree(d_Y_oovv);
+//--
         }
 
         // g_occ
@@ -343,8 +320,42 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            B = d_X_oovv; ldb = nO*nO;
            C = d_r2; ldc = nO*nO;
            cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+/*
+           double * X_oovv = malloc(nO*nO*nV*nV*sizeof(double));
+           lda=nO*nO;
+           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_X_oovv, lda, X_oovv, lda);
+           cudaFree(d_X_oovv);
+
+           double * r2_tmp = malloc(nO*nO*nV*nV*sizeof(double));
+           lda=nO*nO;
+           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, lda, r2_tmp, lda);
+
            for (int j=0 ; j<nV ; ++j) {
              for (int i=0 ; i<nV ; ++i) {
+               double * xx = &(r2_tmp[nO*nO*(i+nV*j)]);
+               const double * yy = &(X_oovv[nO*nO*(j+nV*i)]);
+               for (int k=0 ; k<nO ; ++k) {
+                 for (int l=0 ; l<nO ; ++l) {
+                   xx[l + k*nO] -= yy[k + l*nO];
+                 }
+               }
+             }
+           }
+           free(X_oovv);
+
+           lda=nO*nO;
+           cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2_tmp, lda, d_r2, lda);
+           free(r2_tmp);
+*/
+
+//--
+           for (int i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+           for (int j=0 ; j<nV ; ++j) {
+             for (int 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;
@@ -353,7 +364,13 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                 cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
              }
            }
+           for (int i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
            cudaFree(d_X_oovv);
+//--
+
         }
 
         #pragma omp section
@@ -369,6 +386,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            for (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
              int mbs = nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE;
              for (int gam=iblock ; gam<mbs ; ++gam) {
+               cudaStreamCreate(&(stream[gam]));
+             }
+             for (int gam=iblock ; gam<mbs ; ++gam) {
+               cublasSetStream(handle, stream[gam]);
                alpha = 1.0;
                beta  = 0.0;
                m=nV;  n=nO*nV;  k=cholesky_mo_num;
@@ -377,6 +398,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                C=&(d_X_vovv[nV*nO*nV*(gam-iblock)]);  ldc=nV;
                cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
              }
+             for (int 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;
@@ -387,20 +412,31 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
              cublasDgemm(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 (int i=0 ; i<nV ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
            for (int j=0 ; j<nV ; ++j) {
              for (int i=0 ; i<nV ; ++i) {
-                alpha = 1.0;
+                cublasSetStream(handle, stream[i]);
-                beta =  1.0;
                 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;
                 cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+             }
+             for (int 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;
                 cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
              }
            }
+           for (int i=0 ; i<nV ; ++i) {
+             cudaStreamDestroy(stream[i]);
+           }
+           cublasSetStream(handle, NULL);
            cudaFree(d_Y_oovv);
         }
 
@@ -444,23 +480,351 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
           cudaFree(d_tcc);
           cudaFree(d_tcc2);
 
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 1.0;
+          beta = -1.0;
           for(int gam = 0; gam < nV; gam++){
             for(int bet = 0; bet < nV; bet++){
-               alpha = 1.0;
+               cublasSetStream(handle, stream[bet]);
-               beta = -1.0;
                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;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for(int 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;
-                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+               C = &(d_r2[nO*nO*(bet+nV*gam)]); ldc = nO;
+               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
-
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
           cudaFree(d_X_ovvo);
         }
 
-      }
+        #pragma omp section
+        {
+          double* d_X_oovv;
+          lda = nO*nO;
+          cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));
+
+          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;
+          cublasDgemm(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;
+          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int j=0 ; j<nV ; ++j) {
+            for (int 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;
+               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          double* d_X_vovo;
+          lda = nV*nO;
+          cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(double));
+
+          alpha = 0.0;
+          beta =  1.0;
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int i=0 ; i<nO ; ++i) {
+            for (int 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;
+               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          double* d_Y_oovo;
+          lda = nO*nO;
+          cudaMalloc((void **)&d_Y_oovo, nO*nO*nV*nO * sizeof(double));
+
+          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;
+          cublasDgemm(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;
+          cublasDgemm(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 (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int j=0 ; j<nV ; ++j) {
+            for (int 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;
+               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+            for (int 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;
+               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int 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;
+          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+        }
+
+        #pragma omp section
+        {
+
+          double* d_J1;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+
+          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;
+          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nV*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+
+
+          double* d_X_ovoo;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(double));
+          alpha = 0.0;
+          beta =  1.0;
+           for (int i=0 ; i<nO ; ++i) {
+             cudaStreamCreate(&(stream[i]));
+           }
+          for (int j=0 ; j<nO ; ++j) {
+            for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nO ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+
+          double* d_Y_ovov;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+
+          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;
+          cublasDgemm(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 (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int j=0 ; j<nO ; ++j) {
+            for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          double* d_tmp_cc;
+          lda = cholesky_mo_num;
+          cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(double));
+
+          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;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          double* d_J1_tmp;
+          lda = nV*nO;
+          cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(double));
+
+          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;
+          cublasDgemm(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 (int i=0 ; i<nO ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO*nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          }
+          for (int i=0 ; i<nO ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_J1_tmp);
+
+//---
+          double* d_X_voov;
+          lda = nV*nO;
+          cudaMalloc((void **)&d_X_voov, nV*nO*nO*nV * sizeof(double));
+
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = 0.5;
+          for (int j=0 ; j<nO ; ++j) {
+            for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          alpha = 1.0;
+          beta  = 0.0;
+          for (int j=0 ; j<nO ; ++j) {
+            for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+
+          double* d_Z_ovvo;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_Z_ovvo, nO*nV*nV*nO * sizeof(double));
+
+          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;
+          cublasDgemm(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 (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          for (int i=0 ; i<nO ; ++i) {
+            for (int 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;
+              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
+            }
+          }
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamDestroy(stream[i]);
+          }
+          cublasSetStream(handle, NULL);
+          cudaFree(d_Z_ovvo);
+
+          lda = nO*nV;
+          cublasGetMatrix(nO*nV, nV*nO, sizeof(double), d_J1, lda, J1, lda);
+          cudaFree(d_J1);
+
+        }
+
+        #pragma omp section
+        {
+        }
+
+      } // end sections
 
 
       lda = cholesky_mo_num * nV;
@@ -537,8 +901,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cudaFree(d_B1);
       cudaFree(d_tmp_cc2);
 
-      cudaFree(d_cc_space_v_vo_chol);
-      cudaFree(d_cc_space_v_vv_chol);
       cudaFree(d_tau);
       cudaFree(d_t1);
       cudaFree(d_tmp_cc);
@@ -554,179 +916,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       free(r2_tmp);
 
       cudaFree(d_r2);
+
       cublasDestroy(handle);
   }
 
-  for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i)
-  {
-     r2[i] += cc_space_v_oovv[i];
-  }
-}
-
-
-
-
-
-
-void compute_g_vir_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
-                double* t1,
-                double* t2,
-                double* tau,
-                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_oovv,
-                double* cc_space_v_vvoo,
-                double* cc_space_f_vo,
-                double* H_vv,
-                double* r2)
-{
-      int m,n,k, lda, ldb, ldc;
-      double alpha, beta;
-      double* A;
-      double* B;
-      double* C;
-
-      int ithread = omp_get_thread_num();
-      int igpu = ithread ;
-
-      cudaSetDevice(igpu);
-      cublasHandle_t handle;
-
-      cublasCreate(&handle);
-
-      double* d_tau;
-      double* d_r2;
-      double* d_cc_space_v_vv_chol;
-      double* d_cc_space_v_vo_chol;
-      double* d_cc_space_v_ov_chol;
-      double* d_t1;
-      double* d_t2;
-      double* d_tmp_cc;
-
-
-
-      lda = nO * nO;
-      cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
-      cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
-
-      lda = nO*nO;
-      cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
-      cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
-      lda = nO;
-      cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
-      cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
-      lda = cholesky_mo_num * nV;
-      cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
-
-      lda = cholesky_mo_num * nO;
-      cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
-      cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
-
- // ---
-           double* d_g_vir;
-           cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
-           cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_g_vir, nV);
-
-           double* d_cc_space_f_vo;
-           cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
-           cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
-
-           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;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_cc_space_f_vo);
-
-           double* d_tmp_k;
-           cudaMalloc((void**)&d_tmp_k, cholesky_mo_num*sizeof(double));
-           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;
-           cublasDgemm(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;
-           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_tmp_k);
-
-
-           double* d_tmp_vo;
-           cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(double));
-           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;
-           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-           double* d_tmp_vo2;
-           cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(double));
-           for (int i=0 ; i<nO ; ++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;
-             cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
-           }
-           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;
-           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-           cudaFree(d_tmp_vo2);
-
-           double* d_Y_oovv;
-           cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(double));
-           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;
-           cublasDgemm(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;
-           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
-           for (int j=0 ; j<nV ; ++j) {
-             for (int i=0 ; i<nV ; ++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;
-                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
-
-             }
-           }
-           cudaFree(d_Y_oovv);
-
-           cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, nO*nO, r2, nO*nO);
 }
 
 

devel/ccsd_gpu/gpu.h

@@ -0,0 +1,15 @@
+typedef struct {
+   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_oovv;
+   double* cc_space_v_vvoo;
+   double* cc_space_v_oovo;
+   double* cc_space_v_ovvo;
+   double* cc_space_v_ovoo;
+   double* cc_space_f_vo;
+} gpu_data;
+

devel/ccsd_gpu/gpu_dgemm.c

@@ -0,0 +1,61 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <omp.h>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+
+#define BLOCK_SIZE 16
+
+void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*,
+		double*, double*, int*);
+
+
+
+
+void gpu_dgemm(char transa, char transb, int m, int n, int k, double alpha,
+                  double* A, int lda, double* B, int ldb, double beta, double* C, int ldc)
+{
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    double * d_A;
+    double * d_B;
+    double * d_C;
+    cublasOperation_t ta, tb;
+
+    if (transa == 'N') {
+      cudaMalloc((void**)&d_A, lda*k*sizeof(double));
+      cublasSetMatrix(m, k, sizeof(double), A, lda, d_A, lda);
+      ta = CUBLAS_OP_N;
+    } else {
+      cudaMalloc((void**)&d_A, lda*m*sizeof(double));
+      cublasSetMatrix(k, m, sizeof(double), A, lda, d_A, lda);
+      ta = CUBLAS_OP_T;
+    }
+
+    if (transb == 'N') {
+      cudaMalloc((void**)&d_B, ldb*n*sizeof(double));
+      cublasSetMatrix(k, n, sizeof(double), B, ldb, d_B, ldb);
+      tb = CUBLAS_OP_N;
+    } else {
+      cudaMalloc((void**)&d_B, ldb*k*sizeof(double));
+      cublasSetMatrix(n, k, sizeof(double), B, ldb, d_B, ldb);
+      tb = CUBLAS_OP_T;
+    }
+
+    cudaMalloc((void**)&d_C, ldc*n*sizeof(double));
+    if (beta != 0.) {
+      cublasSetMatrix(m, n, sizeof(double), C, ldc, d_C, ldc);
+    }
+
+    cublasDgemm(handle, ta, tb, m, n, k, &alpha, d_A, lda, d_B, ldb, &beta, d_C, ldc);
+
+    cublasGetMatrix(m, n,  sizeof(double), d_C, ldc, C, ldc);
+
+    cudaFree(d_A);
+    cudaFree(d_B);
+    cudaFree(d_C);
+    cublasDestroy(handle);
+}
+

devel/ccsd_gpu/gpu_init.c

@@ -0,0 +1,105 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <omp.h>
+#include <cublas_v2.h>
+#include <cuda_runtime.h>
+#include "gpu.h"
+
+#define BLOCK_SIZE 16
+
+gpu_data* gpu_init(
+   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_oovv, double* cc_space_v_vvoo,
+   double* cc_space_v_oovo, double* cc_space_v_ovvo,
+   double* cc_space_v_ovoo, double* cc_space_f_vo)
+{
+    int ngpus = 1;
+    cudaGetDeviceCount(&ngpus);
+
+    gpu_data* data = (gpu_data*) malloc (ngpus*sizeof(gpu_data));
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      int lda;
+      int igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+      cublasHandle_t handle;
+
+      cublasCreate(&handle);
+
+      double* d_cc_space_v_oo_chol;
+      lda = cholesky_mo_num * nO;
+      cudaMalloc((void **)&d_cc_space_v_oo_chol, lda * nO * sizeof(double));
+      cublasSetMatrix(cholesky_mo_num*nO, nO, sizeof(double), cc_space_v_oo_chol, lda, d_cc_space_v_oo_chol, lda);
+
+      double* d_cc_space_v_ov_chol;
+      lda = cholesky_mo_num * nO;
+      cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
+      cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
+
+      double* d_cc_space_v_vo_chol;
+      lda = cholesky_mo_num * nV;
+      cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
+      cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(double), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
+
+      double* d_cc_space_v_vv_chol;
+      lda = cholesky_mo_num * nV;
+      cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
+      cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
+
+      double* d_cc_space_v_oooo;
+      cudaMalloc((void**)&d_cc_space_v_oooo, nO*nO*nO*nO*sizeof(double));
+      cublasSetMatrix(nO*nO, nO*nO, sizeof(double), cc_space_v_oooo, nO*nO, d_cc_space_v_oooo, nO*nO);
+
+      double* d_cc_space_v_vooo;
+      cudaMalloc((void**)&d_cc_space_v_vooo, nV*nO*nO*nO*sizeof(double));
+      cublasSetMatrix(nV*nO, nO*nO, sizeof(double), cc_space_v_vooo, nV*nO, d_cc_space_v_vooo, nV*nO);
+
+      double* d_cc_space_v_oovv;
+      cudaMalloc((void**)&d_cc_space_v_oovv, nO*nO*nV*nV*sizeof(double));
+      cublasSetMatrix(nO*nO, nV*nV, sizeof(double), cc_space_v_oovv, nO*nO, d_cc_space_v_oovv, nO*nO);
+
+      double* d_cc_space_v_vvoo;
+      cudaMalloc((void**)&d_cc_space_v_vvoo, nV*nV*nO*nO*sizeof(double));
+      cublasSetMatrix(nV*nV, nO*nO, sizeof(double), cc_space_v_vvoo, nV*nV, d_cc_space_v_vvoo, nV*nV);
+
+      double* d_cc_space_v_oovo;
+      lda = nO*nO;
+      cudaMalloc((void **)&d_cc_space_v_oovo, nO*nO*nV*nO * sizeof(double));
+      cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_oovo, lda, d_cc_space_v_oovo, lda);
+
+      double* d_cc_space_v_ovvo;
+      lda = nO*nV;
+      cudaMalloc((void **)&d_cc_space_v_ovvo, nO*nV*nV*nO * sizeof(double));
+      cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_ovvo, lda, d_cc_space_v_ovvo, lda);
+
+      double* d_cc_space_v_ovoo;
+      lda = nO*nV;
+      cudaMalloc((void **)&d_cc_space_v_ovoo, nO*nV*nO*nO * sizeof(double));
+      cublasSetMatrix(lda, nO*nO, sizeof(double), cc_space_v_ovoo, lda, d_cc_space_v_ovoo, lda);
+
+      double* d_cc_space_f_vo;
+      cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
+      cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
+
+      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_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_ovoo = d_cc_space_v_ovoo;
+      data[igpu].cc_space_f_vo = d_cc_space_f_vo;
+    }
+    return data;
+}
+
+

devel/ccsd_gpu/gpu_module.f90

@@ -4,15 +4,12 @@ module gpu_module
   implicit none
 
   interface
-    subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
+    type(c_ptr) function gpu_init(nO, nV, cholesky_mo_num, &
-          cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
+      cc_space_v_oo_chol, cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
-          cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
+      cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
-          cc_space_f_vo, H_vv, g_occ, r2) bind(C)
+      cc_space_v_oovo, cc_space_v_ovvo, cc_space_v_ovoo, cc_space_f_vo) bind(C)
-        import c_int, c_double
+        import c_int, c_double, c_ptr
-        integer(c_int), value :: nO, nV, cholesky_mo_num
+        integer(c_int), intent(in), value :: nO, nV, cholesky_mo_num
-        real(c_double), intent(in)  :: t1(nO,nV)
-        real(c_double), intent(in)  :: t2(nO,nO,nV,nV)
-        real(c_double), intent(in)  :: tau(nO,nO,nV,nV)
         real(c_double), intent(in)  :: cc_space_v_oo_chol(cholesky_mo_num,nO,nO)
         real(c_double), intent(in)  :: cc_space_v_ov_chol(cholesky_mo_num,nO,nV)
         real(c_double), intent(in)  :: cc_space_v_vo_chol(cholesky_mo_num,nV,nO)
@@ -21,9 +18,23 @@ module gpu_module
         real(c_double), intent(in)  :: cc_space_v_vooo(nV,nO,nO,nO)
         real(c_double), intent(in)  :: cc_space_v_oovv(nO,nO,nV,nV)
         real(c_double), intent(in)  :: cc_space_v_vvoo(nV,nV,nO,nO)
+        real(c_double), intent(in)  :: cc_space_v_oovo(nO,nO,nV,nO)
+        real(c_double), intent(in)  :: cc_space_v_ovvo(nO,nV,nV,nO)
+        real(c_double), intent(in)  :: cc_space_v_ovoo(nO,nV,nO,nO)
         real(c_double), intent(in)  :: cc_space_f_vo(nV,nO)
+    end function
+
+    subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, gpu_data, t1, t2, tau,&
+          H_vv, g_occ, J1, r2) bind(C)
+        import c_int, c_double, c_ptr
+        integer(c_int), intent(in), value :: nO, nV, cholesky_mo_num
+        type(c_ptr), value    :: gpu_data
+        real(c_double), intent(in)  :: t1(nO,nV)
+        real(c_double), intent(in)  :: t2(nO,nO,nV,nV)
+        real(c_double), intent(in)  :: tau(nO,nO,nV,nV)
         real(c_double), intent(in)  :: H_vv(nV,nV)
         real(c_double), intent(in)  :: g_occ(nO,nO)
+        real(c_double)              :: J1(nO,nV,nV,nO)
         real(c_double), intent(out) :: r2(nO,nO,nV,nV)
     end subroutine