/* glpssx02.c (simplex method, rational arithmetic) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * Copyright (C) 2003-2013 Free Software Foundation, Inc. * Written by Andrew Makhorin . * * 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 . ***********************************************************************/ #include "env.h" #include "glpssx.h" static void show_progress(SSX *ssx, int phase) { /* this auxiliary routine displays information about progress of the search */ int i, def = 0; for (i = 1; i <= ssx->m; i++) if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++; xprintf("%s%6d: %s = %22.15g (%d)\n", phase == 1 ? " " : "*", ssx->it_cnt, phase == 1 ? "infsum" : "objval", mpq_get_d(ssx->bbar[0]), def); #if 0 ssx->tm_lag = utime(); #else ssx->tm_lag = xtime(); #endif return; } /*---------------------------------------------------------------------- // ssx_phase_I - find primal feasible solution. // // This routine implements phase I of the primal simplex method. // // On exit the routine returns one of the following codes: // // 0 - feasible solution found; // 1 - problem has no feasible solution; // 2 - iterations limit exceeded; // 3 - time limit exceeded. ----------------------------------------------------------------------*/ int ssx_phase_I(SSX *ssx) { int m = ssx->m; int n = ssx->n; int *type = ssx->type; mpq_t *lb = ssx->lb; mpq_t *ub = ssx->ub; mpq_t *coef = ssx->coef; int *A_ptr = ssx->A_ptr; int *A_ind = ssx->A_ind; mpq_t *A_val = ssx->A_val; int *Q_col = ssx->Q_col; mpq_t *bbar = ssx->bbar; mpq_t *pi = ssx->pi; mpq_t *cbar = ssx->cbar; int *orig_type, orig_dir; mpq_t *orig_lb, *orig_ub, *orig_coef; int i, k, ret; /* save components of the original LP problem, which are changed by the routine */ orig_type = xcalloc(1+m+n, sizeof(int)); orig_lb = xcalloc(1+m+n, sizeof(mpq_t)); orig_ub = xcalloc(1+m+n, sizeof(mpq_t)); orig_coef = xcalloc(1+m+n, sizeof(mpq_t)); for (k = 1; k <= m+n; k++) { orig_type[k] = type[k]; mpq_init(orig_lb[k]); mpq_set(orig_lb[k], lb[k]); mpq_init(orig_ub[k]); mpq_set(orig_ub[k], ub[k]); } orig_dir = ssx->dir; for (k = 0; k <= m+n; k++) { mpq_init(orig_coef[k]); mpq_set(orig_coef[k], coef[k]); } /* build an artificial basic solution, which is primal feasible, and also build an auxiliary objective function to minimize the sum of infeasibilities for the original problem */ ssx->dir = SSX_MIN; for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1); mpq_set_si(bbar[0], 0, 1); for (i = 1; i <= m; i++) { int t; k = Q_col[i]; /* x[k] = xB[i] */ t = type[k]; if (t == SSX_LO || t == SSX_DB || t == SSX_FX) { /* in the original problem x[k] has lower bound */ if (mpq_cmp(bbar[i], lb[k]) < 0) { /* which is violated */ type[k] = SSX_UP; mpq_set(ub[k], lb[k]); mpq_set_si(lb[k], 0, 1); mpq_set_si(coef[k], -1, 1); mpq_add(bbar[0], bbar[0], ub[k]); mpq_sub(bbar[0], bbar[0], bbar[i]); } } if (t == SSX_UP || t == SSX_DB || t == SSX_FX) { /* in the original problem x[k] has upper bound */ if (mpq_cmp(bbar[i], ub[k]) > 0) { /* which is violated */ type[k] = SSX_LO; mpq_set(lb[k], ub[k]); mpq_set_si(ub[k], 0, 1); mpq_set_si(coef[k], +1, 1); mpq_add(bbar[0], bbar[0], bbar[i]); mpq_sub(bbar[0], bbar[0], lb[k]); } } } /* now the initial basic solution should be primal feasible due to changes of bounds of some basic variables, which turned to implicit artifical variables */ /* compute simplex multipliers and reduced costs */ ssx_eval_pi(ssx); ssx_eval_cbar(ssx); /* display initial progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif show_progress(ssx, 1); /* main loop starts here */ for (;;) { /* display current progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif #if 0 if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001) #else if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001) #endif show_progress(ssx, 1); /* we do not need to wait until all artificial variables have left the basis */ if (mpq_sgn(bbar[0]) == 0) { /* the sum of infeasibilities is zero, therefore the current solution is primal feasible for the original problem */ ret = 0; break; } /* check if the iterations limit has been exhausted */ if (ssx->it_lim == 0) { ret = 2; break; } /* check if the time limit has been exhausted */ #if 0 if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg) #else if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg)) #endif { ret = 3; break; } /* choose non-basic variable xN[q] */ ssx_chuzc(ssx); /* if xN[q] cannot be chosen, the sum of infeasibilities is minimal but non-zero; therefore the original problem has no primal feasible solution */ if (ssx->q == 0) { ret = 1; break; } /* compute q-th column of the simplex table */ ssx_eval_col(ssx); /* choose basic variable xB[p] */ ssx_chuzr(ssx); /* the sum of infeasibilities cannot be negative, therefore the auxiliary lp problem cannot have unbounded solution */ xassert(ssx->p != 0); /* update values of basic variables */ ssx_update_bbar(ssx); if (ssx->p > 0) { /* compute p-th row of the inverse inv(B) */ ssx_eval_rho(ssx); /* compute p-th row of the simplex table */ ssx_eval_row(ssx); xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0); /* update simplex multipliers */ ssx_update_pi(ssx); /* update reduced costs of non-basic variables */ ssx_update_cbar(ssx); } /* xB[p] is leaving the basis; if it is implicit artificial variable, the corresponding residual vanishes; therefore bounds of this variable should be restored to the original values */ if (ssx->p > 0) { k = Q_col[ssx->p]; /* x[k] = xB[p] */ if (type[k] != orig_type[k]) { /* x[k] is implicit artificial variable */ type[k] = orig_type[k]; mpq_set(lb[k], orig_lb[k]); mpq_set(ub[k], orig_ub[k]); xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU); ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL); if (type[k] == SSX_FX) ssx->p_stat = SSX_NS; /* nullify the objective coefficient at x[k] */ mpq_set_si(coef[k], 0, 1); /* since coef[k] has been changed, we need to compute new reduced cost of x[k], which it will have in the adjacent basis */ /* the formula d[j] = cN[j] - pi' * N[j] is used (note that the vector pi is not changed, because it depends on objective coefficients at basic variables, but in the adjacent basis, for which the vector pi has been just recomputed, x[k] is non-basic) */ if (k <= m) { /* x[k] is auxiliary variable */ mpq_neg(cbar[ssx->q], pi[k]); } else { /* x[k] is structural variable */ int ptr; mpq_t temp; mpq_init(temp); mpq_set_si(cbar[ssx->q], 0, 1); for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++) { mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]); mpq_add(cbar[ssx->q], cbar[ssx->q], temp); } mpq_clear(temp); } } } /* jump to the adjacent vertex of the polyhedron */ ssx_change_basis(ssx); /* one simplex iteration has been performed */ if (ssx->it_lim > 0) ssx->it_lim--; ssx->it_cnt++; } /* display final progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif show_progress(ssx, 1); /* restore components of the original problem, which were changed by the routine */ for (k = 1; k <= m+n; k++) { type[k] = orig_type[k]; mpq_set(lb[k], orig_lb[k]); mpq_clear(orig_lb[k]); mpq_set(ub[k], orig_ub[k]); mpq_clear(orig_ub[k]); } ssx->dir = orig_dir; for (k = 0; k <= m+n; k++) { mpq_set(coef[k], orig_coef[k]); mpq_clear(orig_coef[k]); } xfree(orig_type); xfree(orig_lb); xfree(orig_ub); xfree(orig_coef); /* return to the calling program */ return ret; } /*---------------------------------------------------------------------- // ssx_phase_II - find optimal solution. // // This routine implements phase II of the primal simplex method. // // On exit the routine returns one of the following codes: // // 0 - optimal solution found; // 1 - problem has unbounded solution; // 2 - iterations limit exceeded; // 3 - time limit exceeded. ----------------------------------------------------------------------*/ int ssx_phase_II(SSX *ssx) { int ret; /* display initial progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif show_progress(ssx, 2); /* main loop starts here */ for (;;) { /* display current progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif #if 0 if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001) #else if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001) #endif show_progress(ssx, 2); /* check if the iterations limit has been exhausted */ if (ssx->it_lim == 0) { ret = 2; break; } /* check if the time limit has been exhausted */ #if 0 if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg) #else if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg)) #endif { ret = 3; break; } /* choose non-basic variable xN[q] */ ssx_chuzc(ssx); /* if xN[q] cannot be chosen, the current basic solution is dual feasible and therefore optimal */ if (ssx->q == 0) { ret = 0; break; } /* compute q-th column of the simplex table */ ssx_eval_col(ssx); /* choose basic variable xB[p] */ ssx_chuzr(ssx); /* if xB[p] cannot be chosen, the problem has no dual feasible solution (i.e. unbounded) */ if (ssx->p == 0) { ret = 1; break; } /* update values of basic variables */ ssx_update_bbar(ssx); if (ssx->p > 0) { /* compute p-th row of the inverse inv(B) */ ssx_eval_rho(ssx); /* compute p-th row of the simplex table */ ssx_eval_row(ssx); xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0); #if 0 /* update simplex multipliers */ ssx_update_pi(ssx); #endif /* update reduced costs of non-basic variables */ ssx_update_cbar(ssx); } /* jump to the adjacent vertex of the polyhedron */ ssx_change_basis(ssx); /* one simplex iteration has been performed */ if (ssx->it_lim > 0) ssx->it_lim--; ssx->it_cnt++; } /* display final progress of the search */ #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ON) #endif show_progress(ssx, 2); /* return to the calling program */ return ret; } /*---------------------------------------------------------------------- // ssx_driver - base driver to exact simplex method. // // This routine is a base driver to a version of the primal simplex // method using exact (bignum) arithmetic. // // On exit the routine returns one of the following codes: // // 0 - optimal solution found; // 1 - problem has no feasible solution; // 2 - problem has unbounded solution; // 3 - iterations limit exceeded (phase I); // 4 - iterations limit exceeded (phase II); // 5 - time limit exceeded (phase I); // 6 - time limit exceeded (phase II); // 7 - initial basis matrix is exactly singular. ----------------------------------------------------------------------*/ int ssx_driver(SSX *ssx) { int m = ssx->m; int *type = ssx->type; mpq_t *lb = ssx->lb; mpq_t *ub = ssx->ub; int *Q_col = ssx->Q_col; mpq_t *bbar = ssx->bbar; int i, k, ret; ssx->tm_beg = xtime(); /* factorize the initial basis matrix */ if (ssx_factorize(ssx)) #if 0 /* 25/XI-2017 */ { xprintf("Initial basis matrix is singular\n"); #else { if (ssx->msg_lev >= GLP_MSG_ERR) xprintf("Initial basis matrix is singular\n"); #endif ret = 7; goto done; } /* compute values of basic variables */ ssx_eval_bbar(ssx); /* check if the initial basic solution is primal feasible */ for (i = 1; i <= m; i++) { int t; k = Q_col[i]; /* x[k] = xB[i] */ t = type[k]; if (t == SSX_LO || t == SSX_DB || t == SSX_FX) { /* x[k] has lower bound */ if (mpq_cmp(bbar[i], lb[k]) < 0) { /* which is violated */ break; } } if (t == SSX_UP || t == SSX_DB || t == SSX_FX) { /* x[k] has upper bound */ if (mpq_cmp(bbar[i], ub[k]) > 0) { /* which is violated */ break; } } } if (i > m) { /* no basic variable violates its bounds */ ret = 0; goto skip; } /* phase I: find primal feasible solution */ ret = ssx_phase_I(ssx); switch (ret) { case 0: ret = 0; break; case 1: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n"); ret = 1; break; case 2: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n"); ret = 3; break; case 3: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n"); ret = 5; break; default: xassert(ret != ret); } /* compute values of basic variables (actually only the objective value needs to be computed) */ ssx_eval_bbar(ssx); skip: /* compute simplex multipliers */ ssx_eval_pi(ssx); /* compute reduced costs of non-basic variables */ ssx_eval_cbar(ssx); /* if phase I failed, do not start phase II */ if (ret != 0) goto done; /* phase II: find optimal solution */ ret = ssx_phase_II(ssx); switch (ret) { case 0: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("OPTIMAL SOLUTION FOUND\n"); ret = 0; break; case 1: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n"); ret = 2; break; case 2: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n"); ret = 4; break; case 3: #if 1 /* 25/XI-2017 */ if (ssx->msg_lev >= GLP_MSG_ALL) #endif xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n"); ret = 6; break; default: xassert(ret != ret); } done: /* decrease the time limit by the spent amount of time */ if (ssx->tm_lim >= 0.0) #if 0 { ssx->tm_lim -= utime() - ssx->tm_beg; #else { ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg); #endif if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0; } return ret; } /* eof */