comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 22:25:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

still debugging

Browse source
This commit is contained in:
Jeff Hammond 2023-10-10 09:40:44 +03:00
parent 0a1b13399c
commit 6de9fbeb57
3 changed files with 320 additions and 87 deletions
Download patch file Download diff file Expand all files Collapse all files

37
src/nwchem.nw
View file

@ -1,28 +1,39 @@
#
# $Id$
#
title "Water geometry optimization in cc-pvdz basis set"
echo
start h2o
geometry autosym
memory global 4000 mb heap 100 mb stack 4000 mb
geometry autosym noprint
O 0.0 0.0 -0.02
H -0.74 0.0 -0.76
H 0.74 0.0 -0.76
end
basis
H library cc-pvdz
O library cc-pvdz
basis noprint
* library cc-pvdz
end
driver
clear
scf
thresh 1e-10
print none
end
scf; print low; end
tce
ccsd
thresh 1e-12
io ga
2eorb
tilesize 40
attilesize 40
2emet 13
end
task scf optimize
set tce:nts F
set tce:ccsd_t2_8_alg 8
task tce energy
set tce:nts T
task tce energy

5
src/tce/ccsd/ccsd_t2.F
View file

@ -324,6 +324,11 @@ c
else if (ccsd_t2_8_alg.eq.7) then
CALL ccsd_t2_8_dgemm_driver(d_t2,k_t2_offset,d_v2,k_v2_offset,
1 d_i0,k_i0_offset,maxh,maxp)
#endif
#ifdef USE_TCE_CUBLAS
else if (ccsd_t2_8_alg.eq.8) then
CALL ccsd_t2_8_cublas(d_t2,k_t2_offset,d_v2,k_v2_offset,
1 d_i0,k_i0_offset,maxh,maxp)
#endif
else
CALL errquit('ccsd_t2_8_alg',ccsd_t2_8_alg,INPUT_ERR)

365
src/tce/ccsd/ccsd_t2_8.F
View file

@ -1,8 +1,8 @@
SUBROUTINE ccsd_t2_8(d_a,k_a_offset,d_b,k_b_offset,d_c,k_c_offset)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -177,10 +177,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -189,7 +189,7 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
#include "tce.fh"
INTEGER d_a,d_b,d_c
INTEGER k_a_offset,k_b_offset,k_c_offset
INTEGER maxh,maxp,dimhhpp,dimpppp,dimtemp
INTEGER maxh,maxp,dimhhpp,dimpppp
INTEGER next,nprocs,count
INTEGER p5b,p6b,p3b,p4b,h1b,h2b
INTEGER p5b_1,p6b_1,h1b_1,h2b_1
@ -215,7 +215,6 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
dimhhpp = maxh*maxh*maxp*maxp
dimpppp = maxp*maxp*maxp*maxp
dimtemp = max(dimpppp,dimhhpp)
e_a=0
e_b=0
@ -399,10 +398,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
SUBROUTINE ccsd_t2_8_spiral(d_a,k_a_offset,d_b,k_b_offset,
1 d_c,k_c_offset)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -428,16 +427,16 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
integer nbh
double precision alpha
external NXTASK
c
c print*,'entering ccsd_t2_8_spiral (energy)'
c
!
! print*,'entering ccsd_t2_8_spiral (energy)'
!
nprocs = ga_nnodes()
count = 0
next = nxtask(nprocs, 1)
c
!
if (.not.ma_push_get(mt_dbl,tile_dim**4,'c',l_c,k_c))
1 call errquit('ccsd_t2_8',9,MA_ERR)
c
!
do p3b = noab+1,noab+nvab
dim3=int_mb(k_range+p3b-1)
spn3=int_mb(k_spin +p3b-1)
@ -446,11 +445,11 @@ c
dim4=int_mb(k_range+p4b-1)
spn4=int_mb(k_spin +p4b-1)
sym4=int_mb(k_sym +p4b-1)
c
!
dim34 = dim3 * dim4
spn34 = spn3 + spn4
sym34 = ieor(sym3,sym4)
c
!
do p5b = noab+1,noab+nvab
dim5=int_mb(k_range+p5b-1)
spn5=int_mb(k_spin +p5b-1)
@ -459,26 +458,26 @@ c
dim6=int_mb(k_range+p6b-1)
spn6=int_mb(k_spin +p6b-1)
sym6=int_mb(k_sym +p6b-1)
c
!
dim56 = dim5 * dim6
spn56 = spn5 + spn6
sym56 = ieor(sym5,sym6)
c
!
dimb = dim34 * dim56
c
!
if ( (dimb.gt.0) .and. (ieor(sym34,sym56).eq.0)
1 .and. (spn34.eq.spn56) ) then
c
!
if (next.eq.count) then
c
!
call tce_restricted_4(p3b,p4b,p5b,p6b,
1 p3b_2,p4b_2,p5b_2,p6b_2)
c
!
if (.not.ma_push_get(mt_dbl,dimb,'bs',l_bs,k_bs))
1 call errquit('ccsd_t2_8',4,MA_ERR)
if (.not.ma_push_get(mt_dbl,dimb,'b',l_b,k_b))
1 call errquit('ccsd_t2_8',5,MA_ERR)
c
!
if(.not.intorb) then
call get_hash_block(d_b,dbl_mb(k_b),dimb,
1 int_mb(k_b_offset),
@ -492,13 +491,13 @@ c
3 (noab+nvab) * (p4b_2 - 1 + (noab+nvab) *
4 (p3b_2 - 1)))),p6b_2,p5b_2,p4b_2,p3b_2)
end if
c
!
call tce_sort_4(dbl_mb(k_b),dbl_mb(k_bs),
1 dim3,dim4,dim5,dim6,2,1,4,3,1.0d0)
c
!
if (.not.ma_pop_stack(l_b))
1 call errquit('ccsd_t2_8',6,MA_ERR)
c
!
do h1b = 1,noab
dim1=int_mb(k_range+h1b-1)
spn1=int_mb(k_spin +h1b-1)
@ -507,48 +506,48 @@ c
dim2=int_mb(k_range+h2b-1)
spn2=int_mb(k_spin +h2b-1)
sym2=int_mb(k_sym +h2b-1)
c
!
dim12 = dim1 * dim2
spn12 = spn1 + spn2
sym12 = ieor(sym1,sym2)
c
!
dima = dim12 * dim56
c
!
if (dima.gt.0) then
c
!
call tce_restricted_4(p5b,p6b,h1b,h2b,
1 p5b_1,p6b_1,h1b_1,h2b_1)
c
!
if (spn34.eq.spn12) then
if (spn56.eq.spn12) then
if ((.not.restricted).or.((spn34+spn12).ne.8)) then
if (ieor(sym34,sym12).eq.0) then
if (ieor(sym56,sym12).eq.0) then
c
!
dimc = dim12 * dim34
c
!
if (.not.ma_push_get(mt_dbl,dimc,'cs',l_cs,k_cs))
1 call errquit('ccsd_t2_8',0,MA_ERR)
c
!
call dfill(dimc,0.0d0,dbl_mb(k_cs),1)
c
!
if (.not.ma_push_get(mt_dbl,dima,'as',l_as,k_as))
1 call errquit('ccsd_t2_8',1,MA_ERR)
if (.not.ma_push_get(mt_dbl,dima,'a',l_a,k_a))
1 call errquit('ccsd_t2_8',2,MA_ERR)
c
!
call get_hash_block(d_a,dbl_mb(k_a),dima,
1 int_mb(k_a_offset),
2 (h2b_1 - 1 + noab * (h1b_1 - 1 + noab *
3 (p6b_1 - noab - 1 + nvab *
4 (p5b_1 - noab - 1)))))
c
!
call tce_sort_4(dbl_mb(k_a),dbl_mb(k_as),
1 dim5,dim6,dim1,dim2,4,3,2,1,1.0d0)
c
!
if (.not.ma_pop_stack(l_a))
1 call errquit('ccsd_t2_8',3,MA_ERR)
c
!
if (p5b .eq. p6b) then
alpha = 1.0d0
else
@ -557,50 +556,50 @@ c
call dgemm('T','N',dim12,dim34,dim56,alpha,
2 dbl_mb(k_as),dim56,dbl_mb(k_bs),dim56,
3 1.0d0,dbl_mb(k_cs),dim12)
c
!
if (.not.ma_pop_stack(l_as))
1 call errquit('ccsd_t2_8',8,MA_ERR)
c
!
call ga_nbwait(nbh) ! wait until previous put of c is gone before overwriting buffer
c
!
call tce_sort_4(dbl_mb(k_cs),dbl_mb(k_c),
1 dim4,dim3,dim2,dim1,2,1,4,3,0.5d0)
c
!
call add_hash_block_nb(d_c,dbl_mb(k_c),dimc,
1 int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab*
2 (p4b-noab-1+nvab*(p3b-noab-1)))),nbh)
c
!
if (.not.ma_pop_stack(l_cs))
1 call errquit('ccsd_t2_8',11,MA_ERR)
c
!
end if
end if
end if
end if
end if
c
!
endif ! dima>0
c
!
end do
end do
c
!
if (.not.ma_pop_stack(l_bs))
1 call errquit('ccsd_t2_8',7,MA_ERR)
c
!
next = NXTASK(nprocs, 1)
end if ! next=count
count = count + 1
c
!
endif ! dimb>0
c
!
end do
end do
end do
end do
c
!
if (.not.ma_pop_stack(l_c))
1 call errquit('ccsd_t2_8',10,MA_ERR)
c
!
next = NXTASK(-nprocs, 1)
call ga_sync()
RETURN
@ -610,10 +609,10 @@ c
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -766,10 +765,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -1142,10 +1141,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -1323,10 +1322,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
C $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
! $Id: ccsd_t2.F 27404 2015-08-24 14:20:43Z jhammond $
! This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
! Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
@ -1356,3 +1355,221 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v (
END
#endif
#ifdef USE_TCE_CUBLAS
SUBROUTINE ccsd_t2_8_cublas(d_a,k_a_offset,
& d_b,k_b_offset,
& d_c,k_c_offset,
& maxh,maxp)
! i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( p3 p4 p5 p6 )_v
use iso_fortran_env
use cudafor
use cublas, only : cublasHandle, CUBLAS_OP_N, cublascreate,
& cublassetstream, cublasdestroy, cublasDgemm_v2
IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
#include "sym.fh"
#include "errquit.fh"
#include "tce.fh"
INTEGER d_a,d_b,d_c
INTEGER k_a_offset,k_b_offset,k_c_offset
INTEGER maxh,maxp,dimhhpp,dimpppp
INTEGER next,nprocs,count
INTEGER p5b,p6b,p3b,p4b,h1b,h2b
INTEGER p5b_1,p6b_1,h1b_1,h2b_1
INTEGER p3b_2,p4b_2,p5b_2,p6b_2
INTEGER dima,dimb,dimc,dim_common,dima_sort,dimb_sort
integer(kind=INT32) :: h21d, p43d, p65d
double precision, allocatable, pinned :: f_a(:), f_b(:), f_c(:)
double precision, allocatable, pinned :: f_d(:)
double precision, allocatable, device :: x_a(:), x_b(:)
double precision, allocatable, device :: x_c(:)
integer :: e_a
double precision alpha
integer p5b_in,p6b_in
INTEGER NXTASK
EXTERNAL NXTASK
integer :: err
integer(kind=cuda_stream_kind) :: stream
type(cublasHandle) :: handle
integer(INT32), parameter :: cu_op_n = CUBLAS_OP_N
err = cudaStreamCreate(stream)
if (err.ne.0) call errquit('cudaStreamCreate',err,UNKNOWN_ERR)
err = cublasCreate(handle)
if (err.ne.0) call errquit('cublasCreate',err,UNKNOWN_ERR)
err = cublasSetStream(handle, stream)
if (err.ne.0) call errquit('cublasSetStream',err,UNKNOWN_ERR)
nprocs = GA_NNODES()
count = 0
next = NXTASK(nprocs, 1)
dimhhpp = maxh*maxh*maxp*maxp
dimpppp = maxp*maxp*maxp*maxp
allocate(f_a(1:dimhhpp), f_b(1:dimpppp),
& f_c(1:dimhhpp), f_d(1:dimhhpp), stat=e_a)
if (e_a.ne.0) call errquit("MA Halloc",2*dimhhpp+1*dimpppp,MA_ERR)
allocate(x_a(1:dimhhpp), x_b(1:dimpppp),
& x_c(1:dimhhpp), stat=e_a)
if (e_a.ne.0) call errquit("MA Galloc",2*dimhhpp+1*dimpppp,MA_ERR)
DO p3b = noab+1,noab+nvab
DO p4b = p3b,noab+nvab
DO h1b = 1,noab
DO h2b = h1b,noab
IF ((.not.restricted).or.
& ( int_mb(k_spin+p3b-1)+int_mb(k_spin+p4b-1)
& +int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-1).ne.8)) THEN
IF (int_mb(k_spin+p3b-1)+int_mb(k_spin+p4b-1) .eq.
& int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-1)) THEN
IF (ieor(int_mb(k_sym+p3b-1),ieor(int_mb(k_sym+p4b-1),
& ieor(int_mb(k_sym+h1b-1),int_mb(k_sym+h2b-1))))
& .eq. ieor(irrep_v,irrep_t)) THEN
IF (next.eq.count) THEN
dima_sort = int_mb(k_range+h1b-1)
& * int_mb(k_range+h2b-1)
dimb_sort = int_mb(k_range+p3b-1)
& * int_mb(k_range+p4b-1)
dimc = int_mb(k_range+p3b-1) * int_mb(k_range+p4b-1)
& * int_mb(k_range+h1b-1) * int_mb(k_range+h2b-1)
CALL DFILL(dimc,0.0d0,f_c,1)
!print*,'cudaMemcpyAsync C (in)'
!x_c(1:dimc) = f_c(1:dimc)
!err = cudaMemcpyAsync(x_c,f_c,dimc,stream)
!if (err.ne.0) then
! call errquit('cudaMemcpyAsync',err,UNKNOWN_ERR)
!endif
CALL DFILL(dimc,0.0d0,f_d,1)
DO p5b_in =ga_nodeid(),ga_nodeid()+nvab-1
p5b=mod(p5b_in,nvab)+noab+1
DO p6b_in=ga_nodeid(),ga_nodeid()+nvab+noab-p5b
p6b=mod(p6b_in,noab+nvab-p5b+1)+p5b
IF (int_mb(k_spin+p5b-1)+int_mb(k_spin+p6b-1) .eq.
& int_mb(k_spin+h1b-1)+int_mb(k_spin+h2b-1)) THEN
IF (ieor(int_mb(k_sym+p5b-1),ieor(int_mb(k_sym+p6b-1),
& ieor(int_mb(k_sym+h1b-1),int_mb(k_sym+h2b-1))))
& .eq. irrep_t) THEN
CALL TCE_RESTRICTED_4(p5b,p6b,h1b,h2b,
& p5b_1,p6b_1,h1b_1,h2b_1)
CALL TCE_RESTRICTED_4(p3b,p4b,p5b,p6b,
& p3b_2,p4b_2,p5b_2,p6b_2)
dim_common = int_mb(k_range+p5b-1)
& * int_mb(k_range+p6b-1)
dima = dim_common * dima_sort
dimb = dim_common * dimb_sort
IF ((dima .gt. 0) .and. (dimb .gt. 0)) THEN
CALL GET_HASH_BLOCK(d_a,f_a,dima,
& int_mb(k_a_offset),(h2b_1-1+noab*(h1b_1-1+noab*
& (p6b_1-noab-1+nvab*(p5b_1-noab-1)))))
if(.not.intorb) then
CALL GET_HASH_BLOCK(d_b,f_b,dimb,
& int_mb(k_b_offset),(p6b_2-1+(noab+nvab)*
& (p5b_2-1+(noab+nvab)*(p4b_2-1+(noab+nvab)*
& (p3b_2-1)))))
else
CALL GET_HASH_BLOCK_I(d_b,f_b,dimb,
& int_mb(k_b_offset),(p6b_2-1+(noab+nvab)*
& (p5b_2-1+(noab+nvab)*(p4b_2-1+(noab+nvab)*
& (p3b_2-1)))),p6b_2,p5b_2,p4b_2,p3b_2)
end if
if (p5b .eq. p6b) then
alpha = 1.0d0
else
alpha = 2.0d0
end if
h21d = int_mb(k_range+h1b-1)*int_mb(k_range+h2b-1)
p43d = int_mb(k_range+p3b-1)*int_mb(k_range+p4b-1)
p65d = int_mb(k_range+p5b-1)*int_mb(k_range+p6b-1)
!print*,'cudaMemcpyAsync A'
!x_a(1:dima) = f_a(1:dima)
!err = cudaMemcpyAsync(x_a,f_a,dima,stream)
!if (err.ne.0) then
! call errquit('cudaMemcpyAsync',err,UNKNOWN_ERR)
!endif
!print*,'cudaMemcpyAsync B'
!x_b(1:dimb) = f_b(1:dimb)
!err = cudaMemcpyAsync(x_b,f_b,dimb,stream)
!if (err.ne.0) then
! call errquit('cudaMemcpyAsync',err,UNKNOWN_ERR)
!endif
!print*,'cublasDgemm_v2',h21d,p43d,p65d,dima,dimb,dimc
!err = cublasDgemm_v2(handle,cu_op_n,cu_op_n,
& ! h21d,p43d,p65d,
& ! 0.5d0*alpha,x_a(1:dima),h21d,
& ! x_b(1:dimb),p65d,
& ! 1.0d0, x_c(1:dimc),h21d)
!if (err.ne.0) then
! call errquit('cublasDgemm_v2',err,UNKNOWN_ERR)
!endif
!err = cudaDeviceSynchronize()
!if (err.ne.0) then
! call errquit('cudaDeviceSync',err,UNKNOWN_ERR)
!endif
block
integer :: ih21d, ip43d, ip65d
ih21d = int_mb(k_range+h1b-1)*int_mb(k_range+h2b-1)
ip43d = int_mb(k_range+p3b-1)*int_mb(k_range+p4b-1)
ip65d = int_mb(k_range+p5b-1)*int_mb(k_range+p6b-1)
print*,'BLAS DGEMM'
call DGEMM('n','n',ih21d,ip43d,ip65d,
& 0.5d0*alpha,f_a,ih21d,
& f_b,ip65d,
& 1.0d0, f_d,ih21d)
end block
block
use cublas, only : dgemm
print*,'CUBLAS DGEMM'
call DGEMM('n','n',h21d,p43d,p65d,
& 0.5d0*alpha,f_a,h21d,
& f_b,p65d,
& 1.0d0, f_c,h21d)
end block
END IF
END IF
END IF
END DO
END DO
!print*,'cudaMemcpyAsync C (out)'
!f_c(1:dimc) = x_c(1:dimc)
!err = cudaMemcpyAsync(f_c,x_c,dimc,stream)
!if (err.ne.0) then
! call errquit('cudaMemcpyAsync',err,UNKNOWN_ERR)
!endif
!err = cudaDeviceSynchronize()
!if (err.ne.0) then
! call errquit('cudaDeviceSync',err,UNKNOWN_ERR)
!endif
f_c(1:dimc) = f_d(1:dimc)
CALL ADD_HASH_BLOCK(d_c,f_c,dimc,
& int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab*
& (p4b-noab-1+nvab*(p3b-noab-1)))))
next = NXTASK(nprocs, 1)
END IF
count = count + 1
END IF
END IF
END IF
END DO
END DO
END DO
END DO
next = NXTASK(-nprocs, 1)
call GA_SYNC()
deallocate(f_a,f_b,f_c,f_d,stat=e_a)
if (e_a.ne.0) call errquit("MA Hfree",e_a,MA_ERR)
deallocate(x_a,x_b,x_c,stat=e_a)
if (e_a.ne.0) call errquit("MA Gfree",e_a,MA_ERR)
err = cublasDestroy(handle)
if (err.ne.0) call errquit('cublasDestroy',err,UNKNOWN_ERR)
err = cudaStreamDestroy(stream)
if (err.ne.0) call errquit('cudaStreamDestroy',err,UNKNOWN_ERR)
END
#endif