anchor.diff:          real(rk), intent(in) :: A(d + 1, num_points, d + 1, num_points)
           real(rk), intent(in) :: b(d + 1, num_points)
   
  -#ifdef PCREO_USE_MKL
  -#ifdef PCREO_SINGLE_PREC
  +#if defined(PCREO_USE_MKL) .and. defined(PCREO_SINGLE_PREC)
           call ssymv('U', num_vars, -1.0_rk, A, num_vars, b, 1, 0.0_rk, c, 1)
  -#else
  +#elif defined(PCREO_USE_MKL) .and. defined(PCREO_DOUBLE_PREC)
           call dsymv('U', num_vars, -1.0_rk, A, num_vars, b, 1, 0.0_rk, c, 1)
  -#endif
   #else
           integer :: i, j, k, l
   
  
positive.code:         real(rk), intent(inout) :: A(d + 1, num_points, d + 1, num_points)
          real(rk), intent(in) :: c, x(d + 1, num_points), y(d + 1, num_points)
  
  #ifdef PCREO_USE_MKL
  #ifdef PCREO_SINGLE_PREC
          call ssyr2('U', num_vars, c, x, 1, y, 1, A, num_vars)
  #else
          call dsyr2('U', num_vars, c, x, 1, y, 1, A, num_vars)
  #endif
  #else
          integer :: i, j, k, l
  
  
negative.code: #else
          call riesz_energy_force(points, energy, force)
  #endif
          call finish_iteration
      end do
  
  #elif defined(PCREO_BFGS)
  
      do
          step_size = quadratic_line_search( &
                  & points, energy, step_direction, step_size)
          call check_step_size
  