Debugging with gdb

23 MAR 2026

“Blowup of an aliased, non-energy-conserving model is God’s way of protecting you from believing a bad simulation.” - J. P. Boyd

Bugs are a natural part of the software development cycle. Be it a segfault, floating point exception, undefined behavior, or something else, bugs are bound to occur. Fixing them takes a combination of experience, intuition, luck, and tooling. Experience and intuition are built up over time. Luck, although useful, is not reliable. Tooling, like gdb and pygdb, exist to accelerate the debugging in tandem with experience, intuition, and luck.

gdb (the GNU debugger) is the standard tool for debugging C, C++, and Fortran code. It is very full featured, allowing stepping through code, inspecting variables, placing breakpoints, and even inspecting register values and assembly. It is also relatively opaque and unwelcoming to beginners. Even so, it is an invaluable part of a software developer’s toolkit.

We will breakdown the basics of gdb, building a simple workflow that covers most use cases, followed by an example.

gdb

Let’s assume that we have an executable main that we are debugging. We can invoke gdb as

gdb --args ./main # pass command line args as needed

and gdb will launch and load debugging symbols from the executable. Once gdb has launched, you will be dropped into an interactive prompt allowing us to run and interact with the program. Running the code is simple: r (or run). The code will run normally, either to completion or error.

When a crash occurs, the most useful first command is bt (or backtrace): This prints the call stack at the point of failure. It tells you how you got there — which functions were called, and in what order. In many cases, the top few frames of the backtrace are enough to identify the source of the bug.

If we need more information we can inspect the stack at the point of the crash directly. Say the crash happened in a file main.cpp on line 50 – gdb would report this as main.cpp:50. We can mark this location – set a break point – with b main.cpp:50. With that mark, we can run the program again with r. This time the code will run until the first time we reach main.cpp:50. Once here we can view the code around this line with l (or list). If we want to inspect the value of a variable we can use p varname (or print varname). This can quickly reveal errors, such as a zero denominator, without having to stick print statements through our code.

Summary

A minimal set of gdb commands that cover most debugging sessions:

These commands form a simple loop: run the program, break where needed, look at the code, and inspect values. In the next section we will look at an example of erroneous code.

Example

The following C++ code implements a classic Newtonian hydrostatic equilibrium solver with a polytropic equation of state:

\[\frac{dP}{dr} = -\frac{GM(r)\rho(r)}{r^2}\] \[P = K\rho^{1 + \frac{1}{n}}\]

These equations state that pressure balances gravity and that the pressure is a simple function of density. More information on hydrostatic equilibrium and polytropes can be found here and here, respectively.

By specifying a central density and integrating this equation outwards until the pressure becomes zero we can build up a simple model of a star. The code snippet below uses a classic 4th order Runge-Kutta integrator (see also Mike Zingale’s text).

#include <cassert>
#include <cmath>
#include <iostream>
#include <vector>

/**
 * @brief Polytropic pressure
 * p = k rho^(1 + 1/n)
 *
 * @param rho The density in cgs
 * @param k The polytropic constant
 * @param n The polytropic index
 * @return pressure
 */
auto pressure(const double rho, const double k, const double n) -> double {
  assert(k > 0 && "The polytropic constant k must be positive.");
  assert((n > 0 && n < 5) && "The polytropic index n must be in (0,5).");
  return k * std::pow(rho, 1.0 + 1.0 / n);
}

/**
 * @brief polytropic density from pressure
 * rho = std::pow( p / k, n/(n + 1))
 *
 * @param p The pressure in cgs
 * @param k The polytropic constant
 * @param n The polytropic index
 * @return density
 */
auto density_from_pressure(const double p, const double k, const double n)
    -> double {
  assert(k > 0 && "The polytropic constant k must be positive.");
  assert((n > 0 && n < 5) && "The polytropic index n must be in (0,5).");
  return std::pow(p / k, n / (n + 1.0));
}

/**
 * @brief Right hand side of hydrostatic equilibrium
 * dP/dr = - rho g
 * with g = G M / r^2
 */
auto rhs(const double m, const double r, const double p, const double k,
         const double n) -> double {
  constexpr double G = 6.674299999999999e-8; // cgs
  const double rho = density_from_pressure(p, k, n);
  return -G * m * rho / (r * r);
}

auto main(int argc, char **argv) -> int {
  // Declare polytrope values
  const double k = 1.0e15;
  const double n = 3.5;
  const double pi = 3.1415926535897931;
  const double fourpi = 4.0 * pi;

  // Set a central density and calculate central pressure
  const double rho_c = 150.0;
  const double p_c = pressure(rho_c, k, n);

  // Some io
  std::printf("# Building a star with:\n");
  std::printf("# k = %e\n", k);
  std::printf("# n = %f\n", n);
  std::printf("# rho_c = %e\n", rho_c);
  std::printf("# p_c = %e\n", p_c);

  // Now we use the classic 4th order Runge Kutta integration method
  // to integrate hydrostatic equilibrium dP/dr = -rho g outwards until
  // P ~ 0. We'll set our innermost radius to be 1km.
  const double r_c = 5e4;
  const double dr = 10.0 * r_c;

  std::vector<double> radius;
  std::vector<double> mass;
  std::vector<double> density;
  std::vector<double> pressure;

  radius.push_back(r_c);
  density.push_back(rho_c);
  pressure.push_back(p_c);
  mass.push_back((fourpi / 3.0) * std::pow(r_c, 3.0) * rho_c);

  // RK4
  int i = 0;
  const double p_threshold = 0.1; // below this we consider p ~ 0
  while (pressure.back() > p_threshold) {
    const double r = radius[i];
    const double p = pressure[i];
    const double rho = density[i];
    const double m = mass[i];

    const double dm = fourpi * rho * r * r * dr;

    const double k1 = dr * rhs(m, r, p, k, n);
    const double k2 = dr * rhs(m + 0.5 * dm, r + 0.5 * dr, p + 0.5 * k1, k, n);
    const double k3 = dr * rhs(m + 0.5 * dm, r + 0.5 * dr, p + 0.5 * k2, k, n);
    const double k4 = dr * rhs(m + dm, r + dr, p + k3, k, n);

    mass.push_back(m + dm);
    const double dp = (1.0 / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4);

    pressure.push_back(p + dp);
    density.push_back(density_from_pressure(pressure.back(), k, n));
    ++i;
  }

  // summary io
  const double M_sun = 1.98840987e+33; // g
  const double R_sun = 6.957e10;       // cm
  std::printf("\n# We have built a star:\n");
  std::printf("# M = %e [g] %f [Msun]\n", mass.back(), mass.back() / M_sun);
  std::printf("# R = %e [cm] %f [Rsun]\n", radius.back(),
              radius.back() / R_sun);
  return 0;
}

The code above can be compiled with g++ -g main.cpp -o main. Notice the -g – this embeds debugging information in the executable. However, the code contains a bug and crashes quickly. Spinning up gdb with gdb --args ./main and running with r we recover the crash. To locate the source we can get the stack backtrace with bt. We will see something like this:

Program received signal SIGABRT, Aborted.
0x00007ffff7989a2c in ?? () from /usr/lib/libc.so.6
#0  0x00007ffff7989a2c in ?? () from /usr/lib/libc.so.6
#1  0x00007ffff792f1a0 in raise () from /usr/lib/libc.so.6
#2  0x00007ffff79165fe in abort () from /usr/lib/libc.so.6
#3  0x00007ffff7c9a41f in std::__glibcxx_assert_fail(char const*, int, char const*, char const*) () from /usr/lib/libstdc++.so.6
#4  0x0000555555555ee6 in std::vector<double, std::allocator<double> >::operator[] (this=0x7fffffffd590, __n=1) at /usr/include/c++/15.2.1/bits/stl_vector.h:1263
#5  0x000055555555564a in main (argc=1, argv=0x7fffffffd758) at main.cpp:87

The first few lines regarding libc can usually be ignored. That is followed by something regarding std::vector and finally a location: main.cpp:86. This is an out of bounds check. In this build, the standard library catches it and aborts with an assertion failure (rather than a raw segmentation fault). Looking at line 87

86  while (pressure.back() > p_threshold) {
87    const double r = radius[i];
88    const double p = pressure[i];

we see that it is accessing the i-th component of the radius vector. We can investigate further:

Add the following:

radius.push_back(r + dr);

To the end of the integration loop, recompile, and the code runs successfully!

# Building a star with:
# k = 1.000000e+15
# n = 3.500000
# rho_c = 1.500000e+02
# p_c = 6.278149e+17

# We have built a star:
# M = 6.527818e+33 [g] 3.282934 [Msun]
# R = 1.166306e+11 [cm] 1.676449 [Rsun]

The hydrostatic equilibrium solver successfully constructs a star of roughly 3 solar masses – the expected result!

Summary

We covered the basics of gdb, building up a simple but powerful workflow for debugging code. With this is hand we applied it to a real world example of buggy code, finding and fixing a segfault in a hydrostatic equilibrium solver.

Debugging with gdb | Brandon’s website