mirror of
https://gitlab.com/scemama/qp_plugins_scemama.git
synced 2024-12-22 20:34:05 +01:00
Compare commits
8 Commits
0c6e6a1ca0
...
a07c149795
Author | SHA1 | Date | |
---|---|---|---|
a07c149795 | |||
0ff20e5992 | |||
c45db49df5 | |||
2df6c19772 | |||
389b217f8a | |||
c57f940eb3 | |||
ed440c16a2 | |||
5cec1b8a0c |
@ -456,159 +456,29 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
|
|||||||
! internal
|
! internal
|
||||||
integer :: u,v,i,j,beta,gam,a,b
|
integer :: u,v,i,j,beta,gam,a,b
|
||||||
double precision :: max_r2_local
|
double precision :: max_r2_local
|
||||||
|
double precision, allocatable :: Y_oovv(:,:,:,:)
|
||||||
call set_multiple_levels_omp(.False.)
|
|
||||||
|
|
||||||
!$omp parallel &
|
|
||||||
!$omp shared(nO,nV,r2,cc_space_v_oovv) &
|
|
||||||
!$omp private(u,v,gam,beta) &
|
|
||||||
!$omp default(none)
|
|
||||||
!$omp do
|
|
||||||
do gam = 1, nV
|
|
||||||
do beta = 1, nV
|
|
||||||
do v = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
r2(u,v,beta,gam) = cc_space_v_oovv(u,v,beta,gam)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do
|
|
||||||
!$omp end parallel
|
|
||||||
|
|
||||||
double precision, allocatable :: A1(:,:,:,:)
|
|
||||||
allocate(A1(nO,nO,nO,nO))
|
|
||||||
call compute_A1_chol(nO,nV,t1,t2,tau,A1)
|
|
||||||
call dgemm('N','N',nO*nO,nV*nV,nO*nO, &
|
|
||||||
1d0, A1, size(A1,1) * size(A1,2), &
|
|
||||||
tau, size(tau,1) * size(tau,2), &
|
|
||||||
1d0, r2, size(r2,1) * size(r2,2))
|
|
||||||
|
|
||||||
deallocate(A1)
|
|
||||||
|
|
||||||
integer :: block_size, iblock, k
|
integer :: block_size, iblock, k
|
||||||
block_size = 16
|
block_size = 16
|
||||||
double precision, dimension(:,:,:), allocatable :: B1, tmp_cc, tmpB1
|
call set_multiple_levels_omp(.False.)
|
||||||
double precision, dimension(:,:), allocatable :: tmp_cc2
|
|
||||||
|
|
||||||
call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,tau,cc_space_v_vo_chol, &
|
call compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,tau, &
|
||||||
cc_space_v_vv_chol, r2)
|
cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
|
||||||
|
cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
|
||||||
|
cc_space_f_vo, H_vv, r2)
|
||||||
|
|
||||||
! allocate(tmp_cc(cholesky_mo_num,nV,nV))
|
! call compute_g_vir_chol_gpu(nO,nV,cholesky_mo_num,t1,t2,t2, &
|
||||||
! call gemm0(nO, nV, cholesky_mo_num, cc_space_v_vo_chol, t1, tmp_cc)
|
! cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
|
||||||
!
|
! cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
|
||||||
!! call dgemm('N','N', cholesky_mo_num*nV, nV, nO, 1.d0, &
|
! cc_space_f_vo, H_vv, r2)
|
||||||
!! cc_space_v_vo_chol, cholesky_mo_num*nV, t1, nO, 0.d0, tmp_cc, cholesky_mo_num*nV)
|
|
||||||
!
|
|
||||||
! call set_multiple_levels_omp(.False.)
|
|
||||||
!
|
|
||||||
! !$OMP PARALLEL PRIVATE(gam, iblock, B1, tmpB1, tmp_cc2, beta, b, a)
|
|
||||||
! allocate(B1(nV,nV,block_size), tmpB1(nV,block_size,nV), tmp_cc2(cholesky_mo_num,nV))
|
|
||||||
! !$OMP DO
|
|
||||||
! do gam = 1, nV
|
|
||||||
!!
|
|
||||||
! do a=1,nV
|
|
||||||
! do k=1,cholesky_mo_num
|
|
||||||
! tmp_cc2(k,a) = cc_space_v_vv_chol(k,a,gam) - tmp_cc(k,a,gam)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!
|
|
||||||
! do iblock = 1, nV, block_size
|
|
||||||
!
|
|
||||||
! call gemm1(iblock-1, nV, cholesky_mo_num, tmp_cc, cc_space_v_vv_chol(1,1,gam), tmpB1)
|
|
||||||
!
|
|
||||||
!! call dgemm('T', 'N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, &
|
|
||||||
!! -1.d0, tmp_cc(1,1,iblock), cholesky_mo_num, &
|
|
||||||
!! cc_space_v_vv_chol(1,1,gam), cholesky_mo_num, &
|
|
||||||
!! 0.d0, tmpB1, nV*block_size)
|
|
||||||
!
|
|
||||||
! call gemm2(iblock-1, nV, cholesky_mo_num, tmp_cc2, cc_space_v_vv_chol, tmpB1)
|
|
||||||
!
|
|
||||||
!! call dgemm('T','N', nV*min(block_size, nV-iblock+1), nV, cholesky_mo_num, 1.d0, &
|
|
||||||
!! cc_space_v_vv_chol(1,1,iblock), cholesky_mo_num, &
|
|
||||||
!! tmp_cc2, cholesky_mo_num, &
|
|
||||||
!! 1.d0, tmpB1, nV*block_size)
|
|
||||||
!
|
|
||||||
! do beta = iblock, min(nV, iblock+block_size-1)
|
|
||||||
! do b = 1, nV
|
|
||||||
! do a = 1, nV
|
|
||||||
! B1(a,b,beta-iblock+1) = tmpB1(a,beta-iblock+1,b)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!
|
|
||||||
! call gemm3(iblock-1, nO, nV, gam-1, tau, B1, r2)
|
|
||||||
!
|
|
||||||
!! call dgemm('N','N',nO*nO,min(block_size, nV-iblock+1),nV*nV, &
|
|
||||||
!! 1d0, tau, nO*nO, &
|
|
||||||
!! B1 , nV*nV, &
|
|
||||||
!! 1d0, r2(1,1,iblock,gam), nO*nO)
|
|
||||||
! enddo
|
|
||||||
!
|
|
||||||
! enddo
|
|
||||||
! !$OMP ENDDO
|
|
||||||
!
|
|
||||||
! deallocate(B1, tmpB1, tmp_cc2)
|
|
||||||
! !$OMP END PARALLEL
|
|
||||||
!
|
|
||||||
! deallocate(tmp_cc)
|
|
||||||
|
|
||||||
|
|
||||||
double precision, allocatable :: X_oovv(:,:,:,:)
|
!---
|
||||||
allocate(X_oovv(nO,nO,nV,nV))
|
|
||||||
!$omp parallel &
|
|
||||||
!$omp shared(nO,nV,t2,X_oovv) &
|
|
||||||
!$omp private(u,v,gam,a) &
|
|
||||||
!$omp default(none)
|
|
||||||
!$omp do
|
|
||||||
do a = 1, nV
|
|
||||||
do gam = 1, nV
|
|
||||||
do v = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
X_oovv(u,v,gam,a) = t2(u,v,gam,a)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do
|
|
||||||
!$omp end parallel
|
|
||||||
|
|
||||||
double precision, allocatable :: g_vir(:,:)
|
|
||||||
allocate(g_vir(nV,nV))
|
|
||||||
call compute_g_vir_chol(nO,nV,t1,t2,H_vv,g_vir)
|
|
||||||
|
|
||||||
double precision, allocatable :: Y_oovv(:,:,:,:)
|
|
||||||
allocate(Y_oovv(nO,nO,nV,nV))
|
|
||||||
|
|
||||||
call dgemm('N','N',nO*nO*nV,nV,nV, &
|
|
||||||
1d0, X_oovv, size(X_oovv,1) * size(X_oovv,2) * size(X_oovv,3), &
|
|
||||||
g_vir, size(g_vir,1), &
|
|
||||||
0d0, Y_oovv, size(Y_oovv,1) * size(Y_oovv,2) * size(Y_oovv,3))
|
|
||||||
deallocate(g_vir)
|
|
||||||
deallocate(X_oovv)
|
|
||||||
|
|
||||||
!$omp parallel &
|
|
||||||
!$omp shared(nO,nV,r2,Y_oovv) &
|
|
||||||
!$omp private(u,v,gam,beta) &
|
|
||||||
!$omp default(none)
|
|
||||||
!$omp do
|
|
||||||
do gam = 1, nV
|
|
||||||
do beta = 1, nV
|
|
||||||
do v = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
r2(u,v,beta,gam) = r2(u,v,beta,gam) + Y_oovv(u,v,beta,gam) + Y_oovv(v,u,gam,beta)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do
|
|
||||||
!$omp end parallel
|
|
||||||
deallocate(Y_oovv)
|
|
||||||
|
|
||||||
double precision, allocatable :: g_occ(:,:)
|
double precision, allocatable :: g_occ(:,:)
|
||||||
allocate(g_occ(nO,nO))
|
allocate(g_occ(nO,nO))
|
||||||
call compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
|
call compute_g_occ_chol(nO,nV,t1,t2,H_oo,g_occ)
|
||||||
|
|
||||||
|
double precision, allocatable :: X_oovv(:,:,:,:)
|
||||||
allocate(X_oovv(nO,nO,nV,nV))
|
allocate(X_oovv(nO,nO,nV,nV))
|
||||||
call dgemm('N','N',nO,nO*nV*nV,nO, &
|
call dgemm('N','N',nO,nO*nV*nV,nO, &
|
||||||
1d0, g_occ , size(g_occ,1), &
|
1d0, g_occ , size(g_occ,1), &
|
||||||
@ -1021,92 +891,6 @@ subroutine compute_r2_space_chol(nO,nV,t1,t2,tau,H_oo,H_vv,H_vo,r2,max_r2)
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
! A1
|
|
||||||
|
|
||||||
subroutine compute_A1_chol(nO,nV,t1,t2,tau,A1)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
integer, intent(in) :: nO,nV
|
|
||||||
double precision, intent(in) :: t1(nO, nV)
|
|
||||||
double precision, intent(in) :: t2(nO, nO, nV, nV)
|
|
||||||
double precision, intent(in) :: tau(nO, nO, nV, nV)
|
|
||||||
double precision, intent(out) :: A1(nO, nO, nO, nO)
|
|
||||||
|
|
||||||
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta
|
|
||||||
|
|
||||||
double precision, allocatable :: X_vooo(:,:,:,:), Y_oooo(:,:,:,:)
|
|
||||||
allocate(X_vooo(nV,nO,nO,nO), Y_oooo(nO,nO,nO,nO))
|
|
||||||
|
|
||||||
! A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
|
|
||||||
!$omp parallel &
|
|
||||||
!$omp shared(nO,nV,A1,cc_space_v_oooo,cc_space_v_ovoo,X_vooo) &
|
|
||||||
!$omp private(u,v,i,j) &
|
|
||||||
!$omp default(none)
|
|
||||||
!$omp do collapse(2)
|
|
||||||
do j = 1, nO
|
|
||||||
do i = 1, nO
|
|
||||||
do v = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
A1(u,v,i,j) = cc_space_v_oooo(u,v,i,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do nowait
|
|
||||||
|
|
||||||
! A1(u,v,i,j) += cc_space_v_ovoo(u,a,i,j) * t1(v,a) &
|
|
||||||
|
|
||||||
!$omp do collapse(2)
|
|
||||||
do j = 1, nO
|
|
||||||
do i = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
do a = 1, nV
|
|
||||||
X_vooo(a,u,i,j) = cc_space_v_ovoo(u,a,i,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do
|
|
||||||
!$omp end parallel
|
|
||||||
|
|
||||||
call dgemm('N','N', nO, nO*nO*nO, nV, &
|
|
||||||
1d0, t1 , size(t1,1), &
|
|
||||||
X_vooo, size(X_vooo,1), &
|
|
||||||
0d0, Y_oooo, size(Y_oooo,1))
|
|
||||||
|
|
||||||
!$omp parallel &
|
|
||||||
!$omp shared(nO,nV,A1,Y_oooo) &
|
|
||||||
!$omp private(u,v,i,j) &
|
|
||||||
!$omp default(none)
|
|
||||||
!$omp do collapse(2)
|
|
||||||
do j = 1, nO
|
|
||||||
do i = 1, nO
|
|
||||||
do v = 1, nO
|
|
||||||
do u = 1, nO
|
|
||||||
A1(u,v,i,j) = A1(u,v,i,j) + Y_oooo(v,u,i,j)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$omp end do
|
|
||||||
!$omp end parallel
|
|
||||||
|
|
||||||
deallocate(X_vooo,Y_oooo)
|
|
||||||
|
|
||||||
! A1(u,v,i,j) += cc_space_v_vooo(a,v,i,j) * t1(u,a)
|
|
||||||
call dgemm('N','N', nO, nO*nO*nO, nV, &
|
|
||||||
1d0, t1 , size(t1,1), &
|
|
||||||
cc_space_v_vooo, size(cc_space_v_vooo,1), &
|
|
||||||
1d0, A1 , size(A1,1))
|
|
||||||
|
|
||||||
! A1(u,v,i,j) += cc_space_v_vvoo(a,b,i,j) * tau(u,v,a,b)
|
|
||||||
call dgemm('N','N', nO*nO, nO*nO, nV*nV, &
|
|
||||||
1d0, tau , size(tau,1) * size(tau,2), &
|
|
||||||
cc_space_v_vvoo, size(cc_space_v_vvoo,1) * size(cc_space_v_vvoo,2), &
|
|
||||||
1d0, A1 , size(A1,1) * size(A1,2))
|
|
||||||
|
|
||||||
end
|
|
||||||
|
|
||||||
! g_occ
|
! g_occ
|
||||||
|
|
||||||
@ -1156,7 +940,7 @@ end
|
|||||||
! g_vir
|
! g_vir
|
||||||
|
|
||||||
subroutine compute_g_vir_chol(nO,nV,t1,t2,H_vv,g_vir)
|
subroutine compute_g_vir_chol(nO,nV,t1,t2,H_vv,g_vir)
|
||||||
|
use gpu_module
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
integer, intent(in) :: nO,nV
|
integer, intent(in) :: nO,nV
|
||||||
@ -1166,44 +950,44 @@ subroutine compute_g_vir_chol(nO,nV,t1,t2,H_vv,g_vir)
|
|||||||
|
|
||||||
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
|
integer :: a,tmp_a,b,k,l,c,d,tmp_c,tmp_d,i,j,u,v, beta, gam
|
||||||
|
|
||||||
call dgemm('N','N',nV,nV,nO, &
|
! do beta = 1, nV
|
||||||
-1d0, cc_space_f_vo , size(cc_space_f_vo,1), &
|
! do a = 1, nV
|
||||||
t1 , size(t1,1), &
|
! g_vir(a,beta) = H_vv(a,beta)
|
||||||
0d0, g_vir, size(g_vir,1))
|
! enddo
|
||||||
|
! enddo
|
||||||
|
|
||||||
|
! call dgemm('N','N',nV,nV,nO, &
|
||||||
|
! -1d0, cc_space_f_vo , size(cc_space_f_vo,1), &
|
||||||
|
! t1 , size(t1,1), &
|
||||||
|
! 1d0, g_vir, size(g_vir,1))
|
||||||
|
|
||||||
double precision, allocatable :: tmp_k(:), tmp_vo(:,:,:), tmp_vo2(:,:,:)
|
double precision, allocatable :: tmp_k(:), tmp_vo(:,:,:), tmp_vo2(:,:,:)
|
||||||
allocate(tmp_k(cholesky_mo_num))
|
! allocate(tmp_k(cholesky_mo_num))
|
||||||
call dgemm('N','N', cholesky_mo_num, 1, nO*nV, 1.d0, &
|
! call dgemm('N','N', cholesky_mo_num, 1, nO*nV, 1.d0, &
|
||||||
cc_space_v_ov_chol, cholesky_mo_num, t1, nO*nV, 0.d0, tmp_k, cholesky_mo_num)
|
! cc_space_v_ov_chol, cholesky_mo_num, t1, nO*nV, 0.d0, tmp_k, cholesky_mo_num)
|
||||||
|
!
|
||||||
|
! call dgemm('T','N', nV*nV, 1, cholesky_mo_num, 2.d0, &
|
||||||
|
! cc_space_v_vv_chol, cholesky_mo_num, tmp_k, cholesky_mo_num, 1.d0, &
|
||||||
|
! g_vir, nV*nV)
|
||||||
|
! deallocate(tmp_k)
|
||||||
|
|
||||||
call dgemm('T','N', nV*nV, 1, cholesky_mo_num, 2.d0, &
|
! allocate(tmp_vo(cholesky_mo_num,nV,nO))
|
||||||
cc_space_v_vv_chol, cholesky_mo_num, tmp_k, cholesky_mo_num, 1.d0, &
|
! call dgemm('N','T',cholesky_mo_num*nV, nO, nV, 1.d0, &
|
||||||
g_vir, nV*nV)
|
! cc_space_v_vv_chol, cholesky_mo_num*nV, t1, nO, 0.d0, tmp_vo, cholesky_mo_num*nV)
|
||||||
deallocate(tmp_k)
|
!
|
||||||
|
! allocate(tmp_vo2(cholesky_mo_num,nO,nV))
|
||||||
allocate(tmp_vo(cholesky_mo_num,nV,nO))
|
! do beta=1,nV
|
||||||
call dgemm('N','T',cholesky_mo_num*nV, nO, nV, 1.d0, &
|
! do i=1,nO
|
||||||
cc_space_v_vv_chol, cholesky_mo_num*nV, t1, nO, 0.d0, tmp_vo, cholesky_mo_num*nV)
|
! do k=1,cholesky_mo_num
|
||||||
|
! tmp_vo2(k,i,beta) = -tmp_vo(k,beta,i)
|
||||||
allocate(tmp_vo2(cholesky_mo_num,nO,nV))
|
! enddo
|
||||||
do beta=1,nV
|
! enddo
|
||||||
do i=1,nO
|
! enddo
|
||||||
do k=1,cholesky_mo_num
|
! deallocate(tmp_vo)
|
||||||
tmp_vo2(k,i,beta) = -tmp_vo(k,beta,i)
|
!
|
||||||
enddo
|
! call dgemm('T','N', nV, nV, nO*cholesky_mo_num, 1.d0, &
|
||||||
enddo
|
! cc_space_v_ov_chol, cholesky_mo_num*nO, &
|
||||||
enddo
|
! tmp_vo2, cholesky_mo_num*nO, 1.d0, g_vir, nV)
|
||||||
deallocate(tmp_vo)
|
|
||||||
|
|
||||||
do beta = 1, nV
|
|
||||||
do a = 1, nV
|
|
||||||
g_vir(a,beta) = g_vir(a,beta) + H_vv(a,beta)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call dgemm('T','N', nV, nV, nO*cholesky_mo_num, 1.d0, &
|
|
||||||
cc_space_v_ov_chol, cholesky_mo_num*nO, &
|
|
||||||
tmp_vo2, cholesky_mo_num*nO, 1.d0, g_vir, nV)
|
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -1,9 +1,10 @@
|
|||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <omp.h>
|
||||||
#include <cublas_v2.h>
|
#include <cublas_v2.h>
|
||||||
#include <cuda_runtime.h>
|
#include <cuda_runtime.h>
|
||||||
|
|
||||||
|
|
||||||
#define BLOCK_SIZE 16
|
#define BLOCK_SIZE 16
|
||||||
|
|
||||||
void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*,
|
void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*,
|
||||||
@ -11,48 +12,299 @@ 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)
|
||||||
|
{
|
||||||
|
cublasHandle_t handle;
|
||||||
|
cublasCreate(&handle);
|
||||||
|
|
||||||
|
double * d_A;
|
||||||
|
double * d_B;
|
||||||
|
double * d_C;
|
||||||
|
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);
|
||||||
|
ta = CUBLAS_OP_N;
|
||||||
|
} else {
|
||||||
|
cudaMalloc((void**)&d_A, lda*m*sizeof(double));
|
||||||
|
cublasSetMatrix(k, m, sizeof(double), A, lda, d_A, lda);
|
||||||
|
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);
|
||||||
|
tb = CUBLAS_OP_N;
|
||||||
|
} else {
|
||||||
|
cudaMalloc((void**)&d_B, ldb*k*sizeof(double));
|
||||||
|
cublasSetMatrix(n, k, sizeof(double), B, ldb, d_B, ldb);
|
||||||
|
tb = CUBLAS_OP_T;
|
||||||
|
}
|
||||||
|
|
||||||
|
cudaMalloc((void**)&d_C, ldc*n*sizeof(double));
|
||||||
|
if (beta != 0.) {
|
||||||
|
cublasSetMatrix(m, n, sizeof(double), C, ldc, d_C, ldc);
|
||||||
|
}
|
||||||
|
|
||||||
|
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);
|
||||||
|
|
||||||
|
cudaFree(d_A);
|
||||||
|
cudaFree(d_B);
|
||||||
|
cudaFree(d_C);
|
||||||
|
cublasDestroy(handle);
|
||||||
|
}
|
||||||
|
|
||||||
void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
|
void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
|
||||||
double* t1,
|
double* t1,
|
||||||
|
double* t2,
|
||||||
double* tau,
|
double* tau,
|
||||||
|
double* cc_space_v_ov_chol,
|
||||||
double* cc_space_v_vo_chol,
|
double* cc_space_v_vo_chol,
|
||||||
double* cc_space_v_vv_chol,
|
double* cc_space_v_vv_chol,
|
||||||
|
double* cc_space_v_oooo,
|
||||||
|
double* cc_space_v_vooo,
|
||||||
|
double* cc_space_v_oovv,
|
||||||
|
double* cc_space_v_vvoo,
|
||||||
|
double* cc_space_f_vo,
|
||||||
|
double* H_vv,
|
||||||
double* r2)
|
double* r2)
|
||||||
{
|
{
|
||||||
|
|
||||||
|
int ngpus;
|
||||||
|
cudaGetDeviceCount(&ngpus);
|
||||||
|
#pragma omp parallel num_threads(ngpus)
|
||||||
|
{
|
||||||
int m,n,k, lda, ldb, ldc;
|
int m,n,k, lda, ldb, ldc;
|
||||||
double alpha, beta;
|
double alpha, beta;
|
||||||
double* A;
|
double* A;
|
||||||
double* B;
|
double* B;
|
||||||
double* C;
|
double* C;
|
||||||
|
|
||||||
|
int ithread = omp_get_thread_num();
|
||||||
|
int igpu = ithread ;
|
||||||
|
|
||||||
|
cudaSetDevice(igpu);
|
||||||
cublasHandle_t handle;
|
cublasHandle_t handle;
|
||||||
|
|
||||||
cublasCreate(&handle);
|
cublasCreate(&handle);
|
||||||
|
|
||||||
double* d_tau;
|
double* d_tau;
|
||||||
|
double* d_r2;
|
||||||
|
double* d_cc_space_v_vv_chol;
|
||||||
|
double* d_cc_space_v_vo_chol;
|
||||||
|
double* d_cc_space_v_ov_chol;
|
||||||
|
double* d_t1;
|
||||||
|
double* d_t2;
|
||||||
|
double* d_tmp_cc;
|
||||||
|
|
||||||
lda = nO * nO;
|
lda = nO * nO;
|
||||||
cudaMalloc((void **)&d_tau, lda * nV * nV * sizeof(double));
|
cudaMalloc((void **)&d_tau, lda * nV * nV * sizeof(double));
|
||||||
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), tau, lda, d_tau, lda);
|
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), tau, lda, d_tau, lda);
|
||||||
|
|
||||||
double* d_r2;
|
|
||||||
lda = nO * nO;
|
lda = nO * nO;
|
||||||
cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
|
cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
|
||||||
|
memset(r2, 0, nO*nO*nV*nV*sizeof(double));
|
||||||
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
|
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
|
||||||
|
|
||||||
double* d_cc_space_v_vv_chol;
|
|
||||||
lda = cholesky_mo_num * nV;
|
lda = cholesky_mo_num * nV;
|
||||||
cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
|
cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
|
||||||
cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
|
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_vo_chol;
|
lda = cholesky_mo_num * nO;
|
||||||
|
cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
|
||||||
|
cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
|
||||||
|
|
||||||
lda = cholesky_mo_num * nV;
|
lda = cholesky_mo_num * nV;
|
||||||
cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
|
cudaMalloc((void **)&d_cc_space_v_vo_chol, lda * nO * sizeof(double));
|
||||||
cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(double), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
|
cublasSetMatrix(cholesky_mo_num*nV, nO, sizeof(double), cc_space_v_vo_chol, lda, d_cc_space_v_vo_chol, lda);
|
||||||
|
|
||||||
double* d_t1;
|
|
||||||
lda = nO;
|
lda = nO;
|
||||||
cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
|
cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
|
||||||
cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
|
cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
|
||||||
|
|
||||||
double* d_tmp_cc;
|
lda = nO*nO;
|
||||||
|
cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
|
||||||
|
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
|
||||||
|
|
||||||
|
#pragma omp sections
|
||||||
|
{
|
||||||
|
|
||||||
|
#pragma omp section
|
||||||
|
{
|
||||||
|
double* d_cc_space_v_vooo;
|
||||||
|
cudaMalloc((void**)&d_cc_space_v_vooo, nV*nO*nO*nO*sizeof(double));
|
||||||
|
cublasSetMatrix(nV*nO, nO*nO, sizeof(double), cc_space_v_vooo, nV*nO, d_cc_space_v_vooo, nV*nO);
|
||||||
|
|
||||||
|
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);
|
||||||
|
cudaFree(d_cc_space_v_vooo);
|
||||||
|
|
||||||
|
double* d_A1;
|
||||||
|
cudaMalloc((void**)&d_A1, nO*nO*nO*nO*sizeof(double));
|
||||||
|
|
||||||
|
double* d_cc_space_v_oooo;
|
||||||
|
cudaMalloc((void**)&d_cc_space_v_oooo, nO*nO*nO*nO*sizeof(double));
|
||||||
|
cublasSetMatrix(nO*nO, nO*nO, sizeof(double), cc_space_v_oooo, nO*nO, d_cc_space_v_oooo, nO*nO);
|
||||||
|
|
||||||
|
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 j=0 ; j<nO ; ++j) {
|
||||||
|
for (int i=0 ; i<nO ; ++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);
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
cudaFree(d_Y_oooo);
|
||||||
|
|
||||||
|
double* d_cc_space_v_vvoo;
|
||||||
|
cudaMalloc((void**)&d_cc_space_v_vvoo, nV*nV*nO*nO*sizeof(double));
|
||||||
|
cublasSetMatrix(nV*nV, nO*nO, sizeof(double), cc_space_v_vvoo, nV*nV, d_cc_space_v_vvoo, nV*nV);
|
||||||
|
|
||||||
|
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);
|
||||||
|
cudaFree(d_cc_space_v_vvoo);
|
||||||
|
|
||||||
|
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));
|
||||||
|
cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_g_vir, nV);
|
||||||
|
|
||||||
|
double* d_cc_space_f_vo;
|
||||||
|
cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
|
||||||
|
cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
|
||||||
|
|
||||||
|
alpha = -1.0;
|
||||||
|
beta = 1.0;
|
||||||
|
m=nV ; n=nV; k=nO;
|
||||||
|
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);
|
||||||
|
cudaFree(d_cc_space_f_vo);
|
||||||
|
|
||||||
|
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) {
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
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 j=0 ; j<nV ; ++j) {
|
||||||
|
for (int i=0 ; i<nV ; ++i) {
|
||||||
|
alpha = 1.0;
|
||||||
|
beta = 1.0;
|
||||||
|
A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
|
||||||
|
B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
|
||||||
|
C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
|
||||||
|
cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
cudaFree(d_Y_oovv);
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
lda = cholesky_mo_num * nV;
|
lda = cholesky_mo_num * nV;
|
||||||
cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(double));
|
cudaMalloc((void **)&d_tmp_cc, lda * nV * sizeof(double));
|
||||||
|
|
||||||
@ -61,8 +313,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
|
|||||||
A = d_cc_space_v_vo_chol; B = d_t1; C = d_tmp_cc;
|
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);
|
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, A, m, B, k, &beta, C, m);
|
||||||
|
|
||||||
#pragma omp parallel
|
|
||||||
{
|
|
||||||
double* d_tmp_cc2;
|
double* d_tmp_cc2;
|
||||||
cudaMalloc((void **)&d_tmp_cc2, cholesky_mo_num*nV*sizeof(double));
|
cudaMalloc((void **)&d_tmp_cc2, cholesky_mo_num*nV*sizeof(double));
|
||||||
|
|
||||||
@ -118,7 +368,6 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
|
|||||||
|
|
||||||
alpha=1.0; beta=1.0;
|
alpha=1.0; beta=1.0;
|
||||||
m=nO*nO; n=mbs; k=nV*nV;
|
m=nO*nO; n=mbs; k=nV*nV;
|
||||||
|
|
||||||
A=d_tau; lda=nO*nO;
|
A=d_tau; lda=nO*nO;
|
||||||
B=d_B1 ; ldb=nV*nV;
|
B=d_B1 ; ldb=nV*nV;
|
||||||
C=&(d_r2[nO*nO*(iblock + nV*gam)]); ldc=nO*nO;
|
C=&(d_r2[nO*nO*(iblock + nV*gam)]); ldc=nO*nO;
|
||||||
@ -126,42 +375,197 @@ void compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo
|
|||||||
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
cudaFree(d_tmpB1);
|
||||||
lda=nO*nO;
|
cudaFree(d_B1);
|
||||||
cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, lda, r2, lda);
|
cudaFree(d_tmp_cc2);
|
||||||
|
|
||||||
cudaFree(d_cc_space_v_vo_chol);
|
cudaFree(d_cc_space_v_vo_chol);
|
||||||
cudaFree(d_cc_space_v_vv_chol);
|
cudaFree(d_cc_space_v_vv_chol);
|
||||||
cudaFree(d_tau);
|
cudaFree(d_tau);
|
||||||
cudaFree(d_t1);
|
cudaFree(d_t1);
|
||||||
cudaFree(d_tmp_cc);
|
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<nO*nO*nV*nV ; ++i) {
|
||||||
|
r2[i] += r2_tmp[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
free(r2_tmp);
|
||||||
|
|
||||||
cudaFree(d_r2);
|
cudaFree(d_r2);
|
||||||
cublasDestroy(handle);
|
cublasDestroy(handle);
|
||||||
|
}
|
||||||
|
|
||||||
|
for (size_t i=0 ; i<nO*nO*nV*nV ; ++i)
|
||||||
|
{
|
||||||
|
r2[i] += cc_space_v_oovv[i];
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
/*
|
|
||||||
|
|
||||||
|
void compute_g_vir_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
|
||||||
|
double* t1,
|
||||||
|
double* t2,
|
||||||
|
double* tau,
|
||||||
|
double* cc_space_v_ov_chol,
|
||||||
|
double* cc_space_v_vo_chol,
|
||||||
|
double* cc_space_v_vv_chol,
|
||||||
|
double* cc_space_v_oooo,
|
||||||
|
double* cc_space_v_vooo,
|
||||||
|
double* cc_space_v_oovv,
|
||||||
|
double* cc_space_v_vvoo,
|
||||||
|
double* cc_space_f_vo,
|
||||||
|
double* H_vv,
|
||||||
|
double* r2)
|
||||||
|
{
|
||||||
|
int m,n,k, lda, ldb, ldc;
|
||||||
|
double alpha, beta;
|
||||||
|
double* A;
|
||||||
|
double* B;
|
||||||
|
double* C;
|
||||||
|
|
||||||
|
int ithread = omp_get_thread_num();
|
||||||
|
int igpu = ithread ;
|
||||||
|
|
||||||
|
cudaSetDevice(igpu);
|
||||||
cublasHandle_t handle;
|
cublasHandle_t handle;
|
||||||
|
|
||||||
cublasCreate(&handle);
|
cublasCreate(&handle);
|
||||||
|
|
||||||
double *d_A, *d_B, *d_C;
|
double* d_tau;
|
||||||
cudaMalloc((void **)&d_A, m * k * sizeof(double));
|
double* d_r2;
|
||||||
cudaMalloc((void **)&d_B, k * n * sizeof(double));
|
double* d_cc_space_v_vv_chol;
|
||||||
cudaMalloc((void **)&d_C, m * n * sizeof(double));
|
double* d_cc_space_v_vo_chol;
|
||||||
|
double* d_cc_space_v_ov_chol;
|
||||||
|
double* d_t1;
|
||||||
|
double* d_t2;
|
||||||
|
double* d_tmp_cc;
|
||||||
|
|
||||||
cublasSetMatrix(m, k, sizeof(double), A, m, d_A, m);
|
|
||||||
cublasSetMatrix(k, n, sizeof(double), B, k, d_B, k);
|
|
||||||
cublasSetMatrix(m, n, sizeof(double), C, m, d_C, m);
|
|
||||||
|
|
||||||
cublasDgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k, &alpha, d_A, m, d_B, k, &beta, d_C, m);
|
|
||||||
|
|
||||||
cublasGetMatrix(m, n, sizeof(double), d_C, m, C, m);
|
lda = nO * nO;
|
||||||
|
cudaMalloc((void **)&d_r2, lda * nV * nV * sizeof(double));
|
||||||
|
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), r2, lda, d_r2, lda);
|
||||||
|
|
||||||
|
lda = nO*nO;
|
||||||
|
cudaMalloc((void **)&d_t2, nO*nO*nV*nV * sizeof(double));
|
||||||
|
cublasSetMatrix(nO*nO, nV*nV, sizeof(double), t2, lda, d_t2, lda);
|
||||||
|
lda = nO;
|
||||||
|
cudaMalloc((void **)&d_t1, nO * nV * sizeof(double));
|
||||||
|
cublasSetMatrix(nO, nV, sizeof(double), t1, lda, d_t1, lda);
|
||||||
|
lda = cholesky_mo_num * nV;
|
||||||
|
cudaMalloc((void **)&d_cc_space_v_vv_chol, lda * nV * sizeof(double));
|
||||||
|
cublasSetMatrix(cholesky_mo_num*nV, nV, sizeof(double), cc_space_v_vv_chol, lda, d_cc_space_v_vv_chol, lda);
|
||||||
|
|
||||||
|
lda = cholesky_mo_num * nO;
|
||||||
|
cudaMalloc((void **)&d_cc_space_v_ov_chol, lda * nV * sizeof(double));
|
||||||
|
cublasSetMatrix(cholesky_mo_num*nO, nV, sizeof(double), cc_space_v_ov_chol, lda, d_cc_space_v_ov_chol, lda);
|
||||||
|
|
||||||
|
// ---
|
||||||
|
double* d_g_vir;
|
||||||
|
cudaMalloc((void**)&d_g_vir, nV*nV*sizeof(double));
|
||||||
|
cublasSetMatrix(nV, nV, sizeof(double), H_vv, nV, d_g_vir, nV);
|
||||||
|
|
||||||
|
double* d_cc_space_f_vo;
|
||||||
|
cudaMalloc((void**)&d_cc_space_f_vo, nV*nO*sizeof(double));
|
||||||
|
cublasSetMatrix(nV, nO, sizeof(double), cc_space_f_vo, nV, d_cc_space_f_vo, nV);
|
||||||
|
|
||||||
|
alpha = -1.0;
|
||||||
|
beta = 1.0;
|
||||||
|
m=nV ; n=nV; k=nO;
|
||||||
|
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);
|
||||||
|
cudaFree(d_cc_space_f_vo);
|
||||||
|
|
||||||
|
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) {
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
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 j=0 ; j<nV ; ++j) {
|
||||||
|
for (int i=0 ; i<nV ; ++i) {
|
||||||
|
alpha = 1.0;
|
||||||
|
beta = 1.0;
|
||||||
|
A = &(d_r2[nO*nO*(i+nV*j)]); lda = nO;
|
||||||
|
B = &(d_Y_oovv[nO*nO*(j+nV*i)]); ldb = nO;
|
||||||
|
C = &(d_r2[nO*nO*(i+nV*j)]); ldc = nO;
|
||||||
|
cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_T, nO, nO, &alpha, A, lda, &beta, B, ldb, C, ldc);
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
cudaFree(d_Y_oovv);
|
||||||
|
|
||||||
|
cublasGetMatrix(nO*nO, nV*nV, sizeof(double), d_r2, nO*nO, r2, nO*nO);
|
||||||
|
}
|
||||||
|
|
||||||
cudaFree(d_A);
|
|
||||||
cudaFree(d_B);
|
|
||||||
cudaFree(d_C);
|
|
||||||
cublasDestroy(handle);
|
|
||||||
|
|
||||||
*/
|
|
||||||
|
@ -4,48 +4,56 @@ module gpu_module
|
|||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
interface
|
interface
|
||||||
subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, &
|
subroutine compute_r2_space_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
|
||||||
t1,tau,cc_space_v_vo_chol,cc_space_v_vv_chol, r2) bind(C)
|
cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
|
||||||
|
cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
|
||||||
|
cc_space_f_vo, H_vv, r2) bind(C)
|
||||||
import c_int, c_double
|
import c_int, c_double
|
||||||
integer(c_int), value :: nO, nV, cholesky_mo_num
|
integer(c_int), value :: nO, nV, cholesky_mo_num
|
||||||
real(c_double), intent(in) :: t1(nO,nV)
|
real(c_double), intent(in) :: t1(nO,nV)
|
||||||
|
real(c_double), intent(in) :: t2(nO,nO,nV,nV)
|
||||||
real(c_double), intent(in) :: tau(nO,nO,nV,nV)
|
real(c_double), intent(in) :: tau(nO,nO,nV,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_ov_chol(cholesky_mo_num,nO,nV)
|
||||||
real(c_double), intent(in) :: cc_space_v_vo_chol(cholesky_mo_num,nV,nO)
|
real(c_double), intent(in) :: cc_space_v_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_vv_chol(cholesky_mo_num,nV,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_oooo(nO,nO,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vooo(nV,nO,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_oovv(nO,nO,nV,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vvoo(nV,nV,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_f_vo(nV,nO)
|
||||||
|
real(c_double), intent(in) :: H_vv(nV,nV)
|
||||||
real(c_double), intent(out) :: r2(nO,nO,nV,nV)
|
real(c_double), intent(out) :: r2(nO,nO,nV,nV)
|
||||||
end subroutine
|
end subroutine
|
||||||
|
|
||||||
subroutine gemm0(nO, nV, cholesky_mo_num, cc_space_v_vo_chol, t1, tmp_cc) bind(C, name="gemm0")
|
subroutine compute_g_vir_chol_gpu(nO,nV,cholesky_mo_num, t1, t2, tau,&
|
||||||
|
cc_space_v_ov_chol, cc_space_v_vo_chol, cc_space_v_vv_chol, &
|
||||||
|
cc_space_v_oooo, cc_space_v_vooo, cc_space_v_oovv, cc_space_v_vvoo, &
|
||||||
|
cc_space_f_vo, H_vv, r2) bind(C)
|
||||||
import c_int, c_double
|
import c_int, c_double
|
||||||
integer(c_int), value :: nO, nV, cholesky_mo_num
|
integer(c_int), value :: nO, nV, cholesky_mo_num
|
||||||
real(c_double) :: cc_space_v_vo_chol(cholesky_mo_num,nV,nO)
|
real(c_double), intent(in) :: t1(nO,nV)
|
||||||
real(c_double) :: t1(nO,nV)
|
real(c_double), intent(in) :: t2(nO,nO,nV,nV)
|
||||||
real(c_double) :: tmp_cc(cholesky_mo_num,nV,nV)
|
real(c_double), intent(in) :: tau(nO,nO,nV,nV)
|
||||||
end subroutine gemm0
|
real(c_double), intent(in) :: cc_space_v_ov_chol(cholesky_mo_num,nO,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vo_chol(cholesky_mo_num,nV,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vv_chol(cholesky_mo_num,nV,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_oooo(nO,nO,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vooo(nV,nO,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_oovv(nO,nO,nV,nV)
|
||||||
|
real(c_double), intent(in) :: cc_space_v_vvoo(nV,nV,nO,nO)
|
||||||
|
real(c_double), intent(in) :: cc_space_f_vo(nV,nO)
|
||||||
|
real(c_double), intent(in) :: H_vv(nV,nV)
|
||||||
|
real(c_double), intent(out) :: r2(nO,nO,nV,nV)
|
||||||
|
end subroutine
|
||||||
|
|
||||||
subroutine gemm1(iblock, nV, cholesky_mo_num, tmp_cc, cc_space_v_vv_chol_, tmpB1) bind(C, name="gemm1")
|
|
||||||
import c_int, c_double
|
|
||||||
integer(c_int), value :: iblock, nV, cholesky_mo_num
|
|
||||||
real(c_double) :: tmp_cc(cholesky_mo_num,nV,nV)
|
|
||||||
real(c_double) :: cc_space_v_vv_chol_(cholesky_mo_num,nV)
|
|
||||||
real(c_double) :: tmpB1(nV,16,nV)
|
|
||||||
end subroutine gemm1
|
|
||||||
|
|
||||||
subroutine gemm2(iblock, nV, cholesky_mo_num, tmp_cc2, cc_space_v_vv_chol, tmpB1) bind(C, name="gemm2")
|
subroutine gpu_dgemm(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc) bind(C)
|
||||||
import c_int, c_double
|
import c_int, c_double, c_char
|
||||||
integer(c_int), value :: iblock, nV, cholesky_mo_num
|
character(c_char), value :: transa, transb
|
||||||
real(c_double) :: tmp_cc2(cholesky_mo_num,nV)
|
integer(c_int), value :: m,n,k,lda,ldb,ldc
|
||||||
real(c_double) :: cc_space_v_vv_chol(cholesky_mo_num,nV,nV)
|
real(c_double), value :: alpha, beta
|
||||||
real(c_double) :: tmpB1(nV,16,nV)
|
real(c_double) :: A(lda,*), B(ldb,*), C(ldc,*)
|
||||||
end subroutine gemm2
|
end subroutine
|
||||||
|
|
||||||
subroutine gemm3(iblock, nO, nV, gam, tau, B1, r2) bind(C, name="gemm3")
|
|
||||||
import c_int, c_double
|
|
||||||
integer(c_int), value :: iblock, nO, nV, gam
|
|
||||||
real(c_double) :: tau(nO,nO,nV,nV)
|
|
||||||
real(c_double) :: B1(nV,nV,*)
|
|
||||||
real(c_double) :: r2(nO,nO,nV,nV)
|
|
||||||
end subroutine gemm3
|
|
||||||
|
|
||||||
end interface
|
end interface
|
||||||
|
|
||||||
|
1679
devel/ccsd_gpu/mo_integrals_cc.irp.f
Normal file
1679
devel/ccsd_gpu/mo_integrals_cc.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user