diff --git a/src/tce/ccsd/GNUmakefile b/src/tce/ccsd/GNUmakefile index 5a0c2edb7f..5d593e51ca 100644 --- a/src/tce/ccsd/GNUmakefile +++ b/src/tce/ccsd/GNUmakefile @@ -13,7 +13,7 @@ LIBRARY = libtce.a USES_BLAS = ccsd_e.F ccsd_t1.F ccsd_t2.F cc2_t1.F cc2_t2.F \ ccsd_1prdm_hh.F ccsd_1prdm_hp.F ccsd_1prdm_ph.F \ ccsd_1prdm_pp.F ccsd_1prdm.F \ - icsd_t1.F icsd_t2.F ccsd_t2_8.F + icsd_t1.F icsd_t2.F ccsd_t2_8.F ccsd_kernels.F LIB_DEFINES = -DDEBUG_PRINT diff --git a/src/tce/ccsd/ccsd_kernels.F b/src/tce/ccsd/ccsd_kernels.F index 0fea703fa1..14d9e32a73 100644 --- a/src/tce/ccsd/ccsd_kernels.F +++ b/src/tce/ccsd/ccsd_kernels.F @@ -1,90 +1,39 @@ ! ! Reference implementation written by Karol Kowalski +! Loops fully converted to DGEMM later. ! subroutine t2_p8(h1d,h2d,p3d,p4d,p5d,p6d, 1 t2sub,v2sub,r2sub,factor) IMPLICIT NONE - integer h1d,h2d,p3d,p4d,p5d,p6d - integer h1,h2,p3,p4,p5,p6 - double precision t2sub(h2d,h1d,p6d,p5d) - double precision v2sub(p6d,p5d,p4d,p3d) - double precision r2sub(h2d,h1d,p4d,p3d) - double precision factor - if ((p5d.lt.8).or.(p6d.lt.8)) then -!$omp parallel do collapse(2) -!$omp& default(shared) schedule(static) -!$omp& private(h1,h2,p3,p4,p5,p6) - do p3=1,p3d - do p4=1,p4d - do h1=1,h1d - do h2=1,h2d - do p5=1,p5d -!$omp simd - do p6=1,p6d - r2sub(h2,h1,p4,p3)=r2sub(h2,h1,p4,p3) - & + factor*t2sub(h2,h1,p6,p5)*v2sub(p6,p5,p4,p3) - enddo -!$omp end simd - enddo - enddo - enddo - enddo - enddo -!$omp end parallel do -! -! All dims are at least 8, so more SIMD optimizations allowed. -! - else ! inner loops at least 8 -!$omp parallel do collapse(2) -!$omp& default(shared) schedule(static) -!$omp& private(h1,h2,p3,p4,p5,p6) - do p3=1,p3d - do p4=1,p4d - do h1=1,h1d - do h2=1,h2d -!dir$ loop count min(8), max(80), avg(32) -!dec$ unroll_and_jam = 8 - do p5=1,p5d -!dir$ loop count min(8), max(80), avg(32) -!dec$ unroll_and_jam = 8 -!$omp simd - do p6=1,p6d - r2sub(h2,h1,p4,p3)=r2sub(h2,h1,p4,p3) - & + factor*t2sub(h2,h1,p6,p5)*v2sub(p6,p5,p4,p3) - enddo -!$omp end simd - enddo - enddo - enddo - enddo - enddo -!$omp end parallel do - endif ! inner loops at least 8 - return + integer, intent(in) :: h1d,h2d,p3d,p4d,p5d,p6d + double precision, intent(in) :: factor + double precision, intent(in) :: t2sub(h2d*h1d,p6d*p5d) + double precision, intent(in) :: v2sub(p6d*p5d,p4d*p3d) + double precision, intent(inout) :: r2sub(h2d*h1d,p4d*p3d) +#if 0 + r2sub = r2sub + factor * matmul(t2sub,v2sub) +#else + call DGEMM('n','n',h2d*h1d,p4d*p3d,p6d*p5d, + & factor,t2sub,h2d*h1d, + & v2sub,p6d*p5d, + & 1.0d0, r2sub,h2d*h1d) +#endif end subroutine t2_p8_x(h1d,h2d,p3d,p4d,p5d,p6d, 1 t2sub,v2sub,r2sub,factor) IMPLICIT NONE - integer h1d,h2d,p3d,p4d,p5d,p6d - integer h1,h2,p3,p4,p5,p6 - double precision t2sub(h2d,h1d,p6d,p5d) - double precision v2sub(p6d,p5d,p4d,p3d) - double precision r2sub(h2d,h1d,p4d,p3d) - double precision factor - do p3=1,p3d - do p4=1,p4d - do h1=1,h1d - do h2=1,h2d - do p5=1,p5d - do p6=1,p6d - r2sub(h2,h1,p4,p3)=r2sub(h2,h1,p4,p3) - & + factor*t2sub(h2,h1,p6,p5)*v2sub(p6,p5,p4,p3) - enddo - enddo - enddo - enddo - enddo - enddo - return + integer, intent(in) :: h1d,h2d,p3d,p4d,p5d,p6d + double precision, intent(in) :: factor + double precision, intent(in) :: t2sub(h2d*h1d,p6d*p5d) + double precision, intent(in) :: v2sub(p6d*p5d,p4d*p3d) + double precision, intent(inout) :: r2sub(h2d*h1d,p4d*p3d) +#if 0 + r2sub = r2sub + factor * matmul(t2sub,v2sub) +#else + call DGEMM('n','n',h2d*h1d,p4d*p3d,p6d*p5d, + & factor,t2sub,h2d*h1d, + & v2sub,p6d*p5d, + & 1.0d0, r2sub,h2d*h1d) +#endif end diff --git a/src/tce/ccsd/ccsd_t2_8.F b/src/tce/ccsd/ccsd_t2_8.F index 4aa10195eb..fbbb18ad38 100644 --- a/src/tce/ccsd/ccsd_t2_8.F +++ b/src/tce/ccsd/ccsd_t2_8.F @@ -195,22 +195,16 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( 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 :: h21d, p43d, p65d #ifdef USE_F90_ALLOCATABLE - double precision, allocatable :: f_a(:) - double precision, allocatable :: f_b(:) - double precision, allocatable :: f_c(:) - double precision, allocatable :: f_t(:) + double precision, allocatable :: f_a(:), f_b(:), f_c(:) #ifdef USE_FASTMEM - !dec$ attributes fastmem :: f_a,f_b,f_c,f_t + !dec$ attributes fastmem :: f_a,f_b,f_c #endif - integer :: e_a,e_b,e_c,e_t #else - integer k_a, l_a - integer k_b, l_b - integer k_c, l_c - integer k_t, l_t - integer e_a,e_b,e_c,e_t + integer :: k_a, l_a, k_b, l_b, k_c, l_c #endif + integer :: e_a,e_b,e_c double precision alpha integer p5b_in,p6b_in INTEGER NXTASK @@ -223,32 +217,21 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( dimpppp = maxp*maxp*maxp*maxp dimtemp = max(dimpppp,dimhhpp) + e_a=0 + e_b=0 + e_c=0 #ifdef USE_F90_ALLOCATABLE allocate(f_a(1:dimhhpp),stat=e_a) allocate(f_b(1:dimpppp),stat=e_b) allocate(f_c(1:dimhhpp),stat=e_c) -# ifndef USE_LOOPS_NOT_DGEMM - allocate(f_t(1:dimtemp),stat=e_t) -# endif #else - e_a=0 if(.not.MA_PUSH_GET(mt_dbl,dimhhpp,"a",l_a,k_a)) e_a=-1 - e_b=0 if(.not.MA_PUSH_GET(mt_dbl,dimpppp,"b",l_b,k_b)) e_b=-1 - e_c=0 if(.not.MA_PUSH_GET(mt_dbl,dimhhpp,"c",l_c,k_c)) e_c=-1 -# ifndef USE_LOOPS_NOT_DGEMM - e_t=0 - if(.not.MA_PUSH_GET(mt_dbl,dimtemp,"t",l_t,k_t)) e_t=-1 -# else - dimtemp=-12345 - e_t=.false. -# endif #endif if (e_a.ne.0) call errquit("MA a",dimhhpp,MA_ERR) if (e_b.ne.0) call errquit("MA b",dimpppp,MA_ERR) if (e_c.ne.0) call errquit("MA c",dimhhpp,MA_ERR) - if (e_t.ne.0) call errquit("MA t",dimtemp,MA_ERR) DO p3b = noab+1,noab+nvab DO p4b = p3b,noab+nvab DO h1b = 1,noab @@ -274,10 +257,10 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( DO p5b = noab+1,noab+nvab DO p6b = p5b,noab+nvab #else - 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 + 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 #endif 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 @@ -293,92 +276,51 @@ C i0 ( p3 p4 h1 h2 )_vt + = 1/2 * Sum ( p5 p6 ) * t ( p5 p6 h1 h2 )_t * v ( dima = dim_common * dima_sort dimb = dim_common * dimb_sort IF ((dima .gt. 0) .and. (dimb .gt. 0)) THEN -#ifdef USE_LOOPS_NOT_DGEMM 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))))) -#else - CALL GET_HASH_BLOCK(d_a,f_t,dima, - & int_mb(k_a_offset),(h2b_1-1+noab*(h1b_1-1+noab* - & (p6b_1-noab-1+nvab*(p5b_1-noab-1))))) - CALL TCE_SORT_4(f_t,f_a, - & int_mb(k_range+p5b-1),int_mb(k_range+p6b-1), - & int_mb(k_range+h1b-1),int_mb(k_range+h2b-1), - & 4,3,2,1,1.0d0) -#endif if(.not.intorb) then -#ifdef USE_LOOPS_NOT_DGEMM 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(d_b,f_t,dimb, - & int_mb(k_b_offset),(p6b_2-1+(noab+nvab)* - & (p5b_2-1+(noab+nvab)*(p4b_2-1+(noab+nvab)* - & (p3b_2-1))))) -#endif else -#ifdef USE_LOOPS_NOT_DGEMM 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) -#else - CALL GET_HASH_BLOCK_I(d_b,f_t,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) -#endif end if -#ifndef USE_LOOPS_NOT_DGEMM - CALL TCE_SORT_4(f_t,f_b, - & int_mb(k_range+p3b-1),int_mb(k_range+p4b-1), - & int_mb(k_range+p5b-1),int_mb(k_range+p6b-1), - & 2,1,4,3,1.0d0) -#endif if (p5b .eq. p6b) then alpha = 1.0d0 else alpha = 2.0d0 end if -#ifdef USE_LOOPS_NOT_DGEMM - call t2_p8(int_mb(k_range+h1b-1), - & int_mb(k_range+h2b-1), - & int_mb(k_range+p3b-1), - & int_mb(k_range+p4b-1), - & int_mb(k_range+p5b-1), - & int_mb(k_range+p6b-1), - & f_a,f_b,f_c, - & 0.5d0*alpha) -#else - CALL DGEMM('T','N',dima_sort,dimb_sort,dim_common, - & alpha,f_a,dim_common,f_b, - & dim_common,1.0d0,f_c,dima_sort) -#endif + 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) + call DGEMM('n','n',h21d,p43d,p65d, + & 0.5d0*alpha,f_a,h21d, + & f_b,p65d, + & 1.0d0, f_c,h21d) END IF END IF END IF END DO END DO -#ifdef USE_LOOPS_NOT_DGEMM 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))))) -#else - CALL TCE_SORT_4(f_c,f_t, - & int_mb(k_range+p4b-1),int_mb(k_range+p3b-1), - & int_mb(k_range+h2b-1),int_mb(k_range+h1b-1), - & 2,1,4,3,0.5d0) - CALL ADD_HASH_BLOCK(d_c,f_t,dimc, - & int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab* - & (p4b-noab-1+nvab*(p3b-noab-1))))) -#endif #else -celse// USE_F90_ALLOCATABLE CALL DFILL(dimc,0.0d0,dbl_mb(k_c),1) +#if 0 DO p5b = noab+1,noab+nvab DO p6b = p5b,noab+nvab +#else + 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 +#endif 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), @@ -393,89 +335,41 @@ celse// USE_F90_ALLOCATABLE dima = dim_common * dima_sort dimb = dim_common * dimb_sort IF ((dima .gt. 0) .and. (dimb .gt. 0)) THEN -#ifdef USE_LOOPS_NOT_DGEMM CALL GET_HASH_BLOCK(d_a,dbl_mb(k_a),dima, & int_mb(k_a_offset),(h2b_1-1+noab*(h1b_1-1+noab* & (p6b_1-noab-1+nvab*(p5b_1-noab-1))))) -#else - CALL GET_HASH_BLOCK(d_a,dbl_mb(k_t),dima, - & int_mb(k_a_offset),(h2b_1-1+noab*(h1b_1-1+noab* - & (p6b_1-noab-1+nvab*(p5b_1-noab-1))))) - CALL TCE_SORT_4(dbl_mb(k_t),dbl_mb(k_a), - & int_mb(k_range+p5b-1),int_mb(k_range+p6b-1), - & int_mb(k_range+h1b-1),int_mb(k_range+h2b-1), - & 4,3,2,1,1.0d0) -#endif if(.not.intorb) then -#ifdef USE_LOOPS_NOT_DGEMM CALL GET_HASH_BLOCK(d_b,dbl_mb(k_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(d_b,dbl_mb(k_t),dimb, - & int_mb(k_b_offset),(p6b_2-1+(noab+nvab)* - & (p5b_2-1+(noab+nvab)*(p4b_2-1+(noab+nvab)* - & (p3b_2-1))))) -#endif else -#ifdef USE_LOOPS_NOT_DGEMM CALL GET_HASH_BLOCK_I(d_b,dbl_mb(k_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) -#else - CALL GET_HASH_BLOCK_I(d_b,dbl_mb(k_t),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) -#endif end if -#ifndef USE_LOOPS_NOT_DGEMM - CALL TCE_SORT_4(dbl_mb(k_t),dbl_mb(k_b), - & int_mb(k_range+p3b-1),int_mb(k_range+p4b-1), - & int_mb(k_range+p5b-1),int_mb(k_range+p6b-1), - & 2,1,4,3,1.0d0) -#endif if (p5b .eq. p6b) then alpha = 1.0d0 else alpha = 2.0d0 end if -#ifdef USE_LOOPS_NOT_DGEMM - call t2_p8(int_mb(k_range+h1b-1), - & int_mb(k_range+h2b-1), - & int_mb(k_range+p3b-1), - & int_mb(k_range+p4b-1), - & int_mb(k_range+p5b-1), - & int_mb(k_range+p6b-1), - & dbl_mb(k_a),dbl_mb(k_b),dbl_mb(k_c), - & 0.5d0*alpha) -#else - CALL DGEMM('T','N',dima_sort,dimb_sort,dim_common, - & alpha,dbl_mb(k_a),dim_common,dbl_mb(k_b), - & dim_common,1.0d0,dbl_mb(k_c),dima_sort) -#endif + 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) + call DGEMM('n','n',h21d,p43d,p65d, + & 0.5d0*alpha,dbl_mb(k_a),h21d, + & dbl_mb(k_b),p65d, + & 1.0d0, dbl_mb(k_c),h21d) END IF END IF END IF END DO END DO -#ifdef USE_LOOPS_NOT_DGEMM CALL ADD_HASH_BLOCK(d_c,dbl_mb(k_c),dimc, & int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab* & (p4b-noab-1+nvab*(p3b-noab-1))))) -#else - CALL TCE_SORT_4(dbl_mb(k_c),dbl_mb(k_t), - & int_mb(k_range+p4b-1),int_mb(k_range+p3b-1), - & int_mb(k_range+h2b-1),int_mb(k_range+h1b-1), - & 2,1,4,3,0.5d0) - CALL ADD_HASH_BLOCK(d_c,dbl_mb(k_t),dimc, - & int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab* - & (p4b-noab-1+nvab*(p3b-noab-1))))) -#endif #endif -cendif// USE_F90_ALLOCATABLE next = NXTASK(nprocs, 1) END IF count = count + 1 @@ -493,22 +387,11 @@ cendif// USE_F90_ALLOCATABLE deallocate(f_a,stat=e_a) deallocate(f_b,stat=e_b) deallocate(f_c,stat=e_c) -# ifndef USE_LOOPS_NOT_DGEMM - deallocate(f_t,stat=e_t) -# endif #else -# ifndef USE_LOOPS_NOT_DGEMM - e_t=0 - if(.not.MA_POP_STACK(l_t)) e_t=-1 -# else - l_t=-12345 - e_t=0 -# endif e_a=0 if(.not.MA_CHOP_STACK(l_a)) e_a=-1 #endif if (e_a.ne.0) call errquit("MA pops a",0,MA_ERR) - if (e_t.ne.0) call errquit("MA pops t",1,MA_ERR) RETURN END diff --git a/src/tce/ccsd/icsd_t2.F b/src/tce/ccsd/icsd_t2.F index f3355e4db4..7e4bb8dd2f 100644 --- a/src/tce/ccsd/icsd_t2.F +++ b/src/tce/ccsd/icsd_t2.F @@ -8159,8 +8159,8 @@ c old way next = NXTASK(-nprocs, 1) c old way call GA_SYNC() RETURN END - SUBROUTINE icsd_t2_8(d_a,k_a_offset,d_b,k_b_offset,d_c,k_c_offset, - &ctx,icounter) + SUBROUTINE icsd_t2_8_old(d_a,k_a_offset,d_b,k_b_offset, + & d_c,k_c_offset,ctx,icounter) C $Id$ C This is a Fortran77 program generated by Tensor Contraction Engine v.1.0 C Copyright (c) Battelle & Pacific Northwest National Laboratory (2002) @@ -9022,3 +9022,239 @@ c c c c + SUBROUTINE icsd_t2_8(d_a,k_a_offset, + & d_b,k_b_offset, + & d_c,k_c_offset, + & ctx,icounter) +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 + 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,dimtemp + 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 :: h21d, p43d, p65d +#ifdef USE_F90_ALLOCATABLE + double precision, allocatable :: f_a(:), f_b(:), f_c(:) +#ifdef USE_FASTMEM + !dec$ attributes fastmem :: f_a,f_b,f_c +#endif +#else + integer :: k_a, l_a, k_b, l_b, k_c, l_c +#endif + integer :: e_a,e_b,e_c + double precision alpha + integer p5b_in,p6b_in + integer ctx,icounter + external nxt_ctx_create, nxt_ctx_destroy, nxt_ctx_next + integer :: p,h + nprocs = GA_NNODES() + count = 0 + !next = NXTASK(nprocs, 1) + call nxt_ctx_next(ctx, icounter, next) + + ! TODO hoist this like ccsd_t2 path + maxp = 0 + do p = noab+1,noab+nvab + maxp = max(maxp,int_mb(k_range+p-1)) + enddo + maxh = 0 + do h = 1,noab + maxh = max(maxh,int_mb(k_range+h-1)) + enddo + + dimhhpp = maxh*maxh*maxp*maxp + dimpppp = maxp*maxp*maxp*maxp + dimtemp = max(dimpppp,dimhhpp) + + e_a=0 + e_b=0 + e_c=0 +#ifdef USE_F90_ALLOCATABLE + allocate(f_a(1:dimhhpp),stat=e_a) + allocate(f_b(1:dimpppp),stat=e_b) + allocate(f_c(1:dimhhpp),stat=e_c) +#else + if(.not.MA_PUSH_GET(mt_dbl,dimhhpp,"a",l_a,k_a)) e_a=-1 + if(.not.MA_PUSH_GET(mt_dbl,dimpppp,"b",l_b,k_b)) e_b=-1 + if(.not.MA_PUSH_GET(mt_dbl,dimhhpp,"c",l_c,k_c)) e_c=-1 +#endif + if (e_a.ne.0) call errquit("MA a",dimhhpp,MA_ERR) + if (e_b.ne.0) call errquit("MA b",dimpppp,MA_ERR) + if (e_c.ne.0) call errquit("MA c",dimhhpp,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) +#ifdef USE_F90_ALLOCATABLE + CALL DFILL(dimc,0.0d0,f_c,1) +#if 0 + DO p5b = noab+1,noab+nvab + DO p6b = p5b,noab+nvab +#else + 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 +#endif + 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) + call DGEMM('n','n',h21d,p43d,p65d, + & 0.5d0*alpha,f_a,h21d, + & f_b,p65d, + & 1.0d0, f_c,h21d) + END IF + END IF + END IF + END DO + END DO + 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))))) +#else + CALL DFILL(dimc,0.0d0,dbl_mb(k_c),1) +#if 0 + DO p5b = noab+1,noab+nvab + DO p6b = p5b,noab+nvab +#else + 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 +#endif + 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,dbl_mb(k_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,dbl_mb(k_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,dbl_mb(k_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) + call DGEMM('n','n',h21d,p43d,p65d, + & 0.5d0*alpha,dbl_mb(k_a),h21d, + & dbl_mb(k_b),p65d, + & 1.0d0, dbl_mb(k_c),h21d) + END IF + END IF + END IF + END DO + END DO + CALL ADD_HASH_BLOCK(d_c,dbl_mb(k_c),dimc, + & int_mb(k_c_offset),(h2b-1+noab*(h1b-1+noab* + & (p4b-noab-1+nvab*(p3b-noab-1))))) +#endif + !next = NXTASK(nprocs, 1) + call nxt_ctx_next(ctx, icounter, next) + 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() + +#ifdef USE_F90_ALLOCATABLE + deallocate(f_a,stat=e_a) + deallocate(f_b,stat=e_b) + deallocate(f_c,stat=e_c) +#else + e_a=0 + if(.not.MA_CHOP_STACK(l_a)) e_a=-1 +#endif + if (e_a.ne.0) call errquit("MA pops a",0,MA_ERR) + RETURN + END + +