rxd_extracellular.cpp

Go to the documentation of this file.

#include <../../nrnconf.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include "grids.h"
#include "rxd.h"
#include <nrnwrap_Python.h>
#include <cmath>
 
#define loc(x, y, z) ((z) + (y) *grid->size_z + (x) *grid->size_z * grid->size_y)
 
static void ecs_refresh_reactions(int);
/*
    Globals
*/
 
/*Store the onset/offset for each threads reaction tasks*/
ReactGridData* threaded_reactions_tasks;
 
extern int NUM_THREADS;
extern pthread_t* Threads;
extern TaskQueue* AllTasks;
extern double* t_ptr;
extern double* states;
 
Reaction* ecs_reactions = NULL;
 
int states_cvode_offset;
 
/*Update the global array of reaction tasks when the number of reactions
 *or threads change.
 *n - the old number of threads - use to free the old threaded_reactions_tasks*/
static void ecs_refresh_reactions(const int n) {
    int k;
    if (threaded_reactions_tasks != NULL) {
        for (k = 0; k < NUM_THREADS; k++) {
            SAFE_FREE(threaded_reactions_tasks[k].onset);
            SAFE_FREE(threaded_reactions_tasks[k].offset);
        }
        SAFE_FREE(threaded_reactions_tasks);
    }
    threaded_reactions_tasks = create_threaded_reactions(n);
}
 
void set_num_threads_3D(const int n) {
    Grid_node* grid;
    if (n != NUM_THREADS) {
        for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
            grid->set_num_threads(n);
        }
    }
    ecs_refresh_reactions(n);
}
 
 
/*Removal all reactions*/
void clear_rates_ecs(void) {
    Reaction *r, *tmp;
    Grid_node* grid;
    ECS_Grid_node* g;
 
    for (r = ecs_reactions; r != NULL; r = tmp) {
        SAFE_FREE(r->species_states);
        if (r->subregion) {
            SAFE_FREE(r->subregion);
        }
 
        tmp = r->next;
        SAFE_FREE(r);
    }
    ecs_reactions = NULL;
 
    ecs_refresh_reactions(NUM_THREADS);
    for (Grid_node* grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        g = dynamic_cast<ECS_Grid_node*>(grid);
        if (g)
            g->clear_multicompartment_reaction();
    }
}
 
/*Create a reaction
 * list_idx - the index for the linked list in the global array Parallel_grids
 * grid_id - the grid id within the linked list
 * ECSReactionRate - the reaction function
 * subregion - either NULL or a boolean array the same size as the grid
 *  indicating if reaction occurs at a specific voxel
 */
Reaction* ecs_create_reaction(int list_idx,
                              int num_species,
                              int num_params,
                              int* species_ids,
                              ECSReactionRate f,
                              unsigned char* subregion,
                              uint64_t* mc3d_start_indices,
                              int mc3d_region_size,
                              double* mc3d_mults) {
    Grid_node* grid;
    Reaction* r;
    int i, j;
 
    r = (Reaction*) malloc(sizeof(Reaction));
    assert(r);
    r->reaction = f;
    /*place reaction on the top of the stack of reactions*/
    r->next = ecs_reactions;
    ecs_reactions = r;
 
    for (grid = Parallel_grids[list_idx], i = 0; grid != NULL; grid = grid->next, i++) {
        /* Assume all species have the same grid */
        if (i == species_ids[0]) {
            if (mc3d_region_size > 0) {
                r->subregion = NULL;
                r->region_size = mc3d_region_size;
                r->mc3d_indices_offsets = (uint64_t*) malloc(sizeof(uint64_t) *
                                                             (num_species + num_params));
                memcpy(r->mc3d_indices_offsets,
                       mc3d_start_indices,
                       sizeof(uint64_t) * (num_species + num_params));
                r->mc3d_mults = (double**) malloc(sizeof(double*) * (num_species + num_params));
                int mult_idx = 0;
                for (int row = 0; row < num_species + num_params; row++) {
                    r->mc3d_mults[row] = (double*) malloc(sizeof(double) * mc3d_region_size);
                    for (int c = 0; c < mc3d_region_size; c++, mult_idx++) {
                        r->mc3d_mults[row][c] = mc3d_mults[mult_idx];
                    }
                }
 
 
            } else if (subregion == NULL) {
                r->subregion = NULL;
                r->region_size = grid->size_x * grid->size_y * grid->size_z;
                r->mc3d_indices_offsets = NULL;
            } else {
                for (r->region_size = 0, j = 0; j < grid->size_x * grid->size_y * grid->size_z; j++)
                    r->region_size += subregion[j];
                r->subregion = subregion;
                r->mc3d_indices_offsets = NULL;
            }
        }
    }
 
    r->num_species_involved = num_species;
    r->num_params_involved = num_params;
    r->species_states = (double**) malloc(sizeof(Grid_node*) * (num_species + num_params));
    assert(r->species_states);
 
    for (i = 0; i < num_species + num_params; i++) {
        /*Assumes grids are the same size (no interpolation)
         *Assumes all the species will be in the same Parallel_grids list*/
        for (grid = Parallel_grids[list_idx], j = 0; grid != NULL; grid = grid->next, j++) {
            if (j == species_ids[i])
                r->species_states[i] = grid->states;
        }
    }
 
    return r;
}
 
/*ics_register_reaction is called from python it creates a reaction with
 * list_idx - the index for the linked list in the global array Parallel_grids
 *  (currently this is always 0)
 * grid_id - the grid id within the linked list - this corresponds to species
 * ECSReactionRate - the reaction function
 */
extern "C" void ics_register_reaction(int list_idx,
                                      int num_species,
                                      int num_params,
                                      int* species_id,
                                      uint64_t* mc3d_start_indices,
                                      int mc3d_region_size,
                                      double* mc3d_mults,
                                      ECSReactionRate f) {
    ecs_create_reaction(list_idx,
                        num_species,
                        num_params,
                        species_id,
                        f,
                        NULL,
                        mc3d_start_indices,
                        mc3d_region_size,
                        mc3d_mults);
    ecs_refresh_reactions(NUM_THREADS);
}
 
/*ecs_register_reaction is called from python it creates a reaction with
 * list_idx - the index for the linked list in the global array Parallel_grids
 *  (currently this is always 0)
 * grid_id - the grid id within the linked list - this corresponds to species
 * ECSReactionRate - the reaction function
 */
extern "C" void ecs_register_reaction(int list_idx,
                                      int num_species,
                                      int num_params,
                                      int* species_id,
                                      ECSReactionRate f) {
    ecs_create_reaction(list_idx, num_species, num_params, species_id, f, NULL, NULL, 0, NULL);
    ecs_refresh_reactions(NUM_THREADS);
}
 
 
/*register_subregion_reaction is called from python it creates a reaction with
 * list_idx - the index for the linked list in the global array Parallel_grids
 *  (currently this is always 0)
 * grid_id - the grid id within the linked list - this corresponds to species
 * my_subregion - a boolean array indicating the voxels where a reaction occurs
 * ECSReactionRate - the reaction function
 */
void ecs_register_subregion_reaction_ecs(int list_idx,
                                         int num_species,
                                         int num_params,
                                         int* species_id,
                                         unsigned char* my_subregion,
                                         ECSReactionRate f) {
    ecs_create_reaction(
        list_idx, num_species, num_params, species_id, f, my_subregion, NULL, 0, NULL);
    ecs_refresh_reactions(NUM_THREADS);
}
 
 
/*create_threaded_reactions is used to generate evenly distribution of reactions
 * across NUM_THREADS
 * returns ReactGridData* which can be used as the global
 *  threaded_reactions_tasks
 */
ReactGridData* create_threaded_reactions(const int n) {
    unsigned int i;
    int load, k;
    int react_count = 0;
    Reaction* react;
 
    for (react = ecs_reactions; react != NULL; react = react->next)
        react_count += react->region_size;
 
    if (react_count == 0)
        return NULL;
 
    /*Divide the number of reactions between the threads*/
    const int tasks_per_thread = (react_count) / n;
    ReactGridData* tasks = (ReactGridData*) calloc(sizeof(ReactGridData), n);
    const int extra = react_count % n;
 
    tasks[0].onset = (ReactSet*) malloc(sizeof(ReactSet));
    tasks[0].onset->reaction = ecs_reactions;
    tasks[0].onset->idx = 0;
 
    for (k = 0, load = 0, react = ecs_reactions; react != NULL; react = react->next) {
        for (i = 0; i < react->region_size; i++) {
            if (!react->subregion || react->subregion[i])
                load++;
            if (load >= tasks_per_thread + (extra > k)) {
                tasks[k].offset = (ReactSet*) malloc(sizeof(ReactSet));
                tasks[k].offset->reaction = react;
                tasks[k].offset->idx = i;
 
                if (++k < n) {
                    tasks[k].onset = (ReactSet*) malloc(sizeof(ReactSet));
                    tasks[k].onset->reaction = react;
                    tasks[k].onset->idx = i + 1;
                    load = 0;
                }
            }
            if (k == n - 1 && react->next == NULL) {
                tasks[k].offset = (ReactSet*) malloc(sizeof(ReactSet));
                tasks[k].offset->reaction = react;
                tasks[k].offset->idx = i;
            }
        }
    }
    return tasks;
}
 
/*ecs_do_reactions takes ReactGridData which defines the set of reactions to be
 * performed. It calculate the reaction based on grid->old_states and updates
 * grid->states
 */
void* ecs_do_reactions(void* dataptr) {
    ReactGridData task = *(ReactGridData*) dataptr;
    unsigned char started = FALSE, stop = FALSE;
    unsigned int i, j, k, n, start_idx, stop_idx, offset_idx;
    double temp, ge_value, val_to_set;
    double dt = *dt_ptr;
    Reaction* react;
 
    double* states_cache = NULL;
    double* params_cache = NULL;
    double* states_cache_dx = NULL;
    double* results_array = NULL;
    double* results_array_dx = NULL;
    double* mc_mults_array = NULL;
    double dx = FLT_EPSILON;
    double pd;
    MAT* jacobian;
    VEC* x;
    VEC* b;
    PERM* pivot;
 
    for (react = ecs_reactions; react != NULL; react = react->next) {
        // TODO: This is bad. Need to refactor
        if (react->mc3d_indices_offsets != NULL) {
            /*find start location for this thread*/
            if (started || react == task.onset->reaction) {
                if (!started) {
                    started = TRUE;
                    start_idx = task.onset->idx;
                } else {
                    start_idx = 0;
                }
 
                if (task.offset->reaction == react) {
                    stop_idx = task.offset->idx;
                    stop = TRUE;
                } else {
                    stop_idx = react->region_size - 1;
                    stop = FALSE;
                }
                if (react->num_species_involved == 0)
                    continue;
                /*allocate data structures*/
                jacobian = m_get(react->num_species_involved, react->num_species_involved);
                b = v_get(react->num_species_involved);
                x = v_get(react->num_species_involved);
                pivot = px_get(jacobian->m);
                states_cache = (double*) malloc(sizeof(double) * react->num_species_involved);
                params_cache = (double*) malloc(sizeof(double) * react->num_params_involved);
                states_cache_dx = (double*) malloc(sizeof(double) * react->num_species_involved);
                results_array = (double*) malloc(react->num_species_involved * sizeof(double));
                results_array_dx = (double*) malloc(react->num_species_involved * sizeof(double));
                mc_mults_array = (double*) malloc(react->num_species_involved * sizeof(double));
                for (i = start_idx; i <= stop_idx; i++) {
                    if (!react->subregion || react->subregion[i]) {
                        // TODO: this assumes grids are the same size/shape
                        //      add interpolation in case they aren't
                        for (j = 0; j < react->num_species_involved; j++) {
                            offset_idx = i + react->mc3d_indices_offsets[j];
                            states_cache[j] = react->species_states[j][offset_idx];
                            states_cache_dx[j] = react->species_states[j][offset_idx];
                            mc_mults_array[j] = react->mc3d_mults[j][i];
                        }
                        MEM_ZERO(results_array, react->num_species_involved * sizeof(double));
                        for (k = 0; j < react->num_species_involved + react->num_params_involved;
                             k++, j++) {
                            offset_idx = i + react->mc3d_indices_offsets[j];
                            params_cache[k] = react->species_states[j][offset_idx];
                            mc_mults_array[k] = react->mc3d_mults[j][i];
                        }
                        react->reaction(states_cache, params_cache, results_array, mc_mults_array);
 
                        for (j = 0; j < react->num_species_involved; j++) {
                            states_cache_dx[j] += dx;
                            MEM_ZERO(results_array_dx,
                                     react->num_species_involved * sizeof(double));
                            react->reaction(states_cache_dx,
                                            params_cache,
                                            results_array_dx,
                                            mc_mults_array);
                            v_set_val(b, j, dt * results_array[j]);
 
                            for (k = 0; k < react->num_species_involved; k++) {
                                pd = (results_array_dx[k] - results_array[k]) / dx;
                                m_set_val(jacobian, k, j, (j == k) - dt * pd);
                            }
                            states_cache_dx[j] -= dx;
                        }
                        // solve for x
                        if (react->num_species_involved == 1) {
                            react->species_states[0][i] += v_get_val(b, 0) /
                                                           m_get_val(jacobian, 0, 0);
                        } else {
                            // find entry in leftmost column with largest absolute value
                            // Pivot
                            for (j = 0; j < react->num_species_involved; j++) {
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    if (abs(m_get_val(jacobian, j, j)) <
                                        abs(m_get_val(jacobian, k, j))) {
                                        for (n = 0; n < react->num_species_involved; n++) {
                                            temp = m_get_val(jacobian, j, n);
                                            m_set_val(jacobian, j, n, m_get_val(jacobian, k, n));
                                            m_set_val(jacobian, k, n, temp);
                                        }
                                    }
                                }
                            }
 
                            for (j = 0; j < react->num_species_involved - 1; j++) {
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    ge_value = m_get_val(jacobian, k, j) /
                                               m_get_val(jacobian, j, j);
                                    for (n = 0; n < react->num_species_involved; n++) {
                                        val_to_set = m_get_val(jacobian, k, n) -
                                                     ge_value * m_get_val(jacobian, j, n);
                                        m_set_val(jacobian, k, n, val_to_set);
                                    }
                                    v_set_val(b, k, v_get_val(b, k) - ge_value * v_get_val(b, j));
                                }
                            }
 
                            for (j = react->num_species_involved - 1; j + 1 > 0; j--) {
                                v_set_val(x, j, v_get_val(b, j));
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    if (k != j) {
                                        v_set_val(x,
                                                  j,
                                                  v_get_val(x, j) -
                                                      m_get_val(jacobian, j, k) * v_get_val(x, k));
                                    }
                                }
                                v_set_val(x, j, v_get_val(x, j) / m_get_val(jacobian, j, j));
                            }
                            for (j = 0; j < react->num_species_involved; j++) {
                                // I think this should be something like
                                // react->species_states[j][mc3d_indices[i]] += v_get_val(x,j);
                                // Since the grid has a uniform discretization, the mc3d_indices
                                // should be the same length. So just need to access the correct
                                // mc3d_indices[i] maybe do two lines?: index =
                                // react->species_indices[j][i] react->species_states[j][index] +=
                                // v_get_val(x,j);
                                offset_idx = i + react->mc3d_indices_offsets[j];
                                react->species_states[j][offset_idx] += v_get_val(x, j);
                            }
                        }
                    }
                }
                m_free(jacobian);
                v_free(b);
                v_free(x);
                px_free(pivot);
 
                SAFE_FREE(states_cache);
                SAFE_FREE(states_cache_dx);
                SAFE_FREE(params_cache);
                SAFE_FREE(results_array);
                SAFE_FREE(results_array_dx);
                SAFE_FREE(mc_mults_array);
 
                if (stop)
                    return NULL;
            }
        } else {
            /*find start location for this thread*/
            if (started || react == task.onset->reaction) {
                if (!started) {
                    started = TRUE;
                    start_idx = task.onset->idx;
                } else {
                    start_idx = 0;
                }
 
                if (task.offset->reaction == react) {
                    stop_idx = task.offset->idx;
                    stop = TRUE;
                } else {
                    stop_idx = react->region_size - 1;
                    stop = FALSE;
                }
                if (react->num_species_involved == 0)
                    continue;
                /*allocate data structures*/
                jacobian = m_get(react->num_species_involved, react->num_species_involved);
                b = v_get(react->num_species_involved);
                x = v_get(react->num_species_involved);
                pivot = px_get(jacobian->m);
                states_cache = (double*) malloc(sizeof(double) * react->num_species_involved);
                params_cache = (double*) malloc(sizeof(double) * react->num_params_involved);
                states_cache_dx = (double*) malloc(sizeof(double) * react->num_species_involved);
                results_array = (double*) malloc(react->num_species_involved * sizeof(double));
                results_array_dx = (double*) malloc(react->num_species_involved * sizeof(double));
                for (i = start_idx; i <= stop_idx; i++) {
                    if (!react->subregion || react->subregion[i]) {
                        // TODO: this assumes grids are the same size/shape
                        //      add interpolation in case they aren't
                        for (j = 0; j < react->num_species_involved; j++) {
                            states_cache[j] = react->species_states[j][i];
                            states_cache_dx[j] = react->species_states[j][i];
                        }
                        MEM_ZERO(results_array, react->num_species_involved * sizeof(double));
                        for (k = 0; j < react->num_species_involved + react->num_params_involved;
                             k++, j++) {
                            params_cache[k] = react->species_states[j][i];
                        }
                        react->reaction(states_cache, params_cache, results_array, NULL);
 
                        for (j = 0; j < react->num_species_involved; j++) {
                            states_cache_dx[j] += dx;
                            MEM_ZERO(results_array_dx,
                                     react->num_species_involved * sizeof(double));
                            react->reaction(states_cache_dx, params_cache, results_array_dx, NULL);
                            v_set_val(b, j, dt * results_array[j]);
 
                            for (k = 0; k < react->num_species_involved; k++) {
                                pd = (results_array_dx[k] - results_array[k]) / dx;
                                m_set_val(jacobian, k, j, (j == k) - dt * pd);
                            }
                            states_cache_dx[j] -= dx;
                        }
                        // solve for x
                        if (react->num_species_involved == 1) {
                            react->species_states[0][i] += v_get_val(b, 0) /
                                                           m_get_val(jacobian, 0, 0);
                        } else {
                            // find entry in leftmost column with largest absolute value
                            // Pivot
                            for (j = 0; j < react->num_species_involved; j++) {
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    if (abs(m_get_val(jacobian, j, j)) <
                                        abs(m_get_val(jacobian, k, j))) {
                                        for (n = 0; n < react->num_species_involved; n++) {
                                            temp = m_get_val(jacobian, j, n);
                                            m_set_val(jacobian, j, n, m_get_val(jacobian, k, n));
                                            m_set_val(jacobian, k, n, temp);
                                        }
                                    }
                                }
                            }
 
                            for (j = 0; j < react->num_species_involved - 1; j++) {
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    ge_value = m_get_val(jacobian, k, j) /
                                               m_get_val(jacobian, j, j);
                                    for (n = 0; n < react->num_species_involved; n++) {
                                        val_to_set = m_get_val(jacobian, k, n) -
                                                     ge_value * m_get_val(jacobian, j, n);
                                        m_set_val(jacobian, k, n, val_to_set);
                                    }
                                    v_set_val(b, k, v_get_val(b, k) - ge_value * v_get_val(b, j));
                                }
                            }
 
                            for (j = react->num_species_involved - 1; j + 1 > 0; j--) {
                                v_set_val(x, j, v_get_val(b, j));
                                for (k = j + 1; k < react->num_species_involved; k++) {
                                    if (k != j) {
                                        v_set_val(x,
                                                  j,
                                                  v_get_val(x, j) -
                                                      m_get_val(jacobian, j, k) * v_get_val(x, k));
                                    }
                                }
                                v_set_val(x, j, v_get_val(x, j) / m_get_val(jacobian, j, j));
                            }
                            for (j = 0; j < react->num_species_involved; j++) {
                                // I think this should be something like
                                // react->species_states[j][mc3d_indices[i]] += v_get_val(x,j);
                                // Since the grid has a uniform discretization, the mc3d_indices
                                // should be the same length. So just need to access the correct
                                // mc3d_indices[i] maybe do two lines?: index =
                                // react->species_indices[j][i] react->species_states[j][index] +=
                                // v_get_val(x,j);
                                react->species_states[j][i] += v_get_val(x, j);
                            }
                        }
                    }
                }
                m_free(jacobian);
                v_free(b);
                v_free(x);
                px_free(pivot);
 
                SAFE_FREE(states_cache);
                SAFE_FREE(states_cache_dx);
                SAFE_FREE(params_cache);
                SAFE_FREE(results_array);
                SAFE_FREE(results_array_dx);
 
                if (stop)
                    return NULL;
            }
        }
    }
    return NULL;
}
 
/* run_threaded_reactions
 * Array ReactGridData tasks length NUM_THREADS and calls ecs_do_reactions and
 * executes each task with a separate thread
 */
static void run_threaded_reactions(ReactGridData* tasks) {
    int k;
    /* launch threads */
    for (k = 0; k < NUM_THREADS - 1; k++) {
        TaskQueue_add_task(AllTasks, &ecs_do_reactions, &tasks[k], NULL);
    }
    /* run one task in the main thread */
    ecs_do_reactions(&tasks[NUM_THREADS - 1]);
 
    /* wait for them to finish */
    TaskQueue_sync(AllTasks);
}
 
void _fadvance_fixed_step_3D(void) {
    Grid_node* grid;
    ECS_Grid_node* g;
    double dt = (*dt_ptr);
 
    /*Currents to broadcast via MPI*/
    /*TODO: Handle multiple grids with one pass*/
    /*Maybe TODO: Should check #currents << #voxels and not the other way round*/
    int id;
 
    if (threaded_reactions_tasks != NULL)
        run_threaded_reactions(threaded_reactions_tasks);
 
    for (id = 0, grid = Parallel_grids[0]; grid != NULL; grid = grid->next, id++) {
        MEM_ZERO(grid->states_cur, sizeof(double) * grid->size_x * grid->size_y * grid->size_z);
        g = dynamic_cast<ECS_Grid_node*>(grid);
        if (g)
            g->do_multicompartment_reactions(NULL);
        grid->do_grid_currents(grid->states_cur, dt, id);
        grid->apply_node_flux3D(dt, NULL);
 
        // grid->volume_setup();
        if (grid->hybrid) {
            grid->hybrid_connections();
        }
        grid->dg_adi();
    }
    /* transfer concentrations */
    scatter_concentrations();
}
 
extern "C" void scatter_concentrations(void) {
    /* transfer concentrations to classic NEURON */
    Grid_node* grid;
 
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid->scatter_grid_concentrations();
    }
}
 
/*****************************************************************************
 *
 * Begin variable step code
 *
 *****************************************************************************/
 
/* count the total number of state variables AND store their offset (passed in) in the cvode vector
 */
int ode_count(const int offset) {
    int count = 0;
    states_cvode_offset = offset;
    Grid_node* grid;
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        count += grid->size_x * grid->size_y * grid->size_z;
    }
    return count;
}
 
 
void ecs_atolscale(double* y) {
    Grid_node* grid;
    ssize_t i;
    int grid_size;
    y += states_cvode_offset;
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_size = grid->size_x * grid->size_y * grid->size_z;
        for (i = 0; i < grid_size; i++) {
            y[i] *= grid->atolscale;
        }
        y += grid_size;
    }
}
 
void _ecs_ode_reinit(double* y) {
    Grid_node* grid;
    ECS_Grid_node* g;
 
    ssize_t i;
    int grid_size;
    double* grid_states;
    y += states_cvode_offset;
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_states = grid->states;
        grid_size = grid->size_x * grid->size_y * grid->size_z;
 
        for (i = 0; i < grid_size; i++) {
            y[i] = grid_states[i];
        }
        y += grid_size;
        g = dynamic_cast<ECS_Grid_node*>(grid);
        if (g)
            g->initialize_multicompartment_reaction();
    }
}
 
 
void _rhs_variable_step_ecs(const double* states, double* ydot) {
    Grid_node* grid;
    ECS_Grid_node* g;
    ssize_t i;
    int grid_size;
    double dt = *dt_ptr;
    double* grid_states;
    double* const orig_ydot = ydot + states_cvode_offset;
    double const* const orig_states = states + states_cvode_offset;
    const unsigned char calculate_rhs = ydot == NULL ? 0 : 1;
    states = orig_states;
    ydot = orig_ydot;
    /* prepare for advance by syncing data with local copy */
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_states = grid->states;
        grid_size = grid->size_x * grid->size_y * grid->size_z;
 
        /* copy the passed in states to local memory (needed to make reaction rates correct) */
        for (i = 0; i < grid_size; i++) {
            grid_states[i] = states[i];
        }
        states += grid_size;
    }
    /* transfer concentrations to classic NEURON states */
    scatter_concentrations();
    if (!calculate_rhs) {
        return;
    }
 
    states = orig_states;
    ydot = orig_ydot;
    /* TODO: reactions contribute to adaptive step-size*/
    if (threaded_reactions_tasks != NULL) {
        run_threaded_reactions(threaded_reactions_tasks);
    }
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_states = grid->states;
        grid_size = grid->size_x * grid->size_y * grid->size_z;
        for (i = 0; i < grid_size; i++) {
            ydot[i] += (grid_states[i] - states[i]) / dt;
            grid_states[i] = states[i];
        }
        states += grid_size;
        ydot += grid_size;
    }
 
    ydot = orig_ydot;
    states = orig_states;
    /* process currents */
    for (i = 0, grid = Parallel_grids[0]; grid != NULL; grid = grid->next, i++) {
        g = dynamic_cast<ECS_Grid_node*>(grid);
        if (g)
            g->do_multicompartment_reactions(ydot);
        grid->do_grid_currents(ydot, 1.0, i);
 
        /*Add node fluxes to the result*/
        grid->apply_node_flux3D(1.0, ydot);
 
        ydot += grid_size;
    }
    ydot = orig_ydot;
 
    /* do the diffusion rates */
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_size = grid->size_x * grid->size_y * grid->size_z;
        grid->variable_step_diffusion(states, ydot);
 
        ydot += grid_size;
        states += grid_size;
    }
}
 
 
void _rhs_variable_step_helper(Grid_node* g, double const* const states, double* ydot) {
    int num_states_x = g->size_x, num_states_y = g->size_y, num_states_z = g->size_z;
    double dx = g->dx, dy = g->dy, dz = g->dz;
    int i, j, k, stop_i, stop_j, stop_k;
    double dc_x = g->dc_x, dc_y = g->dc_y, dc_z = g->dc_z;
 
    double rate_x = dc_x / (dx * dx);
    double rate_y = dc_y / (dy * dy);
    double rate_z = dc_z / (dz * dz);
 
    /*indices*/
    int index, prev_i, prev_j, prev_k, next_i, next_j, next_k;
    double div_x, div_y, div_z;
 
    /* Euler advance x, y, z (all points) */
    stop_i = num_states_x - 1;
    stop_j = num_states_y - 1;
    stop_k = num_states_z - 1;
    if (g->bc->type == NEUMANN) {
        for (i = 0,
            index = 0,
            prev_i = num_states_y * num_states_z,
            next_i = num_states_y * num_states_z;
             i < num_states_x;
             i++) {
            /*Zero flux boundary conditions*/
            div_x = (i == 0 || i == stop_i) ? 2. : 1.;
 
            for (j = 0, prev_j = index + num_states_z, next_j = index + num_states_z;
                 j < num_states_y;
                 j++) {
                div_y = (j == 0 || j == stop_j) ? 2. : 1.;
 
                for (k = 0, prev_k = index + 1, next_k = index + 1; k < num_states_z;
                     k++, index++, prev_i++, next_i++, prev_j++, next_j++) {
                    div_z = (k == 0 || k == stop_k) ? 2. : 1.;
 
                    if (stop_i > 0) {
                        ydot[index] += rate_x *
                                       (states[prev_i] - 2.0 * states[index] + states[next_i]) /
                                       div_x;
                    }
                    if (stop_j > 0) {
                        ydot[index] += rate_y *
                                       (states[prev_j] - 2.0 * states[index] + states[next_j]) /
                                       div_y;
                    }
                    if (stop_k > 0) {
                        ydot[index] += rate_z *
                                       (states[prev_k] - 2.0 * states[index] + states[next_k]) /
                                       div_z;
                    }
                    next_k = (k == stop_k - 1) ? index : index + 2;
                    prev_k = index;
                }
                prev_j = index - num_states_z;
                next_j = (j == stop_j - 1) ? prev_j : index + num_states_z;
            }
            prev_i = index - num_states_y * num_states_z;
            next_i = (i == stop_i - 1) ? prev_i : index + num_states_y * num_states_z;
        }
    } else {
        for (i = 0, index = 0, prev_i = 0, next_i = num_states_y * num_states_z; i < num_states_x;
             i++) {
            for (j = 0, prev_j = index - num_states_z, next_j = index + num_states_z;
                 j < num_states_y;
                 j++) {
                for (k = 0, prev_k = index - 1, next_k = index + 1; k < num_states_z;
                     k++, index++, prev_i++, next_i++, prev_j++, next_j++, next_k++, prev_k++) {
                    if (i == 0 || i == stop_i || j == 0 || j == stop_j || k == 0 || k == stop_k) {
                        // set to zero to prevent currents altering concentrations at the boundary
                        ydot[index] = 0;
                    } else {
                        ydot[index] += rate_x *
                                       (states[prev_i] - 2.0 * states[index] + states[next_i]);
 
                        ydot[index] += rate_y *
                                       (states[prev_j] - 2.0 * states[index] + states[next_j]);
 
                        ydot[index] += rate_z *
                                       (states[prev_k] - 2.0 * states[index] + states[next_k]);
                    }
                }
                prev_j = index - num_states_z;
                next_j = index + num_states_z;
            }
            prev_i = index - num_states_y * num_states_z;
            next_i = index + num_states_y * num_states_z;
        }
    }
    /*
            for (i = 1, index = (num_states_y+1)*num_states_z + 1,
                 prev_i = num_states_z + 1, next_i = (2*num_states_y+1)*num_states_z + 1;
                 i < stop_i; i++) {
                for (j = 1, prev_j = index - num_states_z, next_j = index + num_states_z; j <
       stop_j; j++) {
 
                    for(k = 1, prev_k = index - 1, next_k = index + 1; k < stop_k;
                        k++, index++, prev_i++, next_i++, prev_j++, next_j++, next_k++) {
                        if(index != i*num_states_y*num_states_z + j*num_states_z + k)
                        {
                            fprintf(stderr,"%i \t %i %i %i\n", index,i,j,k);
                        }
                        if(prev_i != (i-1)*num_states_y*num_states_z + j*num_states_z + k)
                        {
                            fprintf(stderr,"prev_i %i %i \t %i %i %i\n",
       prev_i,(i-1)*num_states_y*num_states_z + j*num_states_z + k,i,j,k);
                        }
                        if(next_i != (i+1)*num_states_y*num_states_z + j*num_states_z + k)
                        {
                            fprintf(stderr,"next_i %i %i \t %i %i %i\n",
       next_i,(i+1)*num_states_y*num_states_z + j*num_states_z + k,i,j,k);
                        }
                        if(prev_j != i*num_states_y*num_states_z + (j-1)*num_states_z + k)
                        {
                            fprintf(stderr,"prev_j %i %i \t %i %i %i\n",
       prev_j,i*num_states_y*num_states_z + (j-1)*num_states_z + k,i,j,k);
                        }
                        if(next_j != i*num_states_y*num_states_z + (j+1)*num_states_z + k)
                        {
                            fprintf(stderr,"next_j %i %i \t %i %i %i\n",
       next_j,i*num_states_y*num_states_z + (j+1)*num_states_z + k,i,j,k);
                        }
                        if(prev_k != i*num_states_y*num_states_z + j*num_states_z + k-1)
                        {
                            fprintf(stderr,"prev_k %i %i \t %i %i %i\n",
       prev_k,i*num_states_y*num_states_z + j*num_states_z + k-1,i,j,k);
                        }
                        if(next_k != i*num_states_y*num_states_z + j*num_states_z + k+1)
                        {
                            fprintf(stderr,"next_k %i %i \t %i %i %i\n",
       next_k,i*num_states_y*num_states_z + j*num_states_z + k+1,i,j,k);
                        }
 
                        ydot[index] += rate_x * (states[prev_i] -
                            2.0 * states[index] + states[next_i]);
 
                        ydot[index] += rate_y * (states[prev_j] -
                            2.0 * states[index] + states[next_j]);
 
                        ydot[index] += rate_z * (states[prev_k] -
                            2.0 * states[index] + states[next_k]);
 
                    }
                    index += 2; //skip z boundaries
                    //prev_j = index - num_states_z;
                    //next_j = index + num_states_z;
                }
                index += 2*num_states_z; //skip y boundaries
                //prev_i = index - num_states_y*num_states_z;
                //next_i = index + num_states_y*num_states_z;
            }
        }
        */
}
 
// p1 = b  p2 = states
void ics_ode_solve(double dt, double* RHS, const double* states) {
    Grid_node* grid;
    ssize_t i;
    int grid_size = 0;
    double* grid_states;
    double const* const orig_states = states + states_cvode_offset;
    const unsigned char calculate_rhs = RHS == NULL ? 0 : 1;
    double* const orig_RHS = RHS + states_cvode_offset;
    states = orig_states;
    RHS = orig_RHS;
 
    /* prepare for advance by syncing data with local copy */
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid_states = grid->states;
        grid_size = grid->size_x * grid->size_y * grid->size_z;
 
        /* copy the passed in states to local memory (needed to make reaction rates correct) */
        for (i = 0; i < grid_size; i++) {
            grid_states[i] = states[i];
        }
        states += grid_size;
    }
    /* transfer concentrations to classic NEURON states */
    scatter_concentrations();
    if (!calculate_rhs) {
        return;
    }
 
    states = orig_states;
    RHS = orig_RHS;
 
    /* TODO: reactions contribute to adaptive step-size*/
    if (threaded_reactions_tasks != NULL) {
        run_threaded_reactions(threaded_reactions_tasks);
    }
    /* do the diffusion rates */
    for (grid = Parallel_grids[0]; grid != NULL; grid = grid->next) {
        grid->variable_step_ode_solve(RHS, dt);
        RHS += grid_size;
        states += grid_size;
    }
}
/*****************************************************************************
 *
 * End variable step code
 *
 *****************************************************************************/
 
 
/* solve_dd_clhs_tridiag uses Thomas Algorithm to solves a
 * diagonally dominant tridiagonal matrix (Ax=b), where the
 * triple (lower, diagonal and upper) are repeated to form A.
 * N        -   length of the matrix
 * l_diag   -   lower diagonal
 * diag     -   diagonal
 * u_diag   -   upper diagonal
 * b        -   pointer to the RHS  (length N)
 * lbc_diag -   the lower boundary diagonal (the first diagonal element)
 * lbc_u_diag - the lower boundary upper diagonal (the first upper
 *              diagonal element)
 * ubc_l_diag - the upper boundary lower diagonal (the last lower
 *              diagonal element)
 * ubc_diag -   the upper boundary diagonal (the last diagonal element)
 * The solution (x) is stored in b.
 * c          - scratchpad array, N - 1 doubles long
 */
static int solve_dd_clhs_tridiag(const int N,
                                 const double l_diag,
                                 const double diag,
                                 const double u_diag,
                                 const double lbc_diag,
                                 const double lbc_u_diag,
                                 const double ubc_l_diag,
                                 const double ubc_diag,
                                 double* const b,
                                 double* const c) {
    int i;
 
 
    c[0] = lbc_u_diag / lbc_diag;
    b[0] = b[0] / lbc_diag;
 
    for (i = 1; i < N - 1; i++) {
        c[i] = u_diag / (diag - l_diag * c[i - 1]);
        b[i] = (b[i] - l_diag * b[i - 1]) / (diag - l_diag * c[i - 1]);
    }
    b[N - 1] = (b[N - 1] - ubc_l_diag * b[N - 2]) / (ubc_diag - ubc_l_diag * c[N - 2]);
 
    /*back substitution*/
    for (i = N - 2; i >= 0; i--) {
        b[i] = b[i] - c[i] * b[i + 1];
    }
    return 0;
}
 
/*
static int solve_dd_clhs_tridiag_rev(const int N, const double l_diag, const double diag,
    const double u_diag, const double lbc_diag, const double lbc_u_diag,
    const double ubc_l_diag, const double ubc_diag, double* const d, double* const a)
{
    int i;
 
 
    a[N-2] = ubc_l_diag/ubc_diag;
    d[N-1] = d[N-1]/ubc_diag;
 
    for(i=N-2; i>0; i--)
    {
        a[i-1] = l_diag/(diag - u_diag*a[i]);
        d[i] = (d[i] - u_diag*d[i+1])/(diag - u_diag*a[i]);
    }
    d[0] = (d[0] - lbc_u_diag*d[1])/(lbc_diag - lbc_u_diag*a[0]);
 
    for(i=1; i<N; i++)
    {
        d[i] = d[i] - a[i-1]*d[i-1];
 
    }
    return 0;
}
unsigned char callcount = TRUE;
static int solve_dd_clhs_tridiag(const int N, const double l_diag, const double diag,
    const double u_diag, const double lbc_diag, const double lbc_u_diag,
    const double ubc_l_diag, const double ubc_diag, double* const b, double* const c)
{
    if(callcount)
    {
        callcount = FALSE;
        return solve_dd_clhs_tridiag_for(N, l_diag, diag, u_diag, lbc_diag, lbc_u_diag, ubc_l_diag,
ubc_diag, b, c);
    }
    else
    {
        callcount = TRUE;
        return solve_dd_clhs_tridiag_rev(N, l_diag, diag, u_diag, lbc_diag, lbc_u_diag, ubc_l_diag,
ubc_diag, b, c);
    }
}
*/
 
 
/* dg_adi_x performs the first of 3 steps in DG-ADI
 * g    -   the parameters and state of the grid
 * dt   -   the time step
 * y    -   the index for the y plane
 * z    -   the index for the z plane
 * state    -   where the output of this step is stored
 * scratch  - scratchpad array of doubles, length g->size_x - 1
 */
static void ecs_dg_adi_x(ECS_Grid_node* g,
                         const double dt,
                         const int y,
                         const int z,
                         double const* const state,
                         double* const RHS,
                         double* const scratch) {
    int yp, ym, zp, zm;
    int x;
    double r = g->dc_x * dt / SQ(g->dx);
    double div_y, div_z;
    /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
    if (g->bc->type == DIRICHLET &&
        (y == 0 || z == 0 || y == g->size_y - 1 || z == g->size_z - 1)) {
        for (x = 0; x < g->size_x; x++)
            RHS[x] = g->bc->value;
        return;
    }
 
    if (g->size_y > 1) {
        yp = (y == g->size_y - 1) ? y - 1 : y + 1;
        ym = (y == 0) ? y + 1 : y - 1;
        div_y = (y == 0 || y == g->size_y - 1) ? 2. : 1.;
    } else {
        yp = 0;
        ym = 0;
        div_y = 1;
    }
    if (g->size_z > 1) {
        zp = (z == g->size_z - 1) ? z - 1 : z + 1;
        zm = (z == 0) ? z + 1 : z - 1;
        div_z = (z == 0 || z == g->size_z - 1) ? 2. : 1.;
    } else {
        zp = 0;
        zm = 0;
        div_z = 1;
    }
 
    if (g->bc->type == NEUMANN) {
        /*zero flux boundary condition*/
        RHS[0] = state[IDX(0, y, z)] + g->states_cur[IDX(0, y, z)] +
                 dt *
                     ((g->dc_y / SQ(g->dy)) *
                          (state[IDX(0, yp, z)] - 2. * state[IDX(0, y, z)] + state[IDX(0, ym, z)]) /
                          div_y +
                      (g->dc_z / SQ(g->dz)) *
                          (state[IDX(0, y, zp)] - 2. * state[IDX(0, y, z)] + state[IDX(0, y, zm)]) /
                          div_z);
        if (g->size_x > 1) {
            RHS[0] += dt * (g->dc_x / SQ(g->dx)) * (state[IDX(1, y, z)] - state[IDX(0, y, z)]);
            x = g->size_x - 1;
            RHS[x] =
                state[IDX(x, y, z)] + g->states_cur[IDX(x, y, z)] +
                dt * ((g->dc_x / SQ(g->dx)) * (state[IDX(x - 1, y, z)] - state[IDX(x, y, z)]) +
                      (g->dc_y / SQ(g->dy)) *
                          (state[IDX(x, yp, z)] - 2. * state[IDX(x, y, z)] + state[IDX(x, ym, z)]) /
                          div_y +
                      (g->dc_z / SQ(g->dz)) *
                          (state[IDX(x, y, zp)] - 2. * state[IDX(x, y, z)] + state[IDX(x, y, zm)]) /
                          div_z);
        }
    } else {
        RHS[0] = g->bc->value;
        RHS[g->size_x - 1] = g->bc->value;
    }
    for (x = 1; x < g->size_x - 1; x++) {
#ifndef __PGI
        __builtin_prefetch(&(state[IDX(x + PREFETCH, y, z)]), 0, 1);
        __builtin_prefetch(&(state[IDX(x + PREFETCH, yp, z)]), 0, 0);
        __builtin_prefetch(&(state[IDX(x + PREFETCH, ym, z)]), 0, 0);
#endif
        RHS[x] = state[IDX(x, y, z)] +
                 dt *
                     ((g->dc_x / SQ(g->dx)) *
                          (state[IDX(x + 1, y, z)] - 2. * state[IDX(x, y, z)] +
                           state[IDX(x - 1, y, z)]) /
                          2. +
                      (g->dc_y / SQ(g->dy)) *
                          (state[IDX(x, yp, z)] - 2. * state[IDX(x, y, z)] + state[IDX(x, ym, z)]) /
                          div_y +
                      (g->dc_z / SQ(g->dz)) *
                          (state[IDX(x, y, zp)] - 2. * state[IDX(x, y, z)] + state[IDX(x, y, zm)]) /
                          div_z) +
                 g->states_cur[IDX(x, y, z)];
    }
    if (g->size_x > 1) {
        if (g->bc->type == NEUMANN)
            solve_dd_clhs_tridiag(g->size_x,
                                  -r / 2.0,
                                  1.0 + r,
                                  -r / 2.0,
                                  1.0 + r / 2.0,
                                  -r / 2.0,
                                  -r / 2.0,
                                  1.0 + r / 2.0,
                                  RHS,
                                  scratch);
        else
            solve_dd_clhs_tridiag(
                g->size_x, -r / 2.0, 1.0 + r, -r / 2.0, 1.0, 0, 0, 1.0, RHS, scratch);
    }
}
 
 
/* dg_adi_y performs the second of 3 steps in DG-ADI
 * g    -   the parameters and state of the grid
 * dt   -   the time step
 * x    -   the index for the x plane
 * z    -   the index for the z plane
 * state    -   the values from the first step, which are
 *              overwritten by the output of this step
 * scratch  -   scratchpad array of doubles, length g->size_y - 1
 */
static void ecs_dg_adi_y(ECS_Grid_node* g,
                         double const dt,
                         int const x,
                         int const z,
                         double const* const state,
                         double* const RHS,
                         double* const scratch) {
    int y;
    double r = (g->dc_y * dt / SQ(g->dy));
    /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
    if (g->bc->type == DIRICHLET &&
        (x == 0 || z == 0 || x == g->size_x - 1 || z == g->size_z - 1)) {
        for (y = 0; y < g->size_y; y++)
            RHS[y] = g->bc->value;
        return;
    }
    if (g->size_y == 1) {
        if (g->bc->type == NEUMANN)
            RHS[0] = state[x + z * g->size_x];
        else
            RHS[0] = g->bc->value;
        return;
    }
    if (g->bc->type == NEUMANN) {
        /*zero flux boundary condition*/
        RHS[0] = state[x + z * g->size_x] -
                 (g->dc_y * dt / SQ(g->dy)) *
                     (g->states[IDX(x, 1, z)] - 2.0 * g->states[IDX(x, 0, z)] +
                      g->states[IDX(x, 1, z)]) /
                     4.0;
        y = g->size_y - 1;
        RHS[y] = state[x + (z + y * g->size_z) * g->size_x] -
                 (g->dc_y * dt / SQ(g->dy)) *
                     (g->states[IDX(x, y - 1, z)] - 2. * g->states[IDX(x, y, z)] +
                      g->states[IDX(x, y - 1, z)]) /
                     4.0;
    } else {
        RHS[0] = g->bc->value;
        RHS[g->size_y - 1] = g->bc->value;
    }
    for (y = 1; y < g->size_y - 1; y++) {
#ifndef __PGI
        __builtin_prefetch(&state[x + (z + (y + PREFETCH) * g->size_z) * g->size_x], 0, 0);
        __builtin_prefetch(&(g->states[IDX(x, y + PREFETCH, z)]), 0, 1);
#endif
        RHS[y] = state[x + (z + y * g->size_z) * g->size_x] -
                 (g->dc_y * dt / SQ(g->dy)) *
                     (g->states[IDX(x, y + 1, z)] - 2. * g->states[IDX(x, y, z)] +
                      g->states[IDX(x, y - 1, z)]) /
                     2.0;
    }
    if (g->bc->type == NEUMANN)
        solve_dd_clhs_tridiag(g->size_y,
                              -r / 2.0,
                              1.0 + r,
                              -r / 2.0,
                              1.0 + r / 2.0,
                              -r / 2.0,
                              -r / 2.0,
                              1.0 + r / 2.0,
                              RHS,
                              scratch);
    else
        solve_dd_clhs_tridiag(g->size_y, -r / 2., 1. + r, -r / 2., 1.0, 0, 0, 1.0, RHS, scratch);
}
 
 
/* dg_adi_z performs the final step in DG-ADI
 * g    -   the parameters and state of the grid
 * dt   -   the time step
 * x    -   the index for the x plane
 * y    -   the index for the y plane
 * state    -   the values from the second step, which are
 *              overwritten by the output of this step
 * scratch  -   scratchpad array of doubles, length g->size_z - 1
 */
static void ecs_dg_adi_z(ECS_Grid_node* g,
                         double const dt,
                         int const x,
                         int const y,
                         double const* const state,
                         double* const RHS,
                         double* const scratch) {
    int z;
    double r = g->dc_z * dt / SQ(g->dz);
    /*TODO: Get rid of this by not calling dg_adi when on the boundary for DIRICHLET conditions*/
    if (g->bc->type == DIRICHLET &&
        (x == 0 || y == 0 || x == g->size_x - 1 || y == g->size_y - 1)) {
        for (z = 0; z < g->size_z; z++)
            RHS[z] = g->bc->value;
        return;
    }
 
    if (g->size_z == 1) {
        if (g->bc->type == NEUMANN)
            RHS[0] = state[y + g->size_y * x];
        else
            RHS[0] = g->bc->value;
        return;
    }
 
    if (g->bc->type == NEUMANN) {
        /*zero flux boundary condition*/
        RHS[0] = state[y + g->size_y * (x * g->size_z)] -
                 (g->dc_z * dt / SQ(g->dz)) *
                     (g->states[IDX(x, y, 1)] - 2.0 * g->states[IDX(x, y, 0)] +
                      g->states[IDX(x, y, 1)]) /
                     4.0;
        z = g->size_z - 1;
        RHS[z] = state[y + g->size_y * (x * g->size_z + z)] -
                 (g->dc_z * dt / SQ(g->dz)) *
                     (g->states[IDX(x, y, z - 1)] - 2.0 * g->states[IDX(x, y, z)] +
                      g->states[IDX(x, y, z - 1)]) /
                     4.0;
 
    } else {
        RHS[0] = g->bc->value;
        RHS[g->size_z - 1] = g->bc->value;
    }
    for (z = 1; z < g->size_z - 1; z++) {
        RHS[z] = state[y + g->size_y * (x * g->size_z + z)] -
                 (g->dc_z * dt / SQ(g->dz)) *
                     (g->states[IDX(x, y, z + 1)] - 2. * g->states[IDX(x, y, z)] +
                      g->states[IDX(x, y, z - 1)]) /
                     2.;
    }
 
    if (g->bc->type == NEUMANN)
        solve_dd_clhs_tridiag(g->size_z,
                              -r / 2.0,
                              1.0 + r,
                              -r / 2.0,
                              1.0 + r / 2.0,
                              -r / 2.0,
                              -r / 2.0,
                              1.0 + r / 2.0,
                              RHS,
                              scratch);
    else
        solve_dd_clhs_tridiag(g->size_z, -r / 2., 1. + r, -r / 2., 1.0, 0, 0, 1.0, RHS, scratch);
}
 
static void* ecs_do_dg_adi(void* dataptr) {
    ECSAdiGridData* data = (ECSAdiGridData*) dataptr;
    int start = data->start;
    int stop = data->stop;
    int i, j, k;
    ECSAdiDirection* ecs_adi_dir = data->ecs_adi_dir;
    double dt = *dt_ptr;
    int sizej = data->sizej;
    ECS_Grid_node* g = data->g;
    double* state_in = ecs_adi_dir->states_in;
    double* state_out = ecs_adi_dir->states_out;
    int offset = ecs_adi_dir->line_size;
    double* scratchpad = data->scratchpad;
    void (*ecs_dg_adi_dir)(
        ECS_Grid_node*, double, int, int, double const* const, double* const, double* const) =
        ecs_adi_dir->ecs_dg_adi_dir;
    for (k = start; k < stop; k++) {
        i = k / sizej;
        j = k % sizej;
        ecs_dg_adi_dir(g, dt, i, j, state_in, &state_out[k * offset], scratchpad);
    }
 
    return NULL;
}
 
void ecs_run_threaded_dg_adi(const int i,
                             const int j,
                             ECS_Grid_node* g,
                             ECSAdiDirection* ecs_adi_dir,
                             const int n) {
    int k;
    /* when doing any given direction, the number of tasks is the product of the other two, so
     * multiply everything then divide out the current direction */
    const int tasks_per_thread = (g->size_x * g->size_y * g->size_z / n) / NUM_THREADS;
    const int extra = (g->size_x * g->size_y * g->size_z / n) % NUM_THREADS;
 
    g->ecs_tasks[0].start = 0;
    g->ecs_tasks[0].stop = tasks_per_thread + (extra > 0);
    g->ecs_tasks[0].sizej = j;
    g->ecs_tasks[0].ecs_adi_dir = ecs_adi_dir;
    for (k = 1; k < NUM_THREADS; k++) {
        g->ecs_tasks[k].start = g->ecs_tasks[k - 1].stop;
        g->ecs_tasks[k].stop = g->ecs_tasks[k].start + tasks_per_thread + (extra > k);
        g->ecs_tasks[k].sizej = j;
        g->ecs_tasks[k].ecs_adi_dir = ecs_adi_dir;
    }
    g->ecs_tasks[NUM_THREADS - 1].stop = i * j;
    /* launch threads */
    for (k = 0; k < NUM_THREADS - 1; k++) {
        TaskQueue_add_task(AllTasks, &ecs_do_dg_adi, &(g->ecs_tasks[k]), NULL);
    }
    /* run one task in the main thread */
    ecs_do_dg_adi(&(g->ecs_tasks[NUM_THREADS - 1]));
    /* wait for them to finish */
    TaskQueue_sync(AllTasks);
}
 
void ecs_set_adi_homogeneous(ECS_Grid_node* g) {
    g->ecs_adi_dir_x->ecs_dg_adi_dir = ecs_dg_adi_x;
    g->ecs_adi_dir_y->ecs_dg_adi_dir = ecs_dg_adi_y;
    g->ecs_adi_dir_z->ecs_dg_adi_dir = ecs_dg_adi_z;
}