Increase performance of a spate of 2x2 matrix multiplication 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&gt; bath state
      L_pairs(i) = abs(conjg(Tfd(i)%elements(1, 1))*Tfu(i)%elements(1, 1) &amp;
           - 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
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 &lt; 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.