Reducing compile resources using GSL for an ODE system

Question

I'm using C++ and GSL (Gnu Scientific Library) to integrate a large ODE system with 4503 equations. As the system is stiff, I need to setup a Jacobian for the integration and the resulting matrix has about 20 million entries.

To initialize the ODE system and the Jacobian, GSL requires the use of a functions with specific signatures, that is,

int system(double t, const double y[], double dydt[], void *params);
int jac (double t, const double y[], double *dfdy, double dfdt[], void *params);

Where dydt[] and *dfdy are the arrays for the system and Jacobian, respectively.

I have declared these functions in two separate .cpp files. For the system, I have something like

// ODESystem.cpp

int system(double t, const double y[], double dydt[], void *params) {
    (void)(t);
    double l = *(static_cast<double*>(params));

    dydt[0] = 0;

    // lots of assignment statements like the above

    return GSL_SUCCESS;
}

and for the Jacobian I have something like

// ODEJacobian.cpp

int jacobian(double t, const double y[], double *dfdy, double dfdt[], void *params) {
    (void)(t);
    (void)(y);
    double l = *(static_cast<double*>(params));

    dfdy[0] = 0;

    // lots of assignment statements like the above

    return GSL_SUCCESS;
}

So far, GCC is able to handle ODESystem.cpp quite well, but it takes a very long time (and about 45Gs of ram) to compile ODEJacobian.cpp, (which is, i guess, understandable due to the size of the expression tree for the assignment operations in this monster) but it eventually does so without errors or warnings.

Clang seems unable to cope with the size of the function, and while it exits normally, it issues a warning about "reaching the end of a non void function and finding no return statement", which seems to me as if clang "gives up" compiling somewhere before the Jacobian function ends. For smaller sizes, however (about 5,745,609 assignments) Clang does fine and issues no warnings.

My question is: What can be done in this case to reduce the amount of resources (time and memory) necessary to compile the Jacobian while respecting the constraints imposed by GSL with either (or both) GCC and Clang?

Answer 1

Don't write assignments for this. Instead copy the data into the array.

I would suggest reading the values first from a file into a std::vector<double> vec; at program startup. You can then keep reusing its contents. This also makes it easier to modify the values without recompilation and the compiler won't have to do so much unnecessary work.

Then

std::copy(vec.begin(), vec.end(), dydt);

This requires #include<algorithm> for std::copy and #include<vector> for std::vector.