diff -r d59bea55db9b -r c445c931472f src/glpapi06.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glpapi06.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,806 @@ +/* glpapi06.c (simplex method routines) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics, +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. +* E-mail: . +* +* 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 "glpios.h" +#include "glpnpp.h" +#include "glpspx.h" + +/*********************************************************************** +* NAME +* +* glp_simplex - solve LP problem with the simplex method +* +* SYNOPSIS +* +* int glp_simplex(glp_prob *P, const glp_smcp *parm); +* +* DESCRIPTION +* +* The routine glp_simplex is a driver to the LP solver based on the +* simplex method. This routine retrieves problem data from the +* specified problem object, calls the solver to solve the problem +* instance, and stores results of computations back into the problem +* object. +* +* The simplex solver has a set of control parameters. Values of the +* control parameters can be passed in a structure glp_smcp, which the +* parameter parm points to. +* +* The parameter parm can be specified as NULL, in which case the LP +* solver uses default settings. +* +* RETURNS +* +* 0 The LP problem instance has been successfully solved. This code +* does not necessarily mean that the solver has found optimal +* solution. It only means that the solution process was successful. +* +* GLP_EBADB +* Unable to start the search, because the initial basis specified +* in the problem object is invalid--the number of basic (auxiliary +* and structural) variables is not the same as the number of rows in +* the problem object. +* +* GLP_ESING +* Unable to start the search, because the basis matrix correspodning +* to the initial basis is singular within the working precision. +* +* GLP_ECOND +* Unable to start the search, because the basis matrix correspodning +* to the initial basis is ill-conditioned, i.e. its condition number +* is too large. +* +* GLP_EBOUND +* Unable to start the search, because some double-bounded variables +* have incorrect bounds. +* +* GLP_EFAIL +* The search was prematurely terminated due to the solver failure. +* +* GLP_EOBJLL +* The search was prematurely terminated, because the objective +* function being maximized has reached its lower limit and continues +* decreasing (dual simplex only). +* +* GLP_EOBJUL +* The search was prematurely terminated, because the objective +* function being minimized has reached its upper limit and continues +* increasing (dual simplex only). +* +* GLP_EITLIM +* The search was prematurely terminated, because the simplex +* iteration limit has been exceeded. +* +* GLP_ETMLIM +* The search was prematurely terminated, because the time limit has +* been exceeded. +* +* GLP_ENOPFS +* The LP problem instance has no primal feasible solution (only if +* the LP presolver is used). +* +* GLP_ENODFS +* The LP problem instance has no dual feasible solution (only if the +* LP presolver is used). */ + +static void trivial_lp(glp_prob *P, const glp_smcp *parm) +{ /* solve trivial LP which has empty constraint matrix */ + GLPROW *row; + GLPCOL *col; + int i, j; + double p_infeas, d_infeas, zeta; + P->valid = 0; + P->pbs_stat = P->dbs_stat = GLP_FEAS; + P->obj_val = P->c0; + P->some = 0; + p_infeas = d_infeas = 0.0; + /* make all auxiliary variables basic */ + for (i = 1; i <= P->m; i++) + { row = P->row[i]; + row->stat = GLP_BS; + row->prim = row->dual = 0.0; + /* check primal feasibility */ + if (row->type == GLP_LO || row->type == GLP_DB || + row->type == GLP_FX) + { /* row has lower bound */ + if (row->lb > + parm->tol_bnd) + { P->pbs_stat = GLP_NOFEAS; + if (P->some == 0 && parm->meth != GLP_PRIMAL) + P->some = i; + } + if (p_infeas < + row->lb) + p_infeas = + row->lb; + } + if (row->type == GLP_UP || row->type == GLP_DB || + row->type == GLP_FX) + { /* row has upper bound */ + if (row->ub < - parm->tol_bnd) + { P->pbs_stat = GLP_NOFEAS; + if (P->some == 0 && parm->meth != GLP_PRIMAL) + P->some = i; + } + if (p_infeas < - row->ub) + p_infeas = - row->ub; + } + } + /* determine scale factor for the objective row */ + zeta = 1.0; + for (j = 1; j <= P->n; j++) + { col = P->col[j]; + if (zeta < fabs(col->coef)) zeta = fabs(col->coef); + } + zeta = (P->dir == GLP_MIN ? +1.0 : -1.0) / zeta; + /* make all structural variables non-basic */ + for (j = 1; j <= P->n; j++) + { col = P->col[j]; + if (col->type == GLP_FR) + col->stat = GLP_NF, col->prim = 0.0; + else if (col->type == GLP_LO) +lo: col->stat = GLP_NL, col->prim = col->lb; + else if (col->type == GLP_UP) +up: col->stat = GLP_NU, col->prim = col->ub; + else if (col->type == GLP_DB) + { if (zeta * col->coef > 0.0) + goto lo; + else if (zeta * col->coef < 0.0) + goto up; + else if (fabs(col->lb) <= fabs(col->ub)) + goto lo; + else + goto up; + } + else if (col->type == GLP_FX) + col->stat = GLP_NS, col->prim = col->lb; + col->dual = col->coef; + P->obj_val += col->coef * col->prim; + /* check dual feasibility */ + if (col->type == GLP_FR || col->type == GLP_LO) + { /* column has no upper bound */ + if (zeta * col->dual < - parm->tol_dj) + { P->dbs_stat = GLP_NOFEAS; + if (P->some == 0 && parm->meth == GLP_PRIMAL) + P->some = P->m + j; + } + if (d_infeas < - zeta * col->dual) + d_infeas = - zeta * col->dual; + } + if (col->type == GLP_FR || col->type == GLP_UP) + { /* column has no lower bound */ + if (zeta * col->dual > + parm->tol_dj) + { P->dbs_stat = GLP_NOFEAS; + if (P->some == 0 && parm->meth == GLP_PRIMAL) + P->some = P->m + j; + } + if (d_infeas < + zeta * col->dual) + d_infeas = + zeta * col->dual; + } + } + /* simulate the simplex solver output */ + if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0) + { xprintf("~%6d: obj = %17.9e infeas = %10.3e\n", P->it_cnt, + P->obj_val, parm->meth == GLP_PRIMAL ? p_infeas : d_infeas); + } + if (parm->msg_lev >= GLP_MSG_ALL && parm->out_dly == 0) + { if (P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS) + xprintf("OPTIMAL SOLUTION FOUND\n"); + else if (P->pbs_stat == GLP_NOFEAS) + xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n"); + else if (parm->meth == GLP_PRIMAL) + xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n"); + else + xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n"); + } + return; +} + +static int solve_lp(glp_prob *P, const glp_smcp *parm) +{ /* solve LP directly without using the preprocessor */ + int ret; + if (!glp_bf_exists(P)) + { ret = glp_factorize(P); + if (ret == 0) + ; + else if (ret == GLP_EBADB) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf("glp_simplex: initial basis is invalid\n"); + } + else if (ret == GLP_ESING) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf("glp_simplex: initial basis is singular\n"); + } + else if (ret == GLP_ECOND) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf( + "glp_simplex: initial basis is ill-conditioned\n"); + } + else + xassert(ret != ret); + if (ret != 0) goto done; + } + if (parm->meth == GLP_PRIMAL) + ret = spx_primal(P, parm); + else if (parm->meth == GLP_DUALP) + { ret = spx_dual(P, parm); + if (ret == GLP_EFAIL && P->valid) + ret = spx_primal(P, parm); + } + else if (parm->meth == GLP_DUAL) + ret = spx_dual(P, parm); + else + xassert(parm != parm); +done: return ret; +} + +static int preprocess_and_solve_lp(glp_prob *P, const glp_smcp *parm) +{ /* solve LP using the preprocessor */ + NPP *npp; + glp_prob *lp = NULL; + glp_bfcp bfcp; + int ret; + if (parm->msg_lev >= GLP_MSG_ALL) + xprintf("Preprocessing...\n"); + /* create preprocessor workspace */ + npp = npp_create_wksp(); + /* load original problem into the preprocessor workspace */ + npp_load_prob(npp, P, GLP_OFF, GLP_SOL, GLP_OFF); + /* process LP prior to applying primal/dual simplex method */ + ret = npp_simplex(npp, parm); + if (ret == 0) + ; + else if (ret == GLP_ENOPFS) + { if (parm->msg_lev >= GLP_MSG_ALL) + xprintf("PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION\n"); + } + else if (ret == GLP_ENODFS) + { if (parm->msg_lev >= GLP_MSG_ALL) + xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n"); + } + else + xassert(ret != ret); + if (ret != 0) goto done; + /* build transformed LP */ + lp = glp_create_prob(); + npp_build_prob(npp, lp); + /* if the transformed LP is empty, it has empty solution, which + is optimal */ + if (lp->m == 0 && lp->n == 0) + { lp->pbs_stat = lp->dbs_stat = GLP_FEAS; + lp->obj_val = lp->c0; + if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0) + { xprintf("~%6d: obj = %17.9e infeas = %10.3e\n", P->it_cnt, + lp->obj_val, 0.0); + } + if (parm->msg_lev >= GLP_MSG_ALL) + xprintf("OPTIMAL SOLUTION FOUND BY LP PREPROCESSOR\n"); + goto post; + } + if (parm->msg_lev >= GLP_MSG_ALL) + { xprintf("%d row%s, %d column%s, %d non-zero%s\n", + lp->m, lp->m == 1 ? "" : "s", lp->n, lp->n == 1 ? "" : "s", + lp->nnz, lp->nnz == 1 ? "" : "s"); + } + /* inherit basis factorization control parameters */ + glp_get_bfcp(P, &bfcp); + glp_set_bfcp(lp, &bfcp); + /* scale the transformed problem */ + { ENV *env = get_env_ptr(); + int term_out = env->term_out; + if (!term_out || parm->msg_lev < GLP_MSG_ALL) + env->term_out = GLP_OFF; + else + env->term_out = GLP_ON; + glp_scale_prob(lp, GLP_SF_AUTO); + env->term_out = term_out; + } + /* build advanced initial basis */ + { ENV *env = get_env_ptr(); + int term_out = env->term_out; + if (!term_out || parm->msg_lev < GLP_MSG_ALL) + env->term_out = GLP_OFF; + else + env->term_out = GLP_ON; + glp_adv_basis(lp, 0); + env->term_out = term_out; + } + /* solve the transformed LP */ + lp->it_cnt = P->it_cnt; + ret = solve_lp(lp, parm); + P->it_cnt = lp->it_cnt; + /* only optimal solution can be postprocessed */ + if (!(ret == 0 && lp->pbs_stat == GLP_FEAS && lp->dbs_stat == + GLP_FEAS)) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf("glp_simplex: unable to recover undefined or non-op" + "timal solution\n"); + if (ret == 0) + { if (lp->pbs_stat == GLP_NOFEAS) + ret = GLP_ENOPFS; + else if (lp->dbs_stat == GLP_NOFEAS) + ret = GLP_ENODFS; + else + xassert(lp != lp); + } + goto done; + } +post: /* postprocess solution from the transformed LP */ + npp_postprocess(npp, lp); + /* the transformed LP is no longer needed */ + glp_delete_prob(lp), lp = NULL; + /* store solution to the original problem */ + npp_unload_sol(npp, P); + /* the original LP has been successfully solved */ + ret = 0; +done: /* delete the transformed LP, if it exists */ + if (lp != NULL) glp_delete_prob(lp); + /* delete preprocessor workspace */ + npp_delete_wksp(npp); + return ret; +} + +int glp_simplex(glp_prob *P, const glp_smcp *parm) +{ /* solve LP problem with the simplex method */ + glp_smcp _parm; + int i, j, ret; + /* check problem object */ + if (P == NULL || P->magic != GLP_PROB_MAGIC) + xerror("glp_simplex: P = %p; invalid problem object\n", P); + if (P->tree != NULL && P->tree->reason != 0) + xerror("glp_simplex: operation not allowed\n"); + /* check control parameters */ + if (parm == NULL) + parm = &_parm, glp_init_smcp((glp_smcp *)parm); + if (!(parm->msg_lev == GLP_MSG_OFF || + parm->msg_lev == GLP_MSG_ERR || + parm->msg_lev == GLP_MSG_ON || + parm->msg_lev == GLP_MSG_ALL || + parm->msg_lev == GLP_MSG_DBG)) + xerror("glp_simplex: msg_lev = %d; invalid parameter\n", + parm->msg_lev); + if (!(parm->meth == GLP_PRIMAL || + parm->meth == GLP_DUALP || + parm->meth == GLP_DUAL)) + xerror("glp_simplex: meth = %d; invalid parameter\n", + parm->meth); + if (!(parm->pricing == GLP_PT_STD || + parm->pricing == GLP_PT_PSE)) + xerror("glp_simplex: pricing = %d; invalid parameter\n", + parm->pricing); + if (!(parm->r_test == GLP_RT_STD || + parm->r_test == GLP_RT_HAR)) + xerror("glp_simplex: r_test = %d; invalid parameter\n", + parm->r_test); + if (!(0.0 < parm->tol_bnd && parm->tol_bnd < 1.0)) + xerror("glp_simplex: tol_bnd = %g; invalid parameter\n", + parm->tol_bnd); + if (!(0.0 < parm->tol_dj && parm->tol_dj < 1.0)) + xerror("glp_simplex: tol_dj = %g; invalid parameter\n", + parm->tol_dj); + if (!(0.0 < parm->tol_piv && parm->tol_piv < 1.0)) + xerror("glp_simplex: tol_piv = %g; invalid parameter\n", + parm->tol_piv); + if (parm->it_lim < 0) + xerror("glp_simplex: it_lim = %d; invalid parameter\n", + parm->it_lim); + if (parm->tm_lim < 0) + xerror("glp_simplex: tm_lim = %d; invalid parameter\n", + parm->tm_lim); + if (parm->out_frq < 1) + xerror("glp_simplex: out_frq = %d; invalid parameter\n", + parm->out_frq); + if (parm->out_dly < 0) + xerror("glp_simplex: out_dly = %d; invalid parameter\n", + parm->out_dly); + if (!(parm->presolve == GLP_ON || parm->presolve == GLP_OFF)) + xerror("glp_simplex: presolve = %d; invalid parameter\n", + parm->presolve); + /* basic solution is currently undefined */ + P->pbs_stat = P->dbs_stat = GLP_UNDEF; + P->obj_val = 0.0; + P->some = 0; + /* check bounds of double-bounded variables */ + for (i = 1; i <= P->m; i++) + { GLPROW *row = P->row[i]; + if (row->type == GLP_DB && row->lb >= row->ub) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf("glp_simplex: row %d: lb = %g, ub = %g; incorrec" + "t bounds\n", i, row->lb, row->ub); + ret = GLP_EBOUND; + goto done; + } + } + for (j = 1; j <= P->n; j++) + { GLPCOL *col = P->col[j]; + if (col->type == GLP_DB && col->lb >= col->ub) + { if (parm->msg_lev >= GLP_MSG_ERR) + xprintf("glp_simplex: column %d: lb = %g, ub = %g; incor" + "rect bounds\n", j, col->lb, col->ub); + ret = GLP_EBOUND; + goto done; + } + } + /* solve LP problem */ + if (parm->msg_lev >= GLP_MSG_ALL) + { xprintf("GLPK Simplex Optimizer, v%s\n", glp_version()); + xprintf("%d row%s, %d column%s, %d non-zero%s\n", + P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" : "s", + P->nnz, P->nnz == 1 ? "" : "s"); + } + if (P->nnz == 0) + trivial_lp(P, parm), ret = 0; + else if (!parm->presolve) + ret = solve_lp(P, parm); + else + ret = preprocess_and_solve_lp(P, parm); +done: /* return to the application program */ + return ret; +} + +/*********************************************************************** +* NAME +* +* glp_init_smcp - initialize simplex method control parameters +* +* SYNOPSIS +* +* void glp_init_smcp(glp_smcp *parm); +* +* DESCRIPTION +* +* The routine glp_init_smcp initializes control parameters, which are +* used by the simplex solver, with default values. +* +* Default values of the control parameters are stored in a glp_smcp +* structure, which the parameter parm points to. */ + +void glp_init_smcp(glp_smcp *parm) +{ parm->msg_lev = GLP_MSG_ALL; + parm->meth = GLP_PRIMAL; + parm->pricing = GLP_PT_PSE; + parm->r_test = GLP_RT_HAR; + parm->tol_bnd = 1e-7; + parm->tol_dj = 1e-7; + parm->tol_piv = 1e-10; + parm->obj_ll = -DBL_MAX; + parm->obj_ul = +DBL_MAX; + parm->it_lim = INT_MAX; + parm->tm_lim = INT_MAX; + parm->out_frq = 500; + parm->out_dly = 0; + parm->presolve = GLP_OFF; + return; +} + +/*********************************************************************** +* NAME +* +* glp_get_status - retrieve generic status of basic solution +* +* SYNOPSIS +* +* int glp_get_status(glp_prob *lp); +* +* RETURNS +* +* The routine glp_get_status reports the generic status of the basic +* solution for the specified problem object as follows: +* +* GLP_OPT - solution is optimal; +* GLP_FEAS - solution is feasible; +* GLP_INFEAS - solution is infeasible; +* GLP_NOFEAS - problem has no feasible solution; +* GLP_UNBND - problem has unbounded solution; +* GLP_UNDEF - solution is undefined. */ + +int glp_get_status(glp_prob *lp) +{ int status; + status = glp_get_prim_stat(lp); + switch (status) + { case GLP_FEAS: + switch (glp_get_dual_stat(lp)) + { case GLP_FEAS: + status = GLP_OPT; + break; + case GLP_NOFEAS: + status = GLP_UNBND; + break; + case GLP_UNDEF: + case GLP_INFEAS: + status = status; + break; + default: + xassert(lp != lp); + } + break; + case GLP_UNDEF: + case GLP_INFEAS: + case GLP_NOFEAS: + status = status; + break; + default: + xassert(lp != lp); + } + return status; +} + +/*********************************************************************** +* NAME +* +* glp_get_prim_stat - retrieve status of primal basic solution +* +* SYNOPSIS +* +* int glp_get_prim_stat(glp_prob *lp); +* +* RETURNS +* +* The routine glp_get_prim_stat reports the status of the primal basic +* solution for the specified problem object as follows: +* +* GLP_UNDEF - primal solution is undefined; +* GLP_FEAS - primal solution is feasible; +* GLP_INFEAS - primal solution is infeasible; +* GLP_NOFEAS - no primal feasible solution exists. */ + +int glp_get_prim_stat(glp_prob *lp) +{ int pbs_stat = lp->pbs_stat; + return pbs_stat; +} + +/*********************************************************************** +* NAME +* +* glp_get_dual_stat - retrieve status of dual basic solution +* +* SYNOPSIS +* +* int glp_get_dual_stat(glp_prob *lp); +* +* RETURNS +* +* The routine glp_get_dual_stat reports the status of the dual basic +* solution for the specified problem object as follows: +* +* GLP_UNDEF - dual solution is undefined; +* GLP_FEAS - dual solution is feasible; +* GLP_INFEAS - dual solution is infeasible; +* GLP_NOFEAS - no dual feasible solution exists. */ + +int glp_get_dual_stat(glp_prob *lp) +{ int dbs_stat = lp->dbs_stat; + return dbs_stat; +} + +/*********************************************************************** +* NAME +* +* glp_get_obj_val - retrieve objective value (basic solution) +* +* SYNOPSIS +* +* double glp_get_obj_val(glp_prob *lp); +* +* RETURNS +* +* The routine glp_get_obj_val returns value of the objective function +* for basic solution. */ + +double glp_get_obj_val(glp_prob *lp) +{ /*struct LPXCPS *cps = lp->cps;*/ + double z; + z = lp->obj_val; + /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/ + return z; +} + +/*********************************************************************** +* NAME +* +* glp_get_row_stat - retrieve row status +* +* SYNOPSIS +* +* int glp_get_row_stat(glp_prob *lp, int i); +* +* RETURNS +* +* The routine glp_get_row_stat returns current status assigned to the +* auxiliary variable associated with i-th row as follows: +* +* GLP_BS - basic variable; +* GLP_NL - non-basic variable on its lower bound; +* GLP_NU - non-basic variable on its upper bound; +* GLP_NF - non-basic free (unbounded) variable; +* GLP_NS - non-basic fixed variable. */ + +int glp_get_row_stat(glp_prob *lp, int i) +{ if (!(1 <= i && i <= lp->m)) + xerror("glp_get_row_stat: i = %d; row number out of range\n", + i); + return lp->row[i]->stat; +} + +/*********************************************************************** +* NAME +* +* glp_get_row_prim - retrieve row primal value (basic solution) +* +* SYNOPSIS +* +* double glp_get_row_prim(glp_prob *lp, int i); +* +* RETURNS +* +* The routine glp_get_row_prim returns primal value of the auxiliary +* variable associated with i-th row. */ + +double glp_get_row_prim(glp_prob *lp, int i) +{ /*struct LPXCPS *cps = lp->cps;*/ + double prim; + if (!(1 <= i && i <= lp->m)) + xerror("glp_get_row_prim: i = %d; row number out of range\n", + i); + prim = lp->row[i]->prim; + /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/ + return prim; +} + +/*********************************************************************** +* NAME +* +* glp_get_row_dual - retrieve row dual value (basic solution) +* +* SYNOPSIS +* +* double glp_get_row_dual(glp_prob *lp, int i); +* +* RETURNS +* +* The routine glp_get_row_dual returns dual value (i.e. reduced cost) +* of the auxiliary variable associated with i-th row. */ + +double glp_get_row_dual(glp_prob *lp, int i) +{ /*struct LPXCPS *cps = lp->cps;*/ + double dual; + if (!(1 <= i && i <= lp->m)) + xerror("glp_get_row_dual: i = %d; row number out of range\n", + i); + dual = lp->row[i]->dual; + /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/ + return dual; +} + +/*********************************************************************** +* NAME +* +* glp_get_col_stat - retrieve column status +* +* SYNOPSIS +* +* int glp_get_col_stat(glp_prob *lp, int j); +* +* RETURNS +* +* The routine glp_get_col_stat returns current status assigned to the +* structural variable associated with j-th column as follows: +* +* GLP_BS - basic variable; +* GLP_NL - non-basic variable on its lower bound; +* GLP_NU - non-basic variable on its upper bound; +* GLP_NF - non-basic free (unbounded) variable; +* GLP_NS - non-basic fixed variable. */ + +int glp_get_col_stat(glp_prob *lp, int j) +{ if (!(1 <= j && j <= lp->n)) + xerror("glp_get_col_stat: j = %d; column number out of range\n" + , j); + return lp->col[j]->stat; +} + +/*********************************************************************** +* NAME +* +* glp_get_col_prim - retrieve column primal value (basic solution) +* +* SYNOPSIS +* +* double glp_get_col_prim(glp_prob *lp, int j); +* +* RETURNS +* +* The routine glp_get_col_prim returns primal value of the structural +* variable associated with j-th column. */ + +double glp_get_col_prim(glp_prob *lp, int j) +{ /*struct LPXCPS *cps = lp->cps;*/ + double prim; + if (!(1 <= j && j <= lp->n)) + xerror("glp_get_col_prim: j = %d; column number out of range\n" + , j); + prim = lp->col[j]->prim; + /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/ + return prim; +} + +/*********************************************************************** +* NAME +* +* glp_get_col_dual - retrieve column dual value (basic solution) +* +* SYNOPSIS +* +* double glp_get_col_dual(glp_prob *lp, int j); +* +* RETURNS +* +* The routine glp_get_col_dual returns dual value (i.e. reduced cost) +* of the structural variable associated with j-th column. */ + +double glp_get_col_dual(glp_prob *lp, int j) +{ /*struct LPXCPS *cps = lp->cps;*/ + double dual; + if (!(1 <= j && j <= lp->n)) + xerror("glp_get_col_dual: j = %d; column number out of range\n" + , j); + dual = lp->col[j]->dual; + /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/ + return dual; +} + +/*********************************************************************** +* NAME +* +* glp_get_unbnd_ray - determine variable causing unboundedness +* +* SYNOPSIS +* +* int glp_get_unbnd_ray(glp_prob *lp); +* +* RETURNS +* +* The routine glp_get_unbnd_ray returns the number k of a variable, +* which causes primal or dual unboundedness. If 1 <= k <= m, it is +* k-th auxiliary variable, and if m+1 <= k <= m+n, it is (k-m)-th +* structural variable, where m is the number of rows, n is the number +* of columns in the problem object. If such variable is not defined, +* the routine returns 0. +* +* COMMENTS +* +* If it is not exactly known which version of the simplex solver +* detected unboundedness, i.e. whether the unboundedness is primal or +* dual, it is sufficient to check the status of the variable reported +* with the routine glp_get_row_stat or glp_get_col_stat. If the +* variable is non-basic, the unboundedness is primal, otherwise, if +* the variable is basic, the unboundedness is dual (the latter case +* means that the problem has no primal feasible dolution). */ + +int glp_get_unbnd_ray(glp_prob *lp) +{ int k; + k = lp->some; + xassert(k >= 0); + if (k > lp->m + lp->n) k = 0; + return k; +} + +/* eof */