scemama/qp_plugins_scemama
1
0
Fork 0
mirror of https://gitlab.com/scemama/qp_plugins_scemama.git synced 2024-12-22 04:13:40 +01:00
Code Issues Releases Wiki Activity

Added H_vv

Browse Source
This commit is contained in:
Anthony Scemama 2023-08-04 12:09:07 +02:00
parent ecd9531ee6
commit 98948f1947
4 changed files with 97 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
devel/ccsd_gpu/ccsd_space_orb_sub.irp.f
View File

@ -127,10 +127,10 @@ subroutine run_ccsd_space_orb
! Residue
if (do_ao_cholesky) then
call compute_H_vv_chol(nO,nV,tau_x,H_vv)
call compute_H_vo_chol(nO,nV,t1,H_vo)
call gpu_upload(gpu_data, nO, nV, t1, t2, tau, tau_x, H_vv);
call gpu_upload(gpu_data, nO, nV, t1, t2, tau, tau_x);
call compute_H_oo_chol_gpu(gpu_data,nO,nV,0,H_oo)
call compute_H_vv_chol_gpu(gpu_data,nO,nV,0,H_vv)
call compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
call compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2)

137
devel/ccsd_gpu/ccsd_space_orb_sub_chol.irp.f
View File

@ -290,105 +290,6 @@ subroutine compute_r1_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r1,max_r1)
end
! H_oo
subroutine compute_H_oo_chol(nO,nV,tau_x,H_oo)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: tau_x(nO, nO, nV, nV)
double precision, intent(out) :: H_oo(nO, nO)
integer :: a,b,i,j,u,k
double precision, allocatable :: tau_kau(:,:,:), tmp_vov(:,:,:)
allocate(tau_kau(cholesky_mo_num,nV,nO))
!$omp parallel &
!$omp default(shared) &
!$omp private(i,u,j,k,a,b,tmp_vov)
allocate(tmp_vov(nV,nO,nV) )
!$omp do
do u = 1, nO
do b=1,nV
do j=1,nO
do a=1,nV
tmp_vov(a,j,b) = tau_x(u,j,a,b)
enddo
enddo
enddo
call dgemm('N','T',cholesky_mo_num,nV,nO*nV,1.d0, &
cc_space_v_ov_chol, cholesky_mo_num, tmp_vov, nV, &
0.d0, tau_kau(1,1,u), cholesky_mo_num)
enddo
!$omp end do nowait
deallocate(tmp_vov)
!$omp do
do i = 1, nO
do u = 1, nO
H_oo(u,i) = cc_space_f_oo(u,i)
enddo
enddo
!$omp end do nowait
!$omp barrier
!$omp end parallel
call dgemm('T', 'N', nO, nO, cholesky_mo_num*nV, 1.d0, &
tau_kau, cholesky_mo_num*nV, cc_space_v_vo_chol, cholesky_mo_num*nV, &
1.d0, H_oo, nO)
end
! H_vv
subroutine compute_H_vv_chol(nO,nV,tau_x,H_vv)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: tau_x(nO, nO, nV, nV)
double precision, intent(out) :: H_vv(nV, nV)
integer :: a,b,i,j,u,k, beta
double precision, allocatable :: tau_kia(:,:,:), tmp_oov(:,:,:)
allocate(tau_kia(cholesky_mo_num,nO,nV))
!$omp parallel &
!$omp default(shared) &
!$omp private(i,beta,j,k,a,b,tmp_oov)
allocate(tmp_oov(nO,nO,nV) )
!$omp do
do a = 1, nV
do b=1,nV
do j=1,nO
do i=1,nO
tmp_oov(i,j,b) = tau_x(i,j,a,b)
enddo
enddo
enddo
call dgemm('N','T',cholesky_mo_num,nO,nO*nV,1.d0, &
cc_space_v_ov_chol, cholesky_mo_num, tmp_oov, nO, &
0.d0, tau_kia(1,1,a), cholesky_mo_num)
enddo
!$omp end do nowait
deallocate(tmp_oov)
!$omp do
do beta = 1, nV
do a = 1, nV
H_vv(a,beta) = cc_space_f_vv(a,beta)
enddo
enddo
!$omp end do nowait
!$omp barrier
!$omp end parallel
call dgemm('T', 'N', nV, nV, cholesky_mo_num*nO, -1.d0, &
tau_kia, cholesky_mo_num*nO, cc_space_v_ov_chol, cholesky_mo_num*nO, &
1.d0, H_vv, nV)
end
! H_vo
subroutine compute_H_vo_chol(nO,nV,t1,H_vo)
@ -439,41 +340,3 @@ subroutine compute_H_vo_chol(nO,nV,t1,H_vo)
end
subroutine compute_H_oo_chol2(nO,nV,tau_x,H_oo)
implicit none
integer, intent(in) :: nO,nV
double precision, intent(in) :: tau_x(nO, nO, nV, nV)
double precision, intent(out) :: H_oo(nO, nO)
integer :: a,b,i,j,u,k
double precision, allocatable :: tau_kau(:,:,:), tmp_vov(:,:,:), tmp_ovv(:,:,:)
allocate(tau_kau(cholesky_mo_num,nV,nO))
allocate(tmp_vov(nV,nO,nV) )
allocate(tmp_ovv(nO,nV,nV) )
do u = 1, nO
call dcopy(nO*nV*nV, tau_x(u,1,1,1), nO, tmp_ovv, 1)
print *, u
print *, tmp_ovv
do b=1,nV
do j=1,nO
do a=1,nV
tmp_vov(a,j,b) = tmp_ovv(j,a,b)
enddo
enddo
enddo
call dgemm('N','T',cholesky_mo_num,nV,nO*nV,1.d0, &
cc_space_v_ov_chol, cholesky_mo_num, tmp_vov, nV, &
0.d0, tau_kau(1,1,u), cholesky_mo_num)
enddo
deallocate(tmp_vov)
call dcopy(nO*nO, cc_space_f_oo, 1, H_oo, 1);
call dgemm('T', 'N', nO, nO, cholesky_mo_num*nV, 1.d0, &
tau_kau, cholesky_mo_num*nV, cc_space_v_vo_chol, cholesky_mo_num*nV, &
1.d0, H_oo, nO)
end

93
devel/ccsd_gpu/gpu.c
View File

@ -11,8 +11,7 @@ void gpu_upload(gpu_data* data,
double* t1,
double* t2,
double* tau,
double* tau_x,
double* H_vv)
double* tau_x)
{
int lda;
const int cholesky_mo_num = data->cholesky_mo_num;
@ -40,10 +39,6 @@ void gpu_upload(gpu_data* data,
double* d_t2 = data[igpu].t2;
lda = nO*nO;
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
double* d_H_vv = data[igpu].H_vv;
lda = nV;
cublasSetMatrix(nV, nV, sizeof(double), H_vv, lda, d_H_vv, lda);
}
}
@ -1353,6 +1348,7 @@ void compute_h_oo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_o
{
int ngpus = 1;
if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
igpu = igpu % ngpus;
const int cholesky_mo_num = data[igpu].cholesky_mo_num;
cudaSetDevice(igpu);
@ -1437,3 +1433,88 @@ void compute_h_oo_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_o
cublasDestroy(handle);
}
void compute_h_vv_chol_gpu(gpu_data* data, int nO, int nV, int igpu, double* H_vv)
{
int ngpus = 1;
if (MULTIGPU == 1) cudaGetDeviceCount(&ngpus);
igpu = igpu % ngpus;
const int cholesky_mo_num = data[igpu].cholesky_mo_num;
cudaSetDevice(igpu);
int m,n,k, lda, ldb, ldc;
double alpha, beta;
double* A;
double* B;
double* C;
cudaStream_t stream[nV];
cublasHandle_t handle;
cublasCreate(&handle);
double* d_H_vv = data[igpu].H_vv;
double* d_tau_x = data[igpu].tau_x;
double* d_cc_space_f_vv = data[igpu].cc_space_f_vv;
double* d_cc_space_v_ov_chol = data[igpu].cc_space_v_ov_chol;
double* d_tau_kia;
cudaMalloc((void **)&d_tau_kia, cholesky_mo_num*nO*nV * sizeof(double));
double* d_tmp_oov;
cudaMalloc((void **)&d_tmp_oov, nO*nO*nV * sizeof(double));
alpha = 1.0;
beta = 0.0;
// for (int i=0 ; i<nV ; ++i) {
// cudaStreamCreate(&(stream[i]));
// }
for (int a=0 ; a<nV ; ++a) {
for (int b=0 ; b<nV ; ++b) {
// cublasSetStream(handle, stream[b]);
cublasDcopy(handle, nO*nO, &(d_tau_x[nO*nO*(a+nV*b)]), 1, &(d_tmp_oov[nO*nO*b]), 1);
}
// cudaDeviceSynchronize();
// cublasSetStream(handle, NULL);
alpha = 1.0;
beta = 0.0;
m=cholesky_mo_num; n=nO; k=nO*nV;
A=d_cc_space_v_ov_chol; lda=cholesky_mo_num;
B=d_tmp_oov; ldb=nO;
C=&(d_tau_kia[cholesky_mo_num*nO*a]); ldc=cholesky_mo_num;
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
}
// for (int i=0 ; i<nV ; ++i) {
// cudaStreamDestroy(stream[i]);
// }
cudaFree(d_tmp_oov);
cublasDcopy(handle, nV*nV, d_cc_space_f_vv, 1, d_H_vv, 1);
alpha = -1.0;
beta = 1.0;
m=nV; n=nV; k=cholesky_mo_num*nO;
A=d_tau_kia; lda=cholesky_mo_num*nO;
B=d_cc_space_v_ov_chol; ldb=cholesky_mo_num*nO;
C=d_H_vv; ldc=nV;
cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, m, n, k, &alpha, A, lda, B, ldb, &beta, C, ldc);
cudaFree(d_tau_kia);
// double* H_vv = malloc(nO*nO*sizeof(double));
cublasGetMatrix(nV, nV, sizeof(double), d_H_vv, nV, H_vv, nV);
for (int i=0 ; i<ngpus ; ++i) {
if (i != igpu) {
double* d_H_vv = data[i].H_vv;
cudaSetDevice(i);
cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_H_vv, nV);
}
}
// free(H_vv);
cublasDestroy(handle);
}

10
devel/ccsd_gpu/gpu_module.f90
View File

@ -28,7 +28,7 @@ module gpu_module
real(c_double), intent(in) :: cc_space_f_vv(nV,nV)
end function
subroutine gpu_upload(gpu_data, nO, nV, t1, t2, tau, tau_x, H_vv) bind(C)
subroutine gpu_upload(gpu_data, nO, nV, t1, t2, tau, tau_x) bind(C)
import c_int, c_double, c_ptr
type(c_ptr), value :: gpu_data
integer(c_int), intent(in), value :: nO, nV
@ -36,7 +36,6 @@ module gpu_module
real(c_double), intent(in) :: t2(nO,nO,nV,nV)
real(c_double), intent(in) :: tau(nO,nO,nV,nV)
real(c_double), intent(in) :: tau_x(nO,nO,nV,nV)
real(c_double), intent(in) :: H_vv(nV,nV)
end subroutine
subroutine compute_H_oo_chol_gpu(gpu_data, nO, nV, igpu, H_oo) bind(C)
@ -46,6 +45,13 @@ module gpu_module
real(c_double), intent(out) :: H_oo(nO,nO)
end subroutine
subroutine compute_H_vv_chol_gpu(gpu_data, nO, nV, igpu, H_vv) bind(C)
import c_int, c_double, c_ptr
type(c_ptr), value :: gpu_data
integer(c_int), intent(in), value :: nO, nV, igpu
real(c_double), intent(out) :: H_vv(nO,nO)
end subroutine
subroutine compute_r2_space_chol_gpu(gpu_data, nO, nV, t1, r2, max_r2) bind(C)
import c_int, c_double, c_ptr
type(c_ptr), value :: gpu_data