Comments (2)
Isn't numpy row major? Do you need to flip the indices on the Python version to get the equivalent data access pattern as Fortran?
In any case, I'm afraid that a contracted dimension of 32 is about an order-of-magnitude too small to saturate the compute capability of a modern Intel CPU (K=384 saturates, according to former colleagues in the MKL team), so I'm afraid you are limited by bandwidth limitations and other overheads.
It would be very interesting to perform the CTF comparisons on a very contractions, particularly the more expensive N^6 terms in CCSD. The so-called "four-particle ladder" term (https://github.com/nwchemgit/nwchem/blob/master/src/tce/ccsd/ccsd_t2_8.F) is the most expensive in TCE's CCSD, for a reasonable number of virtual orbitals.
All the terms in CC2 are N^5 (like MP2) and aren't super interesting from a compute perspective, because all the efficient implementations of CC2 and MP2 don't store the two-electron integrals in the fully transformed representation.
Unrelated to CTF, I had some fun porting your Fortran test to use GPUs using NVIDIA StdPar support (maps DO CONCURRENT
and MATMUL
to OpenACC and CUTENSOR, respectively). I'm happy to share and discuss that, but this isn't the right place for it.
from ctf.
Hi, thank you for your comments and suggestions.
I tried to increase the contracted dimension to 384, others 32, the output is
ctf 0.05176 ± 0.00117 s
ctf-symm 0.05084 ± 0.00063 s
I can transpose the array A
before the contraction, the result is similar
ctf 0.0558 ± 0.00089 s
ctf-symm 0.0546 ± 0.00062 s
The revised Python
and Fortran
code (admittedly it is messy) for different dimensions in each index
import numpy as np
import time
import scipy.stats
import warnings
#import timeit
import ctf
warnings.filterwarnings("ignore", category=DeprecationWarning)
# from https://stackoverflow.com/questions/15033511/compute-a-confidence-interval-from-sample-data
def mean_confidence_interval(data, var_name, unit, confidence=0.95):
a = 1.0 * np.array(data)
n = len(a)
m, se = np.mean(a), scipy.stats.sem(a)
h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)
print(var_name, round(m, 5), "\u00B1", round(h, 5), unit )
print('ab,bcde->acd')
# A[a,b] * B[b,c,d,e] = C[a,c,d,e]
na = nc = nd = ne = 32
nb = 384
#nt = 10
n_times = 128
#print('number of common dimension', dim_comm)
print('number of averaged time', n_times)
A = np.random.random((na,nb))
B = np.random.random((nb,nc,nd,ne))
B_symm= np.zeros((nb,nc,nd,ne))
# symmetrize B
for id in range(nd):
for ie in range(ne):
B_symm[:,:,id,ie] = (B[:,:,id,ie] + B[:,:,ie,id]) * 0.5
#N = dim_comm
C_ref = np.einsum('ab,bcde->acde', A, B_symm)
tA = ctf.tensor([na,nb],sym=[ctf.SYM.NS,ctf.SYM.NS])
tAT = ctf.tensor([nb,na],sym=[ctf.SYM.NS,ctf.SYM.NS])
tB2 = ctf.tensor([nb,nc,nd,ne],sym=[ctf.SYM.NS, ctf.SYM.NS, ctf.SYM.NS,ctf.SYM.NS])
tB2_symm = ctf.tensor([nb,nc,nd,ne],sym=[ctf.SYM.NS, ctf.SYM.NS, ctf.SYM.SY,ctf.SYM.NS])
tC = ctf.tensor([na,nc,nd,ne],sym=[ctf.SYM.NS, ctf.SYM.NS, ctf.SYM.NS,ctf.SYM.NS])
tC_symm = ctf.tensor([na,nc,nd,ne],sym=[ctf.SYM.NS, ctf.SYM.NS, ctf.SYM.SY,ctf.SYM.NS])
tA = ctf.from_nparray(A)
tAT = ctf.from_nparray(np.transpose(A))
tB2 = ctf.from_nparray(B_symm)
tB2_symm = ctf.from_nparray(B_symm)
list_time = []
for i in range(n_times):
start_time = time.time()
tC = ctf.einsum('ab,bcde->acde', tA, tB2)
finish_time = time.time()
list_time.append(finish_time - start_time)
mean_confidence_interval(list_time, 'ctf', 's' )
list_time = []
for i in range(n_times):
start_time = time.time()
tC_symm = ctf.einsum('ab,bcde->acde', tA, tB2_symm)
finish_time = time.time()
list_time.append(finish_time - start_time)
mean_confidence_interval(list_time, 'ctf-symm', 's' )
tC_temp = ctf.to_nparray(tC)
print('check ctf',np.allclose(C_ref,tC_temp))
tC_symm_temp = ctf.to_nparray(tC_symm)
print('check ctf symm',np.allclose(C_ref,tC_symm_temp))
list_time = []
for i in range(n_times):
start_time = time.time()
tC = ctf.einsum('ba,bcde->acde', tAT, tB2)
finish_time = time.time()
list_time.append(finish_time - start_time)
mean_confidence_interval(list_time, 'ctf', 's' )
list_time = []
for i in range(n_times):
start_time = time.time()
tC_symm = ctf.einsum('ba,bcde->acde', tAT, tB2_symm)
finish_time = time.time()
list_time.append(finish_time - start_time)
mean_confidence_interval(list_time, 'ctf-symm', 's' )
tC_temp = ctf.to_nparray(tC)
print('check ctf',np.allclose(C_ref,tC_temp))
tC_symm_temp = ctf.to_nparray(tC_symm)
print('check ctf symm',np.allclose(C_ref,tC_symm_temp))
Program test_contraction
!A[a,b] * B[b,c,d,e] = C[a,c,d,e]
integer, parameter :: dp = selected_real_kind(15, 307)
Real( dp ), Dimension( :, : ), Allocatable :: a
Real( dp ), Dimension( :, :, :, : ), Allocatable :: b0, b
Real( dp ), Dimension( :, :, :, : ), Allocatable :: c0, c1, c2
Integer :: na, nb, nc, nd, ne, m, n, m_iter
Integer :: ia, ib, ic, id, ie
Integer :: start, finish, rate
real(dp) :: sum_time
Write( *, * ) 'na, nb, nc, nd, ne ?'
Read( *, * ) na, nb, nc, nd, ne
Allocate( a ( 1:na, 1:nb ) )
Allocate( b0 ( 1:nb, 1:nc, 1:nd, 1:ne ) )
Allocate( b ( 1:nb, 1:nc, 1:nd, 1:ne ) )
Allocate( c0( 1:na, 1:nc, 1:nd, 1:ne ) )
Allocate( c1( 1:na, 1:nc, 1:nd, 1:ne ) )
Allocate( c2( 1:na, 1:nc, 1:nd, 1:ne ) )
Call Random_number( a )
Call Random_number( b0 )
c1 = 0.0_dp
c2 = 0.0_dp
m_iter = 10
write (*,*) 'm_iter average', m_iter
!symmetrize
do id = 1, nd
do ie = 1, ne
b(:, :, id, ie) = ( b0(:, :, id, ie) + b0(:, :, ie, id) )*0.5
end do
end do
! refrence value by nested loops
sum_time = 0.0
do m = 1, m_iter
c0 = 0.0_dp
Call System_clock( start, rate )
do ie = 1, ne
do id = 1, nd
do ic = 1, nc
do ib = 1, nb
do ia = 1, na
c0(ia, ic, id, ie) = c0(ia, ic, id, ie) + a(ia, ib) * b(ib, ic, id, ie)
end do
end do
end do
end do
end do
Call System_clock( finish, rate )
sum_time = sum_time + Real( finish - start, dp ) / rate
end do
Write( *, * ) 'Time for naive loop edcba', sum_time / m_iter
sum_time = 0.0
do m = 1, m_iter
do ie = 1, ne
do id = 1, nd
Call System_clock( start, rate )
Call dgemm( 'N', 'N', na, nc, nb, 1.0_dp, a , Size( a , Dim = 1 ), &
b(:,:,id, ie) , Size( b , Dim = 1 ), &
0.0_dp, c1(:,:,id, ie), Size( c1, Dim = 1 ) )
!c5(:,id)
Call System_clock( finish, rate )
sum_time = sum_time + Real( finish - start, dp ) / rate
end do
end do
end do
Write( *, * ) 'Time for loop dgemm-2', sum_time / m_iter
sum_time = 0.0
do m = 1, m_iter
do ie = 1, ne
do id = 1, ie
Call System_clock( start, rate )
Call dgemm( 'N', 'N', na, nc, nb, 1.0_dp, a , Size( a , Dim = 1 ), &
b(:,:,id, ie) , Size( b , Dim = 1 ), &
0.0_dp, c2(:,:,id, ie), Size( c2, Dim = 1 ) )
!c5(:,id)
Call System_clock( finish, rate )
sum_time = sum_time + Real( finish - start, dp ) / rate
c2(:,:,ie, id) = c2(:,:,id, ie)
end do
end do
end do
Write( *, * ) 'Time for loop dgemm-symmetry', sum_time / m_iter
do ia = 1, na
do ic = 1, nc
do id = 1, nd
do ie = 1, ne
if ( dabs(c1(ia,ic,id,ie) - c0(ia,ic,id,ie)) > 1.e-6 ) then
write (*,*) '!!! c1', ia,ic,id,ie, c0(ia,ic,id,ie), c1(ia,ic,id,ie)
endif
if ( dabs(c2(ia,ic,id,ie) - c0(ia,ic,id,ie)) > 1.e-6 ) then
write (*,*) '!!! c2', ia,ic,id,ie, c2(ia,ic,id,ie), c1(ia,ic,id,ie)
endif
enddo
enddo
enddo
enddo
end
from ctf.
Related Issues (20)
- Efficiency comparisons with einsum and opt-einsum and how to utilize symmetry? HOT 4
- ctf-einsum.py:78: DeprecationWarning: `np.complex` is a deprecated alias for the builtin `complex`. HOT 1
- issues using sparse file io to load tensors HOT 4
- Compile error with GCC 11 HOT 10
- Wrong results when slicing a symmetric sparse tensor in python lib
- understanding performance overheads in CTF HOT 27
- ModuleNotFoundError: No module named 'ctf.core' HOT 1
- segfault executing sparse inner product HOT 18
- oom/memory corruption running an SDDMM (using TTTP specialized routine) HOT 1
- unexpected performance for SpMV operation HOT 1
- [question] setting all nonzeros to a value HOT 2
- Warnings
- make test failure in SVD test HOT 3
- compiling issue when including ctf.hpp
- test_suite failure on Apple HOT 1
- compile issues with undefined references to mkl commands (that appear in the relevant folders) HOT 9
- Scale with Endomorphism
- Set values with low memory footprint
- ctf header error during installation on perlmutter HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from ctf.