Skip to content

Commit

Permalink
added multithread jacobian evaluation, slight improvement (10%)
Browse files Browse the repository at this point in the history
  • Loading branch information
@hhijazi
hhijazi committed Jan 20, 2018
1 parent fa7ebd6 commit bbe2f55
Show file tree
Hide file tree
Showing 4 changed files with 100 additions and 28 deletions.
3 changes: 2 additions & 1 deletion examples/MINLP/Power/ACOPF/ACOPF_main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ using namespace gravity;

int main (int argc, char * argv[])
{
string fname = "../data_sets/Power/nesta_case5_pjm.m", mtype = "ACPOL";
string fname = "../data_sets/Power/nesta_case5_pjm.m", mtype = "ACRECT";
fname = "/Users/hlh/Dropbox/Work/Dev/pglib-opf/sad/pglib_opf_case9241_pegase__sad.m";
DebugOn("argv[0] =" << argv[0] << endl);
string path = argv[0];
int output = 0;
Expand Down
2 changes: 1 addition & 1 deletion include/gravity/constraint.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ namespace gravity {
string _name = "no_name";

public:

unsigned _jac_cstr_idx; /* Firt index of the corresponding non-zero values in the Jacobian */
unsigned _id = -1;
ConstraintType _ctype = leq; /**< Constraint type: leq, geq or eq */
double _rhs = 0;
Expand Down
2 changes: 1 addition & 1 deletion include/gravity/model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ namespace gravity {
size_t _nnz_g = 0; /* Number of non zeros in the Jacobian */
size_t _nnz_h = 0; /* Number of non zeros in the Hessian */
size_t _nnz_g_obj = 0; /* Number of non zeros in the Objective gradient */
vector<double> _jac_vals; /* Jacobian values stored in sparse format */
vector<double> _jac_vals; /* Jacobian values stored in sparse format */
vector<double> _obj_grad_vals; /* Objective gradient values stored in sparse format */
vector<double> _hess_vals; /* Hessian values stored in sparse format */
map<unsigned, param_*> _params; /**< Sorted map pointing to all parameters contained in this model */
Expand Down
121 changes: 96 additions & 25 deletions src/model.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -933,6 +933,7 @@ void Model::fill_in_jac_nnz(int* iRow , int* jCol){
for(auto& c_p :_cons)
{
c = c_p.second.get();
c->_jac_cstr_idx = idx;
auto nb_ins = c->_nb_instances;
for (auto &v_p: c->get_vars()){
v = v_p.second.first.get();
Expand Down Expand Up @@ -960,6 +961,53 @@ void Model::fill_in_jac_nnz(int* iRow , int* jCol){
}
}


void compute_jac(vector<Constraint*>& vec, double* res, unsigned i, unsigned j, bool first_call, vector<double>& jac_vals){
size_t cid = 0;
unique_id vid;
Constraint* c = NULL;
param_* v = NULL;
shared_ptr<func_> dfdx;
auto idx = vec[i]->_jac_cstr_idx;
for (unsigned id = i; id < j; id++) {
c = vec[id];
auto nb_ins = c->_nb_instances;
if (c->is_linear() && !first_call) {
// if (false) {
DebugOff("Linear constraint, using stored jacobian!\n");
for (unsigned i = 0; i<nb_ins; i++) {
for (unsigned j = 0; j<c->get_nb_vars(i); j++) {
res[idx] = jac_vals[idx];
idx++;
}
}
}
else {
for (auto &v_p: c->get_vars()){
v = v_p.second.first.get();
vid = v->_unique_id;
dfdx = c->get_stored_derivative(vid);
for (int inst = 0; inst< nb_ins; inst++){
cid = c->_id+inst;
if (v->_is_vector) {
auto dim = v->get_dim(inst);
for (int j = 0; j<dim; j++) {
res[idx] = dfdx->eval(inst,j);
jac_vals[idx] = res[idx];
idx++;
}
}
else {
res[idx] = dfdx->eval(inst);
jac_vals[idx] = res[idx];
idx++;
}
}
}
}
}
}

/* Fill the nonzero values in the jacobian */
void Model::fill_in_jac(const double* x , double* res, bool new_x){
// if (!_first_call_jac && (!new_x || _type==lin_m)) { /* No need to recompute jacobian for linear models */
Expand All @@ -980,6 +1028,29 @@ void Model::fill_in_jac(const double* x , double* res, bool new_x){
Constraint* c = NULL;
param_* v = NULL;
shared_ptr<func_> dfdx;
vector<Constraint*> cons;
if (_type!=nlin_m) {//Quadratic or Linear
for(auto& c_p: _cons)
{
c = c_p.second.get();
c->_new = false;
cons.push_back(c);
}
unsigned nr_threads = 4;
vector<thread> threads;
/* Split cons into nr_threads parts */
vector<int> limits = bounds(nr_threads, cons.size());
/* Launch all threads in parallel */
for (int i = 0; i < nr_threads; ++i) {
threads.push_back(thread(compute_jac, ref(cons), res, limits[i], limits[i+1], _first_call_jac, ref(_jac_vals)));
}
/* Join the threads with the main thread */
for(auto &t : threads){
t.join();
}
_first_call_jac = false;
return;
}
for(auto& c_p :_cons)
{
c = c_p.second.get();
Expand All @@ -995,7 +1066,7 @@ void Model::fill_in_jac(const double* x , double* res, bool new_x){
}
}
else {
if (_type==nlin_m) {
// if (_type==nlin_m) {
for (auto &v_p: c->get_vars()){
v = v_p.second.first.get();
vid = v->_unique_id;
Expand Down Expand Up @@ -1046,30 +1117,30 @@ void Model::fill_in_jac(const double* x , double* res, bool new_x){
}
}
}
}
else {
for (auto &v_p: c->get_vars()){
v = v_p.second.first.get();
vid = v->_unique_id;
dfdx = c->get_stored_derivative(vid);
for (int inst = 0; inst< nb_ins; inst++){
cid = c->_id+inst;
if (v->_is_vector) {
auto dim = v->get_dim(inst);
for (int j = 0; j<dim; j++) {
res[idx] = dfdx->eval(inst,j);
_jac_vals[idx] = res[idx];
idx++;
}
}
else {
res[idx] = dfdx->eval(inst);
_jac_vals[idx] = res[idx];
idx++;
}
}
}
}
// }
// else {
// for (auto &v_p: c->get_vars()){
// v = v_p.second.first.get();
// vid = v->_unique_id;
// dfdx = c->get_stored_derivative(vid);
// for (int inst = 0; inst< nb_ins; inst++){
// cid = c->_id+inst;
// if (v->_is_vector) {
// auto dim = v->get_dim(inst);
// for (int j = 0; j<dim; j++) {
// res[idx] = dfdx->eval(inst,j);
// _jac_vals[idx] = res[idx];
// idx++;
// }
// }
// else {
// res[idx] = dfdx->eval(inst);
// _jac_vals[idx] = res[idx];
// idx++;
// }
// }
// }
// }
}
}
_first_call_jac = false;
Expand Down

0 comments on commit bbe2f55

Please sign in to comment.