Code Monkey home page Code Monkey logo

Comments (2)

jeffhammond avatar jeffhammond commented on August 17, 2024

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.

Lancashire3000 avatar Lancashire3000 commented on August 17, 2024

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)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo 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.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.