4
\$\begingroup\$

I would like to improve the performance of the piece of code below (in Fortran). It gives good results but the profiling tells me that it is where it spends most of its running time.

Basically, it increments over a time window (j loop) and performs a spate of 2x2 matrix multiplications. I was wondering where I can refactor it to make it much more efficient. For instance, is there a smarter and faster way to initialise the identity matrices?

do j=1,nb_pts_t + 1
   t = t + dt
   ! Compute Zgates for all pairs in up/down states
   call Z_gate (eigen_ener_up, t / dble(2 * CP_seq), Zgate_u)
   call Z_gate (eigen_ener_down, t / dble(2 * CP_seq), Zgate_d)

   ! Initialize identity matrices
   Tfu%elements(1,1) = dcmplx(1.d0, 0.d0)
   Tfu%elements(1,2) = dcmplx(0.d0, 0.d0)
   Tfu%elements(2,1) = dcmplx(0.d0, 0.d0)
   Tfu%elements(2,2) = dcmplx(1.d0, 0.d0)
   Tfd%elements(1,1) = dcmplx(1.d0, 0.d0)
   Tfd%elements(1,2) = dcmplx(0.d0, 0.d0)
   Tfd%elements(2,1) = dcmplx(0.d0, 0.d0)
   Tfd%elements(2,2) = dcmplx(1.d0, 0.d0)

   ! work out transition matrices in up/down states
   ! and compute the decoherence
   do i=1,nb_pairs
      ! rotate to eigenbasis and propagate
      Tu(i)%elements = matmul(Zgate_u(i)%elements, matrot_u(i)%elements)
      ! rotate back to bath basis
      Tu(i)%elements = matmul(matrottrans_u(i)%elements, Tu(i)%elements)
      ! Idem down state
      Td(i)%elements = matmul(Zgate_d(i)%elements, matrot_d(i)%elements)
      Td(i)%elements = matmul(matrottrans_d(i)%elements, Td(i)%elements)

      ! Initialize Tud and Tdu
      Tud(i)%elements = matmul(Tu(i)%elements, Td(i)%elements)
      Tdu(i)%elements = matmul(Td(i)%elements, Tu(i)%elements)

     do k=1,CP_seq
         if(modulo(k, 2) .ne. 0) then
            Tfu(i)%elements = matmul(Tfu(i)%elements, Tud(i)%elements)
            Tfd(i)%elements = matmul(Tfd(i)%elements, Tdu(i)%elements)
         else if(modulo(k, 2) == 0) then
            Tfu(i)%elements = matmul(Tfu(i)%elements, Tdu(i)%elements)
            Tfd(i)%elements = matmul(Tfd(i)%elements, Tud(i)%elements)
         end if
      end do

      ! Decoherence from initial |down-up> bath state
      L_pairs(i) = abs(conjg(Tfd(i)%elements(1, 1))*Tfu(i)%elements(1, 1) &
           - Tfd(i)%elements(1, 2) * Tfu(i)%elements(2, 1))

      ! average over the bath states
      L_pairs(i) = 0.5d0 + 0.5d0 * L_pairs(i)
   end do

   ! Final decay as the product over all pair decays
   L = product(L_pairs) 

   ! write the output
   write(16, fmt)t, L

end do
Improve this question
\$\endgroup\$
3
\$\begingroup\$

  1. Matmul is good for large matrices but not for small matrices. Write explicitly the 2x2 matrix without loops. The compiler will do a great job for you by reordering all operations in the best possible order.

    A(1,1) = B(1,1)*C(1,1) + B(1,2)*C(2,1)
A(2,1) = etc...

  2. If you want performance, avoid using user-defined types and prefer arrays. Define your matrices as:

    complex*16 :: A(2,2,nb_pairs)

    This will greatly enhance you memory accesses. Moreover, with the Intel Fortran compiler, you can add a directive to align your arrays on a 256-bit boundary:

    !DIR$ ATTRIBUTES ALIGN : 32 :: A, B, C

    This will allow your compiler to use AVX or SSE instructions to vectorize your operations (add, mul and load/store).

  3. Also, you can replace your modulo test by checking if the 1st bit of k is set to zero: A modulo is like an integer division (very expensive) but checking the last bit can be done in 1 CPU cycle, and the compiler can do a better optimization in that case. Replace 

    if (modulo(k, 2) .ne. 0)

    with

    if( iand(k,1) == 1)
Improve this answer
\$\endgroup\$
  • \$\begingroup\$ Cool, thanks Anthony for your help ! I will put this in action. \$\endgroup\$ – Roland Nov 20 '14 at 12:21
  • \$\begingroup\$ Just a question since fastest access to fortran arrays is the leftmost indice, wouldn't it be better to define A(nb_pairs, 2, 2) since the loop is mainly over nb_pairs ? \$\endgroup\$ – Roland Nov 20 '14 at 15:13
  • 1
    \$\begingroup\$ Probably not because it will make all the matrix elements of one matrix far away in memory. You will make a much harder pressure on the memory bandwidth as you will change one prefetch stream into 4 prefetch streams (one stream for each matrix element). If you leave it as the last indice, you will be able to fetch all the matrices in a single flow from the RAM to the CPU caches. Anyway, you are memory-bound as you do very few flops per memory access, so you should put the most effort to fetch your data as one continuous stream. Minimizing the number of streams is also better for openMP. \$\endgroup\$ – Anthony Scemama Nov 20 '14 at 16:13
  • 2
    \$\begingroup\$ @Kyle Kanos: (a) I made a benchmark here where I get a 1.25-7.0x acceleration on 2x2 matrices when removing matmul (ifort and gfortran) : github.com/scemama/matmul (b) OK (c) Accessing an array that is inside a type can be slower because the code will have to calculate the shift between the address of the type instance and the array. If this shift is not known at compile-time or not predictible, it adds an overhead. Also, it adds some constraints on the memory addresses : the alignment of a 2x2 matrix may be such that the matrix is always split in 2 cache lines when it could hold in 1. \$\endgroup\$ – Anthony Scemama Jan 14 '15 at 1:43
  • 1
    \$\begingroup\$ @AnthonyScemama: Don't take this the wrong way, but your code is a Fortran programmer's nightmare: the mixing of non-standard constructs and some elements of F90 into F77. I suppose I'm fortunate you're not using do-continue loops. \$\endgroup\$ – Kyle Kanos Jan 14 '15 at 2:23
4
\$\begingroup\$

Leave matmul in
You can test the speed of matmul versus do-loops versus dot_product for any size arrays via the program found in this SO question. You can test this with any compiler and optimization you want, but the results are that any matrix of size < 50 (probably more, but I didn't test it higher than that), the three methods are identical in execution time.

Identity matrix initialization
This can be done with vector assignment and a simple loop

Tfu(i)%elements(:,:) = cmplx(0d0, 0d0) 
do k=1,2
   Tfu(i)%elements(k,k) = cmplx(1d0,0d0)
end do

While Fortran defaults the imaginary component to zero when the 2nd argument is not present (e.g., cmplx(1d0)), it is probably better to be more clear by including it. Note that this won't change the run time, but it will make the declaration more compact.

Time variable in call to Z_gate
AFAIK, when there is math to be done in subroutine calls, Fortran will make a temporary copy of that variable. While it isn't much space for a double precision variable, you should still declare a new variable t_temp=t / dble(2 * CP_seq) and use that in place of your function call, especially since the math is done twice in subsequent calls to Z_gate.

It might also be worth declaring a new variable outside the j-loop

cp_seq_inv = 1d0/dble(cp_seq)

and making use of it as

t_temp = t * 0.5 * cp_seq_inv

As multiplication is faster than division, computing the inverse once (or at least once per call to this subprogram) should shave off some time.

Looping over k
Another thing to consider is that you have do loop increments at your disposal. That is to say, you can write your do k=1,CP_seq loop as

do k=1,CP_seq,2
! odd step
   Tfu(i)%elements = matmul(Tfu(i)%elements, Tud(i)%elements)
   Tfd(i)%elements = matmul(Tfd(i)%elements, Tdu(i)%elements)
! even step
   Tfu(i)%elements = matmul(Tfu(i)%elements, Tdu(i)%elements)
   Tfd(i)%elements = matmul(Tfd(i)%elements, Tud(i)%elements)
end do

Rather than wasting any time testing the two cases, do both in order and increment k by 2.

Improve this answer
\$\endgroup\$
  • 1
    \$\begingroup\$ Great first post! Hope you stick around :). \$\endgroup\$ – Veedrac Jan 12 '15 at 11:28
  • \$\begingroup\$ Thanks a lot Kyle for your useful answer. I didn't work more on that project but I'm about coming back to it in few days time. I'll consider your comments. \$\endgroup\$ – Roland Jan 13 '15 at 10:29
  • 1
    \$\begingroup\$ The suggestion in Anthony Scemamas answer is to explicitly write down the matrix multiplication, this is not covered by your proposed comparison. Writing the multiply explicitly down, might yield a boost in performance in comparison to loops. The more the compiler explicitly knows, the better it can optimize your code, and for a 2x2 matrix it definitely is not worth to leave the matmul in, in my opinion. \$\endgroup\$ – haraldkl Sep 4 '15 at 21:03
  • \$\begingroup\$ @harald: suppose op wants to change the code and use a 3x2 matrix or even some unknown-at-compile-time size. Anthony's method would either not work or would have you rewrite the entire subroutine whereas leaving the matmul in would change nothing. Using explicit loops could be a better option in this case, rather than matmul, but I think the readability is useful here. \$\endgroup\$ – Kyle Kanos Sep 5 '15 at 15:34
  • 1
    \$\begingroup\$ @KyleKanos, readability and performance are often conflicting goals, and OP asked for performance specifically. Often the parts of the code which take most of the time take only a relatively small amount of the overall program structure, and it is totally legit to optimize this on the expense of readability. TANSTAAFL... \$\endgroup\$ – haraldkl Sep 5 '15 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.