#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <omp.h>
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include "gpu.h"

void gpu_upload(gpu_data* data,
                int nO, int nV,
                double* t1,
                double* t2,
                double* tau,
                double* tau_x)
{
    int lda;
    const int cholesky_mo_num = data->cholesky_mo_num;

    int ngpus = 1;
    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);

    #pragma omp parallel num_threads(ngpus)
    {
        int igpu = omp_get_thread_num();
        cudaSetDevice(igpu);

        double* d_tau = data[igpu].tau;
        lda = nO * nO;
        cublasSetMatrix(nO*nO, nV*nV, sizeof(double), tau, lda, d_tau, lda);

        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);

        double* d_t1 = data[igpu].t1;
        lda = nO;
        cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);

        double* d_t2 = data[igpu].t2;
        lda = nO*nO;
        cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
    }

}


void compute_h_oo_chol_gpu(gpu_data* data, int igpu)
{
    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;
    cudaSetDevice(igpu);

    int m,n,k, lda, ldb, ldc;
    double alpha, beta;
    double* A;
    double* B;
    double* C;
    cudaStream_t stream[nV];

    cublasHandle_t handle;
    cublasCreate(&handle);

    double* d_H_oo               = data[igpu].H_oo;
    double* d_tau_x              = data[igpu].tau_x;
    double* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
    double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;

    double* d_tau_kau;
    cudaMalloc((void **)&d_tau_kau, cholesky_mo_num*nV*nO * sizeof(double));

    double* d_tmp_ovv;
    cudaMalloc((void **)&d_tmp_ovv, nO*nV*nV  * sizeof(double));

    double* d_tmp_vov;
    cudaMalloc((void **)&d_tmp_vov, nV*nO*nV  * sizeof(double));

    for (int i=0 ; i<nV ; ++i) {
      cudaStreamCreate(&(stream[i]));
    }
    alpha = 1.0;
    beta = 0.0;
    for (int u=0 ; u<nO ; ++u) {
      cublasDcopy(handle, nO*nV*nV, &(d_tau_x[u]), nO, d_tmp_ovv, 1);
      for (int b=0 ; b<nV ; ++b) {
            cublasSetStream(handle, stream[b]);
            A = &(d_tmp_ovv[nO*nV*b]); lda = nO;
            B = &(d_tmp_ovv[nO*nV*b]); ldb = nO;
            C = &(d_tmp_vov[nV*nO*b]); ldc = nV;
            cublasDgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
      }
      cudaDeviceSynchronize();
      cublasSetStream(handle, NULL);
      alpha = 1.0;
      beta = 0.0;
      m=cholesky_mo_num;  n=nV;  k=nO*nV;
      A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
      B=d_tmp_vov;  ldb=nV;
      C=&(d_tau_kau[cholesky_mo_num*nV*u]);  ldc=cholesky_mo_num;
      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
    }
    for (int i=0 ; i<nV ; ++i) {
      cudaStreamDestroy(stream[i]);
    }

    cudaFree(d_tmp_vov);
    cudaFree(d_tmp_ovv);

    cublasDcopy(handle, nO*nO, d_cc_space_f_oo, 1, d_H_oo, 1);
    alpha = 1.0;
    beta  = 1.0;
    m=nO;  n=nO;  k=cholesky_mo_num*nV;
    A=d_tau_kau;  lda=cholesky_mo_num*nV;
    B=d_cc_space_v_vo_chol;  ldb=cholesky_mo_num*nV;
    C=d_H_oo;  ldc=nO;
    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    cudaFree(d_tau_kau);

    double* H_oo = malloc(nO*nO*sizeof(double));
    cublasGetMatrix(nO, nO, sizeof(double), d_H_oo, nO, H_oo, nO);
    for (int i=0 ; i<ngpus ; ++i) {
       if (i != igpu) {
         double* d_H_oo = data[i].H_oo;
         cudaSetDevice(i);
         cublasSetMatrix(nO, nO, sizeof(double), H_oo, nO, d_H_oo, nO);
       }
    }
    free(H_oo);

    cublasDestroy(handle);

}


void compute_h_vo_chol_gpu(gpu_data* data, int igpu)
{
    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;
    cudaSetDevice(igpu);

    int m,n,k, lda, ldb, ldc;
    double alpha, beta;
    double* A;
    double* B;
    double* C;
    cudaStream_t stream[nV];

    cublasHandle_t handle;
    cublasCreate(&handle);

    double* d_t1                 = data[igpu].t1;
    double* d_H_vo               = data[igpu].H_vo;
    double* d_tau_x              = data[igpu].tau_x;
    double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
    double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;

    cublasDcopy(handle, nV*nO, d_cc_space_f_vo, 1, d_H_vo, 1);

    double* d_tmp_k;
    cudaMalloc((void **)&d_tmp_k, cholesky_mo_num * sizeof(double));

    alpha = 2.0;
    beta  = 0.0;
    m=cholesky_mo_num;  n=1;  k=nO*nV;
    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
    B=d_t1;  ldb=nO*nV;
    C=d_tmp_k;  ldc=cholesky_mo_num;
    cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    alpha = 1.0;
    beta  = 1.0;
    m=nV*nO;  n=1;  k=cholesky_mo_num;
    A=d_cc_space_v_vo_chol;  lda=cholesky_mo_num;
    B=d_tmp_k;  ldb=cholesky_mo_num;
    C=d_H_vo;  ldc=nV*nO;
    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    cudaFree(d_tmp_k);

    double* d_tmp;
    cudaMalloc((void **)&d_tmp, cholesky_mo_num*nO*nO * sizeof(double));

    alpha = 1.0;
    beta  = 0.0;
    m=cholesky_mo_num*nO;  n=nO;  k=nV;
    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
    B=d_t1;  ldb=nO;
    C=d_tmp;  ldc=cholesky_mo_num*nO;
    cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    double* d_tmp2;
    cudaMalloc((void **)&d_tmp2, cholesky_mo_num*nO*nO * sizeof(double));

    for (int i=0 ; i<nV ; ++i) {
      cudaStreamCreate(&(stream[i]));
    }
    for (int i=0 ; i<nO ; ++i) {
      for (int j=0 ; j<nO ; ++j) {
         cublasSetStream(handle, stream[j]);
         cublasDcopy(handle, cholesky_mo_num, &(d_tmp [cholesky_mo_num*(i+nO*j)]), 1,
                                              &(d_tmp2[cholesky_mo_num*(j+nO*i)]), 1);
      }
    }
    for (int i=0 ; i<nV ; ++i) {
      cudaStreamDestroy(stream[i]);
    }
    cublasSetStream(handle, NULL);

    alpha = -1.0;
    beta  = 1.0;
    m=nV;  n=nO;  k=cholesky_mo_num*nO;
    A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
    B=d_tmp2;  ldb=cholesky_mo_num*nO;
    C=d_H_vo;  ldc=nV;
    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    double* H_vo = malloc(nV*nO*sizeof(double));
    cublasGetMatrix(nV, nO, sizeof(double), d_H_vo, nV, H_vo, nV);
    for (int i=0 ; i<ngpus ; ++i) {
       if (i != igpu) {
         double* d_H_vo = data[i].H_vo;
         cudaSetDevice(i);
         cublasSetMatrix(nV, nO, sizeof(double), H_vo, nV, d_H_vo, nV);
       }
    }
    free(H_vo);

    cublasDestroy(handle);

}



void compute_h_vv_chol_gpu(gpu_data* data, int igpu)
{
    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;
    cudaSetDevice(igpu);

    int m,n,k, lda, ldb, ldc;
    double alpha, beta;
    double* A;
    double* B;
    double* C;
    cudaStream_t stream[nV];

    cublasHandle_t handle;
    cublasCreate(&handle);

    double* d_H_vv               = data[igpu].H_vv;
    double* d_tau_x              = data[igpu].tau_x;
    double* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
    double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;

    double* d_tau_kia;
    cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(double));

    double* d_tmp_oov;
    cudaMalloc((void **)&d_tmp_oov, nO*nO*nV  * sizeof(double));

    alpha = 1.0;
    beta = 0.0;
    for (int i=0 ; i<nV ; ++i) {
      cudaStreamCreate(&(stream[i]));
    }
    for (int a=0 ; a<nV ; ++a) {
      for (int b=0 ; b<nV ; ++b) {
         cublasSetStream(handle, stream[b]);
         cublasDcopy(handle, nO*nO, &(d_tau_x[nO*nO*(a+nV*b)]), 1, &(d_tmp_oov[nO*nO*b]), 1);
      }
      cudaDeviceSynchronize();
      cublasSetStream(handle, NULL);
      alpha = 1.0;
      beta = 0.0;
      m=cholesky_mo_num;  n=nO;  k=nO*nV;
      A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
      B=d_tmp_oov;  ldb=nO;
      C=&(d_tau_kia[cholesky_mo_num*nO*a]);  ldc=cholesky_mo_num;
      cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
    }
    for (int i=0 ; i<nV ; ++i) {
      cudaStreamDestroy(stream[i]);
    }
    cudaFree(d_tmp_oov);

    cublasDcopy(handle, nV*nV, d_cc_space_f_vv, 1, d_H_vv, 1);
    alpha = -1.0;
    beta  = 1.0;
    m=nV;  n=nV;  k=cholesky_mo_num*nO;
    A=d_tau_kia;  lda=cholesky_mo_num*nO;
    B=d_cc_space_v_ov_chol;  ldb=cholesky_mo_num*nO;
    C=d_H_vv;  ldc=nV;
    cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

    cudaFree(d_tau_kia);

    double* H_vv = malloc(nV*nV*sizeof(double));
    cublasGetMatrix(nV, nV, sizeof(double), d_H_vv, nV, H_vv, nV);
    for (int i=0 ; i<ngpus ; ++i) {
       if (i != igpu) {
         double* d_H_vv = data[i].H_vv;
         cudaSetDevice(i);
         cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_H_vv, nV);
       }
    }
    free(H_vv);

    cublasDestroy(handle);

}




void compute_r2_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, double* r2, double* max_r2)
{
    const int cholesky_mo_num = data->cholesky_mo_num;

    int ngpus = 1;
    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);

    double* J1 = malloc(nO*nV*nV*nO*sizeof(double));
    double* K1 = malloc(nO*nV*nV*nO*sizeof(double));

    #pragma omp parallel num_threads(ngpus)
    {
      int m,n,k, lda, ldb, ldc;
      double alpha, beta;
      double* A;
      double* B;
      double* C;
      cudaStream_t stream[nV];

      int igpu = omp_get_thread_num();
      cudaSetDevice(igpu);

      cublasHandle_t handle;
      cublasCreate(&handle);

      double* d_r2;
      lda = nO * nO;
      cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
      cudaMemset(d_r2, 0, nO*nO*nV*nV*sizeof(double));
      memset(r2, 0, nO*nO*nV*nV*sizeof(double));

      double* d_cc_space_v_oo_chol = data[igpu].cc_space_v_oo_chol;
      double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
      double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
      double* d_cc_space_v_vv_chol = data[igpu].cc_space_v_vv_chol;
      double* d_cc_space_v_oooo    = data[igpu].cc_space_v_oooo;
      double* d_cc_space_v_vooo    = data[igpu].cc_space_v_vooo;
      double* d_cc_space_v_oovv    = data[igpu].cc_space_v_oovv;
      double* d_cc_space_v_vvoo    = data[igpu].cc_space_v_vvoo;
      double* d_cc_space_v_oovo    = data[igpu].cc_space_v_oovo;
      double* d_cc_space_v_ovvo    = data[igpu].cc_space_v_ovvo;
      double* d_cc_space_v_ovov    = data[igpu].cc_space_v_ovov;
      double* d_cc_space_v_ovoo    = data[igpu].cc_space_v_ovoo;
      double* d_cc_space_f_oo      = data[igpu].cc_space_f_oo;
      double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
      double* d_cc_space_f_vv      = data[igpu].cc_space_f_vv;
      double* d_tau                = data[igpu].tau;
      double* d_t1                 = data[igpu].t1;
      double* d_t2                 = data[igpu].t2;
      double* d_H_oo               = data[igpu].H_oo;
      double* d_H_vv               = data[igpu].H_vv;

      double* d_K1;
      cudaMalloc((void **)&d_K1, nO*nV*nO*nV * sizeof(double));

      #pragma omp sections
      {

        #pragma omp section
        {
          double* d_J1;
          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));

          alpha = 1.0;
          beta =  0.0;
          A = d_cc_space_v_ovvo; lda = nO*nV;
          B = d_cc_space_v_ovvo; ldb = nO*nV;
          C = d_J1; ldc = nO*nV;
          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nV*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);


          double* d_X_ovoo;
          cudaMalloc((void **)&d_X_ovoo, nO*nV*nO*nO * sizeof(double));
          alpha = 0.0;
          beta =  1.0;
          for (int i=0 ; i<nO ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int j=0 ; j<nO ; ++j) {
            for (int i=0 ; i<nO ; ++i) {
              cublasSetStream(handle, stream[i]);
              A = &(d_X_ovoo[nO*nV*(i+nO*j)]); lda = nO;
              B = &(d_cc_space_v_ovoo[nO*nV*(j+nO*i)]); ldb = nO;
              C = &(d_X_ovoo[nO*nV*(i+nO*j)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nO ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);


          double* d_Y_ovov;
          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV*nO;  n=nV;  k=nO;
          A=d_X_ovoo;  lda=nO*nV*nO;
          B=d_t1;  ldb=nO;
          C=d_Y_ovov;  ldc=nO*nV*nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_X_ovoo);

          alpha = 1.0;
          beta =  -1.0;
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int j=0 ; j<nO ; ++j) {
            for (int i=0 ; i<nV ; ++i) {
              cublasSetStream(handle, stream[i]);
              A = &(d_J1[nO*nV*(i+nV*j)]); lda = nO;
              B = &(d_Y_ovov[nO*nV*(j+nO*i)]); ldb = nO;
              C = &(d_J1[nO*nV*(i+nV*j)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_tmp_cc;
          cudaMalloc((void **)&d_tmp_cc, cholesky_mo_num*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=cholesky_mo_num*nV;  n=nO;  k=nV;
          A=d_cc_space_v_vv_chol;  lda=cholesky_mo_num*nV;
          B=d_t1;  ldb=nO;
          C=d_tmp_cc;  ldc=cholesky_mo_num*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          double* d_J1_tmp;
          cudaMalloc((void **)&d_J1_tmp, nV*nO*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nV*nO;  n=nV*nO;  k=cholesky_mo_num;
          A=d_tmp_cc;  lda=cholesky_mo_num;
          B=d_cc_space_v_vo_chol;  ldb=cholesky_mo_num;
          C=d_J1_tmp;  ldc=nV*nO;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_tmp_cc);

          alpha = 1.0;
          beta =  1.0;
          for (int i=0 ; i<nO ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int i=0 ; i<nO ; ++i) {
              cublasSetStream(handle, stream[i]);
              A = &(d_J1[nO*nV*nV*i]); lda = nO*nV;
              B = &(d_J1_tmp[nV*nO*nV*i]); ldb = nV;
              C = &(d_J1[nO*nV*nV*i]); ldc = nO*nV;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO*nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
          }
          for (int i=0 ; i<nO ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_J1_tmp);

          double* d_X_voov;
          cudaMalloc((void **)&d_X_voov, nV*nO*nO*nV * sizeof(double));

          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          alpha = 0.5;
          for (int j=0 ; j<nO ; ++j) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              beta = t1[j+b*nO];
              A = &(d_t2[nO*(j + nO*nV*b)]); lda = nO*nO;
              B = d_t1; ldb = nO;
              C = &(d_Y_ovov[nO*(b+nV*j)]); ldc = nO*nV*nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          alpha = 1.0;
          beta  = 0.0;
          for (int j=0 ; j<nO ; ++j) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              A = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); lda = nV*nV;
              B = &(d_cc_space_v_vvoo[nV*(b+nV*nO*j)]); ldb = nV*nV;
              C = &(d_X_voov[nV*nO*(j+nO*b)]); ldc = nV;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Z_ovvo;
          cudaMalloc((void **)&d_Z_ovvo, nO*nV*nV*nO * sizeof(double));

          alpha = -1.0;
          beta  = 0.0;
          m=nO*nV;  n=nV*nO;  k=nO*nV;
          A=d_Y_ovov;  lda=nO*nV;
          B=d_X_voov;  ldb=nV*nO;
          C=d_Z_ovvo;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_X_voov);
          cudaFree(d_Y_ovov);

          alpha = 1.0;
          beta =  1.0;
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int i=0 ; i<nO ; ++i) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
              B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb=nO*nV;
              C = &(d_J1[nO*nV*(b+nV*i)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }

          double* d_Y_vovo;
          cudaMalloc((void **)&d_Y_vovo, nV*nO*nV*nO * sizeof(double));

          alpha = 1.0;
          beta = -0.5;
          for (int j=0 ; j<nO ; ++j) {
            for (int i=0 ; i<nO ; ++i) {
              cublasSetStream(handle, stream[i]);
              A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
              B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
              C = &(d_Y_vovo[nV*(i+nO*nV*j)]); ldc = nV*nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }

          double* d_X_ovvo;
          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));

          alpha = 1.0;
          beta =  0.0;
          for (int j=0 ; j<nO ; ++j) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
              B = &(d_t2[nO*(j+nO*nV*b)]); ldb = nO*nO;
              C = &(d_X_ovvo[nO*nV*(b+nV*j)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV;  n=nV*nO;  k=nV*nO;
          A=d_X_ovvo;  lda=nO*nV;
          B=d_Y_vovo;  ldb=nV*nO;
          C=d_Z_ovvo;  ldc=nO*nV;
         cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_X_ovvo);
          cudaFree(d_Y_vovo);

          for (int 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) {
              cublasSetStream(handle, stream[b]);
              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
              B = &(d_Z_ovvo[nO*(b+nV*nV*i)]); ldb = nO*nV;
              C = &(d_J1[nO*nV*(b+nV*i)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          cudaFree(d_Z_ovvo);

          lda = nO*nV;
          cublasGetMatrix(nO*nV, nV*nO, sizeof(double), d_J1, lda, J1, lda);
          cudaFree(d_J1);

        }

        #pragma omp section
        {

          alpha = 1.0;
          beta =  0.0;
          A = d_cc_space_v_ovov; lda = nO*nV;
          B = d_cc_space_v_ovov; ldb = nO*nV;
          C = d_K1; ldc = nO*nV;
          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nV, nO*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);

          alpha = -1.0;
          beta  = 1.0;
          m=nO*nV*nO;  n=nV;  k=nO;
          A=d_cc_space_v_ovoo;  lda=nO*nV*nO;
          B=d_t1;  ldb=nO;
          C=d_K1;  ldc=nO*nV*nO;
          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));

          double* d_Y;
          cudaMalloc((void **)&d_Y, nO*nV*nV*nO * sizeof(double));

          for (int 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) {
              cublasSetStream(handle, stream[i]);
              A = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); lda = nV;
              B = &(d_cc_space_v_vvoo[nV*nV*(i+nO*j)]); ldb = nV;
              C = &(d_X[nV*(j+nO*nV*i)]); ldc = nV*nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }

          alpha = 0.5;
          for (int j=0 ; j<nO ; ++j) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              beta = t1[j+b*nO];
              A = &(d_t2[nO*(j+nO*nV*b)]); lda = nO*nO;
              B = d_t1; ldb = nO;
              C = &(d_Y[nO*(b+nV*nV*j)]); ldc = nO*nV;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Z;
          cudaMalloc((void **)&d_Z, nO*nV*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nV*nO;  n=nO*nV;  k=nV*nO;
          A=d_Y;  lda=nO*nV;
          B=d_X;  ldb=nV*nO;
          C=d_Z;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X);
          cudaFree(d_Y);


          double* d_t1v;
          cudaMalloc((void **)&d_t1v, cholesky_mo_num*nO*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=cholesky_mo_num*nO;  n=nO;  k=nV;
          A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num*nO;
          B=d_t1;  ldb=nO;
          C=d_t1v;  ldc=cholesky_mo_num*nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          double* d_K1tmp;
          cudaMalloc((void **)&d_K1tmp, nO*nO*nV*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nO;  n=nV*nV;  k=cholesky_mo_num;
          A=d_t1v;  lda=cholesky_mo_num;
          B=d_cc_space_v_vv_chol;  ldb=cholesky_mo_num;
          C=d_K1tmp;  ldc=nO*nO;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_t1v);


          for (int 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) {
              cublasSetStream(handle, stream[b]);
              A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
              B = &(d_K1tmp[nO*(i+nO*nV*b)]); ldb = nO*nO;
              C = &(d_K1[nO*nV*(i+nO*b)]); ldc = nO;
              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 (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              A = &(d_K1[nO*nV*(i+nO*b)]); lda = nO;
              B = &(d_Z[nO*(b+nV*nV*i)]); ldb = nO*nV;
              C = &(d_K1[nO*nV*(i+nO*b)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }

          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_K1tmp);
          cudaFree(d_Z);

          lda = nO*nV;
          cublasGetMatrix(nO*nV, nO*nV, sizeof(double), d_K1, lda, K1, lda);

        }

        #pragma omp section
        {
           double* d_Y_oooo;
           cudaMalloc((void**)&d_Y_oooo, nO*nO*nO*nO*sizeof(double));

           alpha = 1.0;
           beta =  0.0;
           m=nO ; n=nO*nO*nO; k=nV;
           A = d_t1 ; lda = nO;
           B = d_cc_space_v_vooo ; ldb = nV;
           C = d_Y_oooo; ldc = nO;
           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));

           alpha = 1.0;
           beta =  1.0;
           A = d_cc_space_v_oooo; lda = nO*nO;
           B = d_Y_oooo; ldb = nO*nO;
           C = d_A1; ldc = nO*nO;
           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
           for (int i=0 ; i<nO ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           for (int j=0 ; j<nO ; ++j) {
             for (int i=0 ; i<nO ; ++i) {
                cublasSetStream(handle, stream[i]);
                alpha = 1.0;
                beta =  1.0;
                A = &(d_A1[nO*nO*(i+nO*j)]); lda = nO;
                B = &(d_Y_oooo[nO*nO*(j+nO*i)]); ldb = nO;
                C = &(d_A1[nO*nO*(i+nO*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<nO ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_Y_oooo);

           alpha = 1.0;
           beta =  1.0;
           m=nO*nO ; n=nO*nO; k=nV*nV;
           A = d_tau ; lda = nO*nO;
           B = d_cc_space_v_vvoo ; ldb = nV*nV;
           C = d_A1; ldc = nO*nO;
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

           alpha = 1.0;
           beta =  0.0;
           m=nO*nO ; n=nV*nV; k=nO*nO;
           A = d_A1  ; lda = nO*nO;
           B = d_tau ; ldb = nO*nO;
           C = d_r2; ldc = nO*nO;
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
           cudaFree(d_A1);
        }

        // g_vir
        #pragma omp section
        {
           double* d_g_vir;
           cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
           cublasDcopy(handle, nV*nV, d_H_vv, 1, d_g_vir, 1);

           alpha = -1.0;
           beta =  1.0;
           m=nV ; n=nV; k=nO;
           A = d_cc_space_f_vo  ; lda = nV;
           B = d_t1 ; ldb = nO;
           C = d_g_vir; ldc = nV;
           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));
           alpha = 1.0;
           beta =  0.0;
           m=cholesky_mo_num ; n=1; k=nO*nV;
           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num;
           B = d_t1 ; ldb = nO*nV;
           C = d_tmp_k; ldc = cholesky_mo_num;
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

           alpha = 2.0;
           beta =  1.0;
           m=nV*nV; n=1; k=cholesky_mo_num;
           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num;
           B = d_tmp_k ; ldb = cholesky_mo_num;
           C = d_g_vir; ldc = nV*nV;
           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
           cudaFree(d_tmp_k);

           double* d_tmp_vo;
           cudaMalloc((void**)&d_tmp_vo, cholesky_mo_num*nV*nO*sizeof(double));
           alpha = 1.0;
           beta =  0.0;
           m=cholesky_mo_num*nV ; n=nO; k=nV;
           A = d_cc_space_v_vv_chol; lda = cholesky_mo_num*nV;
           B = d_t1 ; ldb = nO;
           C = d_tmp_vo; ldc = cholesky_mo_num*nV;
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

           double* d_tmp_vo2;
           cudaMalloc((void**)&d_tmp_vo2, cholesky_mo_num*nV*nO*sizeof(double));
           for (int i=0 ; i<nO ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           for (int i=0 ; i<nO ; ++i) {
             cublasSetStream(handle, stream[i]);
             alpha = -1.0;
             beta =  0.0;
             A = &(d_tmp_vo[cholesky_mo_num*nV*i]); lda = cholesky_mo_num;
             B = &(d_tmp_vo[cholesky_mo_num*nV*i]); ldb = cholesky_mo_num;
             C = &(d_tmp_vo2[cholesky_mo_num*i]); ldc = cholesky_mo_num*nO;
             cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, cholesky_mo_num, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
           }
           for (int i=0 ; i<nO ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_tmp_vo);

           alpha = 1.0;
           beta =  1.0;
           m=nV ; n=nV; k=nO*cholesky_mo_num;
           A = d_cc_space_v_ov_chol; lda = cholesky_mo_num*nO;
           B = d_tmp_vo2 ; ldb = cholesky_mo_num*nO;
           C = d_g_vir; ldc = nV;
           cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
           cudaFree(d_tmp_vo2);

           double* d_Y_oovv;
           cudaMalloc((void**)&d_Y_oovv, nO*nO*nV*nV*sizeof(double));
           alpha = 1.0;
           beta =  0.0;
           m=nO*nO*nV ; n=nV; k=nV;
           A = d_t2; lda = nO*nO*nV;
           B = d_g_vir; ldb = nV;
           C = d_Y_oovv; ldc = nO*nO*nV;
           cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
           cudaFree(d_g_vir);

           alpha = 1.0;
           beta =  1.0;
           A = d_r2; lda = nO*nO;
           B = d_Y_oovv; ldb = nO*nO;
           C = d_r2; ldc = nO*nO;
           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);

           for (int i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           for (int j=0 ; j<nV ; ++j) {
             for (int i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                alpha = 1.0;
                beta =  1.0;
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);

             }
           }
           for (int i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_Y_oovv);
        }

        // g_occ
        #pragma omp section
        {
          double* d_g_occ;
          lda = nO;
          cudaMalloc((void **)&d_g_occ, nO*nO * sizeof(double));
          cublasDcopy(handle, nO*nO, d_H_oo, 1, d_g_occ, 1);

          double* d_X;
          cudaMalloc((void **)&d_X, cholesky_mo_num*sizeof(double));

          alpha = 2.0;
          beta  = 0.0;
          m=cholesky_mo_num;  n=nO*nV;
          A=d_cc_space_v_ov_chol;  lda=cholesky_mo_num;
          B=d_t1;  ldb=1;
          C=d_X;  ldc=1;
          cublasDgemv(handle, CUBLAS_OP_N, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);

          alpha = 1.0;
          beta  = 1.0;
          m=cholesky_mo_num;  n=nO*nO;
          A=d_cc_space_v_oo_chol;  lda=cholesky_mo_num;
          B=d_X;  ldb=1;
          C=d_g_occ;  ldc=1;
          cublasDgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X);
          alpha = -1.0;
          beta  = 1.0;
          m=nO*nV;  n=nO*nO;
          A=d_cc_space_v_ovoo;  lda=nO*nV;
          B=d_t1;  ldb=1;
          C=d_g_occ;  ldc=1;
          cublasDgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);

          alpha = 1.0;
          beta  = 1.0;
          m=nO;  n=nO;  k=nV;
          A=d_t1;  lda=nO;
          B=d_cc_space_f_vo;  ldb=nV;
          C=d_g_occ;  ldc=nO;
          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));

          alpha = 1.0;
          beta  = 0.0;
          m=nO;  n=nO*nV*nV;  k=nO;
          A=d_g_occ;  lda=nO;
          B=d_t2;  ldb=nO;
          C=d_X_oovv;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_g_occ);

           alpha = 1.0;
           beta = -1.0;
           A = d_r2; lda = nO*nO;
           B = d_X_oovv; ldb = nO*nO;
           C = d_r2; ldc = nO*nO;
           cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);

           for (int i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           for (int j=0 ; j<nV ; ++j) {
             for (int i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                alpha = 1.0;
                beta = -1.0;
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
           for (int i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_X_oovv);
        }

      } // end sections

      lda = nO*nV;
      cublasSetMatrix(lda, nO*nV, sizeof(double), K1, lda, d_K1, lda);

#define BLOCK_SIZE 16

      #pragma omp sections
      {

        #pragma omp section
        {
           double* d_X_vovv;
           cudaMalloc((void **)&d_X_vovv, nV*nO*nV*BLOCK_SIZE * sizeof(double));

           double* d_Y_oovv;
           cudaMalloc((void **)&d_Y_oovv, nO*nO*nV*nV * sizeof(double));

           for (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
             int mbs = nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE;
             for (int gam=iblock ; gam<mbs ; ++gam) {
               cudaStreamCreate(&(stream[gam]));
             }
             for (int gam=iblock ; gam<mbs ; ++gam) {
               cublasSetStream(handle, stream[gam]);
               alpha = 1.0;
               beta  = 0.0;
               m=nV;  n=nO*nV;  k=cholesky_mo_num;
               A=&(d_cc_space_v_vv_chol[cholesky_mo_num*nV*gam]);  lda=cholesky_mo_num;
               B=d_cc_space_v_ov_chol;  ldb=cholesky_mo_num;
               C=&(d_X_vovv[nV*nO*nV*(gam-iblock)]);  ldc=nV;
               cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
             }
             for (int gam=iblock ; gam<mbs ; ++gam) {
               cudaStreamDestroy(stream[gam]);
             }
             cublasSetStream(handle, NULL);
             mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;
             alpha = 1.0;
             beta  = 0.0;
             m=nO;  n=nO*nV*mbs;  k=nV;
             A=d_t1;  lda=nO;
             B=d_X_vovv;  ldb=nV;
             C=&(d_Y_oovv[nO*nO*nV*iblock]);  ldc=nO;
             cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
           }
           cudaFree(d_X_vovv);
           alpha = 1.0;
           beta =  1.0;
           for (int i=0 ; i<nV ; ++i) {
             cudaStreamCreate(&(stream[i]));
           }
           for (int j=0 ; j<nV ; ++j) {
             for (int i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_Y_oovv[nO*nO*(i+nV*j)]); ldb = nO;
                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
             for (int i=0 ; i<nV ; ++i) {
                cublasSetStream(handle, stream[i]);
                A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
                B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
                C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }
           }
           for (int i=0 ; i<nV ; ++i) {
             cudaStreamDestroy(stream[i]);
           }
           cublasSetStream(handle, NULL);
           cudaFree(d_Y_oovv);
        }

        #pragma omp section
        {
          double* d_tcc2;
          cudaMalloc((void **)&d_tcc2, cholesky_mo_num*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=cholesky_mo_num*nV;  n=nO;  k=nV;
          A=d_cc_space_v_vv_chol;  lda=cholesky_mo_num*nV;
          B=d_t1;  ldb=nO;
          C=d_tcc2;  ldc=cholesky_mo_num*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          double* d_tcc;
          cudaMalloc((void **)&d_tcc, cholesky_mo_num*nO*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=cholesky_mo_num*nO;  n=nV;  k=nO;
          A=d_cc_space_v_oo_chol;  lda=cholesky_mo_num*nO;
          B=d_t1;  ldb=nO;
          C=d_tcc;  ldc=cholesky_mo_num*nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          double* d_X_ovvo;
          cudaMalloc((void **)&d_X_ovvo, nO*nV*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV;  n=nV*nO;  k=cholesky_mo_num;
          A=d_tcc;  lda=cholesky_mo_num;
          B=d_tcc2;  ldb=cholesky_mo_num;
          C=d_X_ovvo;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_tcc);
          cudaFree(d_tcc2);

          for (int 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++){
               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++){
               cublasSetStream(handle, stream[bet]);
               A = &(d_r2[nO*nO*(bet+nV*gam)]); lda = nO;
               B = &(d_X_ovvo[nO*(gam+nV*bet)]); ldb = nO*nV*nV;
               C = &(d_r2[nO*nO*(bet+nV*gam)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_X_ovvo);
        }

        #pragma omp section
        {
          double* d_X_oovv;
          cudaMalloc((void **)&d_X_oovv, nO*nO*nV*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nO*nV;  n=nV;  k=nO;
          A=d_cc_space_v_oovo;  lda=nO*nO*nV;
          B=d_t1;  ldb=nO;
          C=d_X_oovv;  ldc=nO*nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          alpha = 1.0;
          beta = -1.0;
          A = d_r2; lda = nO*nO;
          B = d_X_oovv; ldb = nO*nO;
          C = d_r2; ldc = nO*nO;
          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int j=0 ; j<nV ; ++j) {
            for (int i=0 ; i<nV ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
               B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }

          double* d_X_vovo;
          cudaMalloc((void **)&d_X_vovo, nV*nO*nV*nO * sizeof(double));

          alpha = 0.0;
          beta =  1.0;
          for (int i=0 ; i<nO ; ++i) {
            for (int gam=0 ; gam<nV ; ++gam) {
               cublasSetStream(handle, stream[gam]);
               A = &(d_X_vovo[nV*nO*(gam+nV*i)]); lda = nV;
               B = &(d_cc_space_v_ovvo[nO*nV*(gam+nV*i)]); ldb = nO;
               C = &(d_X_vovo[nV*nO*(gam+nV*i)]); ldc = nV;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Y_oovo;
          cudaMalloc((void **)&d_Y_oovo, nO*nO*nV*nO * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO;  n=nO*nV*nO;  k=nV;
          A=d_t1;  lda=nO;
          B=d_X_vovo;  ldb=nV;
          C=d_Y_oovo;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_vovo);

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nO*nV;  n=nV;  k=nO;
          A=d_Y_oovo;  lda=nO*nO*nV;
          B=d_t1;  ldb=nO;
          C=d_X_oovv;  ldc=nO*nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_Y_oovo);

          alpha = 1.0;
          beta = -1.0;
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int j=0 ; j<nV ; ++j) {
            for (int i=0 ; i<nV ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
               B = &(d_X_oovv[nO*nO*(i+nV*j)]); ldb = nO;
               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
            for (int i=0 ; i<nV ; ++i) {
               cublasSetStream(handle, stream[i]);
               A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
               B = &(d_X_oovv[nO*nO*(j+nV*i)]); ldb = nO;
               C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
               cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_X_oovv);
        }

        #pragma omp section
        {
          alpha = 1.0;
          beta =  1.0;
          A = d_r2; lda = nO*nO;
          B = d_cc_space_v_oovv; ldb = nO*nO;
          C = d_r2; ldc = nO*nO;
          cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO*nO, nV*nV, &alpha, A, lda, &beta, B, ldb, C, ldc);

        }

        #pragma omp section
        {
          double* d_J1;
          lda = nO*nV;
          cudaMalloc((void **)&d_J1, nO*nV*nV*nO * sizeof(double));
          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));

          for (int 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) {
              cublasSetStream(handle, stream[b]);
              A = &(d_J1[nO*nV*(b+nV*i)]); lda = nO;
              B = &(d_K1[nO*nV*(i+nO*b)]); ldb = nO;
              C = &(d_X_ovvo[nO*(b+nV*nV*i)]); ldc = nO*nV;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_J1);

          double* d_Y_voov;
          cudaMalloc((void **)&d_Y_voov, nV*nO*nO*nV * sizeof(double));

          for (int 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) {
              cublasSetStream(handle, stream[g]);
              A = &(d_t2[nO*(v+nO*nV*g)]); lda = nO*nO;
              B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
              C = &(d_Y_voov[nV*nO*(v+nO*g)]); ldc = nV;
              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) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Z_ovov;
          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV;  n=nO*nV;  k=nV*nO;
          A=d_X_ovvo;  lda=nO*nV;
          B=d_Y_voov;  ldb=nV*nO;
          C=d_Z_ovov;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
          cudaFree(d_X_ovvo);
          cudaFree(d_Y_voov);

          for (int 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) {
              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) {
              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_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_Z_ovov);

        }

        #pragma omp section
        {
          double* d_X_ovov;
          cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));

          double* d_Y_ovov;
          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));

          for (int 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) {
              cublasSetStream(handle, stream[b]);
              A = &(d_K1[nO*(a+nV*nO*b)]); lda = nO*nV;
              B = &(d_K1[nO*(a+nV*nO*b)]); ldb = nO*nV;
              C = &(d_X_ovov[nO*(a+nV*nO*b)]); ldc = nO*nV;
              cublasDgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          alpha = 1.0;
          for (int v=0 ; v<nO ; ++v) {
            for (int g=0 ; g<nV ; ++g) {
              cublasSetStream(handle, stream[g]);
              A = &(d_t2[nO*(v+nO*g)]); lda = nO*nO*nV;
              B = &(d_t2[nO*(v+nO*g)]); ldb = nO*nO*nV;
              C = &(d_Y_ovov[nO*nV*(v+nO*g)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Z_ovov;
          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV;  n=nO*nV;  k=nO*nV;
          A=d_X_ovov;  lda=nO*nV;
          B=d_Y_ovov;  ldb=nO*nV;
          C=d_Z_ovov;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_ovov);
          cudaFree(d_Y_ovov);

          for (int 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) {
              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) {
              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_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_Z_ovov);

        }

        #pragma omp section
        {
          double* d_X_ovov;
          cudaMalloc((void **)&d_X_ovov, nO*nV*nO*nV * sizeof(double));

          double* d_Y_ovov;
          cudaMalloc((void **)&d_Y_ovov, nO*nV*nO*nV * sizeof(double));


          for (int 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) {
              cublasSetStream(handle, stream[g]);
              A = &(d_K1[nO*(a+nV*nO*g)]); lda = nO*nV;
              B = &(d_K1[nO*(a+nV*nO*g)]); ldb = nO*nV;
              C = &(d_X_ovov[nO*(g+nV*nO*a)]); ldc = nO*nV;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          alpha = 1.0;
          for (int v=0 ; v<nO ; ++v) {
            for (int b=0 ; b<nV ; ++b) {
              cublasSetStream(handle, stream[b]);
              A = &(d_t2[nO*(v+nO*b)]); lda = nO*nO*nV;
              B = &(d_t2[nO*(v+nO*b)]); ldb = nO*nO*nV;
              C = &(d_Y_ovov[nO*nV*(v+nO*b)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nO, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          double* d_Z_ovov;
          cudaMalloc((void **)&d_Z_ovov, nO*nV*nO*nV * sizeof(double));

          alpha = 1.0;
          beta  = 0.0;
          m=nO*nV;  n=nO*nV;  k=nO*nV;
          A=d_X_ovov;  lda=nO*nV;
          B=d_Y_ovov;  ldb=nO*nV;
          C=d_Z_ovov;  ldc=nO*nV;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_ovov);
          cudaFree(d_Y_ovov);

          for (int 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) {
              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) {
              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_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);
          cudaFree(d_Z_ovov);

        }
      } // end sections

      cudaFree(d_K1);

      double* d_tmp_cc;
      lda = cholesky_mo_num * nV;
      cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(double));

      alpha=1.0; beta=0.0;
      m=cholesky_mo_num*nV;  n=nV;  k=nO;
      A = d_cc_space_v_vo_chol; B = d_t1; C = d_tmp_cc;
      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));

      double* d_B1;
      cudaMalloc((void**)&d_B1, nV*nV*BLOCK_SIZE*sizeof(double));

      double* d_tmpB1;
      cudaMalloc((void**)&d_tmpB1, nV*BLOCK_SIZE*nV*sizeof(double));

      #pragma omp for
      for (int gam=0 ; gam<nV ; ++gam)
      {
         double* d_tmp_cc_ = &(d_tmp_cc[gam*nV*cholesky_mo_num]);
         double* d_cc_space_v_vv_chol_ = &(d_cc_space_v_vv_chol[gam*nV*cholesky_mo_num]);

         alpha = 1.0;
         beta = -1.0;
         A = d_cc_space_v_vv_chol_; lda = cholesky_mo_num;
         B = d_tmp_cc_; ldb = cholesky_mo_num;
         C = d_tmp_cc2 ; ldc = cholesky_mo_num;
         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)
         {
              const int mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;

              alpha=-1.0; beta=0.0;
              m=nV*mbs;  n=nV;  k=cholesky_mo_num;

              A=&(d_tmp_cc[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
              B=d_cc_space_v_vv_chol_; ldb=cholesky_mo_num;
              C=d_tmpB1 ; ldc=nV*BLOCK_SIZE;
              cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, lda, &beta, C, ldc);

              alpha=1.0; beta=1.0;
              m=nV*mbs;  n=nV;  k=cholesky_mo_num;

              A=&(d_cc_space_v_vv_chol[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
              B=d_tmp_cc2; ldb=cholesky_mo_num;
              C=d_tmpB1 ; ldc=nV*BLOCK_SIZE;
              cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, lda, &beta, C, ldc);

              for (size_t bet=iblock ; bet<(nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE) ; ++bet)
              {
                alpha = 1.0;
                beta = 0.0;
                A = &(d_tmpB1[nV*(bet-iblock)]); lda = nV*BLOCK_SIZE;
                B = d_tmpB1; ldb = nV;
                C = &(d_B1[nV*nV*(bet-iblock)]) ; ldc = nV;
                cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, nV, nV, &alpha, A, lda, &beta, B, ldb, C, ldc);
             }

             alpha=1.0; beta=1.0;
             m=nO*nO;  n=mbs;  k=nV*nV;
             A=d_tau; lda=nO*nO;
             B=d_B1 ; ldb=nV*nV;
             C=&(d_r2[nO*nO*(iblock + nV*gam)]); ldc=nO*nO;
             cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

         }
      }
      cudaFree(d_tmpB1);
      cudaFree(d_B1);
      cudaFree(d_tmp_cc2);
      cudaFree(d_tmp_cc);

      double * r2_tmp = malloc(nO*nO*nV*nV*sizeof(double));
      lda=nO*nO;
      cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, lda, r2_tmp, lda);
      #pragma omp critical
      {
         for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
           r2[i] -= r2_tmp[i];
         }
      }
      free(r2_tmp);

      cudaFree(d_r2);

      cublasDestroy(handle);
  }
  free(K1);
  free(J1);

  *max_r2 = 0.;
  for (size_t i=0 ; i<(size_t) nO*nO*nV*nV ; ++i) {
    const double x = r2[i] > 0. ? r2[i] : -r2[i];
    *max_r2 = *max_r2 > x ? *max_r2 : x;
  }

}






void compute_r1_space_chol_gpu(gpu_data* data, int nO, int nV, double* t1, double* r1, double* max_r1)
{
    const int cholesky_mo_num = data->cholesky_mo_num;

    int ngpus = 1;
    if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);

    #pragma omp parallel num_threads(ngpus)
    {
      int m,n,k, lda, ldb, ldc;
      double alpha, beta;
      double* A;
      double* B;
      double* C;
      cudaStream_t stream[nV];

      int igpu = omp_get_thread_num();
      cudaSetDevice(igpu);

      cublasHandle_t handle;
      cublasCreate(&handle);

      double* d_r1;
      lda = nO ;
      cudaMalloc((void **)&d_r1, lda * nV * sizeof(double));
      cudaMemset(d_r1, 0, nO*nV*sizeof(double));
      memset(r1, 0, nO*nV*sizeof(double));

      double* d_cc_space_v_vo_chol = data[igpu].cc_space_v_vo_chol;
      double* d_cc_space_v_vv_chol = data[igpu].cc_space_v_vv_chol;
      double* d_cc_space_v_oovo    = data[igpu].cc_space_v_oovo;
      double* d_cc_space_v_ovov    = data[igpu].cc_space_v_ovov;
      double* d_cc_space_v_voov    = data[igpu].cc_space_v_voov;
      double* d_cc_space_f_ov      = data[igpu].cc_space_f_ov;
      double* d_cc_space_f_vo      = data[igpu].cc_space_f_vo;
      double* d_tau                = data[igpu].tau;
      double* d_t1                 = data[igpu].t1;
      double* d_t2                 = data[igpu].t2;
      double* d_H_oo               = data[igpu].H_oo;
      double* d_H_vo               = data[igpu].H_vo;
      double* d_H_vv               = data[igpu].H_vv;

      #pragma omp sections
      {

        #pragma omp section
        {
          cublasDcopy(handle, nO*nV, d_cc_space_f_ov, 1, d_r1, 1);

          double* d_X_oo;
          cudaMalloc((void **)&d_X_oo, nO*nO * sizeof(double));

          alpha = -2.0;
          beta  = 0.0;
          m=nO;  n=nO;  k=nV;
          A=d_t1;  lda=nO;
          B=d_cc_space_f_vo;  ldb=nV;
          C=d_X_oo;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          alpha = 1.0;
          beta  = 1.0;
          m=nO;  n=nV;  k=nO;
          A=d_X_oo;  lda=nO;
          B=d_t1;  ldb=nO;
          C=d_r1;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_oo);
        }

        #pragma omp section
        {
          alpha = 1.0;
          beta  = 1.0;
          m=nO;  n=nV;  k=nV;
          A=d_t1;  lda=nO;
          B=d_H_vv;  ldb=nV;
          C=d_r1;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
        }

        #pragma omp section
        {
          alpha = -1.0;
          beta  = 1.0;
          m=nO;  n=nV;  k=nO;
          A=d_H_oo;  lda=nO;
          B=d_t1;  ldb=nO;
          C=d_r1;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
        }

        #pragma omp section
        {
          double* d_X_voov;
          cudaMalloc((void **)&d_X_voov, nV* nO* nO* nV * sizeof(double));

          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          alpha = -1.0;
          for (int i=0 ; i<nO ; ++i) {
            for (int bet=0 ; bet<nV ; ++bet) {
              cublasSetStream(handle, stream[bet]);
              beta = t1[i+bet*nO];
              A = &(d_t2[nO*(i+nO*nV*bet)]); lda = nO*nO;
              B = &(d_t1[0]); ldb = nO;
              C = &(d_X_voov[nV*(i+nO*nO*bet)]); ldc = nV*nO;
              cublasDgeam(handle, CUBLAS_OP_T, CUBLAS_OP_T, nV, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          cudaDeviceSynchronize();
          alpha = 1.0;
          beta  = 2.0;
          for (int bet=0 ; bet<nV ; ++bet) {
            cublasSetStream(handle, stream[bet]);
            A = &(d_X_voov[nV*nO*nO*bet]); lda = nV;
            B = &(d_t2[nO*nO*nV*bet]); ldb = nO*nO;
            C = A ; ldc = lda;
            cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nV, nO*nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          alpha = 1.0;
          beta  = 1.0;
          m=nV*nO;  n=nO*nV;
          A=d_X_voov;  lda=nV * nO;
          B=d_H_vo;  ldb=1;
          C=d_r1;  ldc=1;
          cublasDgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_voov);
        }

        #pragma omp section
        {
          double* d_X_ovov;
          cudaMalloc((void **)&d_X_ovov, nO* nV* nO* nV * sizeof(double));
          cublasDcopy(handle, nO*nV*nO*nV, d_cc_space_v_ovov, 1, d_X_ovov, 1);

          for (int 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) {
              cublasSetStream(handle, stream[bet]);
              A = &(d_X_ovov[nO*nV*(u+nO*bet)]); lda = nO;
              B = &(d_cc_space_v_voov[(nV*(u+nO*nO*bet))]); ldb = nV*nO;
              C = A ; ldc = lda;
              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) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          alpha = 1.0;
          beta  = 1.0;
          m=nO*nV;  n=nO*nV;
          A=d_X_ovov;  lda=nO * nV;
          B=d_t1;  ldb=1;
          C=d_r1;  ldc=1;
          cublasDgemv(handle, CUBLAS_OP_T, m, n, &alpha, A, lda, B, ldb, &beta, C, ldc);

          cudaFree(d_X_ovov);
        }


        #pragma omp section
        {
          double* d_T_vvoo;
          cudaMalloc((void **)&d_T_vvoo, nV*nV*nO*nO * sizeof(double));

          alpha = 0.0;
          beta  = 1.0;
          A = d_T_vvoo; lda = nV*nV;
          B = d_tau; ldb = nO*nO;
          C = A ; ldc = lda;
          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));

          double* d_W_vvov_tmp;
          cudaMalloc((void **)&d_W_vvov_tmp, nV*nO*nV*BLOCK_SIZE * sizeof(double));


          for (int iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE) {
            const int mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;

            alpha = 1.0;
            beta  = 0.0;
            m=nV*nO;  n=nV*mbs;  k=cholesky_mo_num;
            A=d_cc_space_v_vo_chol;  lda=cholesky_mo_num;
            B=&(d_cc_space_v_vv_chol[cholesky_mo_num*nV*iblock]);  ldb=cholesky_mo_num;
            C=d_W_vvov_tmp;  ldc=nV*nO;
            cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

            alpha = 2.0;
            beta = -1.0;
            int kk=0;
            for (int i=0 ; i<nV ; ++i) {
              cudaStreamCreate(&(stream[i]));
            }
            for (int i=0 ; i<nO ; ++i) {
              for (int bet=0 ; bet<mbs ; ++bet) {
                cublasSetStream(handle, stream[kk]);
                ++kk;
                if (kk >= nV) kk = 0;
                A = &(d_W_vvov_tmp[nV*(i+nO*nV*bet)]); lda = nV*nO;
                B = &(d_W_vvov_tmp[nV*(i+nO*nV*bet)]); ldb = nV*nO;
                C = &(d_W_vvov[nV*nV*(i+nO*bet)]); ldc = nV;
                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) {
              cudaStreamDestroy(stream[i]);
            }
            cublasSetStream(handle, NULL);

            alpha = 1.0;
            beta  = 1.0;
            m=nO;  n=mbs;  k=nO*nV*nV;
            A=d_T_vvoo;  lda=nV*nV*nO;
            B=d_W_vvov;  ldb=nO*nV*nV;
            C=&(d_r1[nO*iblock]);  ldc=nO;
            cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

          }
          cudaFree(d_W_vvov);
          cudaFree(d_W_vvov_tmp);
          cudaFree(d_T_vvoo);
        }

        #pragma omp section
        {
          double* d_W_oovo;
          cudaMalloc((void **)&d_W_oovo, nO*nO*nV*nO * sizeof(double));

          alpha = 2.0;
          beta  = -1.0;
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamCreate(&(stream[i]));
          }
          for (int u=0 ; u<nO ; ++u) {
            for (int a=0 ; a<nV ; ++a) {
//              cublasSetStream(handle, stream[a]);
              A = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); lda = nO;
              B = &(d_cc_space_v_oovo[nO*nO*(a+nV*u)]); ldb = nO;
              C = &(d_W_oovo[nO*nO*(a+nV*u)]); ldc = nO;
              cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
            }
          }
          for (int i=0 ; i<nV ; ++i) {
            cudaStreamDestroy(stream[i]);
          }
          cublasSetStream(handle, NULL);

          alpha = -1.0;
          beta  = 1.0;
          m=nO;  n=nV;  k=nO*nO*nV;
          A=d_W_oovo;  lda=nO * nO * nV;
          B=d_tau;  ldb=nO * nO * nV;
          C=d_r1;  ldc=nO;
          cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);

        }
      }

      double * r1_tmp = malloc(nO*nV*sizeof(double));
      lda=nO;
      cublasGetMatrix(nO, nV, sizeof(double), d_r1, lda, r1_tmp, lda);
      #pragma omp critical
      {
         for (size_t i=0 ; i<(size_t) nO*nV ; ++i) {
           r1[i] -= r1_tmp[i];
         }
      }
      free(r1_tmp);

      cudaFree(d_r1);

      cublasDestroy(handle);
  }

  *max_r1 = 0.;
  for (size_t i=0 ; i<(size_t) nO*nV ; ++i) {
    const double x = r1[i] > 0. ? r1[i] : -r1[i];
    *max_r1 = *max_r1 > x ? *max_r1 : x;
  }
}