Andrew,
Since you are preparing a new GLPK release, I'm sending you three
patches to consider, that improve the speed of pseudocost
initialization.
The first reuses the same glp_prob in each round to avoid unnecessary
copying, the second introduces a bfd_copy function and the third uses
bfd_copy to avoid recalculating the exact same factorisation in every
simplex call during each initialisation round.
These patches do not affect the execution order of the code, only the
initialisation speed: in all tests I did, I get the exact same number
of nodes and simplex iterations
Best Regards,
Chris Matrakidis
diff --git a/src/glpios09.c b/src/glpios09.c
index d87578c..b9e6b5a 100644
--- a/src/glpios09.c
+++ b/src/glpios09.c
@@ -403,6 +403,8 @@ struct csa
double *up_sum; /* double up_sum[1+n]; */
/* up_sum[j] is the sum of per unit degradations of the objective
over all up_cnt[j] subproblems */
+ glp_prob *work;
+ /* working problem to avoid unnecessary initializations */
};
void *ios_pcost_init(glp_tree *tree)
@@ -418,24 +420,44 @@ void *ios_pcost_init(glp_tree *tree)
{ csa->dn_cnt[j] = csa->up_cnt[j] = 0;
csa->dn_sum[j] = csa->up_sum[j] = 0.0;
}
+ csa->work = NULL;
return csa;
}
-static double eval_degrad(glp_prob *P, int j, double bnd)
+static double eval_degrad(glp_tree *T, int j, double bnd)
{ /* compute degradation of the objective on fixing x[j] at given
value with a limited number of dual simplex iterations */
/* this routine fixes column x[j] at specified value bnd,
solves resulting LP, and returns a lower bound to degradation
of the objective, degrad >= 0 */
- glp_prob *lp;
+ struct csa *csa = T->pcost;
+ glp_prob *lp = csa->work, *P = T->mip;
glp_smcp parm;
- int ret;
+ int i, jjj, ret;
double degrad;
- /* the current basis must be optimal */
- xassert(glp_get_status(P) == GLP_OPT);
- /* create a copy of P */
- lp = glp_create_prob();
- glp_copy_prob(lp, P, 0);
+ /* prepare lp */
+ if (lp == NULL)
+ { /* 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;
+ }
+ else
+ { /* restore original values */
+ for (i = 1; i <= P->m; i++)
+ { lp->row[i]->stat = P->row[i]->stat;
+ lp->row[i]->prim = P->row[i]->prim;
+ lp->row[i]->dual = P->row[i]->dual;
+ }
+ for (jjj = 1; jjj <= P->n; jjj++)
+ { lp->col[jjj]->stat = P->col[jjj]->stat;
+ lp->col[jjj]->prim = P->col[jjj]->prim;
+ lp->col[jjj]->dual = P->col[jjj]->dual;
+ }
+ lp->valid = 0;
+ }
/* fix column x[j] at specified value */
glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
/* try to solve resulting LP */
@@ -478,8 +500,9 @@ static double eval_degrad(glp_prob *P, int j, double bnd)
{ /* the simplex solver failed */
degrad = 0.0;
}
- /* delete the copy of P */
- glp_delete_prob(lp);
+ /* 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);
return degrad;
}
@@ -550,7 +573,7 @@ static double eval_psi(glp_tree *T, int j, int brnch)
if (csa->dn_cnt[j] == 0)
{ /* initialize down pseudocost */
beta = T->mip->col[j]->prim;
- degrad = eval_degrad(T->mip, j, floor(beta));
+ degrad = eval_degrad(T, j, floor(beta));
if (degrad == DBL_MAX)
{ psi = DBL_MAX;
goto done;
@@ -565,7 +588,7 @@ static double eval_psi(glp_tree *T, int j, int brnch)
if (csa->up_cnt[j] == 0)
{ /* initialize up pseudocost */
beta = T->mip->col[j]->prim;
- degrad = eval_degrad(T->mip, j, ceil(beta));
+ degrad = eval_degrad(T, j, ceil(beta));
if (degrad == DBL_MAX)
{ psi = DBL_MAX;
goto done;
@@ -604,9 +627,11 @@ int ios_pcost_branch(glp_tree *T, int *_next)
#endif
int j, jjj, sel;
double beta, psi, d1, d2, d, dmax;
+ struct csa *csa;
/* initialize the working arrays */
if (T->pcost == NULL)
T->pcost = ios_pcost_init(T);
+ csa = T->pcost;
/* nothing has been chosen so far */
jjj = 0, dmax = -1.0;
/* go through the list of branching candidates */
@@ -658,6 +683,11 @@ int ios_pcost_branch(glp_tree *T, int *_next)
jjj = branch_mostf(T, &sel);
}
done: *_next = sel;
+ if (csa->work != NULL)
+ { /* clean up working problem for next round */
+ glp_delete_prob(csa->work);
+ csa->work = NULL;
+ }
return jjj;
}
diff --git a/src/bfd.c b/src/bfd.c
index dece376..89d547a 100644
--- a/src/bfd.c
+++ b/src/bfd.c
@@ -541,4 +541,44 @@ void bfd_delete_it(BFD *bfd)
return;
}
+void bfd_copy(BFD *dst, BFD *src)
+{ /* copy LP basis factorization */
+#ifdef GLP_DEBUG
+ int i;
+#endif
+ xassert(dst != NULL);
+ dst->valid = src->valid;
+ dst->type = src->type;
+ dst->parm = src->parm;
+ dst->upd_cnt = src->upd_cnt;
+ dst->i_norm = src->i_norm;
+#ifdef GLP_DEBUG
+ if (dst->B != NULL)
+ spm_delete_mat(dst->B);
+ dst->B = spm_create_mat(src->B->m, src->B->n);
+ for (i = 1; i <= src->B->m; i++)
+ { SPME *e;
+ for (e = src->B->row[i]; e != NULL; e = e->r_next)
+ spm_new_elem(dst->B, i, e->j, e->val);
+ }
+#endif
+ switch (dst->type)
+ { case 0:
+ break;
+ case 1:
+ if (dst->u.fhvi == NULL)
+ dst->u.fhvi = fhvint_create();
+ fhvint_copy(dst->u.fhvi, src->u.fhvi);
+ break;
+ case 2:
+ if (dst->u.scfi == NULL)
+ dst->u.scfi = scfint_create(src->u.scfi->scf.type);
+ scfint_copy(dst->u.scfi, src->u.scfi);
+ break;
+ default:
+ xassert(dst != dst);
+ }
+ return;
+}
+
/* eof */
diff --git a/src/bfd.h b/src/bfd.h
index 0ef4c02..3f63099 100644
--- a/src/bfd.h
+++ b/src/bfd.h
@@ -102,6 +102,10 @@ int bfd_get_count(BFD *bfd);
void bfd_delete_it(BFD *bfd);
/* delete LP basis factorization */
+#define bfd_copy _glp_bfd_copy
+void bfd_copy(BFD *dst, BFD *src);
+/* copy LP basis factorization */
+
#endif
/* eof */
diff --git a/src/bflib/btfint.c b/src/bflib/btfint.c
index 378d3a8..f4c0123 100644
--- a/src/bflib/btfint.c
+++ b/src/bflib/btfint.c
@@ -404,4 +404,100 @@ void btfint_delete(BTFINT *fi)
return;
}
+void btfint_copy(BTFINT *dst, BTFINT *src)
+{ /* copy interface to BT-factorization */
+ int n, n_max;
+ n = src->btf->n;
+ n_max = dst->n_max;
+ if (n > n_max)
+ { n_max = dst->n_max = src->n_max;
+ if (dst->btf != NULL)
+ { tfree(dst->btf->pp_ind);
+ tfree(dst->btf->pp_inv);
+ tfree(dst->btf->qq_ind);
+ tfree(dst->btf->qq_inv);
+ tfree(dst->btf->beg);
+ tfree(dst->btf->vr_piv);
+ tfree(dst->btf->p1_ind);
+ tfree(dst->btf->p1_inv);
+ tfree(dst->btf->q1_ind);
+ tfree(dst->btf->q1_inv);
+ }
+ else
+ dst->btf = talloc(1, BTF);
+ dst->btf->pp_ind = talloc(1+n_max, int);
+ dst->btf->pp_inv = talloc(1+n_max, int);
+ dst->btf->qq_ind = talloc(1+n_max, int);
+ dst->btf->qq_inv = talloc(1+n_max, int);
+ dst->btf->beg = talloc(1+n_max+1, int);
+ dst->btf->vr_piv = talloc(1+n_max, double);
+ dst->btf->p1_ind = talloc(1+n_max, int);
+ dst->btf->p1_inv = talloc(1+n_max, int);
+ dst->btf->q1_ind = talloc(1+n_max, int);
+ dst->btf->q1_inv = talloc(1+n_max, int);
+ if (dst->sgf != NULL)
+ { tfree(dst->sgf->rs_head);
+ tfree(dst->sgf->rs_prev);
+ tfree(dst->sgf->rs_next);
+ tfree(dst->sgf->cs_head);
+ tfree(dst->sgf->cs_prev);
+ tfree(dst->sgf->cs_next);
+ tfree(dst->sgf->vr_max);
+ tfree(dst->sgf->flag);
+ tfree(dst->sgf->work);
+ }
+ else
+ dst->sgf = talloc(1, SGF);
+ dst->sgf->rs_head = talloc(1+n_max, int);
+ dst->sgf->rs_prev = talloc(1+n_max, int);
+ dst->sgf->rs_next = talloc(1+n_max, int);
+ dst->sgf->cs_head = talloc(1+n_max, int);
+ dst->sgf->cs_prev = talloc(1+n_max, int);
+ dst->sgf->cs_next = talloc(1+n_max, int);
+ dst->sgf->vr_max = talloc(1+n_max, double);
+ dst->sgf->flag = talloc(1+n_max, char);
+ dst->sgf->work = talloc(1+n_max, double);
+ }
+ dst->valid = src->valid;
+ dst->sva_n_max = src->sva_n_max;
+ dst->sva_size = src->sva_size;
+ dst->delta_n0 = src->delta_n0;
+ dst->delta_n = src->delta_n;
+ dst->sgf_piv_tol = src->sgf_piv_tol;
+ dst->sgf_piv_lim = src->sgf_piv_lim;
+ dst->sgf_suhl = src->sgf_suhl;
+ dst->sgf_eps_tol = src->sgf_eps_tol;
+ if (dst->sva == NULL)
+ dst->sva = sva_create_area(src->sva->n_max, src->sva->size);
+ sva_copy_area(dst->sva, src->sva);
+ dst->btf->n = n;
+ dst->btf->num =src->btf->num;
+ dst->btf->ar_ref = src->btf->ar_ref;
+ dst->btf->ac_ref = src->btf->ac_ref;
+ dst->btf->fr_ref = src->btf->fr_ref;
+ dst->btf->fc_ref = src->btf->fc_ref;
+ dst->btf->vr_ref = src->btf->vr_ref;
+ dst->btf->vc_ref = src->btf->vc_ref;
+ memcpy(dst->btf->vr_piv, src->btf->vr_piv,
+ (1+n) * sizeof(double));
+ memcpy(dst->btf->pp_ind, src->btf->pp_ind,(1+n) * sizeof(int));
+ memcpy(dst->btf->pp_inv, src->btf->pp_inv, (1+n) * sizeof(int));
+ memcpy(dst->btf->qq_ind, src->btf->qq_ind, (1+n) * sizeof(int));
+ memcpy(dst->btf->qq_inv, src->btf->qq_inv, (1+n) * sizeof(int));
+ memcpy(dst->btf->beg, src->btf->beg,
+ (1+src->btf->num+1) * sizeof(int));
+ memcpy(dst->btf->p1_ind, src->btf->p1_ind, (1+n) * sizeof(int));
+ memcpy(dst->btf->p1_inv, src->btf->p1_inv, (1+n) * sizeof(int));
+ memcpy(dst->btf->q1_ind, src->btf->q1_ind, (1+n) * sizeof(int));
+ memcpy(dst->btf->q1_inv, src->btf->q1_inv, (1+n) * sizeof(int));
+ dst->btf->sva = dst->sva;
+ dst->sgf->luf = src->sgf->luf;
+ dst->sgf->updat = src->sgf->updat;
+ dst->sgf->piv_tol = src->sgf->piv_tol;
+ dst->sgf->piv_lim = src->sgf->piv_lim;
+ dst->sgf->suhl = src->sgf->suhl;
+ dst->sgf->eps_tol = src->sgf->eps_tol;
+ return;
+}
+
/* eof */
diff --git a/src/bflib/btfint.h b/src/bflib/btfint.h
index 8d0e70e..0bfb9cf 100644
--- a/src/bflib/btfint.h
+++ b/src/bflib/btfint.h
@@ -68,6 +68,10 @@ int btfint_factorize(BTFINT *fi, int n, int (*col)(void *info, int j,
void btfint_delete(BTFINT *fi);
/* delete interface to BT-factorization */
+#define btfint_copy _glp_btfint_copy
+void btfint_copy(BTFINT *dst, BTFINT *src);
+/* copy interface to BT-factorization */
+
#endif
/* eof */
diff --git a/src/bflib/fhvint.c b/src/bflib/fhvint.c
index a21b71c..40852f8 100644
--- a/src/bflib/fhvint.c
+++ b/src/bflib/fhvint.c
@@ -165,4 +165,51 @@ void fhvint_delete(FHVINT *fi)
return;
}
+void fhvint_copy(FHVINT *dst, FHVINT *src)
+{ /* copy interface to FHV-factorization */
+ int nfs_max, n_max, n, nfs;
+ nfs_max = dst->nfs_max;
+ nfs = src->fhv.nfs;
+ n_max = (dst->lufi != NULL) ? dst->lufi->n_max : 0;
+ n = src->lufi->luf->n;
+ if(nfs_max < nfs)
+ { nfs_max = src->nfs_max;
+ if (dst->fhv.hh_ind != NULL)
+ { tfree(dst->fhv.hh_ind);
+ dst->fhv.hh_ind = NULL;
+ }
+ }
+ if(n_max < n)
+ { n_max = src->lufi->n_max;
+ if (dst->fhv.p0_ind != NULL)
+ { tfree(dst->fhv.p0_ind);
+ dst->fhv.p0_ind = NULL;
+ }
+ if (dst->fhv.p0_inv != NULL)
+ { tfree(dst->fhv.p0_inv);
+ dst->fhv.p0_inv = NULL;
+ }
+ }
+ if (dst->fhv.hh_ind == NULL)
+ dst->fhv.hh_ind = talloc(1+nfs_max, int);
+ if (dst->fhv.p0_ind == NULL)
+ dst->fhv.p0_ind = talloc(1+n_max, int);
+ if (dst->fhv.p0_inv == NULL)
+ dst->fhv.p0_inv = talloc(1+n_max, int);
+ if (src->fhv.hh_ind != NULL)
+ memcpy(dst->fhv.hh_ind, src->fhv.hh_ind,
+ (1+nfs) * sizeof(int));
+ if (src->fhv.p0_ind != NULL)
+ memcpy(dst->fhv.p0_ind, src->fhv.p0_ind, (1+n) * sizeof(int));
+ if (src->fhv.p0_inv != NULL)
+ memcpy(dst->fhv.p0_inv, src->fhv.p0_inv, (1+n) * sizeof(int));
+ dst->valid = src->valid;
+ dst->fhv.nfs_max = nfs_max;
+ dst->fhv.nfs = nfs;
+ dst->fhv.hh_ref = src->fhv.hh_ref;
+ dst->nfs_max = nfs_max;
+ lufint_copy(dst->lufi, src->lufi);
+ dst->fhv.luf = dst->lufi->luf;
+ return;
+}
/* eof */
diff --git a/src/bflib/fhvint.h b/src/bflib/fhvint.h
index 000829c..73a78b6 100644
--- a/src/bflib/fhvint.h
+++ b/src/bflib/fhvint.h
@@ -73,6 +73,10 @@ double fhvint_estimate(FHVINT *fi);
void fhvint_delete(FHVINT *fi);
/* delete interface to FHV-factorization */
+#define fhvint_copy _glp_fhvint_copy
+void fhvint_copy(FHVINT *dst, FHVINT *src);
+/* copy interface to FHV-factorization */
+
#endif
/* eof */
diff --git a/src/bflib/lufint.c b/src/bflib/lufint.c
index 7cd0092..07b67f9 100644
--- a/src/bflib/lufint.c
+++ b/src/bflib/lufint.c
@@ -179,4 +179,82 @@ void lufint_delete(LUFINT *fi)
return;
}
+void lufint_copy(LUFINT *dst, LUFINT *src)
+{ /* copy interface to LU-factorization */
+ int n, n_max;
+ n = src->luf->n;
+ n_max = dst->n_max;
+ if (n > n_max)
+ { n_max = dst->n_max = src->n_max;
+ if (dst->luf != NULL)
+ { tfree(dst->luf->vr_piv);
+ tfree(dst->luf->pp_ind);
+ tfree(dst->luf->pp_inv);
+ tfree(dst->luf->qq_ind);
+ tfree(dst->luf->qq_inv);
+ }
+ else
+ dst->luf = talloc(1, LUF);
+ dst->luf->vr_piv = talloc(1+n_max, double);
+ dst->luf->pp_ind = talloc(1+n_max, int);
+ dst->luf->pp_inv = talloc(1+n_max, int);
+ dst->luf->qq_ind = talloc(1+n_max, int);
+ dst->luf->qq_inv = talloc(1+n_max, int);
+ if (dst->sgf != NULL)
+ { tfree(dst->sgf->rs_head);
+ tfree(dst->sgf->rs_prev);
+ tfree(dst->sgf->rs_next);
+ tfree(dst->sgf->cs_head);
+ tfree(dst->sgf->cs_prev);
+ tfree(dst->sgf->cs_next);
+ tfree(dst->sgf->vr_max);
+ tfree(dst->sgf->flag);
+ tfree(dst->sgf->work);
+ }
+ else
+ dst->sgf = talloc(1, SGF);
+ dst->sgf->rs_head = talloc(1+n_max, int);
+ dst->sgf->rs_prev = talloc(1+n_max, int);
+ dst->sgf->rs_next = talloc(1+n_max, int);
+ dst->sgf->cs_head = talloc(1+n_max, int);
+ dst->sgf->cs_prev = talloc(1+n_max, int);
+ dst->sgf->cs_next = talloc(1+n_max, int);
+ dst->sgf->vr_max = talloc(1+n_max, double);
+ dst->sgf->flag = talloc(1+n_max, char);
+ dst->sgf->work = talloc(1+n_max, double);
+ }
+ dst->valid = src->valid;
+ dst->sva_n_max = src->sva_n_max;
+ dst->sva_size = src->sva_size;
+ dst->delta_n0 = src->delta_n0;
+ dst->delta_n = src->delta_n;
+ dst->sgf_updat = src->sgf_updat;
+ dst->sgf_piv_tol = src->sgf_piv_tol;
+ dst->sgf_piv_lim = src->sgf_piv_lim;
+ dst->sgf_suhl = src->sgf_suhl;
+ dst->sgf_eps_tol = src->sgf_eps_tol;
+ if (dst->sva == NULL)
+ dst->sva = sva_create_area(src->sva->n_max, src->sva->size);
+ sva_copy_area(dst->sva, src->sva);
+ dst->luf->n = n;
+ dst->luf->fr_ref = src->luf->fr_ref;
+ dst->luf->fc_ref = src->luf->fc_ref;
+ dst->luf->vr_ref = src->luf->vr_ref;
+ dst->luf->vc_ref = src->luf->vc_ref;
+ memcpy(dst->luf->vr_piv, src->luf->vr_piv,
+ (1+n) * sizeof(double));
+ memcpy(dst->luf->pp_ind, src->luf->pp_ind, (1+n) * sizeof(int));
+ memcpy(dst->luf->pp_inv, src->luf->pp_inv, (1+n) * sizeof(int));
+ memcpy(dst->luf->qq_ind, src->luf->qq_ind, (1+n) * sizeof(int));
+ memcpy(dst->luf->qq_inv, src->luf->qq_inv, (1+n) * sizeof(int));
+ dst->luf->sva = dst->sva;
+ dst->sgf->luf = dst->luf;
+ dst->sgf->updat = src->sgf->updat;
+ dst->sgf->piv_tol = src->sgf->piv_tol;
+ dst->sgf->piv_lim = src->sgf->piv_lim;
+ dst->sgf->suhl = src->sgf->suhl;
+ dst->sgf->eps_tol = src->sgf->eps_tol;
+ return;
+}
+
/* eof */
diff --git a/src/bflib/lufint.h b/src/bflib/lufint.h
index b3ad5b6..f57caa8 100644
--- a/src/bflib/lufint.h
+++ b/src/bflib/lufint.h
@@ -68,6 +68,10 @@ int lufint_factorize(LUFINT *fi, int n, int (*col)(void *info, int j,
void lufint_delete(LUFINT *fi);
/* delete interface to LU-factorization */
+#define lufint_copy _glp_lufint_copy
+void lufint_copy(LUFINT *dst, LUFINT *src);
+/* copy interface to LU-factorization */
+
#endif
/* eof */
diff --git a/src/bflib/scfint.c b/src/bflib/scfint.c
index 06aa8f7..25aed12 100644
--- a/src/bflib/scfint.c
+++ b/src/bflib/scfint.c
@@ -252,4 +252,103 @@ void scfint_delete(SCFINT *fi)
return;
}
+void scfint_copy(SCFINT *dst, SCFINT *src)
+{ /* copy interface to SC-factorization */
+ int n0_max, old_n0_max, nn_max, n0, nn, k, n;
+ xassert(src->scf.type == dst->scf.type);
+ switch (src->scf.type)
+ { case 1:
+ old_n0_max = dst->u.lufi->n_max;
+ lufint_copy(dst->u.lufi, src->u.lufi);
+ n0_max = dst->u.lufi->n_max;
+ n0 = dst->u.lufi->luf->n;
+ dst->scf.sva = dst->u.lufi->sva;
+ dst->scf.a0.luf = dst->u.lufi->luf;
+ break;
+ case 2:
+ old_n0_max = dst->u.btfi->n_max;
+ btfint_copy(dst->u.btfi, src->u.btfi);
+ n0_max = dst->u.btfi->n_max;
+ n0 = dst->u.btfi->btf->n;
+ dst->scf.sva = dst->u.btfi->sva;
+ dst->scf.a0.btf = dst->u.btfi->btf;
+ break;
+ default:
+ xassert(dst != dst);
+ }
+ /* allocate/reallocate arrays, if necessary */
+ if (old_n0_max < n0_max)
+ { if (dst->w1 != NULL)
+ tfree(dst->w1);
+ if (dst->w2 != NULL)
+ tfree(dst->w2);
+ if (dst->w3 != NULL)
+ tfree(dst->w3);
+ dst->w1 = talloc(1+n0_max, double);
+ dst->w2 = talloc(1+n0_max, double);
+ dst->w3 = talloc(1+n0_max, double);
+ }
+ nn_max = dst->scf.nn_max;
+ nn = src->scf.nn;
+ if (nn > nn_max)
+ { nn_max = src->scf.nn_max;
+ if (dst->scf.ifu.f != NULL)
+ tfree(dst->scf.ifu.f);
+ if (dst->scf.ifu.u != NULL)
+ tfree(dst->scf.ifu.u);
+ dst->scf.ifu.f = talloc(nn_max * nn_max, double);
+ dst->scf.ifu.u = talloc(nn_max * nn_max, double);
+ }
+ n = src->scf.ifu.n;
+ if (src->scf.ifu.f != NULL)
+ for (k = 0; k < n; k++)
+ memcpy(dst->scf.ifu.f + k, src->scf.ifu.f + k,
+ n * sizeof(double));
+ if (src->scf.ifu.u != NULL)
+ for (k = 0; k < n; k++)
+ memcpy(dst->scf.ifu.u + k, src->scf.ifu.u + k,
+ n * sizeof(double));
+ if (old_n0_max < n0_max || dst->scf.nn_max != nn_max)
+ { if (dst->scf.pp_ind != NULL)
+ tfree(dst->scf.pp_ind);
+ if (dst->scf.pp_inv != NULL)
+ tfree(dst->scf.pp_inv);
+ if (dst->scf.qq_ind != NULL)
+ tfree(dst->scf.qq_ind);
+ if (dst->scf.qq_inv != NULL)
+ tfree(dst->scf.qq_inv);
+ if (dst->w4 != NULL)
+ tfree(dst->w4);
+ if (dst->w5 != NULL)
+ tfree(dst->w5);
+ dst->scf.pp_ind = talloc(1+n0_max+nn_max, int);
+ dst->scf.pp_inv = talloc(1+n0_max+nn_max, int);
+ dst->scf.qq_ind = talloc(1+n0_max+nn_max, int);
+ dst->scf.qq_inv = talloc(1+n0_max+nn_max, int);
+ dst->w4 = talloc(1+n0_max+nn_max, double);
+ dst->w5 = talloc(1+n0_max+nn_max, double);
+ }
+ k = src->scf.n0 + src->scf.nn;
+ if (src->scf.pp_ind != NULL)
+ memcpy(dst->scf.pp_ind, src->scf.pp_ind, (1+k) * sizeof(int));
+ if (src->scf.pp_inv != NULL)
+ memcpy(dst->scf.pp_inv, src->scf.pp_inv, (1+k) * sizeof(int));
+ if (src->scf.qq_ind != NULL)
+ memcpy(dst->scf.qq_ind, src->scf.qq_ind, (1+k) * sizeof(int));
+ if (src->scf.qq_inv != NULL)
+ memcpy(dst->scf.qq_inv, src->scf.qq_inv, (1+k) * sizeof(int));
+ /* copy SC-factorization */
+ dst->scf.n = src->scf.n;
+ dst->scf.n0 = src->scf.n0;
+ dst->scf.nn_max = nn_max;
+ dst->scf.nn = src->scf.nn;
+ dst->scf.rr_ref = src->scf.rr_ref;
+ dst->scf.ss_ref = src->scf.ss_ref;
+ dst->scf.ifu.n_max = nn_max;
+ dst->scf.ifu.n = src->scf.ifu.n;
+ dst->nn_max = src->nn_max;
+ dst->valid = src->valid;
+ return;
+}
+
/* eof */
diff --git a/src/bflib/scfint.h b/src/bflib/scfint.h
index 3e56355..67ca65d 100644
--- a/src/bflib/scfint.h
+++ b/src/bflib/scfint.h
@@ -84,6 +84,10 @@ double scfint_estimate(SCFINT *fi);
void scfint_delete(SCFINT *fi);
/* delete interface to SC-factorization */
+#define scfint_copy _glp_scfint_copy
+void scfint_copy(SCFINT *dst, SCFINT *src);
+/* copy interface to SC-factorization */
+
#endif
/* eof */
diff --git a/src/bflib/sva.c b/src/bflib/sva.c
index e6a675c..bcf8336 100644
--- a/src/bflib/sva.c
+++ b/src/bflib/sva.c
@@ -569,4 +569,73 @@ void sva_delete_area(SVA *sva)
return;
}
+/***********************************************************************
+* sva_copy_area - copy sparse vector area (SVA)
+*
+* This routine copies a sparse vector area (SVA) into another avoiding
+* new memory allocations if possible */
+
+void sva_copy_area(SVA *dst, SVA *src)
+{ int size, req, delta, k, n, n_max, ptr, len, r_ptr;
+ size = dst->size;
+ r_ptr = src->r_ptr;
+ req = src->m_ptr + src->size + 1 - r_ptr;
+ if (req >= size)
+ { size = dst->size = src->size;
+ tfree(dst->ind);
+ tfree(dst->val);
+ dst->ind = talloc(1+size, int);
+ dst->val = talloc(1+size, double);
+ }
+ delta = size - src->size;
+ n = src->n;
+ n_max = dst->n_max;
+ if (n >= n_max)
+ { n_max = dst->n_max = src->n_max;
+ tfree(dst->ptr);
+ tfree(dst->len);
+ tfree(dst->cap);
+ tfree(dst->prev);
+ tfree(dst->next);
+ dst->ptr = talloc(1+n_max, int);
+ dst->len = talloc(1+n_max, int);
+ dst->cap = talloc(1+n_max, int);
+ dst->prev = talloc(1+n_max, int);
+ dst->next = talloc(1+n_max, int);
+ }
+ dst->n = src->n;
+ dst->m_ptr = src->m_ptr;
+ dst->r_ptr = src->r_ptr + delta;
+ dst->head = src->head;
+ dst->tail = src->tail;
+ dst->talky = src->talky;
+ memcpy(dst->ptr, src->ptr, (1+n) * sizeof(int));
+ memcpy(dst->len, src->len, (1+n) * sizeof(int));
+ memcpy(dst->cap, src->cap, (1+n) * sizeof(int));
+ memcpy(dst->prev, src->prev, (1+n) * sizeof(int));
+ memcpy(dst->next, src->next, (1+n) * sizeof(int));
+ /* use linked list to copy the vectors in the left part */
+ for (k = src->head; k != 0; k = src->next[k])
+ { ptr = src->ptr[k];
+ len = src->len[k];
+ memcpy(dst->ind + ptr, src->ind + ptr, len * sizeof(int));
+ memcpy(dst->val + ptr, src->val + ptr, len * sizeof(double));
+ }
+ /* copy right part and adjust pointers for delta */
+ if (r_ptr <= src->size)
+ { len = src->size + 1 - r_ptr;
+ memcpy(dst->ind + r_ptr + delta, src->ind + r_ptr,
+ len * sizeof(int));
+ memcpy(dst->val + r_ptr + delta, src->val + r_ptr,
+ len * sizeof(double));
+ if (delta != 0)
+ { for (k = 1; k <= n; k++)
+ { if (dst->ptr[k] >= r_ptr)
+ dst->ptr[k] += delta;
+ }
+ }
+ }
+ return;
+}
+
/* eof */
diff --git a/src/bflib/sva.h b/src/bflib/sva.h
index 0eab317..f632a6d 100644
--- a/src/bflib/sva.h
+++ b/src/bflib/sva.h
@@ -156,6 +156,10 @@ void sva_check_area(SVA *sva);
void sva_delete_area(SVA *sva);
/* delete sparse vector area (SVA) */
+#define sva_copy_area _glp_sva_copy_area
+void sva_copy_area(SVA *dst, SVA *src);
+/* copy sparse vector area (SVA) */
+
#endif
/* eof */
diff --git a/src/glpios09.c b/src/glpios09.c
index b9e6b5a..727761c 100644
--- a/src/glpios09.c
+++ b/src/glpios09.c
@@ -24,6 +24,7 @@
#include "env.h"
#include "glpios.h"
+#include "bfd.h"
/***********************************************************************
* NAME
@@ -405,6 +406,14 @@ 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 */
+ int *r_bind;
+ /* copy of row bind values to avoid identical factorizations */
+ int *c_bind;
+ /* copy of column bind values to avoid identical factorizations */
};
void *ios_pcost_init(glp_tree *tree)
@@ -421,6 +430,7 @@ 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;
}
@@ -443,6 +453,20 @@ static double eval_degrad(glp_tree *T, int j, double bnd)
lp = glp_create_prob();
glp_copy_prob(lp, P, 0);
csa->work = lp;
+ /* 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));
+ csa->r_bind = xcalloc(1+P->m, sizeof(int));
+ for (i = 1; i <= P->m; i++)
+ { csa->r_bind[i] = lp->row[i]->bind;
+ csa->head[i] = lp->head[i];
+ }
+ for (jjj = 1; jjj <= P->n; jjj++)
+ csa->c_bind[jjj] = lp->col[jjj]->bind;
}
else
{ /* restore original values */
@@ -450,13 +474,17 @@ static double eval_degrad(glp_tree *T, int j, double bnd)
{ lp->row[i]->stat = P->row[i]->stat;
lp->row[i]->prim = P->row[i]->prim;
lp->row[i]->dual = P->row[i]->dual;
+ lp->row[i]->bind = csa->r_bind[i];
+ lp->head[i] = csa->head[i];
}
for (jjj = 1; jjj <= P->n; jjj++)
{ lp->col[jjj]->stat = P->col[jjj]->stat;
lp->col[jjj]->prim = P->col[jjj]->prim;
lp->col[jjj]->dual = P->col[jjj]->dual;
+ lp->col[jjj]->bind = csa->c_bind[jjj];
}
- lp->valid = 0;
+ bfd_copy(lp->bfd, csa->bfd);
+ lp->valid = 1;
}
/* fix column x[j] at specified value */
glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
@@ -687,6 +715,11 @@ done: *_next = sel;
{ /* clean up working problem for next round */
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);
}
return jjj;
}
_______________________________________________
Help-glpk mailing list
[email protected]
https://lists.gnu.org/mailman/listinfo/help-glpk
