6. Just enough C to read ASYNCH

Code readersPython helps30 minutes

This page teaches the C you need for this codebase, using real code from ASYNCH, and the real bugs listed in Known issues as examples of what goes wrong. If you know Python, the comparisons in italics will help.

6.1 How a C program is built

Python reads your .py files at run time. C is translated to machine code beforehand.

SOURCE advance.c riversys.c … 60 .c files structs.h, … headers: declarations Preprocessor#include, #define, #if Compilereach .c → one .o file Linker+ MPI, HDF5, libpq RESULT asynchthe program libasynch.sothe library make runs these steps again only for the files that changed
How C becomes a program: preprocessor, compiler, linker

  1. Preprocessor. Lines starting with # are text substitutions done before compiling:

    • #include <structs.h> pastes the content of that header file here.

    • #define ASYNCH_SLEEP sleep replaces every ASYNCH_SLEEP with sleep.

    • #if defined(HAVE_POSTGRESQL) ... #endif keeps the code only if the macro is defined (by configure, in config.h). This is how optional features are switched on and off.

  2. Compiler. Each .c file becomes an object file .o, independently of the others. A .h header file contains declarations (“a function named Advance exists and takes these arguments”) so that other .c files can call it.

  3. Linker. Joins the .o files and the libraries (MPI, HDF5, libpq) into the program asynch.

make runs these steps for the files that changed. configure (generated by autotools) finds the libraries and writes the Makefiles. src/Makefile.am lists which .c files are compiled. A .c file that is not listed there is never compiled (that is how issue M-01 was found).

6.2 Types, variables, const

unsigned int dim = link_i->dim;      // integer >= 0, 32 bits
unsigned short model_uid;            // integer 0..65535 (16 bits!), see issue B-06
double h = link_i->h;                // 64-bit floating point (Python's float)
float value;                         // 32-bit float (~7 digits)
bool print_flag = false;             // from <stdbool.h>

Unlike Python, every variable has a fixed type, and integer types have a limited range. Storing 400 000 in an unsigned short silently gives 400 000 mod 65 536 = 6 784.

const double * const global_params means: a pointer that won’t change, to doubles that won’t be changed through it. It is a promise that the function only reads them.

6.3 Pointers and arrays: the most important concept

A pointer is a variable holding a memory address. double *y is “the address of a double”, usually the first of an array of doubles.

double *y_0 = link_i->my->list.tail->y_approx;   // y_0 points to the last state of the link
double q   = y_0[0];                             // first element (discharge)
double s_p = y_0[1];                             // second element

y[i] means “the double located i positions after the address y”. C never checks that i is inside the array.

y the addressof y[0] q y[0] s_p y[1] s_t y[2] s_s y[3] s_precip y[4] V_r y[5] q_b y[6] ? y[7] the 7 states of model 254: the memory reserved for y not part of y: someone else's data C does not check the index: states[7] = 0 silently overwrites the neighbouring memory (bug B-01).
An array of 7 states in memory: y[7] is past its end, in memory that belongs to something else
Reading y[7] from a 7-element array reads whatever happens to be next in memory, and writing there corrupts it. This is exactly issue B-01:

void OutputConstraints_Model256_Hdf5(double* states)   // written for 8 states
{   ...
    if (states[7] < 1e-12)      // model 254 has only states[0..6]: out of bounds!
        states[7] = 0;
}

2-D data is often stored in a 1-D array. The parents’ states in model254:

for (i = 0; i < num_parents; i++)
    ans[0] += y_p[i * dim];         // state 0 of parent i: row i, column 0

y_p is a table with one row per parent, flattened. Row i starts at i * row_length. If the row length used here (dim) differs from the one used when filling the table (max_dim), you read the wrong numbers. That is issue S-01.

&x gives the address of x, and *p gives the value at address p:

int ReadLine(..., unsigned int *flag);  // the function can write into the caller's variable
ReadLine(..., &flag);

This is how C functions “return” several values.

6.4 Structs and ->

A struct groups variables, like a Python class with only attributes:

struct Link {
    unsigned int ID;
    double *params;
    Link **parents;          // array of pointers to the parent links
    Link *child;             // pointer to the downstream link (NULL at the outlet)
    ...
};

link.ID accesses a field of a struct, and link_i->ID accesses a field through a pointer to a struct (the same as (*link_i).ID). You will see -> everywhere: link_i->my->list.tail->y_approx walks from the link to its private data, to its solution list, to the last node, to that node’s state array.

6.5 Function pointers: how models are plugged in

link->differential = &model254;      // in definitions.c (InitRoutines)
...
link_i->differential(t, y, dim, ..., ans);   // in the stepper: calls model254(...)

The solver does not know which model it runs. It calls whatever function the link points to. Like passing a function as an argument in Python. To find what a call really executes, search where the pointer is assigned (grep -n "differential =" src/models/definitions.c).

6.6 Memory: malloc / free

Python frees memory for you. C does not.

double *buf = malloc(n * sizeof(double));   // reserve n doubles (contents are garbage!)
double *z   = calloc(n, sizeof(double));    // same, but zero-filled
buf = realloc(buf, m * sizeof(double));     // resize
free(buf);                                  // give it back, exactly once

The classic errors, all present in ASYNCH:

error

example in ASYNCH

writing past the end of a block (buffer overflow)

B-01

freeing twice (double free)

B-03: fclose twice on the same FILE*

freeing something that was not malloced

B-05: free(&error->abstol) instead of free(error->abstol)

never freeing (leak)

B-07: data_storage

using memory that was never initialised

B-12: malloc + too-short .uini file

These errors often do not crash immediately. They corrupt something that fails later, somewhere else. That is why we use AddressSanitizer (see 09_reproducibility.md): it stops at the exact faulty line.

6.7 switch falls through

switch (model_uid) {
    case 254:
        globals->OutputConstrainsHdf5 = &OutputConstraints_Model254_Hdf5;
                                       // no break: execution CONTINUES into case 256!
    case 256:
        globals->OutputConstrainsHdf5 = &OutputConstraints_Model256_Hdf5;
        break;                         // leave the switch
}

Unlike Python’s match, a C case is only a jump label. Every case needs its own break;. This one missing line is issue B-01.

6.8 Reading files: always check the return value

if (fscanf(initdata, "%lf", y_0_backup + i) == 0)   // B-12: wrong check

fscanf returns the number of items read: 1 on success, 0 if the text is not a number, EOF (−1) at end of file. The correct test is != 1. About 75 calls in ASYNCH do not check at all (B-11).

6.9 MPI in five functions

With mpirun -n 4 asynch x.gbl, 4 independent copies of the program run. Each knows:

MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);   // who am I: 0..3   (global `my_rank` in ASYNCH)
MPI_Comm_size(MPI_COMM_WORLD, &np);        // how many are we: 4 (global `np`)
MPI_Send(buf, n, MPI_DOUBLE, dest, tag, MPI_COMM_WORLD);        // send n doubles to `dest`
MPI_Recv(buf, n, MPI_DOUBLE, src, tag, MPI_COMM_WORLD, &status);  // receive them
MPI_Barrier(MPI_COMM_WORLD);               // wait until everybody reaches this line

Nothing is shared: if process 1 needs a value from process 0, it must be sent. So code often looks like if (my_rank == 0) { read file; send } else { receive }. A consequence: a bug can depend on the number of processes (B-01 only crashes on 1 process, B-02 only shows with several).

6.10 Debugging tools

# a debug build (no optimisation, debug symbols)
../configure CFLAGS="-O0 -g" && make
# run under the debugger, 1 process
gdb --args ./src/asynch ../examples/test.gbl
(gdb) break model254          # stop when the model function is called
(gdb) run
(gdb) print y_i[0]            # inspect a variable
(gdb) bt                      # who called me? (backtrace)
(gdb) next / step / continue

And the quickest tool of all: printf("q=%g at t=%g\n", q, t); followed by a rebuild.

6.11 Anatomy of a model function

void model254(
    double t,                           // current time [min]
    const double * const y_i,           // state of this link, y_i[0..dim-1]
    unsigned int dim,                   // number of states (7)
    const double * const y_p,           // states of the parents, flattened (see 6.3)
    unsigned short num_parents,
    unsigned int max_dim,
    const double * const global_params, // from the .gbl
    const double * const params,        // this link's parameters (after precalculation)
    const double * const forcing_values,// forcing_values[0] = rain [mm/h], [1] = e_pot [mm/month], ...
    const QVSData * const qvs,          // discharge-storage table (dams), unused here
    int state,                          // discontinuity state, unused here
    void* user,                         // free pointer for custom data
    double *ans)                        // OUTPUT: ans[k] = d y_i[k] / dt

Everything a model needs comes in through the arguments, and its only output is ans. That makes model functions the easiest part of ASYNCH to read, test and change.