From 7a116d55120e2bf5178259b24b8a32a45891fc2b Mon Sep 17 00:00:00 2001 From: adam newton Date: Wed, 6 Dec 2023 17:48:23 -0500 Subject: [PATCH 1/2] Moved the memory allocation for rxd reactions to `register_rate`. --- src/nrnpython/rxd.cpp | 336 +++++++++++++++++++----------------------- src/nrnpython/rxd.h | 14 ++ 2 files changed, 163 insertions(+), 187 deletions(-) diff --git a/src/nrnpython/rxd.cpp b/src/nrnpython/rxd.cpp index 5c1117c4a3..d343ec0d2b 100644 --- a/src/nrnpython/rxd.cpp +++ b/src/nrnpython/rxd.cpp @@ -949,6 +949,45 @@ extern "C" void register_rate(int nspecies, } else { react->ecs_state = NULL; } + if (react->num_mult > 0) + react->mc_mult = (double*) malloc(react->num_mult * sizeof(double)); + + if (react->num_ecs_species > 0) { + react->ecs_states_for_reaction = (double*) malloc(react->num_ecs_species * sizeof(double)); + react->ecs_states_for_reaction_dx = (double*) malloc(react->num_ecs_species * + sizeof(double)); + react->ecs_result = (double*) malloc(react->num_ecs_species * sizeof(double)); + react->ecs_result_dx = (double*) malloc(react->num_ecs_species * sizeof(double)); + react->ecsindex = (int*) malloc(react->num_ecs_species * sizeof(int)); + } + if (react->num_ecs_params > 0) + react->ecs_params_for_reaction = (double*) calloc(react->num_ecs_params, sizeof(double)); + + if (_membrane_flux) { + react->flux = (double**) malloc(react->icsN * sizeof(double*)); + for (i = 0; i < react->icsN; i++) + react->flux[i] = (double*) malloc(react->num_regions * sizeof(double)); + } + react->states_for_reaction = (double**) malloc(react->num_species * sizeof(double*)); + react->states_for_reaction_dx = (double**) malloc(react->num_species * sizeof(double*)); + + react->result_array = (double**) malloc(react->num_species * sizeof(double*)); + react->result_array_dx = (double**) malloc(react->num_species * sizeof(double*)); + + for (i = 0; i < react->num_species; i++) { + react->states_for_reaction[i] = (double*) malloc(react->num_regions * sizeof(double)); + react->states_for_reaction_dx[i] = (double*) malloc(react->num_regions * sizeof(double)); + + react->result_array[i] = (double*) malloc(react->num_regions * sizeof(double)); + react->result_array_dx[i] = (double*) malloc(react->num_regions * sizeof(double)); + } + + react->params_for_reaction = (double**) malloc(react->num_params * sizeof(double*)); + for (i = 0; i < react->num_params; i++) + react->params_for_reaction[i] = (double*) malloc(react->num_regions * sizeof(double)); + if (react->num_mult > 0) + react->mc_mult = (double*) malloc(react->num_mult * sizeof(double)); + if (_reactions == NULL) { _reactions = react; react->next = NULL; @@ -990,6 +1029,34 @@ extern "C" void clear_rates() { free(react->state_idx); SAFE_FREE(react->ecs_state); prev = react; + + if (react->num_mult > 0) + free(react->mc_mult); + if (_membrane_flux) { + for (i = 0; i < react->icsN; i++) + free(react->flux[i]); + free(react->flux); + } + if (react->num_ecs_species > 0) { + free(react->ecs_states_for_reaction); + free(react->ecs_states_for_reaction_dx); + free(react->ecs_result); + free(react->ecs_result_dx); + free(react->ecsindex); + } + for (i = 0; i < react->num_species; i++) { + free(react->states_for_reaction[i]); + free(react->result_array[i]); + } + free(react->states_for_reaction); + free(react->result_array); + for (i = 0; i < react->num_params; i++) { + free(react->params_for_reaction[i]); + } + free(react->params_for_reaction); + if (react->num_ecs_params > 0) { + free(react->ecs_params_for_reaction); + } react = react->next; SAFE_FREE(prev); } @@ -1357,95 +1424,67 @@ void _rhs_variable_step(const double* p1, double* p2) { void get_reaction_rates(ICSReactions* react, double* states, double* rates, double* ydot) { int segment; int i, j, k, idx; - double** states_for_reaction = (double**) malloc(react->num_species * sizeof(double*)); - double** params_for_reaction = (double**) malloc(react->num_params * sizeof(double*)); - double** result_array = (double**) malloc(react->num_species * sizeof(double*)); - double* mc_mult = NULL; - double** flux = NULL; - if (react->num_mult > 0) - mc_mult = (double*) malloc(react->num_mult * sizeof(double)); - - double* ecs_states_for_reaction = NULL; - double* ecs_params_for_reaction = NULL; - double* ecs_result = NULL; - int* ecsindex = NULL; double v = 0; - if (react->num_ecs_species > 0) { - ecs_states_for_reaction = (double*) malloc(react->num_ecs_species * sizeof(double)); - ecs_result = (double*) malloc(react->num_ecs_species * sizeof(double)); - } - if (react->num_ecs_params > 0) - ecs_params_for_reaction = (double*) calloc(react->num_ecs_params, sizeof(double)); + for (i = 0; i < react->num_ecs_species; i++) + react->ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]]; if (_membrane_flux) { - flux = (double**) malloc(react->icsN * sizeof(double*)); for (i = 0; i < react->icsN; i++) - flux[i] = (double*) calloc(react->num_regions, sizeof(double)); + bzero(react->flux[i], react->num_regions * sizeof(double)); } - - for (i = 0; i < react->num_species; i++) { - states_for_reaction[i] = (double*) calloc(react->num_regions, sizeof(double)); - result_array[i] = (double*) malloc(react->num_regions * sizeof(double)); - } - for (i = 0; i < react->num_params; i++) - params_for_reaction[i] = (double*) calloc(react->num_regions, sizeof(double)); - ecsindex = (int*) malloc(react->num_ecs_species * sizeof(int)); - for (i = 0; i < react->num_ecs_species; i++) - ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]]; - for (segment = 0; segment < react->num_segments; segment++) { for (i = 0; i < react->num_species; i++) { for (j = 0; j < react->num_regions; j++) { if (react->state_idx[segment][i][j] != SPECIES_ABSENT) { - states_for_reaction[i][j] = states[react->state_idx[segment][i][j]]; + react->states_for_reaction[i][j] = states[react->state_idx[segment][i][j]]; } else { - states_for_reaction[i][j] = NAN; + react->states_for_reaction[i][j] = NAN; } } - memset(result_array[i], 0, react->num_regions * sizeof(double)); + memset(react->result_array[i], 0, react->num_regions * sizeof(double)); } for (k = 0; i < react->num_species + react->num_params; i++, k++) { for (j = 0; j < react->num_regions; j++) { if (react->state_idx[segment][i][j] != SPECIES_ABSENT) { - params_for_reaction[k][j] = states[react->state_idx[segment][i][j]]; + react->params_for_reaction[k][j] = states[react->state_idx[segment][i][j]]; } else { - params_for_reaction[k][j] = NAN; + react->params_for_reaction[k][j] = NAN; } } } for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) { - ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]); + react->ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]); } else { - ecs_states_for_reaction[i] = NAN; + react->ecs_states_for_reaction[i] = NAN; } } for (k = 0; i < react->num_ecs_species + react->num_ecs_params; i++, k++) { if (react->ecs_state[segment][i] != NULL) { - ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]); + react->ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]); } else { - ecs_params_for_reaction[k] = NAN; + react->ecs_params_for_reaction[k] = NAN; } } - memset(ecs_result, 0, react->num_ecs_species * sizeof(double)); + memset(react->ecs_result, 0, react->num_ecs_species * sizeof(double)); for (i = 0; i < react->num_mult; i++) { - mc_mult[i] = react->mc_multiplier[i][segment]; + react->mc_mult[i] = react->mc_multiplier[i][segment]; } if (react->vptrs != NULL) { v = *(react->vptrs[segment]); } - react->reaction(states_for_reaction, - params_for_reaction, - result_array, - mc_mult, - ecs_states_for_reaction, - ecs_params_for_reaction, - ecs_result, - flux, + react->reaction(react->states_for_reaction, + react->params_for_reaction, + react->result_array, + react->mc_mult, + react->ecs_states_for_reaction, + react->ecs_params_for_reaction, + react->ecs_result, + react->flux, v); for (i = 0; i < react->num_species; i++) { @@ -1454,10 +1493,10 @@ void get_reaction_rates(ICSReactions* react, double* states, double* rates, doub if (idx != SPECIES_ABSENT) { if (_membrane_flux && _membrane_lookup[idx] != SPECIES_ABSENT) { _rxd_induced_currents[_membrane_lookup[idx]] -= - _rxd_flux_scale[_membrane_lookup[idx]] * flux[i][j]; + _rxd_flux_scale[_membrane_lookup[idx]] * react->flux[i][j]; } if (rates != NULL) { - rates[idx] += result_array[i][j]; + rates[idx] += react->result_array[i][j]; } } } @@ -1465,36 +1504,12 @@ void get_reaction_rates(ICSReactions* react, double* states, double* rates, doub if (ydot != NULL) { for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) - react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = ecs_result[i]; + react->ecs_grid[i]->all_reaction_states[react->ecsindex[i]++] = + react->ecs_result[i]; // ydot[react->ecs_index[segment][i]] += ecs_result[i]; } } } - /* free allocated memory */ - if (react->num_mult > 0) - free(mc_mult); - if (_membrane_flux) { - for (i = 0; i < react->icsN; i++) - free(flux[i]); - free(flux); - } - if (react->num_ecs_species > 0) { - free(ecs_states_for_reaction); - free(ecs_result); - } - for (i = 0; i < react->num_species; i++) { - free(states_for_reaction[i]); - free(result_array[i]); - } - free(states_for_reaction); - free(result_array); - for (i = 0; i < react->num_params; i++) { - free(params_for_reaction[i]); - } - free(params_for_reaction); - if (react->num_ecs_params > 0) { - free(ecs_params_for_reaction); - } } void solve_reaction(ICSReactions* react, @@ -1511,46 +1526,12 @@ void solve_reaction(ICSReactions* react, auto jacobian = std::make_unique(N, N); auto b = std::make_unique(N); auto x = std::make_unique(N); - - double** states_for_reaction = (double**) malloc(react->num_species * sizeof(double*)); - double** states_for_reaction_dx = (double**) malloc(react->num_species * sizeof(double*)); - double** params_for_reaction = (double**) malloc(react->num_params * sizeof(double*)); - double** result_array = (double**) malloc(react->num_species * sizeof(double*)); - double** result_array_dx = (double**) malloc(react->num_species * sizeof(double*)); - double* mc_mult = NULL; - if (react->num_mult > 0) - mc_mult = (double*) malloc(react->num_mult * sizeof(double)); - - double* ecs_states_for_reaction = NULL; - double* ecs_states_for_reaction_dx = NULL; - double* ecs_params_for_reaction = NULL; - double* ecs_result = NULL; - double* ecs_result_dx = NULL; double v = 0; - int* ecsindex = NULL; if (react->num_ecs_species > 0) { - ecsindex = (int*) malloc(react->num_ecs_species * sizeof(int)); for (i = 0; i < react->num_ecs_species; i++) - ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]]; - ecs_states_for_reaction = (double*) malloc(react->num_ecs_species * sizeof(double)); - ecs_states_for_reaction_dx = (double*) malloc(react->num_ecs_species * sizeof(double)); - ecs_result = (double*) malloc(react->num_ecs_species * sizeof(double)); - ecs_result_dx = (double*) malloc(react->num_ecs_species * sizeof(double)); - } - - if (react->num_ecs_params > 0) - ecs_params_for_reaction = (double*) malloc(react->num_ecs_params * sizeof(double)); - - - for (i = 0; i < react->num_species; i++) { - states_for_reaction[i] = (double*) malloc(react->num_regions * sizeof(double)); - states_for_reaction_dx[i] = (double*) malloc(react->num_regions * sizeof(double)); - result_array[i] = (double*) malloc(react->num_regions * sizeof(double)); - result_array_dx[i] = (double*) malloc(react->num_regions * sizeof(double)); + react->ecsindex[i] = react->ecs_grid[i]->react_offsets[react->ecs_offset_index[i]]; } - for (i = 0; i < react->num_params; i++) - params_for_reaction[i] = (double*) malloc(react->num_regions * sizeof(double)); for (segment = 0; segment < react->num_segments; segment++) { if (react->vptrs != NULL) @@ -1559,22 +1540,22 @@ void solve_reaction(ICSReactions* react, for (i = 0; i < react->num_species; i++) { for (j = 0; j < react->num_regions; j++) { if (react->state_idx[segment][i][j] != SPECIES_ABSENT) { - states_for_reaction[i][j] = states[react->state_idx[segment][i][j]]; - states_for_reaction_dx[i][j] = states_for_reaction[i][j]; + react->states_for_reaction[i][j] = states[react->state_idx[segment][i][j]]; + react->states_for_reaction_dx[i][j] = react->states_for_reaction[i][j]; } else { - states_for_reaction[i][j] = SPECIES_ABSENT; - states_for_reaction_dx[i][j] = states_for_reaction[i][j]; + react->states_for_reaction[i][j] = SPECIES_ABSENT; + react->states_for_reaction_dx[i][j] = react->states_for_reaction[i][j]; } } - memset(result_array[i], 0, react->num_regions * sizeof(double)); - memset(result_array_dx[i], 0, react->num_regions * sizeof(double)); + memset(react->result_array[i], 0, react->num_regions * sizeof(double)); + memset(react->result_array_dx[i], 0, react->num_regions * sizeof(double)); } for (k = 0; i < react->num_species + react->num_params; i++, k++) { for (j = 0; j < react->num_regions; j++) { if (react->state_idx[segment][i][j] != SPECIES_ABSENT) { - params_for_reaction[k][j] = states[react->state_idx[segment][i][j]]; + react->params_for_reaction[k][j] = states[react->state_idx[segment][i][j]]; } else { - params_for_reaction[k][j] = SPECIES_ABSENT; + react->params_for_reaction[k][j] = SPECIES_ABSENT; } } } @@ -1583,35 +1564,35 @@ void solve_reaction(ICSReactions* react, for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) { if (cvode_states != NULL) { - ecs_states_for_reaction[i] = cvode_states[react->ecs_index[segment][i]]; + react->ecs_states_for_reaction[i] = cvode_states[react->ecs_index[segment][i]]; } else { - ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]); + react->ecs_states_for_reaction[i] = *(react->ecs_state[segment][i]); } - ecs_states_for_reaction_dx[i] = ecs_states_for_reaction[i]; + react->ecs_states_for_reaction_dx[i] = react->ecs_states_for_reaction[i]; } } for (k = 0; i < react->num_ecs_species + react->num_ecs_params; i++, k++) { if (react->ecs_state[segment][i] != NULL) { - ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]); + react->ecs_params_for_reaction[k] = *(react->ecs_state[segment][i]); } } if (react->num_ecs_species > 0) { - memset(ecs_result, 0, react->num_ecs_species * sizeof(double)); - memset(ecs_result_dx, 0, react->num_ecs_species * sizeof(double)); + bzero(react->ecs_result, react->num_ecs_species * sizeof(double)); + bzero(react->ecs_result_dx, react->num_ecs_species * sizeof(double)); } for (i = 0; i < react->num_mult; i++) { - mc_mult[i] = react->mc_multiplier[i][segment]; + react->mc_mult[i] = react->mc_multiplier[i][segment]; } - react->reaction(states_for_reaction, - params_for_reaction, - result_array, - mc_mult, - ecs_states_for_reaction, - ecs_params_for_reaction, - ecs_result, + react->reaction(react->states_for_reaction, + react->params_for_reaction, + react->result_array, + react->mc_mult, + react->ecs_states_for_reaction, + react->ecs_params_for_reaction, + react->ecs_result, NULL, v); @@ -1620,24 +1601,24 @@ void solve_reaction(ICSReactions* react, for (j = 0; j < react->num_regions; j++) { if (react->state_idx[segment][i][j] != SPECIES_ABSENT) { if (bval == NULL) - b->elem(idx) = dt * result_array[i][j]; + b->elem(idx) = dt * react->result_array[i][j]; else b->elem(idx) = bval[react->state_idx[segment][i][j]]; // set up the changed states array - states_for_reaction_dx[i][j] += dx; + react->states_for_reaction_dx[i][j] += dx; /* TODO: Handle approximating the Jacobian at a function upper * limit, e.g. acos(1) */ - react->reaction(states_for_reaction_dx, - params_for_reaction, - result_array_dx, - mc_mult, - ecs_states_for_reaction, - ecs_params_for_reaction, - ecs_result_dx, + react->reaction(react->states_for_reaction_dx, + react->params_for_reaction, + react->result_array_dx, + react->mc_mult, + react->ecs_states_for_reaction, + react->ecs_params_for_reaction, + react->ecs_result_dx, NULL, v); @@ -1645,25 +1626,26 @@ void solve_reaction(ICSReactions* react, for (jac_j = 0; jac_j < react->num_regions; jac_j++) { // pd is our Jacobian approximated if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) { - pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) / + pd = (react->result_array_dx[jac_i][jac_j] - + react->result_array[jac_i][jac_j]) / dx; *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd; jac_idx += 1; } - result_array_dx[jac_i][jac_j] = 0; + react->result_array_dx[jac_i][jac_j] = 0; } } for (jac_i = 0; jac_i < react->num_ecs_species; jac_i++) { // pd is our Jacobian approximated if (react->ecs_state[segment][jac_i] != NULL) { - pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx; + pd = (react->ecs_result_dx[jac_i] - react->ecs_result[jac_i]) / dx; *jacobian->mep(jac_idx, idx) = -dt * pd; jac_idx += 1; } - ecs_result_dx[jac_i] = 0; + react->ecs_result_dx[jac_i] = 0; } // reset dx array - states_for_reaction_dx[i][j] -= dx; + react->states_for_reaction_dx[i][j] -= dx; idx++; } } @@ -1673,24 +1655,24 @@ void solve_reaction(ICSReactions* react, for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) { if (bval == NULL) - b->elem(idx) = dt * ecs_result[i]; + b->elem(idx) = dt * react->ecs_result[i]; else b->elem(idx) = cvode_b[react->ecs_index[segment][i]]; // set up the changed states array - ecs_states_for_reaction_dx[i] += dx; + react->ecs_states_for_reaction_dx[i] += dx; /* TODO: Handle approximating the Jacobian at a function upper * limit, e.g. acos(1) */ - react->reaction(states_for_reaction, - params_for_reaction, - result_array_dx, - mc_mult, - ecs_states_for_reaction_dx, - ecs_params_for_reaction, - ecs_result_dx, + react->reaction(react->states_for_reaction, + react->params_for_reaction, + react->result_array_dx, + react->mc_mult, + react->ecs_states_for_reaction_dx, + react->ecs_params_for_reaction, + react->ecs_result_dx, NULL, v); @@ -1698,7 +1680,9 @@ void solve_reaction(ICSReactions* react, for (jac_j = 0; jac_j < react->num_regions; jac_j++) { // pd is our Jacobian approximated if (react->state_idx[segment][jac_i][jac_j] != SPECIES_ABSENT) { - pd = (result_array_dx[jac_i][jac_j] - result_array[jac_i][jac_j]) / dx; + pd = (react->result_array_dx[jac_i][jac_j] - + react->result_array[jac_i][jac_j]) / + dx; *jacobian->mep(jac_idx, idx) = -dt * pd; jac_idx += 1; } @@ -1707,14 +1691,14 @@ void solve_reaction(ICSReactions* react, for (jac_i = 0; jac_i < react->num_ecs_species; jac_i++) { // pd is our Jacobian approximated if (react->ecs_state[segment][jac_i] != NULL) { - pd = (ecs_result_dx[jac_i] - ecs_result[jac_i]) / dx; + pd = (react->ecs_result_dx[jac_i] - react->ecs_result[jac_i]) / dx; *jacobian->mep(jac_idx, idx) = (idx == jac_idx) - dt * pd; jac_idx += 1; } else { *jacobian->mep(idx, idx) = 1.0; } // reset dx array - ecs_states_for_reaction_dx[i] -= dx; + react->ecs_states_for_reaction_dx[i] -= dx; } idx++; } @@ -1734,7 +1718,8 @@ void solve_reaction(ICSReactions* react, } for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) - react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = x->elem(jac_idx++); + react->ecs_grid[i]->all_reaction_states[react->ecsindex[i]++] = x->elem( + jac_idx++); // cvode_b[react->ecs_index[segment][i]] = x->elem(jac_idx++); } } else // fixed-step @@ -1748,34 +1733,11 @@ void solve_reaction(ICSReactions* react, } for (i = 0; i < react->num_ecs_species; i++) { if (react->ecs_state[segment][i] != NULL) - react->ecs_grid[i]->all_reaction_states[ecsindex[i]++] = x->elem(jac_idx++); + react->ecs_grid[i]->all_reaction_states[react->ecsindex[i]++] = x->elem( + jac_idx++); } } } - free(ecsindex); - for (i = 0; i < react->num_species; i++) { - free(states_for_reaction[i]); - free(states_for_reaction_dx[i]); - free(result_array[i]); - free(result_array_dx[i]); - } - for (i = 0; i < react->num_params; i++) - free(params_for_reaction[i]); - if (react->num_mult > 0) - free(mc_mult); - free(states_for_reaction_dx); - free(states_for_reaction); - free(params_for_reaction); - free(result_array); - free(result_array_dx); - if (react->num_ecs_species > 0) { - free(ecs_states_for_reaction); - free(ecs_states_for_reaction_dx); - free(ecs_result); - free(ecs_result_dx); - } - if (react->num_ecs_params > 0) - free(ecs_params_for_reaction); } void do_ics_reactions(double* states, double* b, double* cvode_states, double* cvode_b) { diff --git a/src/nrnpython/rxd.h b/src/nrnpython/rxd.h index 7371bca493..41591d836d 100644 --- a/src/nrnpython/rxd.h +++ b/src/nrnpython/rxd.h @@ -62,6 +62,20 @@ typedef struct ICSReactions { double** mc_multiplier; int* mc_flux_idx; double** vptrs; + double** states_for_reaction; + double** states_for_reaction_dx; + double** params_for_reaction; + double** result_array; + double** result_array_dx; + double* mc_mult; + double** flux; + + double* ecs_states_for_reaction; + double* ecs_states_for_reaction_dx; + double* ecs_params_for_reaction; + double* ecs_result; + double* ecs_result_dx; + int* ecsindex; struct ICSReactions* next; } ICSReactions; From e7f15e570bf0bf383517c19431f8c6f129aca613 Mon Sep 17 00:00:00 2001 From: adam newton Date: Thu, 7 Dec 2023 09:20:07 -0500 Subject: [PATCH 2/2] Replaced bzero with memset for windows build. --- src/nrnpython/rxd.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nrnpython/rxd.cpp b/src/nrnpython/rxd.cpp index d343ec0d2b..e6819ff8f3 100644 --- a/src/nrnpython/rxd.cpp +++ b/src/nrnpython/rxd.cpp @@ -1430,7 +1430,7 @@ void get_reaction_rates(ICSReactions* react, double* states, double* rates, doub if (_membrane_flux) { for (i = 0; i < react->icsN; i++) - bzero(react->flux[i], react->num_regions * sizeof(double)); + memset(react->flux[i], 0, react->num_regions * sizeof(double)); } for (segment = 0; segment < react->num_segments; segment++) { for (i = 0; i < react->num_species; i++) { @@ -1578,8 +1578,8 @@ void solve_reaction(ICSReactions* react, } if (react->num_ecs_species > 0) { - bzero(react->ecs_result, react->num_ecs_species * sizeof(double)); - bzero(react->ecs_result_dx, react->num_ecs_species * sizeof(double)); + memset(react->ecs_result, 0, react->num_ecs_species * sizeof(double)); + memset(react->ecs_result_dx, 0, react->num_ecs_species * sizeof(double)); } for (i = 0; i < react->num_mult; i++) {