Hi Andrew, Leaving aside the issue of the API for the time being, I did some additional experiments in this area.
> To resolve this > issue some time ago I implemented another factorization of the basis > based on Schur complement (please see comments in src/bflib/scf.h), > where B0 = L0 * U0 is not changed, so B = B0 can be easily restored as > in case of "eta file". I did try this, and it seems to work OK. The attached two patches implement it on top of the three original patches I sent (pcost1.patch to pcost3.patch). The first one just makes the pseudocost initialisation use Schur complement updates, to use as a baseline since the pseudocosts calculated may be different now. The second ones introduces a bfd_reset() function that restores the original factorisation and uses this instead of bfd_copy(). Unfortunately, it looks like that the speed gain from using bfd_reset does not offset the overhead of Schur complement updates, so Forrest-Tomlin update with bfd_copy seems to be the preferred option here. > This approach, however, limits > the number of simplex iteration (say, to 100-200), since if > refactorization is needed, to keep B0 = L0 * U0 we should stop. I suppose this means that a new simplex flag may be needed, to stop when re-factorisation is required. Best Regards, Chris Matrakidis
commit 9734f4815fda152a20491a617818badb631adbf6 Author: Chris Matrakidis <[email protected]> Date: Sat Jan 14 16:24:45 2017 +0200 use schur complement update in pcost calculation diff --git a/src/glpios09.c b/src/glpios09.c index 727761c..255e681 100644 --- a/src/glpios09.c +++ b/src/glpios09.c @@ -447,12 +447,19 @@ static double eval_degrad(glp_tree *T, int j, double bnd) double degrad; /* prepare lp */ if (lp == NULL) - { /* the current basis must be optimal */ + { glp_bfcp bfcp; + /* the current basis must be optimal */ xassert(glp_get_status(P) == GLP_OPT); /* create a copy of mip */ lp = glp_create_prob(); glp_copy_prob(lp, P, 0); csa->work = lp; + glp_get_bfcp(lp, &bfcp); + if (bfcp.type & GLP_BF_FT) + { bfcp.type &= ~GLP_BF_FT; + bfcp.type |= GLP_BF_BG; + } + glp_set_bfcp(lp, &bfcp); /* factorize and store */ glp_factorize(lp); csa->bfd = bfd_create_it();
commit d9944e1c142a09e20ee0cf58fa660b3d09206475 Author: Chris Matrakidis <[email protected]> Date: Sat Jan 14 19:36:43 2017 +0200 use bfd_reset in pcost initialisation diff --git a/src/bfd.c b/src/bfd.c index 89d547a..a085298 100644 --- a/src/bfd.c +++ b/src/bfd.c @@ -514,6 +514,16 @@ int bfd_update(BFD *bfd, int j, int len, const int ind[], const double return ret; } +void bfd_reset(BFD *bfd) +{ /* reset LP basis factorization to initial state */ + int ret; + xassert(bfd->valid); + xassert(bfd->type == 2); + scfint_reset(bfd->u.scfi); + bfd->upd_cnt = 0; + return; +} + int bfd_get_count(BFD *bfd) { /* determine factorization update count */ return bfd->upd_cnt; diff --git a/src/bfd.h b/src/bfd.h index 3f63099..d652c55 100644 --- a/src/bfd.h +++ b/src/bfd.h @@ -94,6 +94,10 @@ int bfd_update(BFD *bfd, int j, int len, const int ind[], const double val[]); /* update LP basis factorization */ +#define bfd_reset _glp_bfd_reset +void bfd_reset(BFD *bfd); +/* reset LP basis factorization to initial state */ + #define bfd_get_count _glp_bfd_get_count int bfd_get_count(BFD *bfd); /* determine factorization update count */ diff --git a/src/bflib/scfint.c b/src/bflib/scfint.c index 25aed12..6e572b0 100644 --- a/src/bflib/scfint.c +++ b/src/bflib/scfint.c @@ -135,6 +135,26 @@ int scfint_factorize(SCFINT *fi, int n, int (*col)(void *info, int j, return ret; } +void scfint_reset(SCFINT *fi) +{ /* reset SC-factorization to initial state */ + int k; + int n = fi->scf.n0; + xassert(n > 0); + /* initialize SC-factorization */ + fi->scf.n = n; + fi->scf.nn = 0; + fi->scf.rr_ref = sva_alloc_vecs(fi->scf.sva, fi->scf.ifu.n_max); + fi->scf.ss_ref = sva_alloc_vecs(fi->scf.sva, fi->scf.ifu.n_max); + fi->scf.ifu.n = 0; + for (k = 1; k <= n; k++) + { fi->scf.pp_ind[k] = k; + fi->scf.pp_inv[k] = k; + fi->scf.qq_ind[k] = k; + fi->scf.qq_inv[k] = k; + } + return; +} + int scfint_update(SCFINT *fi, int upd, int j, int len, const int ind[], const double val[]) { /* update SC-factorization after replacing j-th column of A */ diff --git a/src/bflib/scfint.h b/src/bflib/scfint.h index 67ca65d..45142a6 100644 --- a/src/bflib/scfint.h +++ b/src/bflib/scfint.h @@ -63,6 +63,10 @@ int scfint_factorize(SCFINT *fi, int n, int (*col)(void *info, int j, int ind[], double val[]), void *info); /* compute SC-factorization of specified matrix A */ +#define scfint_reset _glp_scfint_reset +void scfint_reset(SCFINT *fi); +/* reset SC-factorization to initial state */ + #define scfint_update _glp_scfint_update int scfint_update(SCFINT *fi, int upd, int j, int len, const int ind[], const double val[]); diff --git a/src/glpios09.c b/src/glpios09.c index 255e681..4d3f5d6 100644 --- a/src/glpios09.c +++ b/src/glpios09.c @@ -406,14 +406,12 @@ struct csa over all up_cnt[j] subproblems */ glp_prob *work; /* working problem to avoid unnecessary initializations */ - BFD *bfd; - /* factorization copy to avoid identical factorizations*/ int *head; - /* basis header copy to avoid identical factorizations */ + /* basis header copy to restore factorization */ int *r_bind; - /* copy of row bind values to avoid identical factorizations */ + /* copy of row bind values to restore factorization */ int *c_bind; - /* copy of column bind values to avoid identical factorizations */ + /* copy of column bind values to restore factorization */ }; void *ios_pcost_init(glp_tree *tree) @@ -430,7 +428,6 @@ void *ios_pcost_init(glp_tree *tree) csa->dn_sum[j] = csa->up_sum[j] = 0.0; } csa->work = NULL; - csa->bfd = NULL; return csa; } @@ -462,8 +459,6 @@ static double eval_degrad(glp_tree *T, int j, double bnd) glp_set_bfcp(lp, &bfcp); /* factorize and store */ glp_factorize(lp); - csa->bfd = bfd_create_it(); - bfd_copy(csa->bfd, lp->bfd); /* copy factorization */ csa->head = xcalloc(1+P->m, sizeof(int)); csa->c_bind = xcalloc(1+P->n, sizeof(int)); @@ -490,7 +485,7 @@ static double eval_degrad(glp_tree *T, int j, double bnd) lp->col[jjj]->dual = P->col[jjj]->dual; lp->col[jjj]->bind = csa->c_bind[jjj]; } - bfd_copy(lp->bfd, csa->bfd); + bfd_reset(lp->bfd); lp->valid = 1; } /* fix column x[j] at specified value */ @@ -723,7 +718,6 @@ done: *_next = sel; glp_delete_prob(csa->work); csa->work = NULL; /* clean up stored factorization for next round */ - bfd_delete_it(csa->bfd); xfree(csa->c_bind); xfree(csa->r_bind); xfree(csa->head);
_______________________________________________ Help-glpk mailing list [email protected] https://lists.gnu.org/mailman/listinfo/help-glpk
