From 0d53d4fcaa2457bb32fe53bfb24b1dbc0a957657 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 21 Aug 2023 12:34:55 +0200
Subject: [PATCH] Fix malloc on GPU

---
 devel/ccsd_gpu/ccsd_space_orb_sub.irp.f |  66 ++-
 devel/ccsd_gpu/gpu.c                    | 583 ++++++++++++++----------
 devel/ccsd_gpu/gpu.h                    |  33 +-
 devel/ccsd_gpu/gpu_dgemm.c              |  42 +-
 devel/ccsd_gpu/gpu_init.c               | 119 +++--
 devel/ccsd_gpu/gpu_module.f90           |  81 +++-
 6 files changed, 617 insertions(+), 307 deletions(-)

diff --git a/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f b/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
index 33ab63b..069244d 100644
--- a/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
+++ b/devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
@@ -94,12 +94,21 @@ subroutine run_ccsd_space_orb
 
   ! Init
   type(c_ptr) :: gpu_data
+  logical :: do_sp = .False.
 
-  gpu_data = gpu_init(nO, nV, cholesky_mo_num, &
+  if (do_sp) then
+    gpu_data = gpu_init_sp(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_voov, 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_ov, cc_space_f_vo, cc_space_f_vv)
+  else
+    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_voov, 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_ov, cc_space_f_vo, cc_space_f_vv)
+  endif
 
   if (.not.do_ao_cholesky) then
     print *, 'ao_choleky is required'
@@ -109,12 +118,20 @@ subroutine run_ccsd_space_orb
   call guess_t1(nO,nV,cc_space_f_o,cc_space_f_v,cc_space_f_ov,t1)
   call guess_t2(nO,nV,cc_space_f_o,cc_space_f_v,cc_space_v_oovv,t2)
 
-  call gpu_upload(gpu_data, nO, nV, t1, t2);
+  if (do_sp) then
+    call gpu_upload_sp(gpu_data, nO, nV, t1, t2);
+  else
+    call gpu_upload(gpu_data, nO, nV, t1, t2);
+  endif
 
   !print*,'hf_energy', hf_energy
   call det_energy(det,uncorr_energy)
   print*,'Det energy', uncorr_energy
-  energy = ccsd_energy_space_gpu(gpu_data)
+  if (do_sp) then
+    energy = ccsd_energy_space_gpu_sp(gpu_data)
+  else
+    energy = ccsd_energy_space_gpu(gpu_data)
+  endif
   print*,'Guess energy', uncorr_energy+energy, energy
 
   nb_iter = 0
@@ -133,18 +150,39 @@ subroutine run_ccsd_space_orb
     ! Residue
     !$OMP PARALLEL SECTIONS
     !$OMP SECTION
-    call compute_H_oo_chol_gpu(gpu_data,0)
+    if (do_sp) then
+      call compute_H_oo_chol_gpu_sp(gpu_data,0)
+    else
+      call compute_H_oo_chol_gpu(gpu_data,0)
+    endif
 
     !$OMP SECTION
-    call compute_H_vo_chol_gpu(gpu_data,1)
+    if (do_sp) then
+      call compute_H_vo_chol_gpu_sp(gpu_data,1)
+    else
+      call compute_H_vo_chol_gpu(gpu_data,1)
+    endif
 
     !$OMP SECTION
-    call compute_H_vv_chol_gpu(gpu_data,2)
+    if (do_sp) then
+      call compute_H_vv_chol_gpu_sp(gpu_data,2)
+    else
+      call compute_H_vv_chol_gpu(gpu_data,2)
+    endif
 
     !$OMP END PARALLEL SECTIONS
 
-    call compute_r1_space_chol_gpu(gpu_data, nO, nV, t1, r1, max_r1)
-    call compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2)
+    if (do_sp) then
+      call compute_r1_space_chol_gpu_sp(gpu_data, nO, nV, t1, r1, max_r1)
+    else
+      call compute_r1_space_chol_gpu(gpu_data, nO, nV, t1, r1, max_r1)
+    endif
+
+    if (do_sp) then
+      call compute_r2_space_chol_gpu_sp(gpu_data, nO, nV, t1, r2, max_r2)
+    else
+      call compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2)
+    endif
 
     max_r = max(max_r1,max_r2)
 
@@ -162,10 +200,18 @@ subroutine run_ccsd_space_orb
       print*,'Unkown cc_method_method: '//cc_update_method
     endif
 
-    call gpu_upload(gpu_data, nO, nV, t1, t2);
+    if (do_sp) then
+      call gpu_upload_sp(gpu_data, nO, nV, t1, t2);
+    else
+      call gpu_upload(gpu_data, nO, nV, t1, t2);
+    endif
 
     ! Energy
-    energy = ccsd_energy_space_gpu(gpu_data)
+    if (do_sp) then
+      energy = ccsd_energy_space_gpu_sp(gpu_data)
+    else
+      energy = ccsd_energy_space_gpu(gpu_data)
+        endif
     write(*,'(A3,I6,A3,F18.12,A3,F16.12,A3,ES10.2,A3,ES10.2,A2)') ' | ',nb_iter,' | ', uncorr_energy+energy,' | ', energy,' | ', max_r1,' | ', max_r2,' |'
 
     nb_iter = nb_iter + 1
diff --git a/devel/ccsd_gpu/gpu.c b/devel/ccsd_gpu/gpu.c
index ab2c8c6..fc118e9 100644
--- a/devel/ccsd_gpu/gpu.c
+++ b/devel/ccsd_gpu/gpu.c
@@ -5,34 +5,41 @@
 #include <cublas_v2.h>
 #include <cuda_runtime.h>
 #include "gpu.h"
+#include "assert.h"
 
 void gpu_upload(gpu_data* data,
                 int nO, int nV,
                 double* t1,
                 double* t2)
 {
-    int lda;
+    size_t lda;
 
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
 
     double * tau = malloc(nO*nO*nV*nV * sizeof(double));
+    assert (tau != NULL);
+
     double * tau_x = malloc(nO*nO*nV*nV * sizeof(double));
+    assert (tau_x != NULL);
 
     #pragma omp parallel num_threads(ngpus)
     {
-        int igpu = omp_get_thread_num();
+        cudaError_t cudaStat = cudaSuccess;
+        size_t igpu = omp_get_thread_num();
         cudaSetDevice(igpu);
 
         double* d_t1 = data[igpu].t1;
         lda = nO;
-        cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
+        cudaStat = cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
+        assert (cudaStat == cudaSuccess);
 
         double* d_t2 = data[igpu].t2;
         lda = nO*nO;
-        cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
+        cudaStat = cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
+        assert (cudaStat == cudaSuccess);
 
-        int lda, ldb, ldc;
+        size_t lda, ldb, ldc;
         double alpha, beta;
         double* A;
         double* B;
@@ -46,16 +53,17 @@ void gpu_upload(gpu_data* data,
 
         double* d_tau_x = data[igpu].tau_x;
         lda = nO * nO;
-        cublasSetMatrix(nO*nO, nV*nV, sizeof(double), tau_x, lda, d_tau_x, lda);
+        cudaStat = cublasSetMatrix(nO*nO, nV*nV, sizeof(double), tau_x, lda, d_tau_x, lda);
+        assert (cudaStat == cudaSuccess);
 
         if (igpu == 0) {
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               beta = t1[j+b*nO];
               A = &(d_t2[nO*(j + nO*nV*b)]); lda = nO*nO;
@@ -67,8 +75,8 @@ void gpu_upload(gpu_data* data,
           cudaDeviceSynchronize();
           alpha = 2.0;
           beta = -1.0;
-          for (int b=0 ; b<nV ; ++b) {
-            for (int a=0 ; a<nV ; ++a) {
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t a=0 ; a<nV ; ++a) {
               cublasSetStream(handle, stream[a]);
               A = &(d_tau[nO*nO*(a + nV*b)]); lda = nO;
               B = &(d_tau[nO*nO*(b + nV*a)]); ldb = nO;
@@ -76,7 +84,7 @@ void gpu_upload(gpu_data* data,
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -103,16 +111,17 @@ void gpu_upload(gpu_data* data,
 
 void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
 {
+    cudaError_t cudaStat;
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
     igpu = igpu % ngpus;
 
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    const int nO = data[igpu].nO;
-    const int nV = data[igpu].nV;
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
     cudaSetDevice(igpu);
 
-    int m,n,k, lda, ldb, ldc;
+    size_t m,n,k, lda, ldb, ldc;
     double alpha, beta;
     double* A;
     double* B;
@@ -129,22 +138,25 @@ void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
     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));
+    cudaStat = cudaMalloc((void **)&d_tau_kau, cholesky_mo_num*nV*nO * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     double* d_tmp_ovv;
-    cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(double));
+    cudaStat = cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     double* d_tmp_vov;
-    cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(double));
+    cudaStat = cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
-    for (int i=0 ; i<nV ; ++i) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamCreate(&(stream[i]));
     }
     alpha = 1.0;
     beta = 0.0;
-    for (int u=0 ; u<nO ; ++u) {
+    for (size_t 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) {
+      for (size_t b=0 ; b<nV ; ++b) {
             cublasSetStream(handle, stream[b]);
             A = &(d_tmp_ovv[nO*nV*b]); lda = nO;
             B = &(d_tmp_ovv[nO*nV*b]); ldb = nO;
@@ -161,7 +173,7 @@ void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
       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) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamDestroy(stream[i]);
     }
 
@@ -180,8 +192,9 @@ void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
     cudaFree(d_tau_kau);
 
     double* H_oo = malloc(nO*nO*sizeof(double));
+    assert (H_oo != NULL);
     cublasGetMatrix(nO, nO, sizeof(double), d_H_oo, nO, H_oo, nO);
-    for (int i=0 ; i<ngpus ; ++i) {
+    for (size_t i=0 ; i<ngpus ; ++i) {
        if (i != igpu) {
          double* d_H_oo = data[i].H_oo;
          cudaSetDevice(i);
@@ -197,16 +210,17 @@ void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
 
 void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
 {
+    cudaError_t cudaStat;
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
     igpu = igpu % ngpus;
 
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    const int nO = data[igpu].nO;
-    const int nV = data[igpu].nV;
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
     cudaSetDevice(igpu);
 
-    int m,n,k, lda, ldb, ldc;
+    size_t m,n,k, lda, ldb, ldc;
     double alpha, beta;
     double* A;
     double* B;
@@ -225,7 +239,8 @@ void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
     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));
+    cudaStat = cudaMalloc((void **)&d_tmp_k, cholesky_mo_num * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     alpha = 2.0;
     beta  = 0.0;
@@ -246,7 +261,8 @@ void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
     cudaFree(d_tmp_k);
 
     double* d_tmp;
-    cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(double));
+    cudaStat = cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     alpha = 1.0;
     beta  = 0.0;
@@ -257,19 +273,20 @@ void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
     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));
+    cudaStat = cudaMalloc((void **)&d_tmp2, cholesky_mo_num*nO*nO * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
-    for (int i=0 ; i<nV ; ++i) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamCreate(&(stream[i]));
     }
-    for (int i=0 ; i<nO ; ++i) {
-      for (int j=0 ; j<nO ; ++j) {
+    for (size_t i=0 ; i<nO ; ++i) {
+      for (size_t 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) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamDestroy(stream[i]);
     }
     cublasSetStream(handle, NULL);
@@ -284,7 +301,7 @@ void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
 
     double* H_vo = malloc(nV*nO*sizeof(double));
     cublasGetMatrix(nV, nO, sizeof(double), d_H_vo, nV, H_vo, nV);
-    for (int i=0 ; i<ngpus ; ++i) {
+    for (size_t i=0 ; i<ngpus ; ++i) {
        if (i != igpu) {
          double* d_H_vo = data[i].H_vo;
          cudaSetDevice(i);
@@ -301,16 +318,17 @@ void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
 
 void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
 {
+    cudaError_t cudaStat;
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
     igpu = igpu % ngpus;
 
-    const int cholesky_mo_num = data[igpu].cholesky_mo_num;
-    const int nO = data[igpu].nO;
-    const int nV = data[igpu].nV;
+    const size_t cholesky_mo_num = data[igpu].cholesky_mo_num;
+    const size_t nO = data[igpu].nO;
+    const size_t nV = data[igpu].nV;
     cudaSetDevice(igpu);
 
-    int m,n,k, lda, ldb, ldc;
+    size_t m,n,k, lda, ldb, ldc;
     double alpha, beta;
     double* A;
     double* B;
@@ -326,18 +344,20 @@ void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
     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));
+    cudaStat = cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     double* d_tmp_oov;
-    cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(double));
+    cudaStat = cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(double));
+    assert(cudaStat == cudaSuccess);
 
     alpha = 1.0;
     beta = 0.0;
-    for (int i=0 ; i<nV ; ++i) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamCreate(&(stream[i]));
     }
-    for (int a=0 ; a<nV ; ++a) {
-      for (int b=0 ; b<nV ; ++b) {
+    for (size_t a=0 ; a<nV ; ++a) {
+      for (size_t 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);
       }
@@ -351,7 +371,7 @@ void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
       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) {
+    for (size_t i=0 ; i<nV ; ++i) {
       cudaStreamDestroy(stream[i]);
     }
     cudaFree(d_tmp_oov);
@@ -369,7 +389,7 @@ void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
 
     double* H_vv = malloc(nV*nV*sizeof(double));
     cublasGetMatrix(nV, nV, sizeof(double), d_H_vv, nV, H_vv, nV);
-    for (int i=0 ; i<ngpus ; ++i) {
+    for (size_t i=0 ; i<ngpus ; ++i) {
        if (i != igpu) {
          double* d_H_vv = data[i].H_vv;
          cudaSetDevice(i);
@@ -387,7 +407,7 @@ void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
 
 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;
+    const size_t cholesky_mo_num = data->cholesky_mo_num;
 
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
@@ -397,14 +417,15 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
     #pragma omp parallel num_threads(ngpus)
     {
-      int m,n,k, lda, ldb, ldc;
+      cudaError_t cudaStat;
+      size_t 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();
+      size_t igpu = omp_get_thread_num();
       cudaSetDevice(igpu);
 
       cublasHandle_t handle;
@@ -412,7 +433,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
       double* d_r2;
       lda = nO * nO;
-      cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cudaMemset(d_r2, 0, nO*nO*nV*nV*sizeof(double));
       memset(r2, 0, nO*nO*nV*nV*sizeof(double));
 
@@ -436,7 +458,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
       double* d_H_vv               = data[igpu].H_vv;
 
       double* d_K1;
-      cudaMalloc((void **)&d_K1, nO*nV*nO*nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_K1, nO*nV*nO*nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       #pragma omp sections
       {
@@ -444,7 +467,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_J1;
-          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta =  0.0;
@@ -455,14 +479,15 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
 
           double* d_X_ovoo;
-          cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
           alpha = 0.0;
           beta =  1.0;
-          for (int i=0 ; i<nO ; ++i) {
+          for (size_t i=0 ; i<nO ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int j=0 ; j<nO ; ++j) {
-            for (int i=0 ; i<nO ; ++i) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_X_ovoo[nO*nV*(i+nO*j)]); lda = nO;
               B = &(d_cc_space_v_ovoo[nO*nV*(j+nO*i)]); ldb = nO;
@@ -470,14 +495,15 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nO ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
 
           double* d_Y_ovov;
-          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -490,11 +516,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
           alpha = 1.0;
           beta =  -1.0;
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int j=0 ; j<nO ; ++j) {
-            for (int i=0 ; i<nV ; ++i) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_J1[nO*nV*(i+nV*j)]); lda = nO;
               B = &(d_Y_ovov[nO*nV*(j+nO*i)]); ldb = nO;
@@ -502,13 +528,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_tmp_cc;
-          cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -519,7 +546,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_J1_tmp;
-          cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -533,31 +561,32 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
           alpha = 1.0;
           beta =  1.0;
-          for (int i=0 ; i<nO ; ++i) {
+          for (size_t i=0 ; i<nO ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int i=0 ; i<nO ; ++i) {
+          for (size_t i=0 ; i<nO ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_J1[nO*nV*nV*i]); lda = nO*nV;
               B = &(d_J1_tmp[nV*nO*nV*i]); ldb = nV;
               C = &(d_J1[nO*nV*nV*i]); ldc = nO*nV;
               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) {
+          for (size_t i=0 ; i<nO ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_J1_tmp);
 
           double* d_X_voov;
-          cudaMalloc((void **)&d_X_voov, nV*nO*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_voov, nV*nO*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t 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) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               beta = t1[j+b*nO];
               A = &(d_t2[nO*(j + nO*nV*b)]); lda = nO*nO;
@@ -568,8 +597,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           }
           alpha = 1.0;
           beta  = 0.0;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); lda = nV*nV;
               B = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); ldb = nV*nV;
@@ -577,13 +606,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Z_ovvo;
-          cudaMalloc((void **)&d_Z_ovvo, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Z_ovvo, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = -1.0;
           beta  = 0.0;
@@ -597,11 +627,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
           alpha = 1.0;
           beta =  1.0;
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int i=0 ; i<nO ; ++i) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
               B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb=nO*nV;
@@ -611,12 +641,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           }
 
           double* d_Y_vovo;
-          cudaMalloc((void **)&d_Y_vovo, nV*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_vovo, nV*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta = -0.5;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int i=0 ; i<nO ; ++i) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
               B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
@@ -626,12 +657,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           }
 
           double* d_X_ovvo;
-          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta =  0.0;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
               B = &(d_t2[nO*(j+nO*nV*b)]); ldb = nO*nO;
@@ -639,7 +671,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -654,13 +686,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_X_ovvo);
           cudaFree(d_Y_vovo);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta =  1.0;
-          for (int i=0 ; i<nO ; ++i) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
               B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb = nO*nV;
@@ -668,7 +700,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -700,18 +732,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_X;
-          cudaMalloc((void **)&d_X, nV*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X, nV*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           double* d_Y;
-          cudaMalloc((void **)&d_Y, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha =-1.0;
           beta =  0.0;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int i=0 ; i<nO ; ++i) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t i=0 ; i<nO ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
               B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
@@ -721,8 +755,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           }
 
           alpha = 0.5;
-          for (int j=0 ; j<nO ; ++j) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t j=0 ; j<nO ; ++j) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               beta = t1[j+b*nO];
               A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
@@ -731,13 +765,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Z;
-          cudaMalloc((void **)&d_Z, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Z, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -752,7 +787,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
 
           double* d_t1v;
-          cudaMalloc((void **)&d_t1v, cholesky_mo_num*nO*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_t1v, cholesky_mo_num*nO*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -763,7 +799,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_K1tmp;
-          cudaMalloc((void **)&d_K1tmp, nO*nO*nV*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_K1tmp, nO*nO*nV*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -775,13 +812,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_t1v);
 
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta =  1.0;
-          for (int i=0 ; i<nO ; ++i) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
               B = &(d_K1tmp[nO*(i+nO*nV*b)]); ldb = nO*nO;
@@ -790,8 +827,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
             }
           }
 
-          for (int i=0 ; i<nO ; ++i) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
               B = &(d_Z[nO*(b+nV*nV*i)]); ldb = nO*nV;
@@ -800,7 +837,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
             }
           }
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -815,7 +852,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
            double* d_Y_oooo;
-           cudaMalloc((void**)&d_Y_oooo, nO*nO*nO*nO*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_Y_oooo, nO*nO*nO*nO*sizeof(double));
+           assert (cudaStat == cudaSuccess);
 
            alpha = 1.0;
            beta =  0.0;
@@ -826,7 +864,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
            double* d_A1;
-           cudaMalloc((void**)&d_A1, nO*nO*nO*nO*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_A1, nO*nO*nO*nO*sizeof(double));
+           assert (cudaStat == cudaSuccess);
 
            alpha = 1.0;
            beta =  1.0;
@@ -834,11 +873,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            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) {
+           for (size_t i=0 ; i<nO ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
-           for (int j=0 ; j<nO ; ++j) {
-             for (int i=0 ; i<nO ; ++i) {
+           for (size_t j=0 ; j<nO ; ++j) {
+             for (size_t i=0 ; i<nO ; ++i) {
                 cublasSetStream(handle, stream[i]);
                 alpha = 1.0;
                 beta =  1.0;
@@ -849,7 +888,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
              }
            }
-           for (int i=0 ; i<nO ; ++i) {
+           for (size_t i=0 ; i<nO ; ++i) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);
@@ -877,7 +916,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
            double* d_g_vir;
-           cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
+           assert (cudaStat == cudaSuccess);
            cublasDcopy(handle, nV*nV, d_H_vv, 1, d_g_vir, 1);
 
            alpha = -1.0;
@@ -889,7 +929,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
            double* d_tmp_k;
-           cudaMalloc((void**)&d_tmp_k, cholesky_mo_num*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_tmp_k, cholesky_mo_num*sizeof(double));
+           assert (cudaStat == cudaSuccess);
            alpha = 1.0;
            beta =  0.0;
            m=cholesky_mo_num ; n=1; k=nO*nV;
@@ -908,7 +949,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            cudaFree(d_tmp_k);
 
            double* d_tmp_vo;
-           cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(double));
+           assert (cudaStat == cudaSuccess);
            alpha = 1.0;
            beta =  0.0;
            m=cholesky_mo_num*nV ; n=nO; k=nV;
@@ -918,11 +960,12 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            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) {
+           cudaStat = cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(double));
+           assert (cudaStat == cudaSuccess);
+           for (size_t i=0 ; i<nO ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
-           for (int i=0 ; i<nO ; ++i) {
+           for (size_t i=0 ; i<nO ; ++i) {
              cublasSetStream(handle, stream[i]);
              alpha = -1.0;
              beta =  0.0;
@@ -931,7 +974,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
              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) {
+           for (size_t i=0 ; i<nO ; ++i) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);
@@ -947,7 +990,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            cudaFree(d_tmp_vo2);
 
            double* d_Y_oovv;
-           cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(double));
+           cudaStat = cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(double));
+           assert (cudaStat == cudaSuccess);
            alpha = 1.0;
            beta =  0.0;
            m=nO*nO*nV ; n=nV; k=nV;
@@ -964,11 +1008,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            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) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
-           for (int j=0 ; j<nV ; ++j) {
-             for (int i=0 ; i<nV ; ++i) {
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
                 cublasSetStream(handle, stream[i]);
                 alpha = 1.0;
                 beta =  1.0;
@@ -979,7 +1023,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
              }
            }
-           for (int i=0 ; i<nV ; ++i) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);
@@ -991,11 +1035,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         {
           double* d_g_occ;
           lda = nO;
-          cudaMalloc((void **)&d_g_occ, nO*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_g_occ, nO*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
           cublasDcopy(handle, nO*nO, d_H_oo, 1, d_g_occ, 1);
 
           double* d_X;
-          cudaMalloc((void **)&d_X, cholesky_mo_num*sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X, cholesky_mo_num*sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 2.0;
           beta  = 0.0;
@@ -1031,7 +1077,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_X_oovv;
-          cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1049,11 +1096,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            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) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
-           for (int j=0 ; j<nV ; ++j) {
-             for (int i=0 ; i<nV ; ++i) {
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
                 cublasSetStream(handle, stream[i]);
                 alpha = 1.0;
                 beta = -1.0;
@@ -1063,7 +1110,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                 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) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);
@@ -1083,17 +1130,19 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
            double* d_X_vovv;
-           cudaMalloc((void **)&d_X_vovv, nV*nO*nV*BLOCK_SIZE * sizeof(double));
+           cudaStat = cudaMalloc((void **)&d_X_vovv, nV*nO*nV*BLOCK_SIZE * sizeof(double));
+           assert (cudaStat == cudaSuccess);
 
            double* d_Y_oovv;
-           cudaMalloc((void **)&d_Y_oovv, nO*nO*nV*nV * sizeof(double));
+           cudaStat = cudaMalloc((void **)&d_Y_oovv, nO*nO*nV*nV * sizeof(double));
+           assert (cudaStat == cudaSuccess);
 
-           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) {
+           for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
+             size_t mbs = nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE;
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
                cudaStreamCreate(&(stream[gam]));
              }
-             for (int gam=iblock ; gam<mbs ; ++gam) {
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
                cublasSetStream(handle, stream[gam]);
                alpha = 1.0;
                beta  = 0.0;
@@ -1103,7 +1152,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                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) {
+             for (size_t gam=iblock ; gam<mbs ; ++gam) {
                cudaStreamDestroy(stream[gam]);
              }
              cublasSetStream(handle, NULL);
@@ -1119,18 +1168,18 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
            cudaFree(d_X_vovv);
            alpha = 1.0;
            beta =  1.0;
-           for (int i=0 ; i<nV ; ++i) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
-           for (int j=0 ; j<nV ; ++j) {
-             for (int i=0 ; i<nV ; ++i) {
+           for (size_t j=0 ; j<nV ; ++j) {
+             for (size_t i=0 ; i<nV ; ++i) {
                 cublasSetStream(handle, stream[i]);
                 A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                 B = &(d_Y_oovv[nO*nO*(i+nV*j)]); ldb = nO;
                 C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                 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) {
+             for (size_t i=0 ; i<nV ; ++i) {
                 cublasSetStream(handle, stream[i]);
                 A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                 B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
@@ -1138,7 +1187,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                 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) {
+           for (size_t i=0 ; i<nV ; ++i) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);
@@ -1148,7 +1197,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_tcc2;
-          cudaMalloc((void **)&d_tcc2, cholesky_mo_num*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_tcc2, cholesky_mo_num*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1159,7 +1209,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_tcc;
-          cudaMalloc((void **)&d_tcc, cholesky_mo_num*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_tcc, cholesky_mo_num*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1170,7 +1221,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
 
           double* d_X_ovvo;
-          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1182,20 +1234,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_tcc);
           cudaFree(d_tcc2);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t 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++){
+          for(size_t gam = 0; gam < nV; gam++){
+            for(size_t bet = 0; bet < nV; bet++){
                cublasSetStream(handle, stream[bet]);
                A = &(d_r2[nO*nO*(bet+nV*gam)]); lda = nO;
                B = &(d_X_ovvo[nO*(bet+nV*gam)]); ldb = nO*nV*nV;
                C = &(d_r2[nO*nO*(bet+nV*gam)]); ldc = nO;
                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++){
+            for(size_t bet = 0; bet < nV; bet++){
                cublasSetStream(handle, stream[bet]);
                A = &(d_r2[nO*nO*(bet+nV*gam)]); lda = nO;
                B = &(d_X_ovvo[nO*(gam+nV*bet)]); ldb = nO*nV*nV;
@@ -1203,7 +1255,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1213,7 +1265,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_X_oovv;
-          cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1229,11 +1282,11 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int j=0 ; j<nV ; ++j) {
-            for (int i=0 ; i<nV ; ++i) {
+          for (size_t j=0 ; j<nV ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
@@ -1243,12 +1296,13 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           }
 
           double* d_X_vovo;
-          cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 0.0;
           beta =  1.0;
-          for (int i=0 ; i<nO ; ++i) {
-            for (int gam=0 ; gam<nV ; ++gam) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t gam=0 ; gam<nV ; ++gam) {
                cublasSetStream(handle, stream[gam]);
                A = &(d_X_vovo[nV*nO*(gam+nV*i)]); lda = nV;
                B = &(d_cc_space_v_ovvo[nO*nV*(gam+nV*i)]); ldb = nO;
@@ -1256,13 +1310,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Y_oovo;
-          cudaMalloc((void **)&d_Y_oovo, nO*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_oovo, nO*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1286,18 +1341,18 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
           alpha = 1.0;
           beta = -1.0;
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int j=0 ; j<nV ; ++j) {
-            for (int i=0 ; i<nV ; ++i) {
+          for (size_t j=0 ; j<nV ; ++j) {
+            for (size_t i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_X_oovv[nO*nO*(i+nV*j)]); ldb = nO;
                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                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) {
+            for (size_t i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
@@ -1305,7 +1360,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1327,19 +1382,21 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         {
           double* d_J1;
           lda = nO*nV;
-          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
           cublasSetMatrix(lda, nV*nO, sizeof(double), J1, lda, d_J1, lda);
 
           double* d_X_ovvo;
-          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta = -0.5;
-          for (int i=0 ; i<nO ; ++i) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
               B = &(d_K1[nO*nV*(i+nO*b)]); ldb = nO;
@@ -1347,22 +1404,23 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_J1);
 
           double* d_Y_voov;
-          cudaMalloc((void **)&d_Y_voov, nV*nO*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_voov, nV*nO*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 2.0;
           beta = -1.0;
-          for (int v=0 ; v<nO ; ++v) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_t2[nO*(v+nO*nV*g)]); lda = nO*nO;
               B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
@@ -1370,13 +1428,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               cublasDgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Z_ovov;
-          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1388,20 +1447,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_X_ovvo);
           cudaFree(d_Y_voov);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta =  1.0;
-          for (int b=0 ; b<nV ; ++b) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
               C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
-            for (int g=0 ; g<nV ; ++g) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
@@ -1409,7 +1468,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1420,18 +1479,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_X_ovov;
-          cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           double* d_Y_ovov;
-          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 0.5;
           beta =  0.0;
-          for (int a=0 ; a<nV ; ++a) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t a=0 ; a<nV ; ++a) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_K1[nO*(a+nV*nO*b)]); lda = nO*nV;
               B = &(d_K1[nO*(a+nV*nO*b)]); ldb = nO*nV;
@@ -1440,8 +1501,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
             }
           }
           alpha = 1.0;
-          for (int v=0 ; v<nO ; ++v) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_t2[nO*(v+nO*g)]); lda = nO*nO*nV;
               B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
@@ -1449,13 +1510,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Z_ovov;
-          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1468,20 +1530,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_X_ovov);
           cudaFree(d_Y_ovov);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta = -1.0;
-          for (int b=0 ; b<nV ; ++b) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
               C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
-            for (int g=0 ; g<nV ; ++g) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
@@ -1489,7 +1551,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1500,19 +1562,21 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_X_ovov;
-          cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           double* d_Y_ovov;
-          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta =  0.0;
-          for (int a=0 ; a<nV ; ++a) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t a=0 ; a<nV ; ++a) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_K1[nO*(a+nV*nO*g)]); lda = nO*nV;
               B = &(d_K1[nO*(a+nV*nO*g)]); ldb = nO*nV;
@@ -1521,8 +1585,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
             }
           }
           alpha = 1.0;
-          for (int v=0 ; v<nO ; ++v) {
-            for (int b=0 ; b<nV ; ++b) {
+          for (size_t v=0 ; v<nO ; ++v) {
+            for (size_t b=0 ; b<nV ; ++b) {
               cublasSetStream(handle, stream[b]);
               A = &(d_t2[nO*(v+nO*b)]); lda = nO*nO*nV;
               B = &(d_t2[nO*(v+nO*b)]); ldb = nO*nO*nV;
@@ -1530,13 +1594,14 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
 
           double* d_Z_ovov;
-          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 1.0;
           beta  = 0.0;
@@ -1549,20 +1614,20 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaFree(d_X_ovov);
           cudaFree(d_Y_ovov);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = 1.0;
           beta = -1.0;
-          for (int b=0 ; b<nV ; ++b) {
-            for (int g=0 ; g<nV ; ++g) {
+          for (size_t b=0 ; b<nV ; ++b) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(g+nV*nO*b)]); ldb = nO*nV;
               C = &(d_r2[nO*nO*(b+nV*g)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
-            for (int g=0 ; g<nV ; ++g) {
+            for (size_t g=0 ; g<nV ; ++g) {
               cublasSetStream(handle, stream[g]);
               A = &(d_r2[nO*nO*(b+nV*g)]); lda = nO;
               B = &(d_Z_ovov[nO*(b+nV*nO*g)]); ldb = nO*nV;
@@ -1570,7 +1635,7 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1583,7 +1648,8 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
       double* d_tmp_cc;
       lda = cholesky_mo_num * nV;
-      cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       alpha=1.0; beta=0.0;
       m=cholesky_mo_num*nV;  n=nV;  k=nO;
@@ -1591,16 +1657,19 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
       cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, m, B, k, &beta, C, m);
 
       double* d_tmp_cc2;
-      cudaMalloc((void **)&d_tmp_cc2, cholesky_mo_num*nV*sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_tmp_cc2, cholesky_mo_num*nV*sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_B1;
-      cudaMalloc((void**)&d_B1, nV*nV*BLOCK_SIZE*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_B1, nV*nV*BLOCK_SIZE*sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_tmpB1;
-      cudaMalloc((void**)&d_tmpB1, nV*BLOCK_SIZE*nV*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_tmpB1, nV*BLOCK_SIZE*nV*sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       #pragma omp for
-      for (int gam=0 ; gam<nV ; ++gam)
+      for (size_t 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]);
@@ -1612,9 +1681,9 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
          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 (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE)
+         for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE)
          {
-              const int mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+              const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
 
               alpha=-1.0; beta=0.0;
               m=nV*mbs;  n=nV;  k=cholesky_mo_num;
@@ -1689,21 +1758,22 @@ void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
 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;
+    const size_t cholesky_mo_num = data->cholesky_mo_num;
 
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
 
     #pragma omp parallel num_threads(ngpus)
     {
-      int m,n,k, lda, ldb, ldc;
+      cudaError_t cudaStat;
+      size_t 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();
+      size_t igpu = omp_get_thread_num();
       cudaSetDevice(igpu);
 
       cublasHandle_t handle;
@@ -1711,7 +1781,8 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
       double* d_r1;
       lda = nO ;
-      cudaMalloc((void **)&d_r1, lda * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_r1, lda * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cudaMemset(d_r1, 0, nO*nV*sizeof(double));
       memset(r1, 0, nO*nV*sizeof(double));
 
@@ -1737,7 +1808,8 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           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));
+          cudaStat = cudaMalloc((void **)&d_X_oo, nO*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = -2.0;
           beta  = 0.0;
@@ -1783,14 +1855,15 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_X_voov;
-          cudaMalloc((void **)&d_X_voov, nV* nO* nO* nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_voov, nV* nO* nO* nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t 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) {
+          for (size_t i=0 ; i<nO ; ++i) {
+            for (size_t bet=0 ; bet<nV ; ++bet) {
               cublasSetStream(handle, stream[bet]);
               beta = t1[i+bet*nO];
               A = &(d_t2[nO*(i+nO*nV*bet)]); lda = nO*nO;
@@ -1802,14 +1875,14 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cudaDeviceSynchronize();
           alpha = 1.0;
           beta  = 2.0;
-          for (int bet=0 ; bet<nV ; ++bet) {
+          for (size_t bet=0 ; bet<nV ; ++bet) {
             cublasSetStream(handle, stream[bet]);
             A = &(d_X_voov[nV*nO*nO*bet]); lda = nV;
             B = &(d_t2[nO*nO*nV*bet]); ldb = nO*nO;
             C = A ; ldc = lda;
             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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1828,16 +1901,17 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_X_ovov;
-          cudaMalloc((void **)&d_X_ovov, nO* nV* nO* nV * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_X_ovov, nO* nV* nO* nV * sizeof(double));
+          assert (cudaStat == cudaSuccess);
           cublasDcopy(handle, nO*nV*nO*nV, d_cc_space_v_ovov, 1, d_X_ovov, 1);
 
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           alpha = -1.0;
           beta  = 2.0;
-          for (int u=0 ; u<nO ; ++u) {
-            for (int bet=0 ; bet<nV ; ++bet) {
+          for (size_t u=0 ; u<nO ; ++u) {
+            for (size_t bet=0 ; bet<nV ; ++bet) {
               cublasSetStream(handle, stream[bet]);
               A = &(d_X_ovov[nO*nV*(u+nO*bet)]); lda = nO;
               B = &(d_cc_space_v_voov[(nV*(u+nO*nO*bet))]); ldb = nV*nO;
@@ -1845,7 +1919,7 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1865,7 +1939,8 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_T_vvoo;
-          cudaMalloc((void **)&d_T_vvoo, nV*nV*nO*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_T_vvoo, nV*nV*nO*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 0.0;
           beta  = 1.0;
@@ -1875,14 +1950,16 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV*nV, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
 
           double* d_W_vvov;
-          cudaMalloc((void **)&d_W_vvov, nV*nV*nO*BLOCK_SIZE * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_W_vvov, nV*nV*nO*BLOCK_SIZE * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           double* d_W_vvov_tmp;
-          cudaMalloc((void **)&d_W_vvov_tmp, nV*nO*nV*BLOCK_SIZE * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_W_vvov_tmp, nV*nO*nV*BLOCK_SIZE * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
 
-          for (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
-            const int mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
+          for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
+            const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
 
             alpha = 1.0;
             beta  = 0.0;
@@ -1894,12 +1971,12 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
 
             alpha = 2.0;
             beta = -1.0;
-            int kk=0;
-            for (int i=0 ; i<nV ; ++i) {
+            size_t kk=0;
+            for (size_t i=0 ; i<nV ; ++i) {
               cudaStreamCreate(&(stream[i]));
             }
-            for (int i=0 ; i<nO ; ++i) {
-              for (int bet=0 ; bet<mbs ; ++bet) {
+            for (size_t i=0 ; i<nO ; ++i) {
+              for (size_t bet=0 ; bet<mbs ; ++bet) {
                 cublasSetStream(handle, stream[kk]);
                 ++kk;
                 if (kk >= nV) kk = 0;
@@ -1909,7 +1986,7 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
                 cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
               }
             }
-            for (int i=0 ; i<nV ; ++i) {
+            for (size_t i=0 ; i<nV ; ++i) {
               cudaStreamDestroy(stream[i]);
             }
             cublasSetStream(handle, NULL);
@@ -1931,15 +2008,16 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
         #pragma omp section
         {
           double* d_W_oovo;
-          cudaMalloc((void **)&d_W_oovo, nO*nO*nV*nO * sizeof(double));
+          cudaStat = cudaMalloc((void **)&d_W_oovo, nO*nO*nV*nO * sizeof(double));
+          assert (cudaStat == cudaSuccess);
 
           alpha = 2.0;
           beta  = -1.0;
-          for (int i=0 ; i<nV ; ++i) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
-          for (int u=0 ; u<nO ; ++u) {
-            for (int a=0 ; a<nV ; ++a) {
+          for (size_t u=0 ; u<nO ; ++u) {
+            for (size_t a=0 ; a<nV ; ++a) {
               cublasSetStream(handle, stream[a]);
               A = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); lda = nO;
               B = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); ldb = nO;
@@ -1947,7 +2025,7 @@ void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, doubl
               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) {
+          for (size_t i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
@@ -1992,15 +2070,16 @@ double ccsd_energy_space_gpu(gpu_data* data)
 {
     double result = 0.0;
 
-    const int nO = data->nO;
-    const int nV = data->nV;
+    const size_t nO = data->nO;
+    const size_t nV = data->nV;
 
     int ngpus = 1;
     if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
 
     #pragma omp parallel num_threads(ngpus)
     {
-      int igpu = omp_get_thread_num();
+    cudaError_t cudaStat;
+      size_t igpu = omp_get_thread_num();
       cudaSetDevice(igpu);
 
       cublasHandle_t handle;
diff --git a/devel/ccsd_gpu/gpu.h b/devel/ccsd_gpu/gpu.h
index f32bc7a..a74c54d 100644
--- a/devel/ccsd_gpu/gpu.h
+++ b/devel/ccsd_gpu/gpu.h
@@ -1,3 +1,5 @@
+#define MULTIGPU 1
+
 typedef struct {
    double* cc_space_v_oo_chol;
    double* cc_space_v_ov_chol;
@@ -28,4 +30,33 @@ typedef struct {
    int cholesky_mo_num;
 } gpu_data;
 
-#define MULTIGPU 1
+typedef struct {
+   float* cc_space_v_oo_chol;
+   float* cc_space_v_ov_chol;
+   float* cc_space_v_vo_chol;
+   float* cc_space_v_vv_chol;
+   float* cc_space_v_oooo;
+   float* cc_space_v_vooo;
+   float* cc_space_v_voov;
+   float* cc_space_v_oovv;
+   float* cc_space_v_vvoo;
+   float* cc_space_v_oovo;
+   float* cc_space_v_ovvo;
+   float* cc_space_v_ovov;
+   float* cc_space_v_ovoo;
+   float* cc_space_f_oo;
+   float* cc_space_f_ov;
+   float* cc_space_f_vo;
+   float* cc_space_f_vv;
+   float* tau;
+   float* tau_x;
+   float* t1;
+   float* t2;
+   float* H_oo;
+   float* H_vo;
+   float* H_vv;
+   int nO;
+   int nV;
+   int cholesky_mo_num;
+} gpu_data_sp;
+
diff --git a/devel/ccsd_gpu/gpu_dgemm.c b/devel/ccsd_gpu/gpu_dgemm.c
index 6530b2b..9ade048 100644
--- a/devel/ccsd_gpu/gpu_dgemm.c
+++ b/devel/ccsd_gpu/gpu_dgemm.c
@@ -4,6 +4,7 @@
 #include <omp.h>
 #include <cublas_v2.h>
 #include <cuda_runtime.h>
+#include <assert.h>
 
 #define BLOCK_SIZE 16
 
@@ -16,6 +17,7 @@ void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, 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)
 {
+    cudaError_t cudaStat = cudaSuccess;
     cublasHandle_t handle;
     cublasCreate(&handle);
 
@@ -25,36 +27,48 @@ void gpu_dgemm(char transa, char transb, int m, int n, int k, double alpha,
     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);
+      cudaStat = cudaMalloc((void**)&d_A, (size_t) lda*k*sizeof(double));
+      assert(cudaStat == cudaSuccess);
+      cudaStat = cublasSetMatrix(m, k, sizeof(double), A, lda, d_A, lda);
+      assert(cudaStat == cudaSuccess);
       ta = CUBLAS_OP_N;
     } else {
-      cudaMalloc((void**)&d_A, lda*m*sizeof(double));
-      cublasSetMatrix(k, m, sizeof(double), A, lda, d_A, lda);
+      cudaStat = cudaMalloc((void**)&d_A, (size_t) lda*m*sizeof(double));
+      assert(cudaStat == cudaSuccess);
+      cudaStat = cublasSetMatrix(k, m, sizeof(double), A, lda, d_A, lda);
+      assert(cudaStat == cudaSuccess);
       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);
+      cudaStat = cudaMalloc((void**)&d_B, (size_t) ldb*n*sizeof(double));
+      assert(cudaStat == cudaSuccess);
+      cudaStat = cublasSetMatrix(k, n, sizeof(double), B, ldb, d_B, ldb);
+      assert(cudaStat == cudaSuccess);
       tb = CUBLAS_OP_N;
     } else {
-      cudaMalloc((void**)&d_B, ldb*k*sizeof(double));
-      cublasSetMatrix(n, k, sizeof(double), B, ldb, d_B, ldb);
+      cudaStat = cudaMalloc((void**)&d_B, (size_t) ldb*k*sizeof(double));
+      assert(cudaStat == cudaSuccess);
+      cudaStat = cublasSetMatrix(n, k, sizeof(double), B, ldb, d_B, ldb);
+      assert(cudaStat == cudaSuccess);
       tb = CUBLAS_OP_T;
     }
 
-    cudaMalloc((void**)&d_C, ldc*n*sizeof(double));
+    cudaStat = cudaMalloc((void**)&d_C, (size_t) ldc*n*sizeof(double));
+    assert(cudaStat == cudaSuccess);
     if (beta != 0.) {
-      cublasSetMatrix(m, n, sizeof(double), C, ldc, d_C, ldc);
+      cudaStat = cublasSetMatrix(m, n, sizeof(double), C, ldc, d_C, ldc);
+      assert(cudaStat == cudaSuccess);
     }
 
-    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);
-
+    cudaStat = cublasDgemm(handle, ta, tb, m, n, k, &alpha, d_A, lda, d_B, ldb, &beta, d_C, ldc);
+    assert(cudaStat == cudaSuccess);
     cudaFree(d_A);
     cudaFree(d_B);
+
+    cudaStat = cublasGetMatrix(m, n,  sizeof(double), d_C, ldc, C, ldc);
+    assert(cudaStat == cudaSuccess);
+
     cudaFree(d_C);
     cublasDestroy(handle);
 }
diff --git a/devel/ccsd_gpu/gpu_init.c b/devel/ccsd_gpu/gpu_init.c
index 0882e95..d23ae48 100644
--- a/devel/ccsd_gpu/gpu_init.c
+++ b/devel/ccsd_gpu/gpu_init.c
@@ -5,6 +5,7 @@
 #include <cublas_v2.h>
 #include <cuda_runtime.h>
 #include "gpu.h"
+#include "assert.h"
 
 gpu_data* gpu_init(
    int nO, int nV, int cholesky_mo_num,
@@ -18,117 +19,143 @@ gpu_data* gpu_init(
    double* cc_space_f_vo, double* cc_space_f_vv)
 {
     int ngpus = 1;
-    cudaGetDeviceCount(&ngpus);
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
 
     gpu_data* data = (gpu_data*) malloc (ngpus*sizeof(gpu_data));
+    assert (data != NULL);
 
     #pragma omp parallel num_threads(ngpus)
     {
-      int lda;
+      cudaError_t cudaStat = cudaSuccess;
+      size_t 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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_oo_chol, lda * nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_oooo, nO*nO*nO*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_vooo, nV*nO*nO*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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_voov;
-      cudaMalloc((void**)&d_cc_space_v_voov, nV*nO*nO*nV*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_voov, nV*nO*nO*nV*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(nV*nO, nO*nV, sizeof(double), cc_space_v_voov, nV*nO, d_cc_space_v_voov, nV*nO);
 
       double* d_cc_space_v_oovv;
-      cudaMalloc((void**)&d_cc_space_v_oovv, nO*nO*nV*nV*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_oovv, nO*nO*nV*nV*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void**)&d_cc_space_v_vvoo, nV*nV*nO*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_oovo, nO*nO*nV*nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovvo, nO*nV*nV*nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_ovvo, lda, d_cc_space_v_ovvo, lda);
 
       double* d_cc_space_v_ovov;
       lda = nO*nV;
-      cudaMalloc((void **)&d_cc_space_v_ovov, nO*nV*nV*nO * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovov, nO*nV*nV*nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(lda, nV*nO, sizeof(double), cc_space_v_ovov, lda, d_cc_space_v_ovov, lda);
 
       double* d_cc_space_v_ovoo;
       lda = nO*nV;
-      cudaMalloc((void **)&d_cc_space_v_ovoo, nO*nV*nO*nO * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_cc_space_v_ovoo, nO*nV*nO*nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(lda, nO*nO, sizeof(double), cc_space_v_ovoo, lda, d_cc_space_v_ovoo, lda);
 
       double* d_cc_space_f_oo;
-      cudaMalloc((void**)&d_cc_space_f_oo, nO*nO*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_oo, nO*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(nO, nO, sizeof(double), cc_space_f_oo, nO, d_cc_space_f_oo, nO);
 
       double* d_cc_space_f_vo;
-      cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_ov, nV*nO*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       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));
+      cudaStat = cudaMalloc((void**)&d_cc_space_f_vv, nV*nV*sizeof(double));
+      assert (cudaStat == cudaSuccess);
       cublasSetMatrix(nV, nV, sizeof(double), cc_space_f_vv, nV, d_cc_space_f_vv, nV);
 
       double* d_tau;
       lda = nO * nO;
-      cudaMalloc((void **)&d_tau, lda * nV * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_tau, lda * nV * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_tau_x;
       lda = nO * nO;
-      cudaMalloc((void **)&d_tau_x, lda * nV * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_tau_x, lda * nV * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_t1;
-      cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_t2;
-      cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_H_oo;
-      cudaMalloc((void **)&d_H_oo, nO * nO * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_H_oo, nO * nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_H_vo;
-      cudaMalloc((void **)&d_H_vo, nV * nO * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_H_vo, nV * nO * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       double* d_H_vv;
-      cudaMalloc((void **)&d_H_vv, nV * nV * sizeof(double));
+      cudaStat = cudaMalloc((void **)&d_H_vv, nV * nV * sizeof(double));
+      assert (cudaStat == cudaSuccess);
 
       data[igpu].cc_space_v_oo_chol = d_cc_space_v_oo_chol;
       data[igpu].cc_space_v_ov_chol = d_cc_space_v_ov_chol;
@@ -164,3 +191,41 @@ gpu_data* gpu_init(
 }
 
 
+void gpu_deinit(gpu_data* data)
+{
+    int ngpus = 1;
+    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
+
+    #pragma omp parallel num_threads(ngpus)
+    {
+      size_t lda;
+      int igpu = omp_get_thread_num();
+      cudaSetDevice(igpu);
+
+      free(data[igpu].cc_space_v_oo_chol);
+      free(data[igpu].cc_space_v_ov_chol);
+      free(data[igpu].cc_space_v_vo_chol);
+      free(data[igpu].cc_space_v_vv_chol);
+      free(data[igpu].cc_space_v_oooo);
+      free(data[igpu].cc_space_v_vooo);
+      free(data[igpu].cc_space_v_voov);
+      free(data[igpu].cc_space_v_oovv);
+      free(data[igpu].cc_space_v_vvoo);
+      free(data[igpu].cc_space_v_oovo);
+      free(data[igpu].cc_space_v_ovvo);
+      free(data[igpu].cc_space_v_ovov);
+      free(data[igpu].cc_space_v_ovoo);
+      free(data[igpu].cc_space_f_oo);
+      free(data[igpu].cc_space_f_ov);
+      free(data[igpu].cc_space_f_vo);
+      free(data[igpu].cc_space_f_vv);
+      free(data[igpu].tau);
+      free(data[igpu].tau_x);
+      free(data[igpu].t1);
+      free(data[igpu].t2);
+      free(data[igpu].H_oo);
+      free(data[igpu].H_vo);
+      free(data[igpu].H_vv);
+    }
+}
+
diff --git a/devel/ccsd_gpu/gpu_module.f90 b/devel/ccsd_gpu/gpu_module.f90
index 2dceffc..abd81f0 100644
--- a/devel/ccsd_gpu/gpu_module.f90
+++ b/devel/ccsd_gpu/gpu_module.f90
@@ -30,6 +30,32 @@ module gpu_module
         real(c_double), intent(in)  :: cc_space_f_vv(nV,nV)
     end function
 
+    type(c_ptr) function gpu_init_sp(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_voov, 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_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)
+        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_voov(nV,nO,nO,nV)
+        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_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
+
     subroutine gpu_upload(gpu_data, nO, nV, t1, t2) bind(C)
         import c_int, c_double, c_ptr
         type(c_ptr), value    :: gpu_data
@@ -38,21 +64,29 @@ module gpu_module
         real(c_double), intent(in) :: t2(nO,nO,nV,nV)
     end subroutine
 
+    subroutine gpu_upload_sp(gpu_data, nO, nV, t1, t2) 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(in) :: t2(nO,nO,nV,nV)
+    end subroutine
+
 
     subroutine compute_H_oo_chol_gpu(gpu_data, igpu) bind(C)
-        import c_int, c_double, c_ptr
+        import c_int, c_ptr
         type(c_ptr), value    :: gpu_data
         integer(c_int), intent(in), value  :: igpu
     end subroutine
 
     subroutine compute_H_vo_chol_gpu(gpu_data, igpu) bind(C)
-        import c_int, c_double, c_ptr
+        import c_int, c_ptr
         type(c_ptr), value    :: gpu_data
         integer(c_int), intent(in), value  :: igpu
     end subroutine
 
     subroutine compute_H_vv_chol_gpu(gpu_data, igpu) bind(C)
-        import c_int, c_double, c_ptr
+        import c_int, c_ptr
         type(c_ptr), value    :: gpu_data
         integer(c_int), intent(in), value  :: igpu
     end subroutine
@@ -81,6 +115,47 @@ module gpu_module
     end function
 
 
+    subroutine compute_H_oo_chol_gpu_sp(gpu_data, igpu) bind(C)
+        import c_int, c_ptr
+        type(c_ptr), value    :: gpu_data
+        integer(c_int), intent(in), value  :: igpu
+    end subroutine
+
+    subroutine compute_H_vo_chol_gpu_sp(gpu_data, igpu) bind(C)
+        import c_int, c_ptr
+        type(c_ptr), value    :: gpu_data
+        integer(c_int), intent(in), value  :: igpu
+    end subroutine
+
+    subroutine compute_H_vv_chol_gpu_sp(gpu_data, igpu) bind(C)
+        import c_int, c_ptr
+        type(c_ptr), value    :: gpu_data
+        integer(c_int), intent(in), value  :: igpu
+    end subroutine
+
+    subroutine compute_r1_space_chol_gpu_sp(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_sp(gpu_data, nO, nV, t1, r2, max_r2) 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) :: r2(nO,nO,nV,nV)
+        real(c_double), intent(out) :: max_r2
+    end subroutine
+
+    double precision function ccsd_energy_space_gpu_sp(gpu_data) bind(C)
+        import c_ptr
+        type(c_ptr), value    :: gpu_data
+    end function
+
     subroutine gpu_dgemm(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) bind(C)
       import c_int, c_double, c_char
       character(c_char), value :: transa, transb