More on GPU

devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f

@@ -462,136 +462,44 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
   block_size = 16
   call set_multiple_levels_omp(.False.)
 
-  call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,tau, &
-      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_f_vo, H_vv, r2)
+  double precision, allocatable :: g_occ(:,:)
+  allocate(g_occ(nO,nO))
+  !$omp parallel &
+  !$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) &
+  !$omp private(i,j,a,u) &
+  !$omp default(none)
+  !$omp do
+  do i = 1, nO
+    do u = 1, nO
+      g_occ(u,i) = H_oo(u,i)
+    enddo
+  enddo
+  !$omp end do
 
-!  call compute_g_vir_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,t2, &
-!      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_f_vo, H_vv, r2)
+  !$omp do
+  do i = 1, nO
+    do j = 1, nO
+      do a = 1, nV
+        do u = 1, nO
+          g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
+        enddo
+      enddo
+    enddo
+  enddo
+  !$omp end do
+  !$omp end parallel
+
+  call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,tau, &
+      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_f_vo, H_vv, g_occ, r2)
 
 
 !---
-  double precision, allocatable :: g_occ(:,:)
-  allocate(g_occ(nO,nO))
-  call compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
-
   double precision, allocatable :: X_oovv(:,:,:,:)
-  allocate(X_oovv(nO,nO,nV,nV))
-  call dgemm('N','N',nO,nO*nV*nV,nO, &
-             1d0, g_occ , size(g_occ,1), &
-                  t2    , size(t2,1), &
-             0d0, X_oovv, size(X_oovv,1))
-  deallocate(g_occ)
-
-  !$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_vovv(:,:,:,:)
-
-  allocate(X_vovv(nV,nO,nV,block_size))
-  allocate(Y_oovv(nO,nO,nV,nV))
-
-  do iblock = 1, nV, block_size
-    do gam = iblock, min(nV, iblock+block_size-1)
-      call dgemm('T','N',nV, nO*nV, cholesky_mo_num, 1.d0, &
-        cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, cc_space_v_ov_chol, &
-        cholesky_mo_num, 0.d0, X_vovv(1,1,1,gam-iblock+1), nV)
-
-    enddo
-    call dgemm('N','N',nO,nO*nV*min(block_size, nV-iblock+1),nV, &
-             1d0, t1    , size(t1,1), &
-                  X_vovv, size(X_vovv,1), &
-             0d0, Y_oovv(1,1,1,iblock), size(Y_oovv,1))
-
-  enddo
-  deallocate(X_vovv)
-
-  !$omp parallel &
-  !$omp shared(nO,nV,r2,Y_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) + Y_oovv(v,u,beta,gam) + Y_oovv(u,v,gam,beta)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-  deallocate(Y_oovv)
-
   double precision, allocatable :: X_ovvo(:,:,:,:)
   double precision, allocatable :: tcc(:,:,:), tcc2(:,:,:)
-  allocate(tcc2(cholesky_mo_num,nV,nO), X_ovvo(nO,nV,nV,nO))
-  allocate(tcc(cholesky_mo_num,nO,nV))
-
-  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, tcc2, cholesky_mo_num*nV)
-
-  call dgemm('N','N', cholesky_mo_num*nO, nV, nO, 1.d0, &
-     cc_space_v_oo_chol, cholesky_mo_num*nO, t1, nO, &
-     0.d0, tcc, cholesky_mo_num*nO)
-
-  call dgemm('T','N', nO*nV, nV*nO, cholesky_mo_num, 1.d0, &
-              tcc, cholesky_mo_num, tcc2, cholesky_mo_num, 0.d0, &
-              X_ovvo, nO*nV)
-
-  deallocate(tcc, tcc2)
-
-  !$omp parallel &
-  !$omp shared(nO,nV,r2,X_ovvo) &
-  !$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_ovvo(u,beta,gam,v)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp do
-  do beta = 1, nV
-    do gam = 1, nV
-      do v = 1, nO
-        do u = 1, nO
-          r2(v,u,gam,beta) = r2(v,u,gam,beta) - X_ovvo(u,beta,gam,v)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-  deallocate(X_ovvo)
-  !-----
 
   allocate(X_oovv(nO,nO,nV,nV))
 
@@ -892,106 +800,7 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
 end
 
 
-! g_occ
 
-subroutine compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
-
-  implicit none
-
-  integer, intent(in)           :: nO,nV
-  double precision, intent(in)  :: t1(nO, nV), H_oo(nO, nO)
-  double precision, intent(in)  :: t2(nO, nO, nV, nV)
-  double precision, intent(out) :: g_occ(nO, nO)
-
-  integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
-
-  call dgemm('N','N',nO,nO,nV, &
-             1d0, t1, size(t1,1), &
-                  cc_space_f_vo, size(cc_space_f_vo,1), &
-             0d0, g_occ, size(g_occ,1))
-
-  !$omp parallel &
-  !$omp shared(nO,nV,g_occ,H_oo, cc_space_v_ovoo,t1) &
-  !$omp private(i,j,a,u) &
-  !$omp default(none)
-  !$omp do
-  do i = 1, nO
-    do u = 1, nO
-      g_occ(u,i) = g_occ(u,i) + H_oo(u,i)
-    enddo
-  enddo
-  !$omp end do
-
-  !$omp do
-  do i = 1, nO
-    do j = 1, nO
-      do a = 1, nV
-        do u = 1, nO
-          g_occ(u,i) = g_occ(u,i) + (2d0 * cc_space_v_ovoo(u,a,i,j) - cc_space_v_ovoo(u,a,j,i)) * t1(j,a)
-        enddo
-      enddo
-    enddo
-  enddo
-  !$omp end do
-  !$omp end parallel
-
-end
-
-! g_vir
-
-subroutine compute_g_vir_chol(nO,nV,t1,t2,H_vv,g_vir)
-  use gpu_module
-  implicit none
-
-  integer, intent(in)           :: nO,nV
-  double precision, intent(in)  :: t1(nO, nV), H_vv(nV, nV)
-  double precision, intent(in)  :: t2(nO, nO, nV, nV)
-  double precision, intent(out) :: g_vir(nV, nV)
-
-  integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
-
-!  do beta = 1, nV
-!    do a = 1, nV
-!      g_vir(a,beta) = H_vv(a,beta)
-!    enddo
-!  enddo
-
-!  call dgemm('N','N',nV,nV,nO, &
-!             -1d0, cc_space_f_vo , size(cc_space_f_vo,1), &
-!                   t1   , size(t1,1), &
-!              1d0, g_vir, size(g_vir,1))
-
-  double precision, allocatable :: tmp_k(:), tmp_vo(:,:,:), tmp_vo2(:,:,:)
-!  allocate(tmp_k(cholesky_mo_num))
-!  call dgemm('N','N', cholesky_mo_num, 1, nO*nV, 1.d0, &
-!    cc_space_v_ov_chol, cholesky_mo_num, t1, nO*nV, 0.d0, tmp_k, cholesky_mo_num)
-!
-!  call dgemm('T','N', nV*nV, 1, cholesky_mo_num, 2.d0, &
-!    cc_space_v_vv_chol, cholesky_mo_num, tmp_k, cholesky_mo_num, 1.d0, &
-!    g_vir, nV*nV)
-!  deallocate(tmp_k)
-
-!  allocate(tmp_vo(cholesky_mo_num,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_vo, cholesky_mo_num*nV)
-!
-!  allocate(tmp_vo2(cholesky_mo_num,nO,nV))
-!  do beta=1,nV
-!    do i=1,nO
-!      do k=1,cholesky_mo_num
-!        tmp_vo2(k,i,beta) = -tmp_vo(k,beta,i)
-!      enddo
-!    enddo
-!  enddo
-!  deallocate(tmp_vo)
-!
-!  call dgemm('T','N', nV, nV, nO*cholesky_mo_num, 1.d0, &
-!     cc_space_v_ov_chol, cholesky_mo_num*nO, &
-!     tmp_vo2, cholesky_mo_num*nO, 1.d0, g_vir, nV)
-
-end
-
-! J1
 
 subroutine compute_J1_chol(nO,nV,t1,t2,v_ovvo,v_ovoo,v_vvoo,J1)
   implicit none

devel/ccsd_gpu/gpu.c

@@ -63,6 +63,7 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                 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,
@@ -72,6 +73,7 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
                 double* cc_space_v_vvoo,
                 double* cc_space_f_vo,
                 double* H_vv,
+                double* g_occ,
                 double* r2)
 {
 
@@ -87,7 +89,7 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
 
       int ithread = omp_get_thread_num();
       int igpu = ithread ;
-
+//igpu=1;
       cudaSetDevice(igpu);
       cublasHandle_t handle;
 
@@ -95,9 +97,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
 
       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_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;
@@ -115,6 +118,10 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       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);
@@ -131,6 +138,10 @@ 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
       {
 
@@ -208,10 +219,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            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;
@@ -219,7 +226,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            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));
@@ -302,6 +308,158 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
            cudaFree(d_Y_oovv);
         }
 
+        // g_occ
+        #pragma omp section
+        {
+          double* d_g_occ;
+          lda = nO;
+          cudaMalloc((void **)&d_g_occ, nO*nO * sizeof(double));
+          cublasSetMatrix(lda, nO, sizeof(double), g_occ, lda, d_g_occ, lda);
+
+          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;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          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;  n=nO*nV*nV;  k=nO;
+          A=d_g_occ;  lda=nO;
+          B=d_t2;  ldb=nO;
+          C=d_X_oovv;  ldc=nO;
+          cublasDgemm(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;
+           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_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);
+             }
+           }
+           cudaFree(d_X_oovv);
+        }
+
+        #pragma omp section
+        {
+           double* d_X_vovv;
+           lda = nV*nO;
+           cudaMalloc((void **)&d_X_vovv, nV*nO*nV*BLOCK_SIZE * sizeof(double));
+
+           double* d_Y_oovv;
+           lda = nO*nO;
+           cudaMalloc((void **)&d_Y_oovv, nO*nO*nV*nV * sizeof(double));
+
+           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) {
+               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;
+               cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+             }
+             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;
+             cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+           }
+           cudaFree(d_X_vovv);
+           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*(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);
+                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);
+             }
+           }
+           cudaFree(d_Y_oovv);
+        }
+
+        #pragma omp section
+        {
+          double* d_tcc2;
+          lda = cholesky_mo_num;
+          cudaMalloc((void **)&d_tcc2, 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_tcc2;  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_tcc;
+          lda = cholesky_mo_num;
+          cudaMalloc((void **)&d_tcc, cholesky_mo_num*nO*nV * sizeof(double));
+
+          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;
+          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
+
+          double* d_X_ovvo;
+          lda = nO*nV;
+          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+
+          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;
+          cublasDgemm(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(int gam = 0; gam < nV; gam++){
+            for(int bet = 0; bet < nV; bet++){
+               alpha = 1.0;
+               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);
+               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);
+            }
+          }
+
+          cudaFree(d_X_ovvo);
+        }
+
       }
 
 
@@ -323,7 +481,7 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cudaMalloc((void**)&d_tmpB1, nV*BLOCK_SIZE*nV*sizeof(double));
 
       #pragma omp for
-      for (size_t gam=0 ; gam<nV ; ++gam)
+      for (int gam=0 ; gam<nV ; ++gam)
       {
          double* d_tmp_cc_ = &(d_tmp_cc[gam*nV*cholesky_mo_num]);
          double* d_cc_space_v_vv_chol_ = &(d_cc_space_v_vv_chol[gam*nV*cholesky_mo_num]);
@@ -335,9 +493,9 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
          C = d_tmp_cc2 ; ldc = cholesky_mo_num;
          cublasDgeam(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)
+         for (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE)
          {
-              const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+              const int mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
 
               alpha=-1.0; beta=0.0;
               m=nV*mbs;  n=nV;  k=cholesky_mo_num;
@@ -389,7 +547,7 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, lda, r2_tmp, lda);
       #pragma omp critical
       {
-         for (size_t i=0 ; i<nO*nO*nV*nV ; ++i) {
+         for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
            r2[i] += r2_tmp[i];
          }
       }
@@ -399,14 +557,17 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
       cublasDestroy(handle);
   }
 
-  for (size_t i=0 ; i<nO*nO*nV*nV ; ++i)
+  for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i)
   {
-           r2[i] += cc_space_v_oovv[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,

devel/ccsd_gpu/gpu_module.f90

@@ -5,14 +5,15 @@ module gpu_module
 
   interface
     subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
-          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_f_vo, H_vv, r2) bind(C)
+          cc_space_f_vo, H_vv, g_occ, r2) bind(C)
         import c_int, c_double
         integer(c_int), 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)
         real(c_double), intent(in)  :: cc_space_v_vv_chol(cholesky_mo_num,nV,nV)
@@ -22,29 +23,10 @@ module gpu_module
         real(c_double), intent(in)  :: cc_space_v_vvoo(nV,nV,nO,nO)
         real(c_double), intent(in)  :: cc_space_f_vo(nV,nO)
         real(c_double), intent(in)  :: H_vv(nV,nV)
+        real(c_double), intent(in)  :: g_occ(nO,nO)
         real(c_double), intent(out) :: r2(nO,nO,nV,nV)
     end subroutine
 
-    subroutine compute_g_vir_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
-          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_f_vo, H_vv, r2) bind(C)
-        import c_int, c_double
-        integer(c_int), 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_ov_chol(cholesky_mo_num,nO,nV)
-        real(c_double), intent(in)  :: cc_space_v_vo_chol(cholesky_mo_num,nV,nO)
-        real(c_double), intent(in)  :: cc_space_v_vv_chol(cholesky_mo_num,nV,nV)
-        real(c_double), intent(in)  :: cc_space_v_oooo(nO,nO,nO,nO)
-        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_f_vo(nV,nO)
-        real(c_double), intent(in)  :: H_vv(nV,nV)
-        real(c_double), intent(out) :: r2(nO,nO,nV,nV)
-    end subroutine
 
 
     subroutine gpu_dgemm(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) bind(C)