From a7e0832daef329fa4592c93e823682c3750a1a5b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 4 Aug 2023 14:07:49 +0200
Subject: [PATCH] Starting r1 on GPU

---
 devel/ccsd_gpu/ccsd_space_orb_sub.irp.f      |  12 +-
 devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f |  94 +--
 devel/ccsd_gpu/gpu.c                         | 703 ++++++++++++-------
 devel/ccsd_gpu/gpu.h                         |   1 +
 devel/ccsd_gpu/gpu_init.c                    |   9 +-
 devel/ccsd_gpu/gpu_module.f90                |  12 +-
 6 files changed, 502 insertions(+), 329 deletions(-)

diff --git a/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f b/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
index d3c6be9..b47e54f 100644
--- a/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
+++ b/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
@@ -121,7 +121,7 @@ subroutine run_ccsd_space_orb
       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_ovov, cc_space_v_ovoo, &
-      cc_space_f_oo, cc_space_f_vo, cc_space_f_vv)
+      cc_space_f_oo, cc_space_f_ov, cc_space_f_vo, cc_space_f_vv)
 
   if (.not.do_ao_cholesky) then
     print *, 'ao_choleky is required'
@@ -132,11 +132,19 @@ subroutine run_ccsd_space_orb
 
     ! Residue
     call gpu_upload(gpu_data, nO, nV, t1, t2, tau, tau_x);
+    !$OMP PARALLEL SECTIONS
+    !$OMP SECTION
     call compute_H_oo_chol_gpu(gpu_data,nO,nV,0,H_oo)
+
+    !$OMP SECTION
     call compute_H_vo_chol_gpu(gpu_data,nO,nV,1,H_vo)
+
+    !$OMP SECTION
     call compute_H_vv_chol_gpu(gpu_data,nO,nV,2,H_vv)
 
-    call compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+    !$OMP END PARALLEL SECTIONS
+
+    call compute_r1_space_chol(gpu_data, nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
     call compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2)
 
     max_r = max(max_r1,max_r2)
diff --git a/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f b/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
index 45b3ea7..f3e3ea8 100644
--- a/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
+++ b/devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
@@ -80,11 +80,12 @@ end
 
 ! R1
 
-subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
-
+subroutine compute_r1_space_chol(gpu_data, nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
+  use gpu_module
   implicit none
 
   ! in
+  type(c_ptr), intent(in) :: gpu_data
   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)
@@ -95,68 +96,35 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
   ! internal
   integer                       :: u,i,j,beta,a,b
 
-  !$omp parallel &
-  !$omp shared(nO,nV,r1,cc_space_f_ov) &
-  !$omp private(u,beta) &
-  !$omp default(none)
-  !$omp do
-  do beta = 1, nV
-    do u = 1, nO
-      r1(u,beta) = cc_space_f_ov(u,beta)
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
 
-  double precision, allocatable :: X_oo(:,:)
-  allocate(X_oo(nO,nO))
-  call dgemm('N','N', nO, nO, nV, &
-             -2d0, t1    , size(t1,1), &
-                   cc_space_f_vo, size(cc_space_f_vo,1), &
-              0d0, X_oo  , size(X_oo,1))
+  call compute_r1_space_chol_gpu(gpu_data, nO, nV, t1, r1, max_r1)
 
-  call dgemm('T','N', nO, nV, nO, &
-             1d0, X_oo, size(X_oo,2), &
-                  t1  , size(t1,1), &
-             1d0, r1  , size(r1,1))
-  deallocate(X_oo)
-
-  call dgemm('N','N', nO, nV, nV, &
-             1d0, t1  , size(t1,1), &
-                  H_vv, size(H_vv,1), &
-             1d0, r1  , size(r1,1))
-
-  call dgemm('N','N', nO, nV, nO, &
-             -1d0, H_oo, size(H_oo,1), &
-                   t1  , size(t1,1), &
-              1d0, r1, size(r1,1))
-
-  double precision, allocatable :: X_voov(:,:,:,:)
-  allocate(X_voov(nV, nO, nO, nV))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,X_voov,t2,t1) &
-  !$omp private(u,beta,i,a) &
-  !$omp default(none)
-  !$omp do
-  do beta = 1, nV
-    do u = 1, nO
-      do i = 1, nO
-        do a = 1, nV
-          X_voov(a,i,u,beta) = 2d0 * t2(i,u,a,beta) - t2(u,i,a,beta) + t1(u,a) * t1(i,beta)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-  call dgemv('T', nV*nO, nO*nV, &
-             1d0, X_voov, size(X_voov,1) * size(X_voov,2), &
-                  H_vo  , 1, &
-             1d0, r1    , 1)
-
-  deallocate(X_voov)
+!  double precision, allocatable :: X_voov(:,:,:,:)
+!  allocate(X_voov(nV, nO, nO, nV))
+!
+!  !$omp parallel &
+!  !$omp shared(nO,nV,X_voov,t2,t1) &
+!  !$omp private(u,beta,i,a) &
+!  !$omp default(none)
+!  !$omp do
+!  do beta = 1, nV
+!    do u = 1, nO
+!      do i = 1, nO
+!        do a = 1, nV
+!          X_voov(a,i,u,beta) = 2d0 * t2(i,u,a,beta) - t2(u,i,a,beta) + t1(u,a) * t1(i,beta)
+!        enddo
+!      enddo
+!    enddo
+!  enddo
+!  !$omp end do
+!  !$omp end parallel
+!
+!  call dgemv('T', nV*nO, nO*nV, &
+!             1d0, X_voov, nV * nO, &
+!                  H_vo  , 1, &
+!             1d0, r1    , 1)
+!
+!  deallocate(X_voov)
 
   double precision, allocatable :: X_ovov(:,:,:,:)
   allocate(X_ovov(nO, nV, nO, nV))
@@ -263,7 +231,7 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
   call dgemm('T','N', nO, nV, nO*nO*nV, &
              -1d0, W_oovo, size(W_oovo,1) * size(W_oovo,2) * size(W_oovo,3), &
                    tau   , size(tau,1) * size(tau,2) * size(tau,3), &
-              1d0, r1    , size(r1,1))
+              1d0, r1    , nO)
 
   deallocate(W_oovo)
 
diff --git a/devel/ccsd_gpu/gpu.c b/devel/ccsd_gpu/gpu.c
index 173782a..47a47c2 100644
--- a/devel/ccsd_gpu/gpu.c
+++ b/devel/ccsd_gpu/gpu.c
@@ -43,6 +43,286 @@ void gpu_upload(gpu_data* data,
 
 }
 
+
+void compute_h_oo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_oo)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
+    cudaSetDevice(igpu);
+
+    int m,n,k, lda, ldb, ldc;
+    double alpha, beta;
+    double* A;
+    double* B;
+    double* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    double* d_H_oo               = data[igpu].H_oo;
+    double* d_tau_x              = data[igpu].tau_x;
+    double* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
+    double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
+    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+
+    double* d_tau_kau;
+    cudaMalloc((void **)&d_tau_kau, cholesky_mo_num*nV*nO * sizeof(double));
+
+    double* d_tmp_ovv;
+    cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(double));
+
+    double* d_tmp_vov;
+    cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(double));
+
+    for (int i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    alpha = 1.0;
+    beta = 0.0;
+    for (int u=0 ; u<nO ; ++u) {
+      cublasDcopy(handle, nO*nV*nV, &(d_tau_x[u]), nO, d_tmp_ovv, 1);
+      for (int 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;
+            cublasDgeam(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;
+      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+    }
+    for (int i=0 ; i<nV ; ++i) {
+      cudaStreamDestroy(stream[i]);
+    }
+
+    cudaFree(d_tmp_vov);
+    cudaFree(d_tmp_ovv);
+
+    cublasDcopy(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;
+    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    cudaFree(d_tau_kau);
+
+//    double* H_oo = malloc(nO*nO*sizeof(double));
+    cublasGetMatrix(nO, nO, sizeof(double), d_H_oo, nO, H_oo, nO);
+    for (int i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         double* d_H_oo = data[i].H_oo;
+         cudaSetDevice(i);
+         cublasSetMatrix(nO, nO, sizeof(double), H_oo, nO, d_H_oo, nO);
+       }
+    }
+//    free(H_oo);
+
+    cublasDestroy(handle);
+
+}
+
+
+void compute_h_vo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_vo)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
+    cudaSetDevice(igpu);
+
+    int m,n,k, lda, ldb, ldc;
+    double alpha, beta;
+    double* A;
+    double* B;
+    double* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    double* d_t1                 = data[igpu].t1;
+    double* d_H_vo               = data[igpu].H_vo;
+    double* d_tau_x              = data[igpu].tau_x;
+    double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+    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;
+
+    cublasDcopy(handle, nV*nO, d_cc_space_f_vo, 1, d_H_vo, 1);
+
+    double* d_tmp_k;
+    cudaMalloc((void **)&d_tmp_k, cholesky_mo_num * sizeof(double));
+
+    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;
+    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;
+    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;
+    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;
+    cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(double));
+
+    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;
+    cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    double* d_tmp2;
+    cudaMalloc((void **)&d_tmp2, cholesky_mo_num*nO*nO * sizeof(double));
+
+    for (int i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    for (int i=0 ; i<nO ; ++i) {
+      for (int j=0 ; j<nO ; ++j) {
+         cublasSetStream(handle, stream[j]);
+         cublasDcopy(handle, cholesky_mo_num, &(d_tmp [cholesky_mo_num*(i+nO*j)]), 1,
+                                              &(d_tmp2[cholesky_mo_num*(j+nO*i)]), 1);
+      }
+    }
+    for (int 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;
+    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+//    double* H_vo = malloc(nO*nO*sizeof(double));
+    cublasGetMatrix(nV, nO, sizeof(double), d_H_vo, nV, H_vo, nV);
+    for (int i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         double* d_H_vo = data[i].H_vo;
+         cudaSetDevice(i);
+         cublasSetMatrix(nV, nO, sizeof(double), H_vo, nV, d_H_vo, nV);
+       }
+    }
+//    free(H_vo);
+
+    cublasDestroy(handle);
+
+}
+
+
+
+void compute_h_vv_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_vv)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+    igpu = igpu % ngpus;
+
+    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
+    cudaSetDevice(igpu);
+
+    int m,n,k, lda, ldb, ldc;
+    double alpha, beta;
+    double* A;
+    double* B;
+    double* C;
+    cudaStream_t stream[nV];
+
+    cublasHandle_t handle;
+    cublasCreate(&handle);
+
+    double* d_H_vv               = data[igpu].H_vv;
+    double* d_tau_x              = data[igpu].tau_x;
+    double* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
+    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+
+    double* d_tau_kia;
+    cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(double));
+
+    double* d_tmp_oov;
+    cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(double));
+
+    alpha = 1.0;
+    beta = 0.0;
+    for (int i=0 ; i<nV ; ++i) {
+      cudaStreamCreate(&(stream[i]));
+    }
+    for (int a=0 ; a<nV ; ++a) {
+      for (int b=0 ; b<nV ; ++b) {
+         cublasSetStream(handle, stream[b]);
+         cublasDcopy(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;
+      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+    }
+    for (int i=0 ; i<nV ; ++i) {
+      cudaStreamDestroy(stream[i]);
+    }
+    cudaFree(d_tmp_oov);
+
+    cublasDcopy(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;
+    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+    cudaFree(d_tau_kia);
+
+//    double* H_vv = malloc(nO*nO*sizeof(double));
+    cublasGetMatrix(nV, nV, sizeof(double), d_H_vv, nV, H_vv, nV);
+    for (int i=0 ; i<ngpus ; ++i) {
+       if (i != igpu) {
+         double* d_H_vv = data[i].H_vv;
+         cudaSetDevice(i);
+         cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_H_vv, nV);
+       }
+    }
+//    free(H_vv);
+
+    cublasDestroy(handle);
+
+}
+
+
+
+
 void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, double* r2, double* max_r2)
 {
     const int cholesky_mo_num = data->cholesky_mo_num;
@@ -1294,7 +1574,6 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
               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;
@@ -1344,281 +1623,183 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 }
 
 
-void compute_h_oo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_oo)
+
+
+
+
+void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, double* r1, double* max_r1)
 {
+    const int cholesky_mo_num = data->cholesky_mo_num;
+
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
-    igpu = igpu % ngpus;
 
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    cudaSetDevice(igpu);
+    #pragma omp parallel num_threads(ngpus)
+    {
+      int m,n,k, lda, ldb, ldc;
+      double alpha, beta;
+      double* A;
+      double* B;
+      double* C;
+      cudaStream_t stream[nV];
 
-    int m,n,k, lda, ldb, ldc;
-    double alpha, beta;
-    double* A;
-    double* B;
-    double* C;
-    cudaStream_t stream[nV];
+      int igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
 
-    cublasHandle_t handle;
-    cublasCreate(&handle);
+      cublasHandle_t handle;
+      cublasCreate(&handle);
 
-    double* d_H_oo               = data[igpu].H_oo;
-    double* d_tau_x              = data[igpu].tau_x;
-    double* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
-    double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
-    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
+      double* d_r1;
+      lda = nO ;
+      cudaMalloc((void **)&d_r1, lda * nV * sizeof(double));
+      cudaMemset(d_r1, 0, nO*nV*sizeof(double));
+      memset(r1, 0, nO*nV*sizeof(double));
 
-    double* d_tau_kau;
-    cudaMalloc((void **)&d_tau_kau, cholesky_mo_num*nV*nO * sizeof(double));
+//    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_ovov    = data[igpu].cc_space_v_ovov;
+//    double* d_cc_space_v_ovoo    = data[igpu].cc_space_v_ovoo;
+//    double* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
+      double* d_cc_space_f_ov      = data[igpu].cc_space_f_ov;
+      double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
+//    double* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
+//    double* d_tau                = data[igpu].tau;
+      double* d_t1                 = data[igpu].t1;
+      double* d_t2                 = data[igpu].t2;
+      double* d_H_oo               = data[igpu].H_oo;
+      double* d_H_vo               = data[igpu].H_vo;
+      double* d_H_vv               = data[igpu].H_vv;
 
-    double* d_tmp_ovv;
-    cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(double));
+      #pragma omp sections
+      {
 
-    double* d_tmp_vov;
-    cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(double));
+        #pragma omp section
+        {
+          cublasDcopy(handle, nO*nV, d_cc_space_f_ov, 1, d_r1, 1);
+
+          double* d_X_oo;
+          cudaMalloc((void **)&d_X_oo, nO*nO * sizeof(double));
+
+          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;
+          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;
+          m=nO;  n=nV;  k=nO;
+          A=d_X_oo;  lda=nO;
+          B=d_t1;  ldb=nO;
+          C=d_r1;  ldc=nO;
+          cublasDgemm(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;
+          cublasDgemm(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;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+        }
+
+        #pragma omp section
+        {
+          double* d_X_voov;
+          cudaMalloc((void **)&d_X_voov, nV* nO* nO* nV * sizeof(double));
+
+          for (int i=0 ; i<nV ; ++i) {
+            cudaStreamCreate(&(stream[i]));
+          }
+          alpha = -1.0;
+          for (int i=0 ; i<nO ; ++i) {
+            for (int 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;
+              cublasDgeam(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 (int 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;
+            cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
+          }
+          for (int 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;
+          cublasDgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          cudaFree(d_X_voov);
+        }
+
+        #pragma omp section
+        {
+        }
 
-    for (int i=0 ; i<nV ; ++i) {
-      cudaStreamCreate(&(stream[i]));
-    }
-    alpha = 1.0;
-    beta = 0.0;
-    for (int u=0 ; u<nO ; ++u) {
-      cublasDcopy(handle, nO*nV*nV, &(d_tau_x[u]), nO, d_tmp_ovv, 1);
-      for (int 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;
-            cublasDgeam(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;
-      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-    }
-    for (int i=0 ; i<nV ; ++i) {
-      cudaStreamDestroy(stream[i]);
-    }
 
-    cudaFree(d_tmp_vov);
-    cudaFree(d_tmp_ovv);
+      double * r1_tmp = malloc(nO*nV*sizeof(double));
+      lda=nO;
+      cublasGetMatrix(nO, nV, sizeof(double), 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);
 
-    cublasDcopy(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;
-    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+      cudaFree(d_r1);
 
-    cudaFree(d_tau_kau);
-
-//    double* H_oo = malloc(nO*nO*sizeof(double));
-    cublasGetMatrix(nO, nO, sizeof(double), d_H_oo, nO, H_oo, nO);
-    for (int i=0 ; i<ngpus ; ++i) {
-       if (i != igpu) {
-         double* d_H_oo = data[i].H_oo;
-         cudaSetDevice(i);
-         cublasSetMatrix(nO, nO, sizeof(double), H_oo, nO, d_H_oo, nO);
-       }
-    }
-//    free(H_oo);
-
-    cublasDestroy(handle);
-
-}
-
-
-
-
-void compute_h_vv_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_vv)
-{
-    int ngpus = 1;
-    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
-    igpu = igpu % ngpus;
-
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    cudaSetDevice(igpu);
-
-    int m,n,k, lda, ldb, ldc;
-    double alpha, beta;
-    double* A;
-    double* B;
-    double* C;
-    cudaStream_t stream[nV];
-
-    cublasHandle_t handle;
-    cublasCreate(&handle);
-
-    double* d_H_vv               = data[igpu].H_vv;
-    double* d_tau_x              = data[igpu].tau_x;
-    double* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
-    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
-
-    double* d_tau_kia;
-    cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(double));
-
-    double* d_tmp_oov;
-    cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(double));
-
-    alpha = 1.0;
-    beta = 0.0;
-    for (int i=0 ; i<nV ; ++i) {
-      cudaStreamCreate(&(stream[i]));
-    }
-    for (int a=0 ; a<nV ; ++a) {
-      for (int b=0 ; b<nV ; ++b) {
-         cublasSetStream(handle, stream[b]);
-         cublasDcopy(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;
-      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-    }
-    for (int i=0 ; i<nV ; ++i) {
-      cudaStreamDestroy(stream[i]);
-    }
-    cudaFree(d_tmp_oov);
-
-    cublasDcopy(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;
-    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-    cudaFree(d_tau_kia);
-
-//    double* H_vv = malloc(nO*nO*sizeof(double));
-    cublasGetMatrix(nV, nV, sizeof(double), d_H_vv, nV, H_vv, nV);
-    for (int i=0 ; i<ngpus ; ++i) {
-       if (i != igpu) {
-         double* d_H_vv = data[i].H_vv;
-         cudaSetDevice(i);
-         cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_H_vv, nV);
-       }
-    }
-//    free(H_vv);
-
-    cublasDestroy(handle);
-
-}
-
-
-
-void compute_h_vo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_vo)
-{
-    int ngpus = 1;
-    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
-    igpu = igpu % ngpus;
-
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    cudaSetDevice(igpu);
-
-    int m,n,k, lda, ldb, ldc;
-    double alpha, beta;
-    double* A;
-    double* B;
-    double* C;
-    cudaStream_t stream[nV];
-
-    cublasHandle_t handle;
-    cublasCreate(&handle);
-
-    double* d_t1                 = data[igpu].t1;
-    double* d_H_vo               = data[igpu].H_vo;
-    double* d_tau_x              = data[igpu].tau_x;
-    double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
-    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;
-
-    cublasDcopy(handle, nV*nO, d_cc_space_f_vo, 1, d_H_vo, 1);
-
-    double* d_tmp_k;
-    cudaMalloc((void **)&d_tmp_k, cholesky_mo_num * sizeof(double));
-
-    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;
-    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;
-    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;
-    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;
-    cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(double));
-
-    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;
-    cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-    double* d_tmp2;
-    cudaMalloc((void **)&d_tmp2, cholesky_mo_num*nO*nO * sizeof(double));
-
-    for (int i=0 ; i<nV ; ++i) {
-      cudaStreamCreate(&(stream[i]));
-    }
-    for (int i=0 ; i<nO ; ++i) {
-      for (int j=0 ; j<nO ; ++j) {
-         cublasSetStream(handle, stream[j]);
-         cublasDcopy(handle, cholesky_mo_num, &(d_tmp [cholesky_mo_num*(i+nO*j)]), 1,
-                                              &(d_tmp2[cholesky_mo_num*(j+nO*i)]), 1);
-      }
-    }
-    for (int 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;
-    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
-
-//    double* H_vo = malloc(nO*nO*sizeof(double));
-    cublasGetMatrix(nV, nO, sizeof(double), d_H_vo, nV, H_vo, nV);
-    for (int i=0 ; i<ngpus ; ++i) {
-       if (i != igpu) {
-         double* d_H_vo = data[i].H_vo;
-         cudaSetDevice(i);
-         cublasSetMatrix(nV, nO, sizeof(double), H_vo, nV, d_H_vo, nV);
-       }
-    }
-//    free(H_vo);
-
-    cublasDestroy(handle);
+      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;
+  }
 }
 
diff --git a/devel/ccsd_gpu/gpu.h b/devel/ccsd_gpu/gpu.h
index 48c6f0f..4fbf657 100644
--- a/devel/ccsd_gpu/gpu.h
+++ b/devel/ccsd_gpu/gpu.h
@@ -12,6 +12,7 @@ typedef struct {
    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;
    double* tau;
diff --git a/devel/ccsd_gpu/gpu_init.c b/devel/ccsd_gpu/gpu_init.c
index 0625bca..7fd776f 100644
--- a/devel/ccsd_gpu/gpu_init.c
+++ b/devel/ccsd_gpu/gpu_init.c
@@ -14,8 +14,8 @@ gpu_data* gpu_init(
    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_vo,
-   double* cc_space_f_vv)
+   double* cc_space_f_oo, double* cc_space_f_ov,
+   double* cc_space_f_vo, double* cc_space_f_vv)
 {
     int ngpus = 1;
     cudaGetDeviceCount(&ngpus);
@@ -95,6 +95,10 @@ gpu_data* gpu_init(
       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);
 
+      double* d_cc_space_f_ov;
+      cudaMalloc((void**)&d_cc_space_f_ov, nV*nO*sizeof(double));
+      cublasSetMatrix(nO, nV, sizeof(double), cc_space_f_ov, nO, d_cc_space_f_ov, nO);
+
       double* d_cc_space_f_vv;
       cudaMalloc((void**)&d_cc_space_f_vv, nV*nV*sizeof(double));
       cublasSetMatrix(nV, nV, sizeof(double), cc_space_f_vv, nV, d_cc_space_f_vv, nV);
@@ -135,6 +139,7 @@ gpu_data* gpu_init(
       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;
diff --git a/devel/ccsd_gpu/gpu_module.f90 b/devel/ccsd_gpu/gpu_module.f90
index 088caba..f1f04c0 100644
--- a/devel/ccsd_gpu/gpu_module.f90
+++ b/devel/ccsd_gpu/gpu_module.f90
@@ -8,7 +8,7 @@ module gpu_module
       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_ovov, cc_space_v_ovoo, &
-      cc_space_f_oo, cc_space_f_vo, cc_space_f_vv) bind(C)
+      cc_space_f_oo, cc_space_f_ov, cc_space_f_vo, cc_space_f_vv) bind(C)
         import c_int, c_double, c_ptr
         integer(c_int), intent(in), value :: nO, nV, cholesky_mo_num
         real(c_double), intent(in)  :: cc_space_v_oo_chol(cholesky_mo_num,nO,nO)
@@ -24,6 +24,7 @@ module gpu_module
         real(c_double), intent(in)  :: cc_space_v_ovov(nO,nV,nO,nV)
         real(c_double), intent(in)  :: cc_space_v_ovoo(nO,nV,nO,nO)
         real(c_double), intent(in)  :: cc_space_f_oo(nO,nO)
+        real(c_double), intent(in)  :: cc_space_f_ov(nO,nV)
         real(c_double), intent(in)  :: cc_space_f_vo(nV,nO)
         real(c_double), intent(in)  :: cc_space_f_vv(nV,nV)
     end function
@@ -59,6 +60,15 @@ module gpu_module
         real(c_double), intent(out) :: H_vv(nO,nO)
     end subroutine
 
+    subroutine compute_r1_space_chol_gpu(gpu_data, nO, nV, t1, r1, max_r1) bind(C)
+        import c_int, c_double, c_ptr
+        type(c_ptr), value    :: gpu_data
+        integer(c_int), intent(in), value  :: nO, nV
+        real(c_double), intent(in)  :: t1(nO,nV)
+        real(c_double), intent(out) :: r1(nO,nO,nV,nV)
+        real(c_double), intent(out) :: max_r1
+    end subroutine
+
     subroutine compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2) bind(C)
         import c_int, c_double, c_ptr
         type(c_ptr), value    :: gpu_data