Re: Patches to improve pseudocost initialisatiion speed
Chris Matrakidis <[email protected]>
| Newsgroups | gmane.comp.gnu.glpk |
|---|---|
| Message-ID | <CAEaJ3S7KU1u=aOEPvZhyAH7-k_XJexSm9i+7XA6AOEa_DYeS_Q@mail.gmail.com> |
Andrew, >> I have a draft patch that introduces an internal API for keeping the >> dual simplex state between calls and adjusting it for new bounds. It >> works fine and speeds up pseudocost initialisation for all >> configurations of the solver (i.e. combinations of USE_AT, EXCL and >> SHIFT). However, I'm not happy with the API design, so I don't >> consider it ready for submission. If you think it may be helpful, I >> will send it. > > Please post your patch to the list with minimal comments for archiving. Here it is, broken down in two pieces. The first patch defines the API and all necessary changes and uses it in pseudocost initialisation. This is on top of my previous pseudocost initialisation patches but not the other simplex patches (but should apply without conflicts). The second patch is a tentative initial effort to update the row wise representation of N. A more worked out version, should probably have reset_nt() in spxnt.c and choose between build_nt() and reset_nt() depending on the number of simplex iterations since the last call. Best Regards, Chris Matrakidis _______________________________________________ Help-glpk mailing list [email protected] https://lists.gnu.org/mailman/listinfo/help-glpk
dual_API1.patch
(application/octet-stream, 15.9 KB)
commit f056c1f9981e12c95a6057a277ab70c4262834d5 Author: Chris Matrakidis <[email protected]> Date: Mon Jan 9 14:05:12 2017 +0200 dual simplex intrnal API for pcost diff --git a/src/glpios09.c b/src/glpios09.c index 727761c..1e2bde3 100644 --- a/src/glpios09.c +++ b/src/glpios09.c @@ -25,6 +25,7 @@ #include "env.h" #include "glpios.h" #include "bfd.h" +#include "simplex.h" /*********************************************************************** * NAME @@ -414,6 +415,8 @@ struct csa /* copy of row bind values to avoid identical factorizations */ int *c_bind; /* copy of column bind values to avoid identical factorizations */ + void *duals; + /* internal state for dual simplex reuse */ }; void *ios_pcost_init(glp_tree *tree) @@ -431,6 +434,7 @@ void *ios_pcost_init(glp_tree *tree) } csa->work = NULL; csa->bfd = NULL; + csa->duals = NULL; return csa; } @@ -488,14 +492,18 @@ static double eval_degrad(glp_tree *T, int j, double bnd) } /* fix column x[j] at specified value */ glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd); + /* initialise solver state if empty or update bound */ + if (csa->duals == NULL) + csa->duals = spy_dual_init(lp, 1); + else + spy_dual_update_col_bnd(lp, j, csa->duals); /* try to solve resulting LP */ glp_init_smcp(&parm); parm.msg_lev = GLP_MSG_OFF; parm.meth = GLP_DUAL; parm.it_lim = 30; parm.out_dly = 1000; - parm.meth = GLP_DUAL; - ret = glp_simplex(lp, &parm); + ret = spy_dual_main(lp, &parm, csa->duals); if (ret == 0 || ret == GLP_EITLIM) { if (glp_get_prim_stat(lp) == GLP_NOFEAS) { /* resulting LP has no primal feasible solution */ @@ -531,6 +539,7 @@ static double eval_degrad(glp_tree *T, int j, double bnd) /* restore column x[j] type and bounds */ glp_set_col_bnds(lp, j, P->col[j]->type, P->col[j]->lb, P->col[j]->ub); + spy_dual_update_col_bnd(lp, j, csa->duals); return degrad; } @@ -721,6 +730,11 @@ done: *_next = sel; xfree(csa->r_bind); xfree(csa->head); } + if (csa->duals !=NULL) + { /* clean up solver state for next time */ + spy_dual_free(csa->duals); + csa->duals = NULL; + } return jjj; } diff --git a/src/simplex/simplex.h b/src/simplex/simplex.h index 4820588..3fa998f 100644 --- a/src/simplex/simplex.h +++ b/src/simplex/simplex.h @@ -34,6 +34,22 @@ int spx_primal(glp_prob *P, const glp_smcp *parm); int spy_dual(glp_prob *P, const glp_smcp *parm); /* driver to dual simplex method */ +#define spy_dual_init _glp_spy_dual_init +struct csa *spy_dual_init(glp_prob *P, int delta); +/* initialise dual simplex state to allow reuse */ + +#define spy_dual_main _glp_spy_dual_main +int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa); +/* dual simplex using previously initialised state */ + +#define spy_dual_free _glp_spy_dual_free +void spy_dual_free(struct csa *csa); +/* free dual simplex state */ + +#define spy_dual_update_col_bnd _glp_spy_dual_update_col_bnd +void spy_dual_update_col_bnd(glp_prob *P, int j, struct csa *csa); +/* update column bounds in dual simplex state */ + #endif /* eof */ diff --git a/src/simplex/spydual.c b/src/simplex/spydual.c index e41baf2..5e1c828 100644 --- a/src/simplex/spydual.c +++ b/src/simplex/spydual.c @@ -146,6 +146,12 @@ struct csa #endif int p_stat, d_stat; /* primal and dual solution statuses */ + int *map; /* int map[1+P->m+P->n]; */ + /* mapping of original LP variables to working LP ones */ + int *daeh; /* int map[1+n]; */ + /* array to hold inverse of the working basis header */ + double *delta; /* double delta[1+P->n]; */ + /* if SHIFT is set, holds shift amount (NULL if not used) */ /*--------------------------------------------------------------*/ /* control parameters (see struct glp_smcp) */ int msg_lev; @@ -1770,23 +1776,49 @@ fini: int spy_dual(glp_prob *P, const glp_smcp *parm) { /* driver to dual simplex method */ - struct csa csa_, *csa = &csa_; - SPXLP lp; -#if USE_AT - SPXAT at; -#else - SPXNT nt; -#endif - SPYSE se; - int ret, *map, *daeh; - /* build working LP and its initial basis */ + struct csa *csa; + int ret; + csa = spy_dual_init(P, 0); + ret = spy_dual_main(P, parm, csa); + spy_dual_free(csa); + /* return to calling program */ + return ret; +} + +/*********************************************************************** +* spy_dual_init - initialise dual simplex state to allow reuse +* +* This routine initialises the dual simplex state and returns a pointer +* to the csa structure that can be used in subsequent calls to dual +* simplex functions. */ + +static double * spy_init_delta(glp_prob *P, int *map) +{ int m = P->m; + int n = P->n; + int k; + double *delta = talloc(1+m+n, double); + for (k = 1; k <= n; k++) + { GLPCOL *col = P->col[k]; + if (map[m+k] == 0 || col->type == GLP_FR) + delta[k] = 0; + else if (map[m+k] > 0) + delta[k] = col->lb / col->sjj; + else + delta[k] = col->ub / col->sjj; + } + return delta; +} + +struct csa *spy_dual_init(glp_prob *P, int use_delta) +{ /* initialise dual simplex state to allow reuse */ + struct csa *csa; + csa = talloc(1, struct csa); memset(csa, 0, sizeof(struct csa)); - csa->lp = &lp; + csa->lp = talloc(1, SPXLP); spx_init_lp(csa->lp, P, EXCL); spx_alloc_lp(csa->lp); - map = talloc(1+P->m+P->n, int); - spx_build_lp(csa->lp, P, EXCL, SHIFT, map); - spx_build_basis(csa->lp, P, map); + csa->map = talloc(1+P->m+P->n, int); + spx_build_lp(csa->lp, P, EXCL, SHIFT, csa->map); switch (P->dir) { case GLP_MIN: csa->dir = +1; @@ -1798,46 +1830,13 @@ int spy_dual(glp_prob *P, const glp_smcp *parm) xassert(P != P); } csa->orig_b = talloc(1+csa->lp->m, double); - memcpy(csa->orig_b, csa->lp->b, (1+csa->lp->m) * sizeof(double)); #if PERTURB csa->orig_c = talloc(1+csa->lp->n, double); - memcpy(csa->orig_c, csa->lp->c, (1+csa->lp->n) * sizeof(double)); #endif csa->orig_l = talloc(1+csa->lp->n, double); - memcpy(csa->orig_l, csa->lp->l, (1+csa->lp->n) * sizeof(double)); csa->orig_u = talloc(1+csa->lp->n, double); - memcpy(csa->orig_u, csa->lp->u, (1+csa->lp->n) * sizeof(double)); -#if USE_AT - /* build matrix A in row-wise format */ - csa->at = &at; - csa->nt = NULL; - spx_alloc_at(csa->lp, csa->at); - spx_build_at(csa->lp, csa->at); -#else - /* build matrix N in row-wise format for initial basis */ - csa->at = NULL; - csa->nt = &nt; - spx_alloc_nt(csa->lp, csa->nt); - spx_init_nt(csa->lp, csa->nt); - spx_build_nt(csa->lp, csa->nt); -#endif - /* allocate and initialize working components */ - csa->phase = 0; csa->beta = talloc(1+csa->lp->m, double); - csa->beta_st = 0; csa->d = talloc(1+csa->lp->n-csa->lp->m, double); - csa->d_st = 0; - switch (parm->pricing) - { case GLP_PT_STD: - csa->se = NULL; - break; - case GLP_PT_PSE: - csa->se = &se; - spy_alloc_se(csa->lp, csa->se); - break; - default: - xassert(parm != parm); - } #if 0 /* 30/III-2016 */ csa->list = talloc(1+csa->lp->m, int); csa->trow = talloc(1+csa->lp->n-csa->lp->m, double); @@ -1856,6 +1855,66 @@ int spy_dual(glp_prob *P, const glp_smcp *parm) fvs_alloc_vec(&csa->wrow, csa->lp->n-csa->lp->m); fvs_alloc_vec(&csa->wcol, csa->lp->m); #endif + csa->daeh = talloc(1+csa->lp->n, int); + if (SHIFT && use_delta != 0) + csa->delta = spy_init_delta(P, csa->map); + return csa; +} + +/*********************************************************************** +* spy_dual_main - dual simplex using previously initialised state +* +* This routine is the main body of the dual simplex method. The csa +* structure is assumed to have been initialised by a call to +* spy_dual_init, however repeated calls to spy_dual_main can be done +* reusing the same csa structure. +* +* On exit the return values are as in spy_dual described earlier */ + +int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa) +{ /* dual simplex using previously initialised state */ + int ret; + /* build working LP initial basis */ + spx_build_basis(csa->lp, P, csa->map); + memcpy(csa->orig_b, csa->lp->b, (1+csa->lp->m) * sizeof(double)); +#if PERTURB + memcpy(csa->orig_c, csa->lp->c, (1+csa->lp->n) * sizeof(double)); +#endif + memcpy(csa->orig_l, csa->lp->l, (1+csa->lp->n) * sizeof(double)); + memcpy(csa->orig_u, csa->lp->u, (1+csa->lp->n) * sizeof(double)); +#if USE_AT + /* build matrix A in row-wise format */ + if (csa->at == NULL) + { csa->at = talloc(1, SPXAT); + spx_alloc_at(csa->lp, csa->at); + spx_build_at(csa->lp, csa->at); + } +#else + /* build matrix N in row-wise format for initial basis */ + if (csa->nt == NULL) + { csa->nt = talloc(1, SPXNT); + spx_alloc_nt(csa->lp, csa->nt); + spx_init_nt(csa->lp, csa->nt); + } + spx_build_nt(csa->lp, csa->nt); +#endif + /* initialize working components */ + csa->phase = 0; + csa->beta_st = 0; + csa->d_st = 0; + switch (parm->pricing) + { case GLP_PT_STD: + break; + case GLP_PT_PSE: + if (csa->se == NULL) + { csa->se = talloc(1, SPYSE); + spy_alloc_se(csa->lp, csa->se); + } + csa->se->valid = 0; + break; + default: + xassert(parm != parm); + } /* initialize control parameters */ csa->msg_lev = parm->msg_lev; csa->dualp = (parm->meth == GLP_DUALP); @@ -1876,7 +1935,8 @@ int spy_dual(glp_prob *P, const glp_smcp *parm) case GLP_RT_HAR: break; case GLP_RT_FLIP: - csa->bp = talloc(1+csa->lp->n-csa->lp->m, SPYBP); + if (csa->bp == NULL) + csa->bp = talloc(1+csa->lp->n-csa->lp->m, SPYBP); break; default: xassert(parm != parm); @@ -1927,16 +1987,14 @@ int spy_dual(glp_prob *P, const glp_smcp *parm) goto skip; /* convert working LP basis to original LP basis and store it to * problem object */ - daeh = talloc(1+csa->lp->n, int); - spx_store_basis(csa->lp, P, map, daeh); + spx_store_basis(csa->lp, P, csa->map, csa->daeh); /* compute simplex multipliers for final basic solution found by * the solver */ spx_eval_pi(csa->lp, csa->work); /* convert working LP solution to original LP solution and store * it to problem object */ - spx_store_sol(csa->lp, P, SHIFT, map, daeh, csa->beta, csa->work, - csa->d); - tfree(daeh); + spx_store_sol(csa->lp, P, SHIFT, csa->map, csa->daeh, csa->beta, + csa->work, csa->d); /* save simplex iteration count */ P->it_cnt = csa->it_cnt; /* report auxiliary/structural variable causing unboundedness */ @@ -1949,16 +2007,26 @@ int spy_dual(glp_prob *P, const glp_smcp *parm) xassert(1 <= k && k <= csa->lp->n); /* convert to number of original variable */ for (kk = 1; kk <= P->m + P->n; kk++) - { if (abs(map[kk]) == k) + { if (abs(csa->map[kk]) == k) { P->some = kk; break; } } xassert(P->some != 0); } -skip: /* deallocate working objects and arrays */ +skip: /* return to calling program */ + return ret >= 0 ? ret : GLP_EFAIL; +} + +/*********************************************************************** +* spy_dual_free - free dual simplex state +* +* This routine clears and deallocates the csa structure. */ + +void spy_dual_free(struct csa *csa) +{ /* free dual simplex state */ spx_free_lp(csa->lp); - tfree(map); + tfree(csa->map); tfree(csa->orig_b); #if PERTURB tfree(csa->orig_c); @@ -1995,8 +2063,117 @@ skip: /* deallocate working objects and arrays */ fvs_free_vec(&csa->wrow); fvs_free_vec(&csa->wcol); #endif - /* return to calling program */ - return ret >= 0 ? ret : GLP_EFAIL; + tfree(csa->daeh); + if (csa->delta != NULL) + tfree(csa->delta); + tfree(csa->lp); + if (csa->at != NULL) + tfree(csa->at); + if (csa->nt != NULL) + tfree(csa->nt); + if (csa->se != NULL) + tfree(csa->se); + tfree(csa); + return; +} + +/*********************************************************************** +* spy_dual_update_col_bnd - update column bounds in dual simplex state +* +* This routine updates the csa structure for a changed bound in a +* structural variable of the original LP so that spy_dual_main can be +* called again without needing a new call to csa_dual_init. */ + +void spy_dual_update_col_bnd(glp_prob *P, int j, struct csa *csa) +{ /* update column bounds in dual simplex state */ + SPXLP *lp = csa->lp; + int *map = csa->map; + int m = lp->m; + int n = lp->n; + double *l = lp->l; + double *u = lp->u; + int k; + double delta; + GLPCOL *col = P->col[j]; + GLPAIJ *aij; + k = abs(map[m+j]); + xassert(k != 0); + /* update bounds of variable x[k] in working LP */ + switch (col->type) + { case GLP_FR: + l[k] = -DBL_MAX, u[k] = +DBL_MAX; + break; + case GLP_LO: + l[k] = col->lb / col->sjj, u[k] = +DBL_MAX; + break; + case GLP_UP: + l[k] = -DBL_MAX, u[k] = col->ub / col->sjj; + break; + case GLP_DB: + l[k] = col->lb / col->sjj, u[k] = col->ub / col->sjj; + xassert(l[k] != u[k]); + break; + case GLP_FX: + l[k] = u[k] = col->lb / col->sjj; + break; + default: + xassert(col != col); + } + if (SHIFT) + { /* shift bounds of variable x[k] */ + if (l[k] == -DBL_MAX && u[k] == +DBL_MAX) + { /* x[k] is unbounded variable */ + delta = 0.0; + } + else if (l[k] != -DBL_MAX && u[k] == +DBL_MAX) + { /* shift lower bound to zero */ + delta = l[k]; + l[k] = 0.0; + } + else if (l[k] == -DBL_MAX && u[k] != +DBL_MAX) + { /* shift upper bound to zero */ + map[m+j] = -k; + delta = u[k]; + u[k] = 0.0; + } + else if (l[k] != u[k]) + { /* x[k] is double bounded variable */ + if (fabs(l[k]) <= fabs(u[k])) + { /* shift lower bound to zero */ + delta = l[k]; + l[k] = 0.0, u[k] -= delta; + } + else + { /* shift upper bound to zero */ + map[m+j] = -k; + delta = u[k]; + l[k] -= delta, u[k] = 0.0; + } + xassert(l[k] != u[k]); + } + else + { /* shift fixed value to zero */ + delta = l[k]; + l[k] = u[k] = 0.0; + } + /* update all constraints and the objective function of + * working LP */ + if (delta != csa->delta[j]) + { int ptr, end; + int *A_ptr = lp->A_ptr; + int *A_ind = lp->A_ind; + double *A_val = lp->A_val; + double *b = lp->b; + double *c = lp->c; + ptr = A_ptr[k]; + end = A_ptr[k+1]; + for (; ptr < end; ptr++) + b[A_ind[ptr]] -= A_val[ptr] * (delta - csa->delta[j]); + c[0] += c[k] * (delta - csa->delta[j]); + csa->delta[j] = delta; + } + } + return; } /* eof */
dual_API2.patch
(application/octet-stream, 2.3 KB)
commit ce09e6791c74be80afa1102234e4dfe4d139cc1a Author: Chris Matrakidis <[email protected]> Date: Mon Jan 9 14:10:51 2017 +0200 reset row wise representation of N matrix if re-entering dual simplex diff --git a/src/simplex/spydual.c b/src/simplex/spydual.c index 5e1c828..d1d2493 100644 --- a/src/simplex/spydual.c +++ b/src/simplex/spydual.c @@ -1871,9 +1871,51 @@ struct csa *spy_dual_init(glp_prob *P, int use_delta) * * On exit the return values are as in spy_dual described earlier */ +static void reset_nt(glp_prob *P, struct csa *csa) +{ int i, j, k, jj; + SPXLP *lp = csa->lp; + int *map = csa->map; + int *head = lp->head; + int m = P->m; + int n = P->n; + int status; + /* remove columns that have changed */ + jj = 0; + for (i = 1; i <= m+n; i++) + { if ((k = map[i]) == 0) + continue; + if (k < 0) + k = -k; + status = (i <= m) ? P->row[i]->stat : P->col[i-m]->stat; + if (status == GLP_BS) + continue; + jj++; + if (head[m+jj] != k) + spx_nt_del_col(lp, csa->nt, jj, head[m+jj]); + } + /* restore original columns */ + jj = 0; + for (i = 1; i <= m+n; i++) + { if ((k = map[i]) == 0) + continue; + if (k < 0) + k = -k; + status = (i <= m) ? P->row[i]->stat : P->col[i-m]->stat; + if (status == GLP_BS) + continue; + jj++; + if (head[m+jj] != k) + spx_nt_add_col(lp, csa->nt, jj, k); + } + return; +} + int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa) { /* dual simplex using previously initialised state */ int ret; + /* reset row wise representation of N matrix if re-entering */ + if (csa->nt != NULL) + reset_nt(P, csa); /* build working LP initial basis */ spx_build_basis(csa->lp, P, csa->map); memcpy(csa->orig_b, csa->lp->b, (1+csa->lp->m) * sizeof(double)); @@ -1895,8 +1937,8 @@ int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa) { csa->nt = talloc(1, SPXNT); spx_alloc_nt(csa->lp, csa->nt); spx_init_nt(csa->lp, csa->nt); + spx_build_nt(csa->lp, csa->nt); } - spx_build_nt(csa->lp, csa->nt); #endif /* initialize working components */ csa->phase = 0;