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