anchor.diff:          points = new_points
           energy = new_energy
           force = new_force
  -
  -#endif
  -
  -        iteration_count = iteration_count + 1
  -        cur_time = current_time()
  -        if (cur_time - last_print_time >= print_time) then
  -            call print_optimization_status
  -            last_print_time = cur_time
  +#ifdef PCREO_TRACK_ANGLE
  +        call matrix_multiply_42(new_step_direction, inv_hess, force)
  +        if (iteration_count > 0) then
  +            step_angle = sum(step_direction * new_step_direction) / &
  +                    & (norm2(step_direction) * norm2(new_step_direction))
           end if
  -        if (cur_time - last_save_time >= save_time) then
  -            call save_point_file(points, iteration_count)
  -            last_save_time = cur_time
  -        end if
  +        step_direction = new_step_direction
  +#else
  +        call matrix_multiply_42(step_direction, inv_hess, force)
  +#endif
  +        call finish_iteration
       end do
  +
  +#endif
  +
   
   contains
   
  
positive.code:     call print_table_header
      call print_optimization_status
  
      do
  #if defined(PCREO_GRAD_DESC)
  
          step_size = quadratic_line_search(points, energy, force, step_size)
          call check_step_size
          points = points + step_size * force
  
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
  