calculate_log_der_matrix Subroutine

      subroutine calculate_log_der_matrix(step_, number_of_channels_,          &
         t_matrix_minus_, t_matrix_, t_matrix_plus_, r_matrix_, r_matrix_plus_,&
         log_der_matrix_)
         !! calculates the log-derivative matrix from
         !! \begin{equation}
         !! {Y}_{\rm N} = \frac{1}{h} \Biggl(\Bigl(\frac{1}{2}\mathbf{I}-{T}_{\rm{N}+1}\Bigr)
         !! \Bigl(\mathbf{I}-{T}_{\rm{N}+1}\Bigr)^{-1} {R}_{\rm{N}+1} -
         !! \Bigl(\frac{1}{2}\mathbf{I}-{T}_{\rm{N}-1}\Bigr)
         !! \Bigl(\mathbf{I}-\mathbf{T}_{\rm{N}-1}\Bigr)^{-1}\mathbf{R}_{\rm{N}}^{-1}
         !! \Biggr)\Bigl(\mathbf{I}-{T}_{\rm{N}}\Bigr) 
         !! \end{equation}
         !! called by numerov at the end of the propagation
         !---------------------------------------------------------------------!
         real(dp), intent(in) :: step_
            !! propagator step
         integer(int32), intent(in) :: number_of_channels_
            !! number of scattering channels in the block
         real(dp), intent(in) :: t_matrix_minus_(number_of_channels_, number_of_channels_)
            !! T-matrix at R_{max - 1}
         real(dp), intent(in) :: t_matrix_(number_of_channels_, number_of_channels_)
            !! T-matrix at R_{max}
         real(dp), intent(in) :: t_matrix_plus_(number_of_channels_, number_of_channels_)
            !! T-matrix at R_{max + 1} 
         real(dp), intent(in) :: r_matrix_(number_of_channels_, number_of_channels_)
            !! R-matrix at R_{max} 
         real(dp), intent(in) :: r_matrix_plus_(number_of_channels_, number_of_channels_)
            !! R-matrix at R_{max + 1}
         real(dp), intent(inout) :: log_der_matrix_(number_of_channels_, number_of_channels_)
            !! log-derivative matrix
         !---------------------------------------------------------------------!
         integer(int32) :: channel_index_1_
         real(dp) :: matrix_a_(number_of_channels_,number_of_channels_),       &
            matrix_b_(number_of_channels_,number_of_channels_),                &
            matrix_c_(number_of_channels_,number_of_channels_),                &
            matrix_d_(number_of_channels_,number_of_channels_),                &
            matrix_e_(number_of_channels_,number_of_channels_),                &
            matrix_ab_(number_of_channels_,number_of_channels_),               &
            matrix_cd_(number_of_channels_,number_of_channels_),               &
            left_matrix_(number_of_channels_,number_of_channels_),             &
            right_matrix_(number_of_channels_,number_of_channels_),            &
            working_r_matrix_(number_of_channels_,number_of_channels_),        &
            matrix_difference_(number_of_channels_,number_of_channels_)
         !---------------------------------------------------------------------!
         log_der_matrix_ = 0
         !---------------------------------------------------------------------!
         ! First bracket: (1/2 I - T_{N+1})
         !---------------------------------------------------------------------!
         matrix_a_ = - t_matrix_plus_
         call add_scalar_to_diagonal(matrix_a_, 0.5_dp)
         !---------------------------------------------------------------------!
         ! Second bracket: (I - T_{N+1})^{-1}
         !---------------------------------------------------------------------!
         matrix_b_ = - t_matrix_plus_
         call add_scalar_to_diagonal(matrix_b_, 1.0_dp)
         call invert_symmetric_matrix(matrix_b_)
         call fill_symmetric_matrix(matrix_b_, "u")
         !---------------------------------------------------------------------!
         ! Third bracket: (1/2 I - T_{N-1})
         !---------------------------------------------------------------------!
         matrix_c_ = - t_matrix_minus_
         call add_scalar_to_diagonal(matrix_c_, 0.5_dp)
         !---------------------------------------------------------------------!
         ! Fourth bracket: (I - T_{N-1})^{-1}
         !---------------------------------------------------------------------!
         matrix_d_ = - t_matrix_minus_
         call add_scalar_to_diagonal(matrix_d_, 1.0_dp)
         call invert_symmetric_matrix(matrix_d_)
         call fill_symmetric_matrix(matrix_d_, "u")
         !---------------------------------------------------------------------!
         ! Last bracket: (I - T_{N})
         !---------------------------------------------------------------------!
         matrix_e_ = - t_matrix_
         call add_scalar_to_diagonal(matrix_e_, 1.0_dp)
         !---------------------------------------------------------------------!
         ! Copy R_{N} to another matrix (r_matrix_ is protected as intent(in))
         !---------------------------------------------------------------------!
         working_r_matrix_ = r_matrix_
         call invert_symmetric_matrix(working_r_matrix_)
         call fill_symmetric_matrix(working_r_matrix_, "u")
         !---------------------------------------------------------------------!
         ! The first term in the large bracket:
         ! (1/2 I - T_{N+1}) (I - T_{N+1})^{-1} R_{N+1}
         !---------------------------------------------------------------------!
         CALL DGEMM('N','N',number_of_channels_,number_of_channels_,           &
            number_of_channels_,1.0_dp,matrix_a_,number_of_channels_,matrix_b_,&
            number_of_channels_,0.0_dp,matrix_ab_,number_of_channels_)
         CALL DGEMM('N','N',number_of_channels_,number_of_channels_,           &
            number_of_channels_,1.0_dp,matrix_ab_,number_of_channels_,         &
            r_matrix_plus_,number_of_channels_,0.0_dp,left_matrix_,number_of_channels_)
         !---------------------------------------------------------------------!
         ! The second term in the large bracket:
         ! (1/2 I - T_{N-1}) (I - T_{N-1})^{-1} R_{N}^{-1}
         !---------------------------------------------------------------------!
         CALL DGEMM('N','N',number_of_channels_,number_of_channels_,           &
            number_of_channels_,1.0_dp,matrix_c_,number_of_channels_,matrix_d_,&
            number_of_channels_,0.0_dp,matrix_cd_,number_of_channels_)
         CALL DGEMM('N','N',number_of_channels_,number_of_channels_,           &
            number_of_channels_,1.0_dp,matrix_cd_,number_of_channels_,         &
            working_r_matrix_,number_of_channels_,0.0_dp,right_matrix_,        &
            number_of_channels_)
         !---------------------------------------------------------------------!
         ! Substract the two terms in the large bracket
         !---------------------------------------------------------------------!
         matrix_difference_ = left_matrix_ - right_matrix_
         !---------------------------------------------------------------------!
         ! Multiply the product by (I - T_{N})
         !---------------------------------------------------------------------!
         CALL DGEMM('N','N',number_of_channels_,number_of_channels_,           &
            number_of_channels_,1.0_dp/step_,matrix_difference_,               &
            number_of_channels_,matrix_e_,number_of_channels_,0.0_dp,          &
            log_der_matrix_,number_of_channels_)
         !---------------------------------------------------------------------!
      end subroutine calculate_log_der_matrix