Multi threading exact exchange

[tkoskela]

Collaborators in Bordeaux are developing this and could use some help. Not sure how much is in the main branch.

[davidbowler]

We want the f-exx developed by @lionelalexandre

[tkoskela]

Notes from meeting with @lionelalexandre

• Exact exchange runs with a small number O(1) of atoms per core.
  • We want to speed up calculation per atom with OpenMP
  • Can't be divided to more MPI ranks because you run out of atoms per rank
• Code is in new branch f-exx-opt
• Main part is exx_kernel_default.f90
  • looks like matrix_multiply
    • Lot of code is identical
• Three schemes to compute exact exchange
  • Mainly looking at scheme 1 -- m_kern_exx_cri
  • Most of time spent calling exx_phi_on_grid
  • Comms outside of the loop
  • Look at scheme 3 if time
• Main hotspots
  • Poisson's equation - exx_poisson
    • Does a local fft calculation
  • Accumulation in Phi and Omega - exx_evalpao
  • Matrix multiplication exx_matmult
  • MPI communication should not be significant

Tests

• Lionel has added two tests to testsuite. test_004 and test_005
  • 004 is scheme 1 from before
    • You can modify EXX.GridSpacing to adjust the amount of work. Lower is more work.
    • Test has 6 atoms
      • Running on 1 or 2 MPI processes is meaningful, 4 may not be
      • To focus on OpenMP performance, best to run on 1 MPI process
  • 005 is scheme 3 from before
    • using an analytic way to compute exact exchange
    • Optimizing this has less priority

Tasks

• Profile test_004 in f-exx-opt in serial, understand hotspots
• Add OpenMP parallel directives at appropriate loop level
• Test thread scaling with test_004
• Merge code in f-exx-opt into develop

[ilectra]

Some useful reading:
Linear–scaling implementation of exact exchange using the density matrix: method and accuracy

[tkoskela]

exx_poisson

• Can we parallelise the loop nest where FFTs are being called?

  CONQUEST-release/src/exx_kernel_default.f90

  Line 1212 in 41cf7e2

  call exx_v_on_grid(inode,extent,work_in_3d,work_out_3d,r_int, &

exx_evalpao

• Look at redundant computations being done in the loop that calls

  CONQUEST-release/src/exx_evalpao.f90

  Line 187 in 41cf7e2

  call evaluate_pao(i_dummy,spec,l1,acz,m1,x,y,z,pao_val,exx_cartesian)

  
  for example, r is recomputed for every combination of quantum numbers when x,y and z don't change. Check the other places where evaluate_pao is called.

exx_matmult

in

CONQUEST-release/src/exx_kernel_default.f90

Lines 1256 to 1266 in 41cf7e2

	do r = 1, 2*extent+1
	do s = 1, 2*extent+1
	do t = 1, 2*extent+1
	exx_mat_elem = exx_mat_elem &
	+ phi_i(t,s,r,nsf1) &
	* Ome_kj(t,s,r,nsf3,nsf2) * dv
	end do
	end do
	end do

• Try !$omp parallel do collapse(3) reduction(+:exx_mat_elem) on loop over r
• Try replacing loops with a sum over element-wise product of phi and omega

[tkoskela]

The kpart loop in

CONQUEST-release/src/exx_kernel_default.f90

Line 442 in 41cf7e2

do kpart = 1,mult(S_X_SX)%ahalo%np_in_halo ! Main loop

duplicates large parts of the main loop in

CONQUEST-release/src/multiply_module.f90

Line 217 in 41cf7e2

main_loop: do kpart = 1,a_b_c%ahalo%np_in_halo

We need to understand how to reduce the duplication. I would like to see if we can call mat_mult from exx_kernel_default instead of replicating the code without making too disruptive modifications to mat_mult. Keeping in mind the changes in #290

[lionelalexandre]

exx_poisson : that would be a very difficult task and I am not sure this is possible (using MPI)

[lionelalexandre]

exx_evalpao : actually r is not needed anymore (some vestige of previous versions). I don't see clearly how to avoid the 6 nested loops for evalpao but maybe there is something to do.

[lionelalexandre]

The kpart loop in : (maybe to clarify) as you noticed this part of the code has been inspired by the mat_mult kernel. The idea behind was to get access locally to matrix elements which were managed by remote nodes. At that time, this simplest way to duplicate mat_mult fo this purpose. Now, I think it is possible to use mat_mult 'as if' to get access to the matrix elements, to store them, and then to compute EXX. We may need to discuss about this.

[tkoskela]

Let's close this issue as completed by #324 and start a new one for the follow-up work.