anchor.diff:  #error "Intel MKL does not support quad-precision arithmetic."
   #endif
   
  +#ifdef PCREO_GRAD_DESC
  +#undef PCREO_BFGS
  +#else
  +#define PCREO_BFGS
  +#endif
  +
   
   
   module constants !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
positive.code:         ! Multiply inverse hessian by negative gradient to obtain step direction
          call matrix_multiply_42(step_direction, inv_hess, force)
  #endif
          step_size = quadratic_line_search(points, energy, &
                  & step_direction, step_size)
          if (step_size == 0.0_rk) then
              call print_optimization_status
              call save_point_file(points, iteration_count)
              write(*,*) "Convergence has been achieved (up to numerical&
                      & round-off error). Exiting."
              stop
          end if
          new_points = points + step_size * step_direction
          call constrain_points(new_points)
          delta_points = new_points - points
          call riesz_energy_force(new_points, new_energy, new_force)
          delta_gradient = force - new_force
          call bfgs_update_inverse_hessian( &
                  & inv_hess, delta_points, delta_gradient)
          points = new_points
          energy = new_energy
          force = new_force
  
negative.code:     real(rk) :: last_print_time, last_save_time, cur_time
      integer :: iteration_count
  
  #if defined(PCREO_GRAD_DESC) .and. defined(PCREO_TRACK_ANGLE)
      real(rk) :: step_angle
      real(rk) :: new_force(d + 1, num_points)
  #elif defined(PCREO_BFGS)
      real(rk) :: new_points(d + 1, num_points)
      real(rk) :: new_energy, new_force(d + 1, num_points)
      real(rk) :: step_direction(d + 1, num_points)
      real(rk) :: inv_hess(d + 1, num_points, d + 1, num_points)
      real(rk), dimension(d + 1, num_points) :: delta_points, delta_gradient
  #ifdef PCREO_TRACK_ANGLE
      real(rk) :: step_angle
  