From 6de9fbeb57e91a3491121ea7d625f28f9f019139 Mon Sep 17 00:00:00 2001 From: Jeff Hammond Date: Tue, 10 Oct 2023 09:40:44 +0300 Subject: [PATCH] still debugging --- src/nwchem.nw | 37 ++-- src/tce/ccsd/ccsd_t2.F | 5 + src/tce/ccsd/ccsd_t2_8.F | 365 +++++++++++++++++++++++++++++++-------- 3 files changed, 320 insertions(+), 87 deletions(-) diff --git a/src/nwchem.nw b/src/nwchem.nw index 77b2b04ee9..9e338432fd 100644 --- a/src/nwchem.nw +++ b/src/nwchem.nw @@ -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 diff --git a/src/tce/ccsd/ccsd_t2.F b/src/tce/ccsd/ccsd_t2.F index 774b614377..e2c6b0ea5a 100644 --- a/src/tce/ccsd/ccsd_t2.F +++ b/src/tce/ccsd/ccsd_t2.F @@ -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) diff --git a/src/tce/ccsd/ccsd_t2_8.F b/src/tce/ccsd/ccsd_t2_8.F index 8acd5c2d17..63d9727942 100644 --- a/src/tce/ccsd/ccsd_t2_8.F +++ b/src/tce/ccsd/ccsd_t2_8.F @@ -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 + +