diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/api/cpxbas.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/api/cpxbas.c | 269 |
1 files changed, 269 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/api/cpxbas.c b/test/monniaux/glpk-4.65/src/api/cpxbas.c new file mode 100644 index 00000000..e1c656a7 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/api/cpxbas.c @@ -0,0 +1,269 @@ +/* cpxbas.c (construct Bixby's initial LP basis) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2008-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: <mao@gnu.org>. +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see <http://www.gnu.org/licenses/>. +***********************************************************************/ + +#include "env.h" +#include "prob.h" + +struct var +{ /* structural variable */ + int j; + /* ordinal number */ + double q; + /* penalty value */ +}; + +static int CDECL fcmp(const void *ptr1, const void *ptr2) +{ /* this routine is passed to the qsort() function */ + struct var *col1 = (void *)ptr1, *col2 = (void *)ptr2; + if (col1->q < col2->q) return -1; + if (col1->q > col2->q) return +1; + return 0; +} + +static int get_column(glp_prob *lp, int j, int ind[], double val[]) +{ /* Bixby's algorithm assumes that the constraint matrix is scaled + such that the maximum absolute value in every non-zero row and + column is 1 */ + int k, len; + double big; + len = glp_get_mat_col(lp, j, ind, val); + big = 0.0; + for (k = 1; k <= len; k++) + if (big < fabs(val[k])) big = fabs(val[k]); + if (big == 0.0) big = 1.0; + for (k = 1; k <= len; k++) val[k] /= big; + return len; +} + +static void cpx_basis(glp_prob *lp) +{ /* main routine */ + struct var *C, *C2, *C3, *C4; + int m, n, i, j, jk, k, l, ll, t, n2, n3, n4, type, len, *I, *r, + *ind; + double alpha, gamma, cmax, temp, *v, *val; + xprintf("Constructing initial basis...\n"); + /* determine the number of rows and columns */ + m = glp_get_num_rows(lp); + n = glp_get_num_cols(lp); + /* allocate working arrays */ + C = xcalloc(1+n, sizeof(struct var)); + I = xcalloc(1+m, sizeof(int)); + r = xcalloc(1+m, sizeof(int)); + v = xcalloc(1+m, sizeof(double)); + ind = xcalloc(1+m, sizeof(int)); + val = xcalloc(1+m, sizeof(double)); + /* make all auxiliary variables non-basic */ + for (i = 1; i <= m; i++) + { if (glp_get_row_type(lp, i) != GLP_DB) + glp_set_row_stat(lp, i, GLP_NS); + else if (fabs(glp_get_row_lb(lp, i)) <= + fabs(glp_get_row_ub(lp, i))) + glp_set_row_stat(lp, i, GLP_NL); + else + glp_set_row_stat(lp, i, GLP_NU); + } + /* make all structural variables non-basic */ + for (j = 1; j <= n; j++) + { if (glp_get_col_type(lp, j) != GLP_DB) + glp_set_col_stat(lp, j, GLP_NS); + else if (fabs(glp_get_col_lb(lp, j)) <= + fabs(glp_get_col_ub(lp, j))) + glp_set_col_stat(lp, j, GLP_NL); + else + glp_set_col_stat(lp, j, GLP_NU); + } + /* C2 is a set of free structural variables */ + n2 = 0, C2 = C + 0; + for (j = 1; j <= n; j++) + { type = glp_get_col_type(lp, j); + if (type == GLP_FR) + { n2++; + C2[n2].j = j; + C2[n2].q = 0.0; + } + } + /* C3 is a set of structural variables having excatly one (lower + or upper) bound */ + n3 = 0, C3 = C2 + n2; + for (j = 1; j <= n; j++) + { type = glp_get_col_type(lp, j); + if (type == GLP_LO) + { n3++; + C3[n3].j = j; + C3[n3].q = + glp_get_col_lb(lp, j); + } + else if (type == GLP_UP) + { n3++; + C3[n3].j = j; + C3[n3].q = - glp_get_col_ub(lp, j); + } + } + /* C4 is a set of structural variables having both (lower and + upper) bounds */ + n4 = 0, C4 = C3 + n3; + for (j = 1; j <= n; j++) + { type = glp_get_col_type(lp, j); + if (type == GLP_DB) + { n4++; + C4[n4].j = j; + C4[n4].q = glp_get_col_lb(lp, j) - glp_get_col_ub(lp, j); + } + } + /* compute gamma = max{|c[j]|: 1 <= j <= n} */ + gamma = 0.0; + for (j = 1; j <= n; j++) + { temp = fabs(glp_get_obj_coef(lp, j)); + if (gamma < temp) gamma = temp; + } + /* compute cmax */ + cmax = (gamma == 0.0 ? 1.0 : 1000.0 * gamma); + /* compute final penalty for all structural variables within sets + C2, C3, and C4 */ + switch (glp_get_obj_dir(lp)) + { case GLP_MIN: temp = +1.0; break; + case GLP_MAX: temp = -1.0; break; + default: xassert(lp != lp); + } + for (k = 1; k <= n2+n3+n4; k++) + { j = C[k].j; + C[k].q += (temp * glp_get_obj_coef(lp, j)) / cmax; + } + /* sort structural variables within C2, C3, and C4 in ascending + order of penalty value */ + qsort(C2+1, n2, sizeof(struct var), fcmp); + for (k = 1; k < n2; k++) xassert(C2[k].q <= C2[k+1].q); + qsort(C3+1, n3, sizeof(struct var), fcmp); + for (k = 1; k < n3; k++) xassert(C3[k].q <= C3[k+1].q); + qsort(C4+1, n4, sizeof(struct var), fcmp); + for (k = 1; k < n4; k++) xassert(C4[k].q <= C4[k+1].q); + /*** STEP 1 ***/ + for (i = 1; i <= m; i++) + { type = glp_get_row_type(lp, i); + if (type != GLP_FX) + { /* row i is either free or inequality constraint */ + glp_set_row_stat(lp, i, GLP_BS); + I[i] = 1; + r[i] = 1; + } + else + { /* row i is equality constraint */ + I[i] = 0; + r[i] = 0; + } + v[i] = +DBL_MAX; + } + /*** STEP 2 ***/ + for (k = 1; k <= n2+n3+n4; k++) + { jk = C[k].j; + len = get_column(lp, jk, ind, val); + /* let alpha = max{|A[l,jk]|: r[l] = 0} and let l' be such + that alpha = |A[l',jk]| */ + alpha = 0.0, ll = 0; + for (t = 1; t <= len; t++) + { l = ind[t]; + if (r[l] == 0 && alpha < fabs(val[t])) + alpha = fabs(val[t]), ll = l; + } + if (alpha >= 0.99) + { /* B := B union {jk} */ + glp_set_col_stat(lp, jk, GLP_BS); + I[ll] = 1; + v[ll] = alpha; + /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */ + for (t = 1; t <= len; t++) + { l = ind[t]; + if (val[t] != 0.0) r[l]++; + } + /* continue to the next k */ + continue; + } + /* if |A[l,jk]| > 0.01 * v[l] for some l, continue to the + next k */ + for (t = 1; t <= len; t++) + { l = ind[t]; + if (fabs(val[t]) > 0.01 * v[l]) break; + } + if (t <= len) continue; + /* otherwise, let alpha = max{|A[l,jk]|: I[l] = 0} and let l' + be such that alpha = |A[l',jk]| */ + alpha = 0.0, ll = 0; + for (t = 1; t <= len; t++) + { l = ind[t]; + if (I[l] == 0 && alpha < fabs(val[t])) + alpha = fabs(val[t]), ll = l; + } + /* if alpha = 0, continue to the next k */ + if (alpha == 0.0) continue; + /* B := B union {jk} */ + glp_set_col_stat(lp, jk, GLP_BS); + I[ll] = 1; + v[ll] = alpha; + /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */ + for (t = 1; t <= len; t++) + { l = ind[t]; + if (val[t] != 0.0) r[l]++; + } + } + /*** STEP 3 ***/ + /* add an artificial variable (auxiliary variable for equality + constraint) to cover each remaining uncovered row */ + for (i = 1; i <= m; i++) + if (I[i] == 0) glp_set_row_stat(lp, i, GLP_BS); + /* free working arrays */ + xfree(C); + xfree(I); + xfree(r); + xfree(v); + xfree(ind); + xfree(val); + return; +} + +/*********************************************************************** +* NAME +* +* glp_cpx_basis - construct Bixby's initial LP basis +* +* SYNOPSIS +* +* void glp_cpx_basis(glp_prob *lp); +* +* DESCRIPTION +* +* The routine glp_cpx_basis constructs an advanced initial basis for +* the specified problem object. +* +* The routine is based on Bixby's algorithm described in the paper: +* +* Robert E. Bixby. Implementing the Simplex Method: The Initial Basis. +* ORSA Journal on Computing, Vol. 4, No. 3, 1992, pp. 267-84. */ + +void glp_cpx_basis(glp_prob *lp) +{ if (lp->m == 0 || lp->n == 0) + glp_std_basis(lp); + else + cpx_basis(lp); + return; +} + +/* eof */ |