lemon-project-template-glpk: deps/glpk/src/glpnpp03.c annotate

lemon-project-template-glpk

annotate deps/glpk/src/glpnpp03.c @ 11:4fc6ad2fb8a6

Test GLPK in src/main.cc

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 21:43:29 +0100
parents
children

rev	line source
alpar@9	1 /* glpnpp03.c */
alpar@9	2
alpar@9	3 /***********************************************************************
alpar@9	4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9	5 *
alpar@9	6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9	7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9	8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9	9 * E-mail: <mao@gnu.org>.
alpar@9	10 *
alpar@9	11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9	12 * under the terms of the GNU General Public License as published by
alpar@9	13 * the Free Software Foundation, either version 3 of the License, or
alpar@9	14 * (at your option) any later version.
alpar@9	15 *
alpar@9	16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9	17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9	18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9	19 * License for more details.
alpar@9	20 *
alpar@9	21 * You should have received a copy of the GNU General Public License
alpar@9	22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9	23 ***********************************************************************/
alpar@9	24
alpar@9	25 #include "glpnpp.h"
alpar@9	26
alpar@9	27 /***********************************************************************
alpar@9	28 * NAME
alpar@9	29 *
alpar@9	30 * npp_empty_row - process empty row
alpar@9	31 *
alpar@9	32 * SYNOPSIS
alpar@9	33 *
alpar@9	34 * #include "glpnpp.h"
alpar@9	35 * int npp_empty_row(NPP npp, NPPROW p);
alpar@9	36 *
alpar@9	37 * DESCRIPTION
alpar@9	38 *
alpar@9	39 * The routine npp_empty_row processes row p, which is empty, i.e.
alpar@9	40 * coefficients at all columns in this row are zero:
alpar@9	41 *
alpar@9	42 * L[p] <= sum 0 x[j] <= U[p], (1)
alpar@9	43 *
alpar@9	44 * where L[p] <= U[p].
alpar@9	45 *
alpar@9	46 * RETURNS
alpar@9	47 *
alpar@9	48 * 0 - success;
alpar@9	49 *
alpar@9	50 * 1 - problem has no primal feasible solution.
alpar@9	51 *
alpar@9	52 * PROBLEM TRANSFORMATION
alpar@9	53 *
alpar@9	54 * If the following conditions hold:
alpar@9	55 *
alpar@9	56 * L[p] <= +eps, U[p] >= -eps, (2)
alpar@9	57 *
alpar@9	58 * where eps is an absolute tolerance for row value, the row p is
alpar@9	59 * redundant. In this case it can be replaced by equivalent redundant
alpar@9	60 * row, which is free (unbounded), and then removed from the problem.
alpar@9	61 * Otherwise, the row p is infeasible and, thus, the problem has no
alpar@9	62 * primal feasible solution.
alpar@9	63 *
alpar@9	64 * RECOVERING BASIC SOLUTION
alpar@9	65 *
alpar@9	66 * See the routine npp_free_row.
alpar@9	67 *
alpar@9	68 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	69 *
alpar@9	70 * See the routine npp_free_row.
alpar@9	71 *
alpar@9	72 * RECOVERING MIP SOLUTION
alpar@9	73 *
alpar@9	74 * None needed. */
alpar@9	75
alpar@9	76 int npp_empty_row(NPP npp, NPPROW p)
alpar@9	77 { /* process empty row */
alpar@9	78 double eps = 1e-3;
alpar@9	79 /* the row must be empty */
alpar@9	80 xassert(p->ptr == NULL);
alpar@9	81 /* check primal feasibility */
alpar@9	82 if (p->lb > +eps \|\| p->ub < -eps)
alpar@9	83 return 1;
alpar@9	84 /* replace the row by equivalent free (unbounded) row */
alpar@9	85 p->lb = -DBL_MAX, p->ub = +DBL_MAX;
alpar@9	86 /* and process it */
alpar@9	87 npp_free_row(npp, p);
alpar@9	88 return 0;
alpar@9	89 }
alpar@9	90
alpar@9	91 /***********************************************************************
alpar@9	92 * NAME
alpar@9	93 *
alpar@9	94 * npp_empty_col - process empty column
alpar@9	95 *
alpar@9	96 * SYNOPSIS
alpar@9	97 *
alpar@9	98 * #include "glpnpp.h"
alpar@9	99 * int npp_empty_col(NPP npp, NPPCOL q);
alpar@9	100 *
alpar@9	101 * DESCRIPTION
alpar@9	102 *
alpar@9	103 * The routine npp_empty_col processes column q:
alpar@9	104 *
alpar@9	105 * l[q] <= x[q] <= u[q], (1)
alpar@9	106 *
alpar@9	107 * where l[q] <= u[q], which is empty, i.e. has zero coefficients in
alpar@9	108 * all constraint rows.
alpar@9	109 *
alpar@9	110 * RETURNS
alpar@9	111 *
alpar@9	112 * 0 - success;
alpar@9	113 *
alpar@9	114 * 1 - problem has no dual feasible solution.
alpar@9	115 *
alpar@9	116 * PROBLEM TRANSFORMATION
alpar@9	117 *
alpar@9	118 * The row of the dual system corresponding to the empty column is the
alpar@9	119 * following:
alpar@9	120 *
alpar@9	121 * sum 0 pi[i] + lambda[q] = c[q], (2)
alpar@9	122 * i
alpar@9	123 *
alpar@9	124 * from which it follows that:
alpar@9	125 *
alpar@9	126 * lambda[q] = c[q]. (3)
alpar@9	127 *
alpar@9	128 * If the following condition holds:
alpar@9	129 *
alpar@9	130 * c[q] < - eps, (4)
alpar@9	131 *
alpar@9	132 * where eps is an absolute tolerance for column multiplier, the lower
alpar@9	133 * column bound l[q] must be active to provide dual feasibility (note
alpar@9	134 * that being preprocessed the problem is always minimization). In this
alpar@9	135 * case the column can be fixed on its lower bound and removed from the
alpar@9	136 * problem (if the column is integral, its bounds are also assumed to
alpar@9	137 * be integral). And if the column has no lower bound (l[q] = -oo), the
alpar@9	138 * problem has no dual feasible solution.
alpar@9	139 *
alpar@9	140 * If the following condition holds:
alpar@9	141 *
alpar@9	142 * c[q] > + eps, (5)
alpar@9	143 *
alpar@9	144 * the upper column bound u[q] must be active to provide dual
alpar@9	145 * feasibility. In this case the column can be fixed on its upper bound
alpar@9	146 * and removed from the problem. And if the column has no upper bound
alpar@9	147 * (u[q] = +oo), the problem has no dual feasible solution.
alpar@9	148 *
alpar@9	149 * Finally, if the following condition holds:
alpar@9	150 *
alpar@9	151 * - eps <= c[q] <= +eps, (6)
alpar@9	152 *
alpar@9	153 * dual feasibility does not depend on a particular value of column q.
alpar@9	154 * In this case the column can be fixed either on its lower bound (if
alpar@9	155 * l[q] > -oo) or on its upper bound (if u[q] < +oo) or at zero (if the
alpar@9	156 * column is unbounded) and then removed from the problem.
alpar@9	157 *
alpar@9	158 * RECOVERING BASIC SOLUTION
alpar@9	159 *
alpar@9	160 * See the routine npp_fixed_col. Having been recovered the column
alpar@9	161 * is assigned status GLP_NS. However, if actually it is not fixed
alpar@9	162 * (l[q] < u[q]), its status should be changed to GLP_NL, GLP_NU, or
alpar@9	163 * GLP_NF depending on which bound it was fixed on transformation stage.
alpar@9	164 *
alpar@9	165 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	166 *
alpar@9	167 * See the routine npp_fixed_col.
alpar@9	168 *
alpar@9	169 * RECOVERING MIP SOLUTION
alpar@9	170 *
alpar@9	171 * See the routine npp_fixed_col. */
alpar@9	172
alpar@9	173 struct empty_col
alpar@9	174 { /* empty column */
alpar@9	175 int q;
alpar@9	176 /* column reference number */
alpar@9	177 char stat;
alpar@9	178 /* status in basic solution */
alpar@9	179 };
alpar@9	180
alpar@9	181 static int rcv_empty_col(NPP npp, void info);
alpar@9	182
alpar@9	183 int npp_empty_col(NPP npp, NPPCOL q)
alpar@9	184 { /* process empty column */
alpar@9	185 struct empty_col *info;
alpar@9	186 double eps = 1e-3;
alpar@9	187 /* the column must be empty */
alpar@9	188 xassert(q->ptr == NULL);
alpar@9	189 /* check dual feasibility */
alpar@9	190 if (q->coef > +eps && q->lb == -DBL_MAX)
alpar@9	191 return 1;
alpar@9	192 if (q->coef < -eps && q->ub == +DBL_MAX)
alpar@9	193 return 1;
alpar@9	194 /* create transformation stack entry */
alpar@9	195 info = npp_push_tse(npp,
alpar@9	196 rcv_empty_col, sizeof(struct empty_col));
alpar@9	197 info->q = q->j;
alpar@9	198 /* fix the column */
alpar@9	199 if (q->lb == -DBL_MAX && q->ub == +DBL_MAX)
alpar@9	200 { /* free column */
alpar@9	201 info->stat = GLP_NF;
alpar@9	202 q->lb = q->ub = 0.0;
alpar@9	203 }
alpar@9	204 else if (q->ub == +DBL_MAX)
alpar@9	205 lo: { /* column with lower bound */
alpar@9	206 info->stat = GLP_NL;
alpar@9	207 q->ub = q->lb;
alpar@9	208 }
alpar@9	209 else if (q->lb == -DBL_MAX)
alpar@9	210 up: { /* column with upper bound */
alpar@9	211 info->stat = GLP_NU;
alpar@9	212 q->lb = q->ub;
alpar@9	213 }
alpar@9	214 else if (q->lb != q->ub)
alpar@9	215 { /* double-bounded column */
alpar@9	216 if (q->coef >= +DBL_EPSILON) goto lo;
alpar@9	217 if (q->coef <= -DBL_EPSILON) goto up;
alpar@9	218 if (fabs(q->lb) <= fabs(q->ub)) goto lo; else goto up;
alpar@9	219 }
alpar@9	220 else
alpar@9	221 { /* fixed column */
alpar@9	222 info->stat = GLP_NS;
alpar@9	223 }
alpar@9	224 /* process fixed column */
alpar@9	225 npp_fixed_col(npp, q);
alpar@9	226 return 0;
alpar@9	227 }
alpar@9	228
alpar@9	229 static int rcv_empty_col(NPP npp, void _info)
alpar@9	230 { /* recover empty column */
alpar@9	231 struct empty_col *info = _info;
alpar@9	232 if (npp->sol == GLP_SOL)
alpar@9	233 npp->c_stat[info->q] = info->stat;
alpar@9	234 return 0;
alpar@9	235 }
alpar@9	236
alpar@9	237 /***********************************************************************
alpar@9	238 * NAME
alpar@9	239 *
alpar@9	240 * npp_implied_value - process implied column value
alpar@9	241 *
alpar@9	242 * SYNOPSIS
alpar@9	243 *
alpar@9	244 * #include "glpnpp.h"
alpar@9	245 * int npp_implied_value(NPP npp, NPPCOL q, double s);
alpar@9	246 *
alpar@9	247 * DESCRIPTION
alpar@9	248 *
alpar@9	249 * For column q:
alpar@9	250 *
alpar@9	251 * l[q] <= x[q] <= u[q], (1)
alpar@9	252 *
alpar@9	253 * where l[q] < u[q], the routine npp_implied_value processes its
alpar@9	254 * implied value s[q]. If this implied value satisfies to the current
alpar@9	255 * column bounds and integrality condition, the routine fixes column q
alpar@9	256 * at the given point. Note that the column is kept in the problem in
alpar@9	257 * any case.
alpar@9	258 *
alpar@9	259 * RETURNS
alpar@9	260 *
alpar@9	261 * 0 - column has been fixed;
alpar@9	262 *
alpar@9	263 * 1 - implied value violates to current column bounds;
alpar@9	264 *
alpar@9	265 * 2 - implied value violates integrality condition.
alpar@9	266 *
alpar@9	267 * ALGORITHM
alpar@9	268 *
alpar@9	269 * Implied column value s[q] satisfies to the current column bounds if
alpar@9	270 * the following condition holds:
alpar@9	271 *
alpar@9	272 * l[q] - eps <= s[q] <= u[q] + eps, (2)
alpar@9	273 *
alpar@9	274 * where eps is an absolute tolerance for column value. If the column
alpar@9	275 * is integral, the following condition also must hold:
alpar@9	276 *
alpar@9	277 * \|s[q] - floor(s[q]+0.5)\| <= eps, (3)
alpar@9	278 *
alpar@9	279 * where floor(s[q]+0.5) is the nearest integer to s[q].
alpar@9	280 *
alpar@9	281 * If both condition (2) and (3) are satisfied, the column can be fixed
alpar@9	282 * at the value s[q], or, if it is integral, at floor(s[q]+0.5).
alpar@9	283 * Otherwise, if s[q] violates (2) or (3), the problem has no feasible
alpar@9	284 * solution.
alpar@9	285 *
alpar@9	286 * Note: If s[q] is close to l[q] or u[q], it seems to be reasonable to
alpar@9	287 * fix the column at its lower or upper bound, resp. rather than at the
alpar@9	288 * implied value. */
alpar@9	289
alpar@9	290 int npp_implied_value(NPP npp, NPPCOL q, double s)
alpar@9	291 { /* process implied column value */
alpar@9	292 double eps, nint;
alpar@9	293 xassert(npp == npp);
alpar@9	294 /* column must not be fixed */
alpar@9	295 xassert(q->lb < q->ub);
alpar@9	296 /* check integrality */
alpar@9	297 if (q->is_int)
alpar@9	298 { nint = floor(s + 0.5);
alpar@9	299 if (fabs(s - nint) <= 1e-5)
alpar@9	300 s = nint;
alpar@9	301 else
alpar@9	302 return 2;
alpar@9	303 }
alpar@9	304 /* check current column lower bound */
alpar@9	305 if (q->lb != -DBL_MAX)
alpar@9	306 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
alpar@9	307 if (s < q->lb - eps) return 1;
alpar@9	308 /* if s[q] is close to l[q], fix column at its lower bound
alpar@9	309 rather than at the implied value */
alpar@9	310 if (s < q->lb + 1e-3 * eps)
alpar@9	311 { q->ub = q->lb;
alpar@9	312 return 0;
alpar@9	313 }
alpar@9	314 }
alpar@9	315 /* check current column upper bound */
alpar@9	316 if (q->ub != +DBL_MAX)
alpar@9	317 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
alpar@9	318 if (s > q->ub + eps) return 1;
alpar@9	319 /* if s[q] is close to u[q], fix column at its upper bound
alpar@9	320 rather than at the implied value */
alpar@9	321 if (s > q->ub - 1e-3 * eps)
alpar@9	322 { q->lb = q->ub;
alpar@9	323 return 0;
alpar@9	324 }
alpar@9	325 }
alpar@9	326 /* fix column at the implied value */
alpar@9	327 q->lb = q->ub = s;
alpar@9	328 return 0;
alpar@9	329 }
alpar@9	330
alpar@9	331 /***********************************************************************
alpar@9	332 * NAME
alpar@9	333 *
alpar@9	334 * npp_eq_singlet - process row singleton (equality constraint)
alpar@9	335 *
alpar@9	336 * SYNOPSIS
alpar@9	337 *
alpar@9	338 * #include "glpnpp.h"
alpar@9	339 * int npp_eq_singlet(NPP npp, NPPROW p);
alpar@9	340 *
alpar@9	341 * DESCRIPTION
alpar@9	342 *
alpar@9	343 * The routine npp_eq_singlet processes row p, which is equiality
alpar@9	344 * constraint having the only non-zero coefficient:
alpar@9	345 *
alpar@9	346 * a[p,q] x[q] = b. (1)
alpar@9	347 *
alpar@9	348 * RETURNS
alpar@9	349 *
alpar@9	350 * 0 - success;
alpar@9	351 *
alpar@9	352 * 1 - problem has no primal feasible solution;
alpar@9	353 *
alpar@9	354 * 2 - problem has no integer feasible solution.
alpar@9	355 *
alpar@9	356 * PROBLEM TRANSFORMATION
alpar@9	357 *
alpar@9	358 * The equality constraint defines implied value of column q:
alpar@9	359 *
alpar@9	360 * x[q] = s[q] = b / a[p,q]. (2)
alpar@9	361 *
alpar@9	362 * If the implied value s[q] satisfies to the column bounds (see the
alpar@9	363 * routine npp_implied_value), the column can be fixed at s[q] and
alpar@9	364 * removed from the problem. In this case row p becomes redundant, so
alpar@9	365 * it can be replaced by equivalent free row and also removed from the
alpar@9	366 * problem.
alpar@9	367 *
alpar@9	368 * Note that the routine removes from the problem only row p. Column q
alpar@9	369 * becomes fixed, however, it is kept in the problem.
alpar@9	370 *
alpar@9	371 * RECOVERING BASIC SOLUTION
alpar@9	372 *
alpar@9	373 * In solution to the original problem row p is assigned status GLP_NS
alpar@9	374 * (active equality constraint), and column q is assigned status GLP_BS
alpar@9	375 * (basic column).
alpar@9	376 *
alpar@9	377 * Multiplier for row p can be computed as follows. In the dual system
alpar@9	378 * of the original problem column q corresponds to the following row:
alpar@9	379 *
alpar@9	380 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>
alpar@9	381 * i
alpar@9	382 *
alpar@9	383 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q].
alpar@9	384 * i!=p
alpar@9	385 *
alpar@9	386 * Therefore:
alpar@9	387 *
alpar@9	388 * 1
alpar@9	389 * pi[p] = ------ (c[q] - lambda[q] - sum a[i,q] pi[i]), (3)
alpar@9	390 * a[p,q] i!=q
alpar@9	391 *
alpar@9	392 * where lambda[q] = 0 (since column[q] is basic), and pi[i] for all
alpar@9	393 * i != p are known in solution to the transformed problem.
alpar@9	394 *
alpar@9	395 * Value of column q in solution to the original problem is assigned
alpar@9	396 * its implied value s[q].
alpar@9	397 *
alpar@9	398 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	399 *
alpar@9	400 * Multiplier for row p is computed with formula (3). Value of column
alpar@9	401 * q is assigned its implied value s[q].
alpar@9	402 *
alpar@9	403 * RECOVERING MIP SOLUTION
alpar@9	404 *
alpar@9	405 * Value of column q is assigned its implied value s[q]. */
alpar@9	406
alpar@9	407 struct eq_singlet
alpar@9	408 { /* row singleton (equality constraint) */
alpar@9	409 int p;
alpar@9	410 /* row reference number */
alpar@9	411 int q;
alpar@9	412 /* column reference number */
alpar@9	413 double apq;
alpar@9	414 /* constraint coefficient a[p,q] */
alpar@9	415 double c;
alpar@9	416 /* objective coefficient at x[q] */
alpar@9	417 NPPLFE *ptr;
alpar@9	418 /* list of non-zero coefficients a[i,q], i != p */
alpar@9	419 };
alpar@9	420
alpar@9	421 static int rcv_eq_singlet(NPP npp, void info);
alpar@9	422
alpar@9	423 int npp_eq_singlet(NPP npp, NPPROW p)
alpar@9	424 { /* process row singleton (equality constraint) */
alpar@9	425 struct eq_singlet *info;
alpar@9	426 NPPCOL *q;
alpar@9	427 NPPAIJ *aij;
alpar@9	428 NPPLFE *lfe;
alpar@9	429 int ret;
alpar@9	430 double s;
alpar@9	431 /* the row must be singleton equality constraint */
alpar@9	432 xassert(p->lb == p->ub);
alpar@9	433 xassert(p->ptr != NULL && p->ptr->r_next == NULL);
alpar@9	434 /* compute and process implied column value */
alpar@9	435 aij = p->ptr;
alpar@9	436 q = aij->col;
alpar@9	437 s = p->lb / aij->val;
alpar@9	438 ret = npp_implied_value(npp, q, s);
alpar@9	439 xassert(0 <= ret && ret <= 2);
alpar@9	440 if (ret != 0) return ret;
alpar@9	441 /* create transformation stack entry */
alpar@9	442 info = npp_push_tse(npp,
alpar@9	443 rcv_eq_singlet, sizeof(struct eq_singlet));
alpar@9	444 info->p = p->i;
alpar@9	445 info->q = q->j;
alpar@9	446 info->apq = aij->val;
alpar@9	447 info->c = q->coef;
alpar@9	448 info->ptr = NULL;
alpar@9	449 /* save column coefficients a[i,q], i != p (not needed for MIP
alpar@9	450 solution) */
alpar@9	451 if (npp->sol != GLP_MIP)
alpar@9	452 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9	453 { if (aij->row == p) continue; /* skip a[p,q] */
alpar@9	454 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9	455 lfe->ref = aij->row->i;
alpar@9	456 lfe->val = aij->val;
alpar@9	457 lfe->next = info->ptr;
alpar@9	458 info->ptr = lfe;
alpar@9	459 }
alpar@9	460 }
alpar@9	461 /* remove the row from the problem */
alpar@9	462 npp_del_row(npp, p);
alpar@9	463 return 0;
alpar@9	464 }
alpar@9	465
alpar@9	466 static int rcv_eq_singlet(NPP npp, void _info)
alpar@9	467 { /* recover row singleton (equality constraint) */
alpar@9	468 struct eq_singlet *info = _info;
alpar@9	469 NPPLFE *lfe;
alpar@9	470 double temp;
alpar@9	471 if (npp->sol == GLP_SOL)
alpar@9	472 { /* column q must be already recovered as GLP_NS */
alpar@9	473 if (npp->c_stat[info->q] != GLP_NS)
alpar@9	474 { npp_error();
alpar@9	475 return 1;
alpar@9	476 }
alpar@9	477 npp->r_stat[info->p] = GLP_NS;
alpar@9	478 npp->c_stat[info->q] = GLP_BS;
alpar@9	479 }
alpar@9	480 if (npp->sol != GLP_MIP)
alpar@9	481 { /* compute multiplier for row p with formula (3) */
alpar@9	482 temp = info->c;
alpar@9	483 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
alpar@9	484 temp -= lfe->val * npp->r_pi[lfe->ref];
alpar@9	485 npp->r_pi[info->p] = temp / info->apq;
alpar@9	486 }
alpar@9	487 return 0;
alpar@9	488 }
alpar@9	489
alpar@9	490 /***********************************************************************
alpar@9	491 * NAME
alpar@9	492 *
alpar@9	493 * npp_implied_lower - process implied column lower bound
alpar@9	494 *
alpar@9	495 * SYNOPSIS
alpar@9	496 *
alpar@9	497 * #include "glpnpp.h"
alpar@9	498 * int npp_implied_lower(NPP npp, NPPCOL q, double l);
alpar@9	499 *
alpar@9	500 * DESCRIPTION
alpar@9	501 *
alpar@9	502 * For column q:
alpar@9	503 *
alpar@9	504 * l[q] <= x[q] <= u[q], (1)
alpar@9	505 *
alpar@9	506 * where l[q] < u[q], the routine npp_implied_lower processes its
alpar@9	507 * implied lower bound l'[q]. As the result the current column lower
alpar@9	508 * bound may increase. Note that the column is kept in the problem in
alpar@9	509 * any case.
alpar@9	510 *
alpar@9	511 * RETURNS
alpar@9	512 *
alpar@9	513 * 0 - current column lower bound has not changed;
alpar@9	514 *
alpar@9	515 * 1 - current column lower bound has changed, but not significantly;
alpar@9	516 *
alpar@9	517 * 2 - current column lower bound has significantly changed;
alpar@9	518 *
alpar@9	519 * 3 - column has been fixed on its upper bound;
alpar@9	520 *
alpar@9	521 * 4 - implied lower bound violates current column upper bound.
alpar@9	522 *
alpar@9	523 * ALGORITHM
alpar@9	524 *
alpar@9	525 * If column q is integral, before processing its implied lower bound
alpar@9	526 * should be rounded up:
alpar@9	527 *
alpar@9	528 * ( floor(l'[q]+0.5), if \|l'[q] - floor(l'[q]+0.5)\| <= eps
alpar@9	529 * l'[q] := < (2)
alpar@9	530 * ( ceil(l'[q]), otherwise
alpar@9	531 *
alpar@9	532 * where floor(l'[q]+0.5) is the nearest integer to l'[q], ceil(l'[q])
alpar@9	533 * is smallest integer not less than l'[q], and eps is an absolute
alpar@9	534 * tolerance for column value.
alpar@9	535 *
alpar@9	536 * Processing implied column lower bound l'[q] includes the following
alpar@9	537 * cases:
alpar@9	538 *
alpar@9	539 * 1) if l'[q] < l[q] + eps, implied lower bound is redundant;
alpar@9	540 *
alpar@9	541 * 2) if l[q] + eps <= l[q] <= u[q] + eps, current column lower bound
alpar@9	542 * l[q] can be strengthened by replacing it with l'[q]. If in this
alpar@9	543 * case new column lower bound becomes close to current column upper
alpar@9	544 * bound u[q], the column can be fixed on its upper bound;
alpar@9	545 *
alpar@9	546 * 3) if l'[q] > u[q] + eps, implied lower bound violates current
alpar@9	547 * column upper bound u[q], in which case the problem has no primal
alpar@9	548 * feasible solution. */
alpar@9	549
alpar@9	550 int npp_implied_lower(NPP npp, NPPCOL q, double l)
alpar@9	551 { /* process implied column lower bound */
alpar@9	552 int ret;
alpar@9	553 double eps, nint;
alpar@9	554 xassert(npp == npp);
alpar@9	555 /* column must not be fixed */
alpar@9	556 xassert(q->lb < q->ub);
alpar@9	557 /* implied lower bound must be finite */
alpar@9	558 xassert(l != -DBL_MAX);
alpar@9	559 /* if column is integral, round up l'[q] */
alpar@9	560 if (q->is_int)
alpar@9	561 { nint = floor(l + 0.5);
alpar@9	562 if (fabs(l - nint) <= 1e-5)
alpar@9	563 l = nint;
alpar@9	564 else
alpar@9	565 l = ceil(l);
alpar@9	566 }
alpar@9	567 /* check current column lower bound */
alpar@9	568 if (q->lb != -DBL_MAX)
alpar@9	569 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->lb));
alpar@9	570 if (l < q->lb + eps)
alpar@9	571 { ret = 0; /* redundant */
alpar@9	572 goto done;
alpar@9	573 }
alpar@9	574 }
alpar@9	575 /* check current column upper bound */
alpar@9	576 if (q->ub != +DBL_MAX)
alpar@9	577 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->ub));
alpar@9	578 if (l > q->ub + eps)
alpar@9	579 { ret = 4; /* infeasible */
alpar@9	580 goto done;
alpar@9	581 }
alpar@9	582 /* if l'[q] is close to u[q], fix column at its upper bound */
alpar@9	583 if (l > q->ub - 1e-3 * eps)
alpar@9	584 { q->lb = q->ub;
alpar@9	585 ret = 3; /* fixed */
alpar@9	586 goto done;
alpar@9	587 }
alpar@9	588 }
alpar@9	589 /* check if column lower bound changes significantly */
alpar@9	590 if (q->lb == -DBL_MAX)
alpar@9	591 ret = 2; /* significantly */
alpar@9	592 else if (q->is_int && l > q->lb + 0.5)
alpar@9	593 ret = 2; /* significantly */
alpar@9	594 else if (l > q->lb + 0.30 * (1.0 + fabs(q->lb)))
alpar@9	595 ret = 2; /* significantly */
alpar@9	596 else
alpar@9	597 ret = 1; /* not significantly */
alpar@9	598 /* set new column lower bound */
alpar@9	599 q->lb = l;
alpar@9	600 done: return ret;
alpar@9	601 }
alpar@9	602
alpar@9	603 /***********************************************************************
alpar@9	604 * NAME
alpar@9	605 *
alpar@9	606 * npp_implied_upper - process implied column upper bound
alpar@9	607 *
alpar@9	608 * SYNOPSIS
alpar@9	609 *
alpar@9	610 * #include "glpnpp.h"
alpar@9	611 * int npp_implied_upper(NPP npp, NPPCOL q, double u);
alpar@9	612 *
alpar@9	613 * DESCRIPTION
alpar@9	614 *
alpar@9	615 * For column q:
alpar@9	616 *
alpar@9	617 * l[q] <= x[q] <= u[q], (1)
alpar@9	618 *
alpar@9	619 * where l[q] < u[q], the routine npp_implied_upper processes its
alpar@9	620 * implied upper bound u'[q]. As the result the current column upper
alpar@9	621 * bound may decrease. Note that the column is kept in the problem in
alpar@9	622 * any case.
alpar@9	623 *
alpar@9	624 * RETURNS
alpar@9	625 *
alpar@9	626 * 0 - current column upper bound has not changed;
alpar@9	627 *
alpar@9	628 * 1 - current column upper bound has changed, but not significantly;
alpar@9	629 *
alpar@9	630 * 2 - current column upper bound has significantly changed;
alpar@9	631 *
alpar@9	632 * 3 - column has been fixed on its lower bound;
alpar@9	633 *
alpar@9	634 * 4 - implied upper bound violates current column lower bound.
alpar@9	635 *
alpar@9	636 * ALGORITHM
alpar@9	637 *
alpar@9	638 * If column q is integral, before processing its implied upper bound
alpar@9	639 * should be rounded down:
alpar@9	640 *
alpar@9	641 * ( floor(u'[q]+0.5), if \|u'[q] - floor(l'[q]+0.5)\| <= eps
alpar@9	642 * u'[q] := < (2)
alpar@9	643 * ( floor(l'[q]), otherwise
alpar@9	644 *
alpar@9	645 * where floor(u'[q]+0.5) is the nearest integer to u'[q],
alpar@9	646 * floor(u'[q]) is largest integer not greater than u'[q], and eps is
alpar@9	647 * an absolute tolerance for column value.
alpar@9	648 *
alpar@9	649 * Processing implied column upper bound u'[q] includes the following
alpar@9	650 * cases:
alpar@9	651 *
alpar@9	652 * 1) if u'[q] > u[q] - eps, implied upper bound is redundant;
alpar@9	653 *
alpar@9	654 * 2) if l[q] - eps <= u[q] <= u[q] - eps, current column upper bound
alpar@9	655 * u[q] can be strengthened by replacing it with u'[q]. If in this
alpar@9	656 * case new column upper bound becomes close to current column lower
alpar@9	657 * bound, the column can be fixed on its lower bound;
alpar@9	658 *
alpar@9	659 * 3) if u'[q] < l[q] - eps, implied upper bound violates current
alpar@9	660 * column lower bound l[q], in which case the problem has no primal
alpar@9	661 * feasible solution. */
alpar@9	662
alpar@9	663 int npp_implied_upper(NPP npp, NPPCOL q, double u)
alpar@9	664 { int ret;
alpar@9	665 double eps, nint;
alpar@9	666 xassert(npp == npp);
alpar@9	667 /* column must not be fixed */
alpar@9	668 xassert(q->lb < q->ub);
alpar@9	669 /* implied upper bound must be finite */
alpar@9	670 xassert(u != +DBL_MAX);
alpar@9	671 /* if column is integral, round down u'[q] */
alpar@9	672 if (q->is_int)
alpar@9	673 { nint = floor(u + 0.5);
alpar@9	674 if (fabs(u - nint) <= 1e-5)
alpar@9	675 u = nint;
alpar@9	676 else
alpar@9	677 u = floor(u);
alpar@9	678 }
alpar@9	679 /* check current column upper bound */
alpar@9	680 if (q->ub != +DBL_MAX)
alpar@9	681 { eps = (q->is_int ? 1e-3 : 1e-3 + 1e-6 * fabs(q->ub));
alpar@9	682 if (u > q->ub - eps)
alpar@9	683 { ret = 0; /* redundant */
alpar@9	684 goto done;
alpar@9	685 }
alpar@9	686 }
alpar@9	687 /* check current column lower bound */
alpar@9	688 if (q->lb != -DBL_MAX)
alpar@9	689 { eps = (q->is_int ? 1e-5 : 1e-5 + 1e-8 * fabs(q->lb));
alpar@9	690 if (u < q->lb - eps)
alpar@9	691 { ret = 4; /* infeasible */
alpar@9	692 goto done;
alpar@9	693 }
alpar@9	694 /* if u'[q] is close to l[q], fix column at its lower bound */
alpar@9	695 if (u < q->lb + 1e-3 * eps)
alpar@9	696 { q->ub = q->lb;
alpar@9	697 ret = 3; /* fixed */
alpar@9	698 goto done;
alpar@9	699 }
alpar@9	700 }
alpar@9	701 /* check if column upper bound changes significantly */
alpar@9	702 if (q->ub == +DBL_MAX)
alpar@9	703 ret = 2; /* significantly */
alpar@9	704 else if (q->is_int && u < q->ub - 0.5)
alpar@9	705 ret = 2; /* significantly */
alpar@9	706 else if (u < q->ub - 0.30 * (1.0 + fabs(q->ub)))
alpar@9	707 ret = 2; /* significantly */
alpar@9	708 else
alpar@9	709 ret = 1; /* not significantly */
alpar@9	710 /* set new column upper bound */
alpar@9	711 q->ub = u;
alpar@9	712 done: return ret;
alpar@9	713 }
alpar@9	714
alpar@9	715 /***********************************************************************
alpar@9	716 * NAME
alpar@9	717 *
alpar@9	718 * npp_ineq_singlet - process row singleton (inequality constraint)
alpar@9	719 *
alpar@9	720 * SYNOPSIS
alpar@9	721 *
alpar@9	722 * #include "glpnpp.h"
alpar@9	723 * int npp_ineq_singlet(NPP npp, NPPROW p);
alpar@9	724 *
alpar@9	725 * DESCRIPTION
alpar@9	726 *
alpar@9	727 * The routine npp_ineq_singlet processes row p, which is inequality
alpar@9	728 * constraint having the only non-zero coefficient:
alpar@9	729 *
alpar@9	730 * L[p] <= a[p,q] * x[q] <= U[p], (1)
alpar@9	731 *
alpar@9	732 * where L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
alpar@9	733 *
alpar@9	734 * RETURNS
alpar@9	735 *
alpar@9	736 * 0 - current column bounds have not changed;
alpar@9	737 *
alpar@9	738 * 1 - current column bounds have changed, but not significantly;
alpar@9	739 *
alpar@9	740 * 2 - current column bounds have significantly changed;
alpar@9	741 *
alpar@9	742 * 3 - column has been fixed on its lower or upper bound;
alpar@9	743 *
alpar@9	744 * 4 - problem has no primal feasible solution.
alpar@9	745 *
alpar@9	746 * PROBLEM TRANSFORMATION
alpar@9	747 *
alpar@9	748 * Inequality constraint (1) defines implied bounds of column q:
alpar@9	749 *
alpar@9	750 * ( L[p] / a[p,q], if a[p,q] > 0
alpar@9	751 * l'[q] = < (2)
alpar@9	752 * ( U[p] / a[p,q], if a[p,q] < 0
alpar@9	753 *
alpar@9	754 * ( U[p] / a[p,q], if a[p,q] > 0
alpar@9	755 * u'[q] = < (3)
alpar@9	756 * ( L[p] / a[p,q], if a[p,q] < 0
alpar@9	757 *
alpar@9	758 * If these implied bounds do not violate current bounds of column q:
alpar@9	759 *
alpar@9	760 * l[q] <= x[q] <= u[q], (4)
alpar@9	761 *
alpar@9	762 * they can be used to strengthen the current column bounds:
alpar@9	763 *
alpar@9	764 * l[q] := max(l[q], l'[q]), (5)
alpar@9	765 *
alpar@9	766 * u[q] := min(u[q], u'[q]). (6)
alpar@9	767 *
alpar@9	768 * (See the routines npp_implied_lower and npp_implied_upper.)
alpar@9	769 *
alpar@9	770 * Once bounds of row p (1) have been carried over column q, the row
alpar@9	771 * becomes redundant, so it can be replaced by equivalent free row and
alpar@9	772 * removed from the problem.
alpar@9	773 *
alpar@9	774 * Note that the routine removes from the problem only row p. Column q,
alpar@9	775 * even it has been fixed, is kept in the problem.
alpar@9	776 *
alpar@9	777 * RECOVERING BASIC SOLUTION
alpar@9	778 *
alpar@9	779 * Note that the row in the dual system corresponding to column q is
alpar@9	780 * the following:
alpar@9	781 *
alpar@9	782 * sum a[i,q] pi[i] + lambda[q] = c[q] ==>
alpar@9	783 * i
alpar@9	784 * (7)
alpar@9	785 * sum a[i,q] pi[i] + a[p,q] pi[p] + lambda[q] = c[q],
alpar@9	786 * i!=p
alpar@9	787 *
alpar@9	788 * where pi[i] for all i != p are known in solution to the transformed
alpar@9	789 * problem. Row p does not exist in the transformed problem, so it has
alpar@9	790 * zero multiplier there. This allows computing multiplier for column q
alpar@9	791 * in solution to the transformed problem:
alpar@9	792 *
alpar@9	793 * lambda~[q] = c[q] - sum a[i,q] pi[i]. (8)
alpar@9	794 * i!=p
alpar@9	795 *
alpar@9	796 * Let in solution to the transformed problem column q be non-basic
alpar@9	797 * with lower bound active (GLP_NL, lambda~[q] >= 0), and this lower
alpar@9	798 * bound be implied one l'[q]. From the original problem's standpoint
alpar@9	799 * this then means that actually the original column lower bound l[q]
alpar@9	800 * is inactive, and active is that row bound L[p] or U[p] that defines
alpar@9	801 * the implied bound l'[q] (2). In this case in solution to the
alpar@9	802 * original problem column q is assigned status GLP_BS while row p is
alpar@9	803 * assigned status GLP_NL (if a[p,q] > 0) or GLP_NU (if a[p,q] < 0).
alpar@9	804 * Since now column q is basic, its multiplier lambda[q] is zero. This
alpar@9	805 * allows using (7) and (8) to find multiplier for row p in solution to
alpar@9	806 * the original problem:
alpar@9	807 *
alpar@9	808 * 1
alpar@9	809 * pi[p] = ------ (c[q] - sum a[i,q] pi[i]) = lambda~[q] / a[p,q] (9)
alpar@9	810 * a[p,q] i!=p
alpar@9	811 *
alpar@9	812 * Now let in solution to the transformed problem column q be non-basic
alpar@9	813 * with upper bound active (GLP_NU, lambda~[q] <= 0), and this upper
alpar@9	814 * bound be implied one u'[q]. As in the previous case this then means
alpar@9	815 * that from the original problem's standpoint actually the original
alpar@9	816 * column upper bound u[q] is inactive, and active is that row bound
alpar@9	817 * L[p] or U[p] that defines the implied bound u'[q] (3). In this case
alpar@9	818 * in solution to the original problem column q is assigned status
alpar@9	819 * GLP_BS, row p is assigned status GLP_NU (if a[p,q] > 0) or GLP_NL
alpar@9	820 * (if a[p,q] < 0), and its multiplier is computed with formula (9).
alpar@9	821 *
alpar@9	822 * Strengthening bounds of column q according to (5) and (6) may make
alpar@9	823 * it fixed. Thus, if in solution to the transformed problem column q is
alpar@9	824 * non-basic and fixed (GLP_NS), we can suppose that if lambda~[q] > 0,
alpar@9	825 * column q has active lower bound (GLP_NL), and if lambda~[q] < 0,
alpar@9	826 * column q has active upper bound (GLP_NU), reducing this case to two
alpar@9	827 * previous ones. If, however, lambda~[q] is close to zero or
alpar@9	828 * corresponding bound of row p does not exist (this may happen if
alpar@9	829 * lambda~[q] has wrong sign due to round-off errors, in which case it
alpar@9	830 * is expected to be close to zero, since solution is assumed to be dual
alpar@9	831 * feasible), column q can be assigned status GLP_BS (basic), and row p
alpar@9	832 * can be made active on its existing bound. In the latter case row
alpar@9	833 * multiplier pi[p] computed with formula (9) will be also close to
alpar@9	834 * zero, and dual feasibility will be kept.
alpar@9	835 *
alpar@9	836 * In all other cases, namely, if in solution to the transformed
alpar@9	837 * problem column q is basic (GLP_BS), or non-basic with original lower
alpar@9	838 * bound l[q] active (GLP_NL), or non-basic with original upper bound
alpar@9	839 * u[q] active (GLP_NU), constraint (1) is inactive. So in solution to
alpar@9	840 * the original problem status of column q remains unchanged, row p is
alpar@9	841 * assigned status GLP_BS, and its multiplier pi[p] is assigned zero
alpar@9	842 * value.
alpar@9	843 *
alpar@9	844 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	845 *
alpar@9	846 * First, value of multiplier for column q in solution to the original
alpar@9	847 * problem is computed with formula (8). If lambda~[q] > 0 and column q
alpar@9	848 * has implied lower bound, or if lambda~[q] < 0 and column q has
alpar@9	849 * implied upper bound, this means that from the original problem's
alpar@9	850 * standpoint actually row p has corresponding active bound, in which
alpar@9	851 * case its multiplier pi[p] is computed with formula (9). In other
alpar@9	852 * cases, when the sign of lambda~[q] corresponds to original bound of
alpar@9	853 * column q, or when lambda~[q] =~ 0, value of row multiplier pi[p] is
alpar@9	854 * assigned zero value.
alpar@9	855 *
alpar@9	856 * RECOVERING MIP SOLUTION
alpar@9	857 *
alpar@9	858 * None needed. */
alpar@9	859
alpar@9	860 struct ineq_singlet
alpar@9	861 { /* row singleton (inequality constraint) */
alpar@9	862 int p;
alpar@9	863 /* row reference number */
alpar@9	864 int q;
alpar@9	865 /* column reference number */
alpar@9	866 double apq;
alpar@9	867 /* constraint coefficient a[p,q] */
alpar@9	868 double c;
alpar@9	869 /* objective coefficient at x[q] */
alpar@9	870 double lb;
alpar@9	871 /* row lower bound */
alpar@9	872 double ub;
alpar@9	873 /* row upper bound */
alpar@9	874 char lb_changed;
alpar@9	875 /* this flag is set if column lower bound was changed */
alpar@9	876 char ub_changed;
alpar@9	877 /* this flag is set if column upper bound was changed */
alpar@9	878 NPPLFE *ptr;
alpar@9	879 /* list of non-zero coefficients a[i,q], i != p */
alpar@9	880 };
alpar@9	881
alpar@9	882 static int rcv_ineq_singlet(NPP npp, void info);
alpar@9	883
alpar@9	884 int npp_ineq_singlet(NPP npp, NPPROW p)
alpar@9	885 { /* process row singleton (inequality constraint) */
alpar@9	886 struct ineq_singlet *info;
alpar@9	887 NPPCOL *q;
alpar@9	888 NPPAIJ apq, aij;
alpar@9	889 NPPLFE *lfe;
alpar@9	890 int lb_changed, ub_changed;
alpar@9	891 double ll, uu;
alpar@9	892 /* the row must be singleton inequality constraint */
alpar@9	893 xassert(p->lb != -DBL_MAX \|\| p->ub != +DBL_MAX);
alpar@9	894 xassert(p->lb < p->ub);
alpar@9	895 xassert(p->ptr != NULL && p->ptr->r_next == NULL);
alpar@9	896 /* compute implied column bounds */
alpar@9	897 apq = p->ptr;
alpar@9	898 q = apq->col;
alpar@9	899 xassert(q->lb < q->ub);
alpar@9	900 if (apq->val > 0.0)
alpar@9	901 { ll = (p->lb == -DBL_MAX ? -DBL_MAX : p->lb / apq->val);
alpar@9	902 uu = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub / apq->val);
alpar@9	903 }
alpar@9	904 else
alpar@9	905 { ll = (p->ub == +DBL_MAX ? -DBL_MAX : p->ub / apq->val);
alpar@9	906 uu = (p->lb == -DBL_MAX ? +DBL_MAX : p->lb / apq->val);
alpar@9	907 }
alpar@9	908 /* process implied column lower bound */
alpar@9	909 if (ll == -DBL_MAX)
alpar@9	910 lb_changed = 0;
alpar@9	911 else
alpar@9	912 { lb_changed = npp_implied_lower(npp, q, ll);
alpar@9	913 xassert(0 <= lb_changed && lb_changed <= 4);
alpar@9	914 if (lb_changed == 4) return 4; /* infeasible */
alpar@9	915 }
alpar@9	916 /* process implied column upper bound */
alpar@9	917 if (uu == +DBL_MAX)
alpar@9	918 ub_changed = 0;
alpar@9	919 else if (lb_changed == 3)
alpar@9	920 { /* column was fixed on its upper bound due to l'[q] = u[q] */
alpar@9	921 /* note that L[p] < U[p], so l'[q] = u[q] < u'[q] */
alpar@9	922 ub_changed = 0;
alpar@9	923 }
alpar@9	924 else
alpar@9	925 { ub_changed = npp_implied_upper(npp, q, uu);
alpar@9	926 xassert(0 <= ub_changed && ub_changed <= 4);
alpar@9	927 if (ub_changed == 4) return 4; /* infeasible */
alpar@9	928 }
alpar@9	929 /* if neither lower nor upper column bound was changed, the row
alpar@9	930 is originally redundant and can be replaced by free row */
alpar@9	931 if (!lb_changed && !ub_changed)
alpar@9	932 { p->lb = -DBL_MAX, p->ub = +DBL_MAX;
alpar@9	933 npp_free_row(npp, p);
alpar@9	934 return 0;
alpar@9	935 }
alpar@9	936 /* create transformation stack entry */
alpar@9	937 info = npp_push_tse(npp,
alpar@9	938 rcv_ineq_singlet, sizeof(struct ineq_singlet));
alpar@9	939 info->p = p->i;
alpar@9	940 info->q = q->j;
alpar@9	941 info->apq = apq->val;
alpar@9	942 info->c = q->coef;
alpar@9	943 info->lb = p->lb;
alpar@9	944 info->ub = p->ub;
alpar@9	945 info->lb_changed = (char)lb_changed;
alpar@9	946 info->ub_changed = (char)ub_changed;
alpar@9	947 info->ptr = NULL;
alpar@9	948 /* save column coefficients a[i,q], i != p (not needed for MIP
alpar@9	949 solution) */
alpar@9	950 if (npp->sol != GLP_MIP)
alpar@9	951 { for (aij = q->ptr; aij != NULL; aij = aij->c_next)
alpar@9	952 { if (aij == apq) continue; /* skip a[p,q] */
alpar@9	953 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9	954 lfe->ref = aij->row->i;
alpar@9	955 lfe->val = aij->val;
alpar@9	956 lfe->next = info->ptr;
alpar@9	957 info->ptr = lfe;
alpar@9	958 }
alpar@9	959 }
alpar@9	960 /* remove the row from the problem */
alpar@9	961 npp_del_row(npp, p);
alpar@9	962 return lb_changed >= ub_changed ? lb_changed : ub_changed;
alpar@9	963 }
alpar@9	964
alpar@9	965 static int rcv_ineq_singlet(NPP npp, void _info)
alpar@9	966 { /* recover row singleton (inequality constraint) */
alpar@9	967 struct ineq_singlet *info = _info;
alpar@9	968 NPPLFE *lfe;
alpar@9	969 double lambda;
alpar@9	970 if (npp->sol == GLP_MIP) goto done;
alpar@9	971 /* compute lambda~[q] in solution to the transformed problem
alpar@9	972 with formula (8) */
alpar@9	973 lambda = info->c;
alpar@9	974 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
alpar@9	975 lambda -= lfe->val * npp->r_pi[lfe->ref];
alpar@9	976 if (npp->sol == GLP_SOL)
alpar@9	977 { /* recover basic solution */
alpar@9	978 if (npp->c_stat[info->q] == GLP_BS)
alpar@9	979 { /* column q is basic, so row p is inactive */
alpar@9	980 npp->r_stat[info->p] = GLP_BS;
alpar@9	981 npp->r_pi[info->p] = 0.0;
alpar@9	982 }
alpar@9	983 else if (npp->c_stat[info->q] == GLP_NL)
alpar@9	984 nl: { /* column q is non-basic with lower bound active */
alpar@9	985 if (info->lb_changed)
alpar@9	986 { /* it is implied bound, so actually row p is active
alpar@9	987 while column q is basic */
alpar@9	988 npp->r_stat[info->p] =
alpar@9	989 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
alpar@9	990 npp->c_stat[info->q] = GLP_BS;
alpar@9	991 npp->r_pi[info->p] = lambda / info->apq;
alpar@9	992 }
alpar@9	993 else
alpar@9	994 { /* it is original bound, so row p is inactive */
alpar@9	995 npp->r_stat[info->p] = GLP_BS;
alpar@9	996 npp->r_pi[info->p] = 0.0;
alpar@9	997 }
alpar@9	998 }
alpar@9	999 else if (npp->c_stat[info->q] == GLP_NU)
alpar@9	1000 nu: { /* column q is non-basic with upper bound active */
alpar@9	1001 if (info->ub_changed)
alpar@9	1002 { /* it is implied bound, so actually row p is active
alpar@9	1003 while column q is basic */
alpar@9	1004 npp->r_stat[info->p] =
alpar@9	1005 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
alpar@9	1006 npp->c_stat[info->q] = GLP_BS;
alpar@9	1007 npp->r_pi[info->p] = lambda / info->apq;
alpar@9	1008 }
alpar@9	1009 else
alpar@9	1010 { /* it is original bound, so row p is inactive */
alpar@9	1011 npp->r_stat[info->p] = GLP_BS;
alpar@9	1012 npp->r_pi[info->p] = 0.0;
alpar@9	1013 }
alpar@9	1014 }
alpar@9	1015 else if (npp->c_stat[info->q] == GLP_NS)
alpar@9	1016 { /* column q is non-basic and fixed; note, however, that in
alpar@9	1017 in the original problem it is non-fixed */
alpar@9	1018 if (lambda > +1e-7)
alpar@9	1019 { if (info->apq > 0.0 && info->lb != -DBL_MAX \|\|
alpar@9	1020 info->apq < 0.0 && info->ub != +DBL_MAX \|\|
alpar@9	1021 !info->lb_changed)
alpar@9	1022 { /* either corresponding bound of row p exists or
alpar@9	1023 column q remains non-basic with its original lower
alpar@9	1024 bound active */
alpar@9	1025 npp->c_stat[info->q] = GLP_NL;
alpar@9	1026 goto nl;
alpar@9	1027 }
alpar@9	1028 }
alpar@9	1029 if (lambda < -1e-7)
alpar@9	1030 { if (info->apq > 0.0 && info->ub != +DBL_MAX \|\|
alpar@9	1031 info->apq < 0.0 && info->lb != -DBL_MAX \|\|
alpar@9	1032 !info->ub_changed)
alpar@9	1033 { /* either corresponding bound of row p exists or
alpar@9	1034 column q remains non-basic with its original upper
alpar@9	1035 bound active */
alpar@9	1036 npp->c_stat[info->q] = GLP_NU;
alpar@9	1037 goto nu;
alpar@9	1038 }
alpar@9	1039 }
alpar@9	1040 /* either lambda~[q] is close to zero, or corresponding
alpar@9	1041 bound of row p does not exist, because lambda~[q] has
alpar@9	1042 wrong sign due to round-off errors; in the latter case
alpar@9	1043 lambda~[q] is also assumed to be close to zero; so, we
alpar@9	1044 can make row p active on its existing bound and column q
alpar@9	1045 basic; pi[p] will have wrong sign, but it also will be
alpar@9	1046 close to zero (rarus casus of dual degeneracy) */
alpar@9	1047 if (info->lb != -DBL_MAX && info->ub == +DBL_MAX)
alpar@9	1048 { /* row lower bound exists, but upper bound doesn't */
alpar@9	1049 npp->r_stat[info->p] = GLP_NL;
alpar@9	1050 }
alpar@9	1051 else if (info->lb == -DBL_MAX && info->ub != +DBL_MAX)
alpar@9	1052 { /* row upper bound exists, but lower bound doesn't */
alpar@9	1053 npp->r_stat[info->p] = GLP_NU;
alpar@9	1054 }
alpar@9	1055 else if (info->lb != -DBL_MAX && info->ub != +DBL_MAX)
alpar@9	1056 { /* both row lower and upper bounds exist */
alpar@9	1057 /* to choose proper active row bound we should not use
alpar@9	1058 lambda~[q], because its value being close to zero is
alpar@9	1059 unreliable; so we choose that bound which provides
alpar@9	1060 primal feasibility for original constraint (1) */
alpar@9	1061 if (info->apq * npp->c_value[info->q] <=
alpar@9	1062 0.5 * (info->lb + info->ub))
alpar@9	1063 npp->r_stat[info->p] = GLP_NL;
alpar@9	1064 else
alpar@9	1065 npp->r_stat[info->p] = GLP_NU;
alpar@9	1066 }
alpar@9	1067 else
alpar@9	1068 { npp_error();
alpar@9	1069 return 1;
alpar@9	1070 }
alpar@9	1071 npp->c_stat[info->q] = GLP_BS;
alpar@9	1072 npp->r_pi[info->p] = lambda / info->apq;
alpar@9	1073 }
alpar@9	1074 else
alpar@9	1075 { npp_error();
alpar@9	1076 return 1;
alpar@9	1077 }
alpar@9	1078 }
alpar@9	1079 if (npp->sol == GLP_IPT)
alpar@9	1080 { /* recover interior-point solution */
alpar@9	1081 if (lambda > +DBL_EPSILON && info->lb_changed \|\|
alpar@9	1082 lambda < -DBL_EPSILON && info->ub_changed)
alpar@9	1083 { /* actually row p has corresponding active bound */
alpar@9	1084 npp->r_pi[info->p] = lambda / info->apq;
alpar@9	1085 }
alpar@9	1086 else
alpar@9	1087 { /* either bounds of column q are both inactive or its
alpar@9	1088 original bound is active */
alpar@9	1089 npp->r_pi[info->p] = 0.0;
alpar@9	1090 }
alpar@9	1091 }
alpar@9	1092 done: return 0;
alpar@9	1093 }
alpar@9	1094
alpar@9	1095 /***********************************************************************
alpar@9	1096 * NAME
alpar@9	1097 *
alpar@9	1098 * npp_implied_slack - process column singleton (implied slack variable)
alpar@9	1099 *
alpar@9	1100 * SYNOPSIS
alpar@9	1101 *
alpar@9	1102 * #include "glpnpp.h"
alpar@9	1103 * void npp_implied_slack(NPP npp, NPPCOL q);
alpar@9	1104 *
alpar@9	1105 * DESCRIPTION
alpar@9	1106 *
alpar@9	1107 * The routine npp_implied_slack processes column q:
alpar@9	1108 *
alpar@9	1109 * l[q] <= x[q] <= u[q], (1)
alpar@9	1110 *
alpar@9	1111 * where l[q] < u[q], having the only non-zero coefficient in row p,
alpar@9	1112 * which is equality constraint:
alpar@9	1113 *
alpar@9	1114 * sum a[p,j] x[j] + a[p,q] x[q] = b. (2)
alpar@9	1115 * j!=q
alpar@9	1116 *
alpar@9	1117 * PROBLEM TRANSFORMATION
alpar@9	1118 *
alpar@9	1119 * (If x[q] is integral, this transformation must not be used.)
alpar@9	1120 *
alpar@9	1121 * The term a[p,q] x[q] in constraint (2) can be considered as a slack
alpar@9	1122 * variable that allows to carry bounds of column q over row p and then
alpar@9	1123 * remove column q from the problem.
alpar@9	1124 *
alpar@9	1125 * Constraint (2) can be written as follows:
alpar@9	1126 *
alpar@9	1127 * sum a[p,j] x[j] = b - a[p,q] x[q]. (3)
alpar@9	1128 * j!=q
alpar@9	1129 *
alpar@9	1130 * According to (1) constraint (3) is equivalent to the following
alpar@9	1131 * inequality constraint:
alpar@9	1132 *
alpar@9	1133 * L[p] <= sum a[p,j] x[j] <= U[p], (4)
alpar@9	1134 * j!=q
alpar@9	1135 *
alpar@9	1136 * where
alpar@9	1137 *
alpar@9	1138 * ( b - a[p,q] u[q], if a[p,q] > 0
alpar@9	1139 * L[p] = < (5)
alpar@9	1140 * ( b - a[p,q] l[q], if a[p,q] < 0
alpar@9	1141 *
alpar@9	1142 * ( b - a[p,q] l[q], if a[p,q] > 0
alpar@9	1143 * U[p] = < (6)
alpar@9	1144 * ( b - a[p,q] u[q], if a[p,q] < 0
alpar@9	1145 *
alpar@9	1146 * From (2) it follows that:
alpar@9	1147 *
alpar@9	1148 * 1
alpar@9	1149 * x[q] = ------ (b - sum a[p,j] x[j]). (7)
alpar@9	1150 * a[p,q] j!=q
alpar@9	1151 *
alpar@9	1152 * In order to eliminate x[q] from the objective row we substitute it
alpar@9	1153 * from (6) to that row:
alpar@9	1154 *
alpar@9	1155 * z = sum c[j] x[j] + c[q] x[q] + c[0] =
alpar@9	1156 * j!=q
alpar@9	1157 * 1
alpar@9	1158 * = sum c[j] x[j] + c[q] [------ (b - sum a[p,j] x[j])] + c0 =
alpar@9	1159 * j!=q a[p,q] j!=q
alpar@9	1160 *
alpar@9	1161 * = sum c~[j] x[j] + c~[0],
alpar@9	1162 * j!=q
alpar@9	1163 * a[p,j] b
alpar@9	1164 * c~[j] = c[j] - c[q] ------, c~0 = c0 - c[q] ------ (8)
alpar@9	1165 * a[p,q] a[p,q]
alpar@9	1166 *
alpar@9	1167 * are values of objective coefficients and constant term, resp., in
alpar@9	1168 * the transformed problem.
alpar@9	1169 *
alpar@9	1170 * Note that column q is column singleton, so in the dual system of the
alpar@9	1171 * original problem it corresponds to the following row singleton:
alpar@9	1172 *
alpar@9	1173 * a[p,q] pi[p] + lambda[q] = c[q]. (9)
alpar@9	1174 *
alpar@9	1175 * In the transformed problem row (9) would be the following:
alpar@9	1176 *
alpar@9	1177 * a[p,q] pi~[p] + lambda[q] = c~[q] = 0. (10)
alpar@9	1178 *
alpar@9	1179 * Subtracting (10) from (9) we have:
alpar@9	1180 *
alpar@9	1181 * a[p,q] (pi[p] - pi~[p]) = c[q]
alpar@9	1182 *
alpar@9	1183 * that gives the following formula to compute multiplier for row p in
alpar@9	1184 * solution to the original problem using its value in solution to the
alpar@9	1185 * transformed problem:
alpar@9	1186 *
alpar@9	1187 * pi[p] = pi~[p] + c[q] / a[p,q]. (11)
alpar@9	1188 *
alpar@9	1189 * RECOVERING BASIC SOLUTION
alpar@9	1190 *
alpar@9	1191 * Status of column q in solution to the original problem is defined
alpar@9	1192 * by status of row p in solution to the transformed problem and the
alpar@9	1193 * sign of coefficient a[p,q] in the original inequality constraint (2)
alpar@9	1194 * as follows:
alpar@9	1195 *
alpar@9	1196 * +-----------------------+---------+--------------------+
alpar@9	1197 * \| Status of row p \| Sign of \| Status of column q \|
alpar@9	1198 * \| (transformed problem) \| a[p,q] \| (original problem) \|
alpar@9	1199 * +-----------------------+---------+--------------------+
alpar@9	1200 * \| GLP_BS \| + / - \| GLP_BS \|
alpar@9	1201 * \| GLP_NL \| + \| GLP_NU \|
alpar@9	1202 * \| GLP_NL \| - \| GLP_NL \|
alpar@9	1203 * \| GLP_NU \| + \| GLP_NL \|
alpar@9	1204 * \| GLP_NU \| - \| GLP_NU \|
alpar@9	1205 * \| GLP_NF \| + / - \| GLP_NF \|
alpar@9	1206 * +-----------------------+---------+--------------------+
alpar@9	1207 *
alpar@9	1208 * Value of column q is computed with formula (7). Since originally row
alpar@9	1209 * p is equality constraint, its status is assigned GLP_NS, and value of
alpar@9	1210 * its multiplier pi[p] is computed with formula (11).
alpar@9	1211 *
alpar@9	1212 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	1213 *
alpar@9	1214 * Value of column q is computed with formula (7). Row multiplier value
alpar@9	1215 * pi[p] is computed with formula (11).
alpar@9	1216 *
alpar@9	1217 * RECOVERING MIP SOLUTION
alpar@9	1218 *
alpar@9	1219 * Value of column q is computed with formula (7). */
alpar@9	1220
alpar@9	1221 struct implied_slack
alpar@9	1222 { /* column singleton (implied slack variable) */
alpar@9	1223 int p;
alpar@9	1224 /* row reference number */
alpar@9	1225 int q;
alpar@9	1226 /* column reference number */
alpar@9	1227 double apq;
alpar@9	1228 /* constraint coefficient a[p,q] */
alpar@9	1229 double b;
alpar@9	1230 /* right-hand side of original equality constraint */
alpar@9	1231 double c;
alpar@9	1232 /* original objective coefficient at x[q] */
alpar@9	1233 NPPLFE *ptr;
alpar@9	1234 /* list of non-zero coefficients a[p,j], j != q */
alpar@9	1235 };
alpar@9	1236
alpar@9	1237 static int rcv_implied_slack(NPP npp, void info);
alpar@9	1238
alpar@9	1239 void npp_implied_slack(NPP npp, NPPCOL q)
alpar@9	1240 { /* process column singleton (implied slack variable) */
alpar@9	1241 struct implied_slack *info;
alpar@9	1242 NPPROW *p;
alpar@9	1243 NPPAIJ *aij;
alpar@9	1244 NPPLFE *lfe;
alpar@9	1245 /* the column must be non-integral non-fixed singleton */
alpar@9	1246 xassert(!q->is_int);
alpar@9	1247 xassert(q->lb < q->ub);
alpar@9	1248 xassert(q->ptr != NULL && q->ptr->c_next == NULL);
alpar@9	1249 /* corresponding row must be equality constraint */
alpar@9	1250 aij = q->ptr;
alpar@9	1251 p = aij->row;
alpar@9	1252 xassert(p->lb == p->ub);
alpar@9	1253 /* create transformation stack entry */
alpar@9	1254 info = npp_push_tse(npp,
alpar@9	1255 rcv_implied_slack, sizeof(struct implied_slack));
alpar@9	1256 info->p = p->i;
alpar@9	1257 info->q = q->j;
alpar@9	1258 info->apq = aij->val;
alpar@9	1259 info->b = p->lb;
alpar@9	1260 info->c = q->coef;
alpar@9	1261 info->ptr = NULL;
alpar@9	1262 /* save row coefficients a[p,j], j != q, and substitute x[q]
alpar@9	1263 into the objective row */
alpar@9	1264 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
alpar@9	1265 { if (aij->col == q) continue; /* skip a[p,q] */
alpar@9	1266 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9	1267 lfe->ref = aij->col->j;
alpar@9	1268 lfe->val = aij->val;
alpar@9	1269 lfe->next = info->ptr;
alpar@9	1270 info->ptr = lfe;
alpar@9	1271 aij->col->coef -= info->c * (aij->val / info->apq);
alpar@9	1272 }
alpar@9	1273 npp->c0 += info->c * (info->b / info->apq);
alpar@9	1274 /* compute new row bounds */
alpar@9	1275 if (info->apq > 0.0)
alpar@9	1276 { p->lb = (q->ub == +DBL_MAX ?
alpar@9	1277 -DBL_MAX : info->b - info->apq * q->ub);
alpar@9	1278 p->ub = (q->lb == -DBL_MAX ?
alpar@9	1279 +DBL_MAX : info->b - info->apq * q->lb);
alpar@9	1280 }
alpar@9	1281 else
alpar@9	1282 { p->lb = (q->lb == -DBL_MAX ?
alpar@9	1283 -DBL_MAX : info->b - info->apq * q->lb);
alpar@9	1284 p->ub = (q->ub == +DBL_MAX ?
alpar@9	1285 +DBL_MAX : info->b - info->apq * q->ub);
alpar@9	1286 }
alpar@9	1287 /* remove the column from the problem */
alpar@9	1288 npp_del_col(npp, q);
alpar@9	1289 return;
alpar@9	1290 }
alpar@9	1291
alpar@9	1292 static int rcv_implied_slack(NPP npp, void _info)
alpar@9	1293 { /* recover column singleton (implied slack variable) */
alpar@9	1294 struct implied_slack *info = _info;
alpar@9	1295 NPPLFE *lfe;
alpar@9	1296 double temp;
alpar@9	1297 if (npp->sol == GLP_SOL)
alpar@9	1298 { /* assign statuses to row p and column q */
alpar@9	1299 if (npp->r_stat[info->p] == GLP_BS \|\|
alpar@9	1300 npp->r_stat[info->p] == GLP_NF)
alpar@9	1301 npp->c_stat[info->q] = npp->r_stat[info->p];
alpar@9	1302 else if (npp->r_stat[info->p] == GLP_NL)
alpar@9	1303 npp->c_stat[info->q] =
alpar@9	1304 (char)(info->apq > 0.0 ? GLP_NU : GLP_NL);
alpar@9	1305 else if (npp->r_stat[info->p] == GLP_NU)
alpar@9	1306 npp->c_stat[info->q] =
alpar@9	1307 (char)(info->apq > 0.0 ? GLP_NL : GLP_NU);
alpar@9	1308 else
alpar@9	1309 { npp_error();
alpar@9	1310 return 1;
alpar@9	1311 }
alpar@9	1312 npp->r_stat[info->p] = GLP_NS;
alpar@9	1313 }
alpar@9	1314 if (npp->sol != GLP_MIP)
alpar@9	1315 { /* compute multiplier for row p */
alpar@9	1316 npp->r_pi[info->p] += info->c / info->apq;
alpar@9	1317 }
alpar@9	1318 /* compute value of column q */
alpar@9	1319 temp = info->b;
alpar@9	1320 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
alpar@9	1321 temp -= lfe->val * npp->c_value[lfe->ref];
alpar@9	1322 npp->c_value[info->q] = temp / info->apq;
alpar@9	1323 return 0;
alpar@9	1324 }
alpar@9	1325
alpar@9	1326 /***********************************************************************
alpar@9	1327 * NAME
alpar@9	1328 *
alpar@9	1329 * npp_implied_free - process column singleton (implied free variable)
alpar@9	1330 *
alpar@9	1331 * SYNOPSIS
alpar@9	1332 *
alpar@9	1333 * #include "glpnpp.h"
alpar@9	1334 * int npp_implied_free(NPP npp, NPPCOL q);
alpar@9	1335 *
alpar@9	1336 * DESCRIPTION
alpar@9	1337 *
alpar@9	1338 * The routine npp_implied_free processes column q:
alpar@9	1339 *
alpar@9	1340 * l[q] <= x[q] <= u[q], (1)
alpar@9	1341 *
alpar@9	1342 * having non-zero coefficient in the only row p, which is inequality
alpar@9	1343 * constraint:
alpar@9	1344 *
alpar@9	1345 * L[p] <= sum a[p,j] x[j] + a[p,q] x[q] <= U[p], (2)
alpar@9	1346 * j!=q
alpar@9	1347 *
alpar@9	1348 * where l[q] < u[q], L[p] < U[p], L[p] > -oo and/or U[p] < +oo.
alpar@9	1349 *
alpar@9	1350 * RETURNS
alpar@9	1351 *
alpar@9	1352 * 0 - success;
alpar@9	1353 *
alpar@9	1354 * 1 - column lower and/or upper bound(s) can be active;
alpar@9	1355 *
alpar@9	1356 * 2 - problem has no dual feasible solution.
alpar@9	1357 *
alpar@9	1358 * PROBLEM TRANSFORMATION
alpar@9	1359 *
alpar@9	1360 * Constraint (2) can be written as follows:
alpar@9	1361 *
alpar@9	1362 * L[p] - sum a[p,j] x[j] <= a[p,q] x[q] <= U[p] - sum a[p,j] x[j],
alpar@9	1363 * j!=q j!=q
alpar@9	1364 *
alpar@9	1365 * from which it follows that:
alpar@9	1366 *
alpar@9	1367 * alfa <= a[p,q] x[q] <= beta, (3)
alpar@9	1368 *
alpar@9	1369 * where
alpar@9	1370 *
alpar@9	1371 * alfa = inf(L[p] - sum a[p,j] x[j]) =
alpar@9	1372 * j!=q
alpar@9	1373 *
alpar@9	1374 * = L[p] - sup sum a[p,j] x[j] = (4)
alpar@9	1375 * j!=q
alpar@9	1376 *
alpar@9	1377 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j],
alpar@9	1378 * j in Jp j in Jn
alpar@9	1379 *
alpar@9	1380 * beta = sup(L[p] - sum a[p,j] x[j]) =
alpar@9	1381 * j!=q
alpar@9	1382 *
alpar@9	1383 * = L[p] - inf sum a[p,j] x[j] = (5)
alpar@9	1384 * j!=q
alpar@9	1385 *
alpar@9	1386 * = L[p] - sum a[p,j] l[j] - sum a[p,j] u[j],
alpar@9	1387 * j in Jp j in Jn
alpar@9	1388 *
alpar@9	1389 * Jp = {j != q: a[p,j] > 0}, Jn = {j != q: a[p,j] < 0}. (6)
alpar@9	1390 *
alpar@9	1391 * Inequality (3) defines implied bounds of variable x[q]:
alpar@9	1392 *
alpar@9	1393 * l'[q] <= x[q] <= u'[q], (7)
alpar@9	1394 *
alpar@9	1395 * where
alpar@9	1396 *
alpar@9	1397 * ( alfa / a[p,q], if a[p,q] > 0
alpar@9	1398 * l'[q] = < (8a)
alpar@9	1399 * ( beta / a[p,q], if a[p,q] < 0
alpar@9	1400 *
alpar@9	1401 * ( beta / a[p,q], if a[p,q] > 0
alpar@9	1402 * u'[q] = < (8b)
alpar@9	1403 * ( alfa / a[p,q], if a[p,q] < 0
alpar@9	1404 *
alpar@9	1405 * Thus, if l'[q] > l[q] - eps and u'[q] < u[q] + eps, where eps is
alpar@9	1406 * an absolute tolerance for column value, column bounds (1) cannot be
alpar@9	1407 * active, in which case column q can be replaced by equivalent free
alpar@9	1408 * (unbounded) column.
alpar@9	1409 *
alpar@9	1410 * Note that column q is column singleton, so in the dual system of the
alpar@9	1411 * original problem it corresponds to the following row singleton:
alpar@9	1412 *
alpar@9	1413 * a[p,q] pi[p] + lambda[q] = c[q], (9)
alpar@9	1414 *
alpar@9	1415 * from which it follows that:
alpar@9	1416 *
alpar@9	1417 * pi[p] = (c[q] - lambda[q]) / a[p,q]. (10)
alpar@9	1418 *
alpar@9	1419 * Let x[q] be implied free (unbounded) variable. Then column q can be
alpar@9	1420 * only basic, so its multiplier lambda[q] is equal to zero, and from
alpar@9	1421 * (10) we have:
alpar@9	1422 *
alpar@9	1423 * pi[p] = c[q] / a[p,q]. (11)
alpar@9	1424 *
alpar@9	1425 * There are possible three cases:
alpar@9	1426 *
alpar@9	1427 * 1) pi[p] < -eps, where eps is an absolute tolerance for row
alpar@9	1428 * multiplier. In this case, to provide dual feasibility of the
alpar@9	1429 * original problem, row p must be active on its lower bound, and
alpar@9	1430 * if its lower bound does not exist (L[p] = -oo), the problem has
alpar@9	1431 * no dual feasible solution;
alpar@9	1432 *
alpar@9	1433 * 2) pi[p] > +eps. In this case row p must be active on its upper
alpar@9	1434 * bound, and if its upper bound does not exist (U[p] = +oo), the
alpar@9	1435 * problem has no dual feasible solution;
alpar@9	1436 *
alpar@9	1437 * 3) -eps <= pi[p] <= +eps. In this case any (either lower or upper)
alpar@9	1438 * bound of row p can be active, because this does not affect dual
alpar@9	1439 * feasibility.
alpar@9	1440 *
alpar@9	1441 * Thus, in all three cases original inequality constraint (2) can be
alpar@9	1442 * replaced by equality constraint, where the right-hand side is either
alpar@9	1443 * lower or upper bound of row p, and bounds of column q can be removed
alpar@9	1444 * that makes it free (unbounded). (May note that this transformation
alpar@9	1445 * can be followed by transformation "Column singleton (implied slack
alpar@9	1446 * variable)" performed by the routine npp_implied_slack.)
alpar@9	1447 *
alpar@9	1448 * RECOVERING BASIC SOLUTION
alpar@9	1449 *
alpar@9	1450 * Status of row p in solution to the original problem is determined
alpar@9	1451 * by its status in solution to the transformed problem and its bound,
alpar@9	1452 * which was choosen to be active:
alpar@9	1453 *
alpar@9	1454 * +-----------------------+--------+--------------------+
alpar@9	1455 * \| Status of row p \| Active \| Status of row p \|
alpar@9	1456 * \| (transformed problem) \| bound \| (original problem) \|
alpar@9	1457 * +-----------------------+--------+--------------------+
alpar@9	1458 * \| GLP_BS \| L[p] \| GLP_BS \|
alpar@9	1459 * \| GLP_BS \| U[p] \| GLP_BS \|
alpar@9	1460 * \| GLP_NS \| L[p] \| GLP_NL \|
alpar@9	1461 * \| GLP_NS \| U[p] \| GLP_NU \|
alpar@9	1462 * +-----------------------+--------+--------------------+
alpar@9	1463 *
alpar@9	1464 * Value of row multiplier pi[p] (as well as value of column q) in
alpar@9	1465 * solution to the original problem is the same as in solution to the
alpar@9	1466 * transformed problem.
alpar@9	1467 *
alpar@9	1468 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	1469 *
alpar@9	1470 * Value of row multiplier pi[p] in solution to the original problem is
alpar@9	1471 * the same as in solution to the transformed problem.
alpar@9	1472 *
alpar@9	1473 * RECOVERING MIP SOLUTION
alpar@9	1474 *
alpar@9	1475 * None needed. */
alpar@9	1476
alpar@9	1477 struct implied_free
alpar@9	1478 { /* column singleton (implied free variable) */
alpar@9	1479 int p;
alpar@9	1480 /* row reference number */
alpar@9	1481 char stat;
alpar@9	1482 /* row status:
alpar@9	1483 GLP_NL - active constraint on lower bound
alpar@9	1484 GLP_NU - active constraint on upper bound */
alpar@9	1485 };
alpar@9	1486
alpar@9	1487 static int rcv_implied_free(NPP npp, void info);
alpar@9	1488
alpar@9	1489 int npp_implied_free(NPP npp, NPPCOL q)
alpar@9	1490 { /* process column singleton (implied free variable) */
alpar@9	1491 struct implied_free *info;
alpar@9	1492 NPPROW *p;
alpar@9	1493 NPPAIJ apq, aij;
alpar@9	1494 double alfa, beta, l, u, pi, eps;
alpar@9	1495 /* the column must be non-fixed singleton */
alpar@9	1496 xassert(q->lb < q->ub);
alpar@9	1497 xassert(q->ptr != NULL && q->ptr->c_next == NULL);
alpar@9	1498 /* corresponding row must be inequality constraint */
alpar@9	1499 apq = q->ptr;
alpar@9	1500 p = apq->row;
alpar@9	1501 xassert(p->lb != -DBL_MAX \|\| p->ub != +DBL_MAX);
alpar@9	1502 xassert(p->lb < p->ub);
alpar@9	1503 /* compute alfa */
alpar@9	1504 alfa = p->lb;
alpar@9	1505 if (alfa != -DBL_MAX)
alpar@9	1506 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
alpar@9	1507 { if (aij == apq) continue; /* skip a[p,q] */
alpar@9	1508 if (aij->val > 0.0)
alpar@9	1509 { if (aij->col->ub == +DBL_MAX)
alpar@9	1510 { alfa = -DBL_MAX;
alpar@9	1511 break;
alpar@9	1512 }
alpar@9	1513 alfa -= aij->val * aij->col->ub;
alpar@9	1514 }
alpar@9	1515 else /* < 0.0 */
alpar@9	1516 { if (aij->col->lb == -DBL_MAX)
alpar@9	1517 { alfa = -DBL_MAX;
alpar@9	1518 break;
alpar@9	1519 }
alpar@9	1520 alfa -= aij->val * aij->col->lb;
alpar@9	1521 }
alpar@9	1522 }
alpar@9	1523 }
alpar@9	1524 /* compute beta */
alpar@9	1525 beta = p->ub;
alpar@9	1526 if (beta != +DBL_MAX)
alpar@9	1527 { for (aij = p->ptr; aij != NULL; aij = aij->r_next)
alpar@9	1528 { if (aij == apq) continue; /* skip a[p,q] */
alpar@9	1529 if (aij->val > 0.0)
alpar@9	1530 { if (aij->col->lb == -DBL_MAX)
alpar@9	1531 { beta = +DBL_MAX;
alpar@9	1532 break;
alpar@9	1533 }
alpar@9	1534 beta -= aij->val * aij->col->lb;
alpar@9	1535 }
alpar@9	1536 else /* < 0.0 */
alpar@9	1537 { if (aij->col->ub == +DBL_MAX)
alpar@9	1538 { beta = +DBL_MAX;
alpar@9	1539 break;
alpar@9	1540 }
alpar@9	1541 beta -= aij->val * aij->col->ub;
alpar@9	1542 }
alpar@9	1543 }
alpar@9	1544 }
alpar@9	1545 /* compute implied column lower bound l'[q] */
alpar@9	1546 if (apq->val > 0.0)
alpar@9	1547 l = (alfa == -DBL_MAX ? -DBL_MAX : alfa / apq->val);
alpar@9	1548 else /* < 0.0 */
alpar@9	1549 l = (beta == +DBL_MAX ? -DBL_MAX : beta / apq->val);
alpar@9	1550 /* compute implied column upper bound u'[q] */
alpar@9	1551 if (apq->val > 0.0)
alpar@9	1552 u = (beta == +DBL_MAX ? +DBL_MAX : beta / apq->val);
alpar@9	1553 else
alpar@9	1554 u = (alfa == -DBL_MAX ? +DBL_MAX : alfa / apq->val);
alpar@9	1555 /* check if column lower bound l[q] can be active */
alpar@9	1556 if (q->lb != -DBL_MAX)
alpar@9	1557 { eps = 1e-9 + 1e-12 * fabs(q->lb);
alpar@9	1558 if (l < q->lb - eps) return 1; /* yes, it can */
alpar@9	1559 }
alpar@9	1560 /* check if column upper bound u[q] can be active */
alpar@9	1561 if (q->ub != +DBL_MAX)
alpar@9	1562 { eps = 1e-9 + 1e-12 * fabs(q->ub);
alpar@9	1563 if (u > q->ub + eps) return 1; /* yes, it can */
alpar@9	1564 }
alpar@9	1565 /* okay; make column q free (unbounded) */
alpar@9	1566 q->lb = -DBL_MAX, q->ub = +DBL_MAX;
alpar@9	1567 /* create transformation stack entry */
alpar@9	1568 info = npp_push_tse(npp,
alpar@9	1569 rcv_implied_free, sizeof(struct implied_free));
alpar@9	1570 info->p = p->i;
alpar@9	1571 info->stat = -1;
alpar@9	1572 /* compute row multiplier pi[p] */
alpar@9	1573 pi = q->coef / apq->val;
alpar@9	1574 /* check dual feasibility for row p */
alpar@9	1575 if (pi > +DBL_EPSILON)
alpar@9	1576 { /* lower bound L[p] must be active */
alpar@9	1577 if (p->lb != -DBL_MAX)
alpar@9	1578 nl: { info->stat = GLP_NL;
alpar@9	1579 p->ub = p->lb;
alpar@9	1580 }
alpar@9	1581 else
alpar@9	1582 { if (pi > +1e-5) return 2; /* dual infeasibility */
alpar@9	1583 /* take a chance on U[p] */
alpar@9	1584 xassert(p->ub != +DBL_MAX);
alpar@9	1585 goto nu;
alpar@9	1586 }
alpar@9	1587 }
alpar@9	1588 else if (pi < -DBL_EPSILON)
alpar@9	1589 { /* upper bound U[p] must be active */
alpar@9	1590 if (p->ub != +DBL_MAX)
alpar@9	1591 nu: { info->stat = GLP_NU;
alpar@9	1592 p->lb = p->ub;
alpar@9	1593 }
alpar@9	1594 else
alpar@9	1595 { if (pi < -1e-5) return 2; /* dual infeasibility */
alpar@9	1596 /* take a chance on L[p] */
alpar@9	1597 xassert(p->lb != -DBL_MAX);
alpar@9	1598 goto nl;
alpar@9	1599 }
alpar@9	1600 }
alpar@9	1601 else
alpar@9	1602 { /* any bound (either L[p] or U[p]) can be made active */
alpar@9	1603 if (p->ub == +DBL_MAX)
alpar@9	1604 { xassert(p->lb != -DBL_MAX);
alpar@9	1605 goto nl;
alpar@9	1606 }
alpar@9	1607 if (p->lb == -DBL_MAX)
alpar@9	1608 { xassert(p->ub != +DBL_MAX);
alpar@9	1609 goto nu;
alpar@9	1610 }
alpar@9	1611 if (fabs(p->lb) <= fabs(p->ub)) goto nl; else goto nu;
alpar@9	1612 }
alpar@9	1613 return 0;
alpar@9	1614 }
alpar@9	1615
alpar@9	1616 static int rcv_implied_free(NPP npp, void _info)
alpar@9	1617 { /* recover column singleton (implied free variable) */
alpar@9	1618 struct implied_free *info = _info;
alpar@9	1619 if (npp->sol == GLP_SOL)
alpar@9	1620 { if (npp->r_stat[info->p] == GLP_BS)
alpar@9	1621 npp->r_stat[info->p] = GLP_BS;
alpar@9	1622 else if (npp->r_stat[info->p] == GLP_NS)
alpar@9	1623 { xassert(info->stat == GLP_NL \|\| info->stat == GLP_NU);
alpar@9	1624 npp->r_stat[info->p] = info->stat;
alpar@9	1625 }
alpar@9	1626 else
alpar@9	1627 { npp_error();
alpar@9	1628 return 1;
alpar@9	1629 }
alpar@9	1630 }
alpar@9	1631 return 0;
alpar@9	1632 }
alpar@9	1633
alpar@9	1634 /***********************************************************************
alpar@9	1635 * NAME
alpar@9	1636 *
alpar@9	1637 * npp_eq_doublet - process row doubleton (equality constraint)
alpar@9	1638 *
alpar@9	1639 * SYNOPSIS
alpar@9	1640 *
alpar@9	1641 * #include "glpnpp.h"
alpar@9	1642 * NPPCOL npp_eq_doublet(NPP npp, NPPROW *p);
alpar@9	1643 *
alpar@9	1644 * DESCRIPTION
alpar@9	1645 *
alpar@9	1646 * The routine npp_eq_doublet processes row p, which is equality
alpar@9	1647 * constraint having exactly two non-zero coefficients:
alpar@9	1648 *
alpar@9	1649 * a[p,q] x[q] + a[p,r] x[r] = b. (1)
alpar@9	1650 *
alpar@9	1651 * As the result of processing one of columns q or r is eliminated from
alpar@9	1652 * all other rows and, thus, becomes column singleton of type "implied
alpar@9	1653 * slack variable". Row p is not changed and along with column q and r
alpar@9	1654 * remains in the problem.
alpar@9	1655 *
alpar@9	1656 * RETURNS
alpar@9	1657 *
alpar@9	1658 * The routine npp_eq_doublet returns pointer to the descriptor of that
alpar@9	1659 * column q or r which has been eliminated. If, due to some reason, the
alpar@9	1660 * elimination was not performed, the routine returns NULL.
alpar@9	1661 *
alpar@9	1662 * PROBLEM TRANSFORMATION
alpar@9	1663 *
alpar@9	1664 * First, we decide which column q or r will be eliminated. Let it be
alpar@9	1665 * column q. Consider i-th constraint row, where column q has non-zero
alpar@9	1666 * coefficient a[i,q] != 0:
alpar@9	1667 *
alpar@9	1668 * L[i] <= sum a[i,j] x[j] <= U[i]. (2)
alpar@9	1669 * j
alpar@9	1670 *
alpar@9	1671 * In order to eliminate column q from row (2) we subtract from it row
alpar@9	1672 * (1) multiplied by gamma[i] = a[i,q] / a[p,q], i.e. we replace in the
alpar@9	1673 * transformed problem row (2) by its linear combination with row (1).
alpar@9	1674 * This transformation changes only coefficients in columns q and r,
alpar@9	1675 * and bounds of row i as follows:
alpar@9	1676 *
alpar@9	1677 * a~[i,q] = a[i,q] - gamma[i] a[p,q] = 0, (3)
alpar@9	1678 *
alpar@9	1679 * a~[i,r] = a[i,r] - gamma[i] a[p,r], (4)
alpar@9	1680 *
alpar@9	1681 * L~[i] = L[i] - gamma[i] b, (5)
alpar@9	1682 *
alpar@9	1683 * U~[i] = U[i] - gamma[i] b. (6)
alpar@9	1684 *
alpar@9	1685 * RECOVERING BASIC SOLUTION
alpar@9	1686 *
alpar@9	1687 * The transformation of the primal system of the original problem:
alpar@9	1688 *
alpar@9	1689 * L <= A x <= U (7)
alpar@9	1690 *
alpar@9	1691 * is equivalent to multiplying from the left a transformation matrix F
alpar@9	1692 * by components of this primal system, which in the transformed problem
alpar@9	1693 * becomes the following:
alpar@9	1694 *
alpar@9	1695 * F L <= F A x <= F U ==> L~ <= A~x <= U~. (8)
alpar@9	1696 *
alpar@9	1697 * The matrix F has the following structure:
alpar@9	1698 *
alpar@9	1699 * ( 1 -gamma[1] )
alpar@9	1700 * ( )
alpar@9	1701 * ( 1 -gamma[2] )
alpar@9	1702 * ( )
alpar@9	1703 * ( ... ... )
alpar@9	1704 * ( )
alpar@9	1705 * F = ( 1 -gamma[p-1] ) (9)
alpar@9	1706 * ( )
alpar@9	1707 * ( 1 )
alpar@9	1708 * ( )
alpar@9	1709 * ( -gamma[p+1] 1 )
alpar@9	1710 * ( )
alpar@9	1711 * ( ... ... )
alpar@9	1712 *
alpar@9	1713 * where its column containing elements -gamma[i] corresponds to row p
alpar@9	1714 * of the primal system.
alpar@9	1715 *
alpar@9	1716 * From (8) it follows that the dual system of the original problem:
alpar@9	1717 *
alpar@9	1718 * A'pi + lambda = c, (10)
alpar@9	1719 *
alpar@9	1720 * in the transformed problem becomes the following:
alpar@9	1721 *
alpar@9	1722 * A'F'inv(F')pi + lambda = c ==> (A~)'pi~ + lambda = c, (11)
alpar@9	1723 *
alpar@9	1724 * where:
alpar@9	1725 *
alpar@9	1726 * pi~ = inv(F')pi (12)
alpar@9	1727 *
alpar@9	1728 * is the vector of row multipliers in the transformed problem. Thus:
alpar@9	1729 *
alpar@9	1730 * pi = F'pi~. (13)
alpar@9	1731 *
alpar@9	1732 * Therefore, as it follows from (13), value of multiplier for row p in
alpar@9	1733 * solution to the original problem can be computed as follows:
alpar@9	1734 *
alpar@9	1735 * pi[p] = pi~[p] - sum gamma[i] pi~[i], (14)
alpar@9	1736 * i
alpar@9	1737 *
alpar@9	1738 * where pi~[i] = pi[i] is multiplier for row i (i != p).
alpar@9	1739 *
alpar@9	1740 * Note that the statuses of all rows and columns are not changed.
alpar@9	1741 *
alpar@9	1742 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	1743 *
alpar@9	1744 * Multiplier for row p in solution to the original problem is computed
alpar@9	1745 * with formula (14).
alpar@9	1746 *
alpar@9	1747 * RECOVERING MIP SOLUTION
alpar@9	1748 *
alpar@9	1749 * None needed. */
alpar@9	1750
alpar@9	1751 struct eq_doublet
alpar@9	1752 { /* row doubleton (equality constraint) */
alpar@9	1753 int p;
alpar@9	1754 /* row reference number */
alpar@9	1755 double apq;
alpar@9	1756 /* constraint coefficient a[p,q] */
alpar@9	1757 NPPLFE *ptr;
alpar@9	1758 /* list of non-zero coefficients a[i,q], i != p */
alpar@9	1759 };
alpar@9	1760
alpar@9	1761 static int rcv_eq_doublet(NPP npp, void info);
alpar@9	1762
alpar@9	1763 NPPCOL npp_eq_doublet(NPP npp, NPPROW *p)
alpar@9	1764 { /* process row doubleton (equality constraint) */
alpar@9	1765 struct eq_doublet *info;
alpar@9	1766 NPPROW *i;
alpar@9	1767 NPPCOL q, r;
alpar@9	1768 NPPAIJ apq, apr, aiq, air, *next;
alpar@9	1769 NPPLFE *lfe;
alpar@9	1770 double gamma;
alpar@9	1771 /* the row must be doubleton equality constraint */
alpar@9	1772 xassert(p->lb == p->ub);
alpar@9	1773 xassert(p->ptr != NULL && p->ptr->r_next != NULL &&
alpar@9	1774 p->ptr->r_next->r_next == NULL);
alpar@9	1775 /* choose column to be eliminated */
alpar@9	1776 { NPPAIJ a1, a2;
alpar@9	1777 a1 = p->ptr, a2 = a1->r_next;
alpar@9	1778 if (fabs(a2->val) < 0.001 * fabs(a1->val))
alpar@9	1779 { /* only first column can be eliminated, because second one
alpar@9	1780 has too small constraint coefficient */
alpar@9	1781 apq = a1, apr = a2;
alpar@9	1782 }
alpar@9	1783 else if (fabs(a1->val) < 0.001 * fabs(a2->val))
alpar@9	1784 { /* only second column can be eliminated, because first one
alpar@9	1785 has too small constraint coefficient */
alpar@9	1786 apq = a2, apr = a1;
alpar@9	1787 }
alpar@9	1788 else
alpar@9	1789 { /* both columns are appropriate; choose that one which is
alpar@9	1790 shorter to minimize fill-in */
alpar@9	1791 if (npp_col_nnz(npp, a1->col) <= npp_col_nnz(npp, a2->col))
alpar@9	1792 { /* first column is shorter */
alpar@9	1793 apq = a1, apr = a2;
alpar@9	1794 }
alpar@9	1795 else
alpar@9	1796 { /* second column is shorter */
alpar@9	1797 apq = a2, apr = a1;
alpar@9	1798 }
alpar@9	1799 }
alpar@9	1800 }
alpar@9	1801 /* now columns q and r have been chosen */
alpar@9	1802 q = apq->col, r = apr->col;
alpar@9	1803 /* create transformation stack entry */
alpar@9	1804 info = npp_push_tse(npp,
alpar@9	1805 rcv_eq_doublet, sizeof(struct eq_doublet));
alpar@9	1806 info->p = p->i;
alpar@9	1807 info->apq = apq->val;
alpar@9	1808 info->ptr = NULL;
alpar@9	1809 /* transform each row i (i != p), where a[i,q] != 0, to eliminate
alpar@9	1810 column q */
alpar@9	1811 for (aiq = q->ptr; aiq != NULL; aiq = next)
alpar@9	1812 { next = aiq->c_next;
alpar@9	1813 if (aiq == apq) continue; /* skip row p */
alpar@9	1814 i = aiq->row; /* row i to be transformed */
alpar@9	1815 /* save constraint coefficient a[i,q] */
alpar@9	1816 if (npp->sol != GLP_MIP)
alpar@9	1817 { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9	1818 lfe->ref = i->i;
alpar@9	1819 lfe->val = aiq->val;
alpar@9	1820 lfe->next = info->ptr;
alpar@9	1821 info->ptr = lfe;
alpar@9	1822 }
alpar@9	1823 /* find coefficient a[i,r] in row i */
alpar@9	1824 for (air = i->ptr; air != NULL; air = air->r_next)
alpar@9	1825 if (air->col == r) break;
alpar@9	1826 /* if a[i,r] does not exist, create a[i,r] = 0 */
alpar@9	1827 if (air == NULL)
alpar@9	1828 air = npp_add_aij(npp, i, r, 0.0);
alpar@9	1829 /* compute gamma[i] = a[i,q] / a[p,q] */
alpar@9	1830 gamma = aiq->val / apq->val;
alpar@9	1831 /* (row i) := (row i) - gamma[i] * (row p); see (3)-(6) */
alpar@9	1832 /* new a[i,q] is exact zero due to elimnation; remove it from
alpar@9	1833 row i */
alpar@9	1834 npp_del_aij(npp, aiq);
alpar@9	1835 /* compute new a[i,r] */
alpar@9	1836 air->val -= gamma * apr->val;
alpar@9	1837 /* if new a[i,r] is close to zero due to numeric cancelation,
alpar@9	1838 remove it from row i */
alpar@9	1839 if (fabs(air->val) <= 1e-10)
alpar@9	1840 npp_del_aij(npp, air);
alpar@9	1841 /* compute new lower and upper bounds of row i */
alpar@9	1842 if (i->lb == i->ub)
alpar@9	1843 i->lb = i->ub = (i->lb - gamma * p->lb);
alpar@9	1844 else
alpar@9	1845 { if (i->lb != -DBL_MAX)
alpar@9	1846 i->lb -= gamma * p->lb;
alpar@9	1847 if (i->ub != +DBL_MAX)
alpar@9	1848 i->ub -= gamma * p->lb;
alpar@9	1849 }
alpar@9	1850 }
alpar@9	1851 return q;
alpar@9	1852 }
alpar@9	1853
alpar@9	1854 static int rcv_eq_doublet(NPP npp, void _info)
alpar@9	1855 { /* recover row doubleton (equality constraint) */
alpar@9	1856 struct eq_doublet *info = _info;
alpar@9	1857 NPPLFE *lfe;
alpar@9	1858 double gamma, temp;
alpar@9	1859 /* we assume that processing row p is followed by processing
alpar@9	1860 column q as singleton of type "implied slack variable", in
alpar@9	1861 which case row p must always be active equality constraint */
alpar@9	1862 if (npp->sol == GLP_SOL)
alpar@9	1863 { if (npp->r_stat[info->p] != GLP_NS)
alpar@9	1864 { npp_error();
alpar@9	1865 return 1;
alpar@9	1866 }
alpar@9	1867 }
alpar@9	1868 if (npp->sol != GLP_MIP)
alpar@9	1869 { /* compute value of multiplier for row p; see (14) */
alpar@9	1870 temp = npp->r_pi[info->p];
alpar@9	1871 for (lfe = info->ptr; lfe != NULL; lfe = lfe->next)
alpar@9	1872 { gamma = lfe->val / info->apq; /* a[i,q] / a[p,q] */
alpar@9	1873 temp -= gamma * npp->r_pi[lfe->ref];
alpar@9	1874 }
alpar@9	1875 npp->r_pi[info->p] = temp;
alpar@9	1876 }
alpar@9	1877 return 0;
alpar@9	1878 }
alpar@9	1879
alpar@9	1880 /***********************************************************************
alpar@9	1881 * NAME
alpar@9	1882 *
alpar@9	1883 * npp_forcing_row - process forcing row
alpar@9	1884 *
alpar@9	1885 * SYNOPSIS
alpar@9	1886 *
alpar@9	1887 * #include "glpnpp.h"
alpar@9	1888 * int npp_forcing_row(NPP npp, NPPROW p, int at);
alpar@9	1889 *
alpar@9	1890 * DESCRIPTION
alpar@9	1891 *
alpar@9	1892 * The routine npp_forcing row processes row p of general format:
alpar@9	1893 *
alpar@9	1894 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
alpar@9	1895 * j
alpar@9	1896 *
alpar@9	1897 * l[j] <= x[j] <= u[j], (2)
alpar@9	1898 *
alpar@9	1899 * where L[p] <= U[p] and l[j] < u[j] for all a[p,j] != 0. It is also
alpar@9	1900 * assumed that:
alpar@9	1901 *
alpar@9	1902 * 1) if at = 0 then \|L[p] - U'[p]\| <= eps, where U'[p] is implied
alpar@9	1903 * row upper bound (see below), eps is an absolute tolerance for row
alpar@9	1904 * value;
alpar@9	1905 *
alpar@9	1906 * 2) if at = 1 then \|U[p] - L'[p]\| <= eps, where L'[p] is implied
alpar@9	1907 * row lower bound (see below).
alpar@9	1908 *
alpar@9	1909 * RETURNS
alpar@9	1910 *
alpar@9	1911 * 0 - success;
alpar@9	1912 *
alpar@9	1913 * 1 - cannot fix columns due to too small constraint coefficients.
alpar@9	1914 *
alpar@9	1915 * PROBLEM TRANSFORMATION
alpar@9	1916 *
alpar@9	1917 * Implied lower and upper bounds of row (1) are determined by bounds
alpar@9	1918 * of corresponding columns (variables) as follows:
alpar@9	1919 *
alpar@9	1920 * L'[p] = inf sum a[p,j] x[j] =
alpar@9	1921 * j
alpar@9	1922 * (3)
alpar@9	1923 * = sum a[p,j] l[j] + sum a[p,j] u[j],
alpar@9	1924 * j in Jp j in Jn
alpar@9	1925 *
alpar@9	1926 * U'[p] = sup sum a[p,j] x[j] =
alpar@9	1927 * (4)
alpar@9	1928 * = sum a[p,j] u[j] + sum a[p,j] l[j],
alpar@9	1929 * j in Jp j in Jn
alpar@9	1930 *
alpar@9	1931 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
alpar@9	1932 *
alpar@9	1933 * If L[p] =~ U'[p] (at = 0), solution can be primal feasible only when
alpar@9	1934 * all variables take their boundary values as defined by (4):
alpar@9	1935 *
alpar@9	1936 * ( u[j], if j in Jp
alpar@9	1937 * x[j] = < (6)
alpar@9	1938 * ( l[j], if j in Jn
alpar@9	1939 *
alpar@9	1940 * Similarly, if U[p] =~ L'[p] (at = 1), solution can be primal feasible
alpar@9	1941 * only when all variables take their boundary values as defined by (3):
alpar@9	1942 *
alpar@9	1943 * ( l[j], if j in Jp
alpar@9	1944 * x[j] = < (7)
alpar@9	1945 * ( u[j], if j in Jn
alpar@9	1946 *
alpar@9	1947 * Condition (6) or (7) allows fixing all columns (variables x[j])
alpar@9	1948 * in row (1) on their bounds and then removing them from the problem
alpar@9	1949 * (see the routine npp_fixed_col). Due to this row p becomes redundant,
alpar@9	1950 * so it can be replaced by equivalent free (unbounded) row and also
alpar@9	1951 * removed from the problem (see the routine npp_free_row).
alpar@9	1952 *
alpar@9	1953 * 1. To apply this transformation row (1) should not have coefficients
alpar@9	1954 * whose magnitude is too small, i.e. all a[p,j] should satisfy to
alpar@9	1955 * the following condition:
alpar@9	1956 *
alpar@9	1957 * \|a[p,j]\| >= eps * max(1, \|a[p,k]\|), (8)
alpar@9	1958 * k
alpar@9	1959 * where eps is a relative tolerance for constraint coefficients.
alpar@9	1960 * Otherwise, fixing columns may be numerically unreliable and may
alpar@9	1961 * lead to wrong solution.
alpar@9	1962 *
alpar@9	1963 * 2. The routine fixes columns and remove bounds of row p, however,
alpar@9	1964 * it does not remove the row and columns from the problem.
alpar@9	1965 *
alpar@9	1966 * RECOVERING BASIC SOLUTION
alpar@9	1967 *
alpar@9	1968 * In the transformed problem row p being inactive constraint is
alpar@9	1969 * assigned status GLP_BS (as the result of transformation of free
alpar@9	1970 * row), and all columns in this row are assigned status GLP_NS (as the
alpar@9	1971 * result of transformation of fixed columns).
alpar@9	1972 *
alpar@9	1973 * Note that in the dual system of the transformed (as well as original)
alpar@9	1974 * problem every column j in row p corresponds to the following row:
alpar@9	1975 *
alpar@9	1976 * sum a[i,j] pi[i] + a[p,j] pi[p] + lambda[j] = c[j], (9)
alpar@9	1977 * i!=p
alpar@9	1978 *
alpar@9	1979 * from which it follows that:
alpar@9	1980 *
alpar@9	1981 * lambda[j] = c[j] - sum a[i,j] pi[i] - a[p,j] pi[p]. (10)
alpar@9	1982 * i!=p
alpar@9	1983 *
alpar@9	1984 * In the transformed problem values of all multipliers pi[i] are known
alpar@9	1985 * (including pi[i], whose value is zero, since row p is inactive).
alpar@9	1986 * Thus, using formula (10) it is possible to compute values of
alpar@9	1987 * multipliers lambda[j] for all columns in row p.
alpar@9	1988 *
alpar@9	1989 * Note also that in the original problem all columns in row p are
alpar@9	1990 * bounded, not fixed. So status GLP_NS assigned to every such column
alpar@9	1991 * must be changed to GLP_NL or GLP_NU depending on which bound the
alpar@9	1992 * corresponding column has been fixed. This status change may lead to
alpar@9	1993 * dual feasibility violation for solution of the original problem,
alpar@9	1994 * because now column multipliers must satisfy to the following
alpar@9	1995 * condition:
alpar@9	1996 *
alpar@9	1997 * ( >= 0, if status of column j is GLP_NL,
alpar@9	1998 * lambda[j] < (11)
alpar@9	1999 * ( <= 0, if status of column j is GLP_NU.
alpar@9	2000 *
alpar@9	2001 * If this condition holds, solution to the original problem is the
alpar@9	2002 * same as to the transformed problem. Otherwise, we have to perform
alpar@9	2003 * one degenerate pivoting step of the primal simplex method to obtain
alpar@9	2004 * dual feasible (hence, optimal) solution to the original problem as
alpar@9	2005 * follows. If, on problem transformation, row p was made active on its
alpar@9	2006 * lower bound (case at = 0), we change its status to GLP_NL (or GLP_NS)
alpar@9	2007 * and start increasing its multiplier pi[p]. Otherwise, if row p was
alpar@9	2008 * made active on its upper bound (case at = 1), we change its status
alpar@9	2009 * to GLP_NU (or GLP_NS) and start decreasing pi[p]. From (10) it
alpar@9	2010 * follows that:
alpar@9	2011 *
alpar@9	2012 * delta lambda[j] = - a[p,j] * delta pi[p] = - a[p,j] pi[p]. (12)
alpar@9	2013 *
alpar@9	2014 * Simple analysis of formulae (3)-(5) shows that changing pi[p] in the
alpar@9	2015 * specified direction causes increasing lambda[j] for every column j
alpar@9	2016 * assigned status GLP_NL (delta lambda[j] > 0) and decreasing lambda[j]
alpar@9	2017 * for every column j assigned status GLP_NU (delta lambda[j] < 0). It
alpar@9	2018 * is understood that once the last lambda[q], which violates condition
alpar@9	2019 * (11), has reached zero, multipliers lambda[j] for all columns get
alpar@9	2020 * valid signs. Such column q can be determined as follows. Let d[j] be
alpar@9	2021 * initial value of lambda[j] (i.e. reduced cost of column j) in the
alpar@9	2022 * transformed problem computed with formula (10) when pi[p] = 0. Then
alpar@9	2023 * lambda[j] = d[j] + delta lambda[j], and from (12) it follows that
alpar@9	2024 * lambda[j] becomes zero if:
alpar@9	2025 *
alpar@9	2026 * delta lambda[j] = - a[p,j] pi[p] = - d[j] ==>
alpar@9	2027 * (13)
alpar@9	2028 * pi[p] = d[j] / a[p,j].
alpar@9	2029 *
alpar@9	2030 * Therefore, the last column q, for which lambda[q] becomes zero, can
alpar@9	2031 * be determined from the following condition:
alpar@9	2032 *
alpar@9	2033 * \|d[q] / a[p,q]\| = max \|pi[p]\| = max \|d[j] / a[p,j]\|, (14)
alpar@9	2034 * j in D j in D
alpar@9	2035 *
alpar@9	2036 * where D is a set of columns j whose, reduced costs d[j] have invalid
alpar@9	2037 * signs, i.e. violate condition (11). (Thus, if D is empty, solution
alpar@9	2038 * to the original problem is the same as solution to the transformed
alpar@9	2039 * problem, and no correction is needed as was noticed above.) In
alpar@9	2040 * solution to the original problem column q is assigned status GLP_BS,
alpar@9	2041 * since it replaces column of auxiliary variable of row p (becoming
alpar@9	2042 * active) in the basis, and multiplier for row p is assigned its new
alpar@9	2043 * value, which is pi[p] = d[q] / a[p,q]. Note that due to primal
alpar@9	2044 * degeneracy values of all columns having non-zero coefficients in row
alpar@9	2045 * p remain unchanged.
alpar@9	2046 *
alpar@9	2047 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	2048 *
alpar@9	2049 * Value of multiplier pi[p] in solution to the original problem is
alpar@9	2050 * corrected in the same way as for basic solution. Values of all
alpar@9	2051 * columns having non-zero coefficients in row p remain unchanged.
alpar@9	2052 *
alpar@9	2053 * RECOVERING MIP SOLUTION
alpar@9	2054 *
alpar@9	2055 * None needed. */
alpar@9	2056
alpar@9	2057 struct forcing_col
alpar@9	2058 { /* column fixed on its bound by forcing row */
alpar@9	2059 int j;
alpar@9	2060 /* column reference number */
alpar@9	2061 char stat;
alpar@9	2062 /* original column status:
alpar@9	2063 GLP_NL - fixed on lower bound
alpar@9	2064 GLP_NU - fixed on upper bound */
alpar@9	2065 double a;
alpar@9	2066 /* constraint coefficient a[p,j] */
alpar@9	2067 double c;
alpar@9	2068 /* objective coefficient c[j] */
alpar@9	2069 NPPLFE *ptr;
alpar@9	2070 /* list of non-zero coefficients a[i,j], i != p */
alpar@9	2071 struct forcing_col *next;
alpar@9	2072 /* pointer to another column fixed by forcing row */
alpar@9	2073 };
alpar@9	2074
alpar@9	2075 struct forcing_row
alpar@9	2076 { /* forcing row */
alpar@9	2077 int p;
alpar@9	2078 /* row reference number */
alpar@9	2079 char stat;
alpar@9	2080 /* status assigned to the row if it becomes active:
alpar@9	2081 GLP_NS - active equality constraint
alpar@9	2082 GLP_NL - inequality constraint with lower bound active
alpar@9	2083 GLP_NU - inequality constraint with upper bound active */
alpar@9	2084 struct forcing_col *ptr;
alpar@9	2085 /* list of all columns having non-zero constraint coefficient
alpar@9	2086 a[p,j] in the forcing row */
alpar@9	2087 };
alpar@9	2088
alpar@9	2089 static int rcv_forcing_row(NPP npp, void info);
alpar@9	2090
alpar@9	2091 int npp_forcing_row(NPP npp, NPPROW p, int at)
alpar@9	2092 { /* process forcing row */
alpar@9	2093 struct forcing_row *info;
alpar@9	2094 struct forcing_col *col = NULL;
alpar@9	2095 NPPCOL *j;
alpar@9	2096 NPPAIJ apj, aij;
alpar@9	2097 NPPLFE *lfe;
alpar@9	2098 double big;
alpar@9	2099 xassert(at == 0 \|\| at == 1);
alpar@9	2100 /* determine maximal magnitude of the row coefficients */
alpar@9	2101 big = 1.0;
alpar@9	2102 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2103 if (big < fabs(apj->val)) big = fabs(apj->val);
alpar@9	2104 /* if there are too small coefficients in the row, transformation
alpar@9	2105 should not be applied */
alpar@9	2106 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2107 if (fabs(apj->val) < 1e-7 * big) return 1;
alpar@9	2108 /* create transformation stack entry */
alpar@9	2109 info = npp_push_tse(npp,
alpar@9	2110 rcv_forcing_row, sizeof(struct forcing_row));
alpar@9	2111 info->p = p->i;
alpar@9	2112 if (p->lb == p->ub)
alpar@9	2113 { /* equality constraint */
alpar@9	2114 info->stat = GLP_NS;
alpar@9	2115 }
alpar@9	2116 else if (at == 0)
alpar@9	2117 { /* inequality constraint; case L[p] = U'[p] */
alpar@9	2118 info->stat = GLP_NL;
alpar@9	2119 xassert(p->lb != -DBL_MAX);
alpar@9	2120 }
alpar@9	2121 else /* at == 1 */
alpar@9	2122 { /* inequality constraint; case U[p] = L'[p] */
alpar@9	2123 info->stat = GLP_NU;
alpar@9	2124 xassert(p->ub != +DBL_MAX);
alpar@9	2125 }
alpar@9	2126 info->ptr = NULL;
alpar@9	2127 /* scan the forcing row, fix columns at corresponding bounds, and
alpar@9	2128 save column information (the latter is not needed for MIP) */
alpar@9	2129 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2130 { /* column j has non-zero coefficient in the forcing row */
alpar@9	2131 j = apj->col;
alpar@9	2132 /* it must be non-fixed */
alpar@9	2133 xassert(j->lb < j->ub);
alpar@9	2134 /* allocate stack entry to save column information */
alpar@9	2135 if (npp->sol != GLP_MIP)
alpar@9	2136 { col = dmp_get_atom(npp->stack, sizeof(struct forcing_col));
alpar@9	2137 col->j = j->j;
alpar@9	2138 col->stat = -1; /* will be set below */
alpar@9	2139 col->a = apj->val;
alpar@9	2140 col->c = j->coef;
alpar@9	2141 col->ptr = NULL;
alpar@9	2142 col->next = info->ptr;
alpar@9	2143 info->ptr = col;
alpar@9	2144 }
alpar@9	2145 /* fix column j */
alpar@9	2146 if (at == 0 && apj->val < 0.0 \|\| at != 0 && apj->val > 0.0)
alpar@9	2147 { /* at its lower bound */
alpar@9	2148 if (npp->sol != GLP_MIP)
alpar@9	2149 col->stat = GLP_NL;
alpar@9	2150 xassert(j->lb != -DBL_MAX);
alpar@9	2151 j->ub = j->lb;
alpar@9	2152 }
alpar@9	2153 else
alpar@9	2154 { /* at its upper bound */
alpar@9	2155 if (npp->sol != GLP_MIP)
alpar@9	2156 col->stat = GLP_NU;
alpar@9	2157 xassert(j->ub != +DBL_MAX);
alpar@9	2158 j->lb = j->ub;
alpar@9	2159 }
alpar@9	2160 /* save column coefficients a[i,j], i != p */
alpar@9	2161 if (npp->sol != GLP_MIP)
alpar@9	2162 { for (aij = j->ptr; aij != NULL; aij = aij->c_next)
alpar@9	2163 { if (aij == apj) continue; /* skip a[p,j] */
alpar@9	2164 lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE));
alpar@9	2165 lfe->ref = aij->row->i;
alpar@9	2166 lfe->val = aij->val;
alpar@9	2167 lfe->next = col->ptr;
alpar@9	2168 col->ptr = lfe;
alpar@9	2169 }
alpar@9	2170 }
alpar@9	2171 }
alpar@9	2172 /* make the row free (unbounded) */
alpar@9	2173 p->lb = -DBL_MAX, p->ub = +DBL_MAX;
alpar@9	2174 return 0;
alpar@9	2175 }
alpar@9	2176
alpar@9	2177 static int rcv_forcing_row(NPP npp, void _info)
alpar@9	2178 { /* recover forcing row */
alpar@9	2179 struct forcing_row *info = _info;
alpar@9	2180 struct forcing_col col, piv;
alpar@9	2181 NPPLFE *lfe;
alpar@9	2182 double d, big, temp;
alpar@9	2183 if (npp->sol == GLP_MIP) goto done;
alpar@9	2184 /* initially solution to the original problem is the same as
alpar@9	2185 to the transformed problem, where row p is inactive constraint
alpar@9	2186 with pi[p] = 0, and all columns are non-basic */
alpar@9	2187 if (npp->sol == GLP_SOL)
alpar@9	2188 { if (npp->r_stat[info->p] != GLP_BS)
alpar@9	2189 { npp_error();
alpar@9	2190 return 1;
alpar@9	2191 }
alpar@9	2192 for (col = info->ptr; col != NULL; col = col->next)
alpar@9	2193 { if (npp->c_stat[col->j] != GLP_NS)
alpar@9	2194 { npp_error();
alpar@9	2195 return 1;
alpar@9	2196 }
alpar@9	2197 npp->c_stat[col->j] = col->stat; /* original status */
alpar@9	2198 }
alpar@9	2199 }
alpar@9	2200 /* compute reduced costs d[j] for all columns with formula (10)
alpar@9	2201 and store them in col.c instead objective coefficients */
alpar@9	2202 for (col = info->ptr; col != NULL; col = col->next)
alpar@9	2203 { d = col->c;
alpar@9	2204 for (lfe = col->ptr; lfe != NULL; lfe = lfe->next)
alpar@9	2205 d -= lfe->val * npp->r_pi[lfe->ref];
alpar@9	2206 col->c = d;
alpar@9	2207 }
alpar@9	2208 /* consider columns j, whose multipliers lambda[j] has wrong
alpar@9	2209 sign in solution to the transformed problem (where lambda[j] =
alpar@9	2210 d[j]), and choose column q, whose multipler lambda[q] reaches
alpar@9	2211 zero last on changing row multiplier pi[p]; see (14) */
alpar@9	2212 piv = NULL, big = 0.0;
alpar@9	2213 for (col = info->ptr; col != NULL; col = col->next)
alpar@9	2214 { d = col->c; /* d[j] */
alpar@9	2215 temp = fabs(d / col->a);
alpar@9	2216 if (col->stat == GLP_NL)
alpar@9	2217 { /* column j has active lower bound */
alpar@9	2218 if (d < 0.0 && big < temp)
alpar@9	2219 piv = col, big = temp;
alpar@9	2220 }
alpar@9	2221 else if (col->stat == GLP_NU)
alpar@9	2222 { /* column j has active upper bound */
alpar@9	2223 if (d > 0.0 && big < temp)
alpar@9	2224 piv = col, big = temp;
alpar@9	2225 }
alpar@9	2226 else
alpar@9	2227 { npp_error();
alpar@9	2228 return 1;
alpar@9	2229 }
alpar@9	2230 }
alpar@9	2231 /* if column q does not exist, no correction is needed */
alpar@9	2232 if (piv != NULL)
alpar@9	2233 { /* correct solution; row p becomes active constraint while
alpar@9	2234 column q becomes basic */
alpar@9	2235 if (npp->sol == GLP_SOL)
alpar@9	2236 { npp->r_stat[info->p] = info->stat;
alpar@9	2237 npp->c_stat[piv->j] = GLP_BS;
alpar@9	2238 }
alpar@9	2239 /* assign new value to row multiplier pi[p] = d[p] / a[p,q] */
alpar@9	2240 npp->r_pi[info->p] = piv->c / piv->a;
alpar@9	2241 }
alpar@9	2242 done: return 0;
alpar@9	2243 }
alpar@9	2244
alpar@9	2245 /***********************************************************************
alpar@9	2246 * NAME
alpar@9	2247 *
alpar@9	2248 * npp_analyze_row - perform general row analysis
alpar@9	2249 *
alpar@9	2250 * SYNOPSIS
alpar@9	2251 *
alpar@9	2252 * #include "glpnpp.h"
alpar@9	2253 * int npp_analyze_row(NPP npp, NPPROW p);
alpar@9	2254 *
alpar@9	2255 * DESCRIPTION
alpar@9	2256 *
alpar@9	2257 * The routine npp_analyze_row performs analysis of row p of general
alpar@9	2258 * format:
alpar@9	2259 *
alpar@9	2260 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
alpar@9	2261 * j
alpar@9	2262 *
alpar@9	2263 * l[j] <= x[j] <= u[j], (2)
alpar@9	2264 *
alpar@9	2265 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0.
alpar@9	2266 *
alpar@9	2267 * RETURNS
alpar@9	2268 *
alpar@9	2269 * 0x?0 - row lower bound does not exist or is redundant;
alpar@9	2270 *
alpar@9	2271 * 0x?1 - row lower bound can be active;
alpar@9	2272 *
alpar@9	2273 * 0x?2 - row lower bound is a forcing bound;
alpar@9	2274 *
alpar@9	2275 * 0x0? - row upper bound does not exist or is redundant;
alpar@9	2276 *
alpar@9	2277 * 0x1? - row upper bound can be active;
alpar@9	2278 *
alpar@9	2279 * 0x2? - row upper bound is a forcing bound;
alpar@9	2280 *
alpar@9	2281 * 0x33 - row bounds are inconsistent with column bounds.
alpar@9	2282 *
alpar@9	2283 * ALGORITHM
alpar@9	2284 *
alpar@9	2285 * Analysis of row (1) is based on analysis of its implied lower and
alpar@9	2286 * upper bounds, which are determined by bounds of corresponding columns
alpar@9	2287 * (variables) as follows:
alpar@9	2288 *
alpar@9	2289 * L'[p] = inf sum a[p,j] x[j] =
alpar@9	2290 * j
alpar@9	2291 * (3)
alpar@9	2292 * = sum a[p,j] l[j] + sum a[p,j] u[j],
alpar@9	2293 * j in Jp j in Jn
alpar@9	2294 *
alpar@9	2295 * U'[p] = sup sum a[p,j] x[j] =
alpar@9	2296 * (4)
alpar@9	2297 * = sum a[p,j] u[j] + sum a[p,j] l[j],
alpar@9	2298 * j in Jp j in Jn
alpar@9	2299 *
alpar@9	2300 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
alpar@9	2301 *
alpar@9	2302 * (Note that bounds of all columns in row p are assumed to be correct,
alpar@9	2303 * so L'[p] <= U'[p].)
alpar@9	2304 *
alpar@9	2305 * Analysis of row lower bound L[p] includes the following cases:
alpar@9	2306 *
alpar@9	2307 * 1) if L[p] > U'[p] + eps, where eps is an absolute tolerance for row
alpar@9	2308 * value, row lower bound L[p] and implied row upper bound U'[p] are
alpar@9	2309 * inconsistent, ergo, the problem has no primal feasible solution;
alpar@9	2310 *
alpar@9	2311 * 2) if U'[p] - eps <= L[p] <= U'[p] + eps, i.e. if L[p] =~ U'[p],
alpar@9	2312 * the row is a forcing row on its lower bound (see description of
alpar@9	2313 * the routine npp_forcing_row);
alpar@9	2314 *
alpar@9	2315 * 3) if L[p] > L'[p] + eps, row lower bound L[p] can be active (this
alpar@9	2316 * conclusion does not account other rows in the problem);
alpar@9	2317 *
alpar@9	2318 * 4) if L[p] <= L'[p] + eps, row lower bound L[p] cannot be active, so
alpar@9	2319 * it is redundant and can be removed (replaced by -oo).
alpar@9	2320 *
alpar@9	2321 * Analysis of row upper bound U[p] is performed in a similar way and
alpar@9	2322 * includes the following cases:
alpar@9	2323 *
alpar@9	2324 * 1) if U[p] < L'[p] - eps, row upper bound U[p] and implied row lower
alpar@9	2325 * bound L'[p] are inconsistent, ergo the problem has no primal
alpar@9	2326 * feasible solution;
alpar@9	2327 *
alpar@9	2328 * 2) if L'[p] - eps <= U[p] <= L'[p] + eps, i.e. if U[p] =~ L'[p],
alpar@9	2329 * the row is a forcing row on its upper bound (see description of
alpar@9	2330 * the routine npp_forcing_row);
alpar@9	2331 *
alpar@9	2332 * 3) if U[p] < U'[p] - eps, row upper bound U[p] can be active (this
alpar@9	2333 * conclusion does not account other rows in the problem);
alpar@9	2334 *
alpar@9	2335 * 4) if U[p] >= U'[p] - eps, row upper bound U[p] cannot be active, so
alpar@9	2336 * it is redundant and can be removed (replaced by +oo). */
alpar@9	2337
alpar@9	2338 int npp_analyze_row(NPP npp, NPPROW p)
alpar@9	2339 { /* perform general row analysis */
alpar@9	2340 NPPAIJ *aij;
alpar@9	2341 int ret = 0x00;
alpar@9	2342 double l, u, eps;
alpar@9	2343 xassert(npp == npp);
alpar@9	2344 /* compute implied lower bound L'[p]; see (3) */
alpar@9	2345 l = 0.0;
alpar@9	2346 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
alpar@9	2347 { if (aij->val > 0.0)
alpar@9	2348 { if (aij->col->lb == -DBL_MAX)
alpar@9	2349 { l = -DBL_MAX;
alpar@9	2350 break;
alpar@9	2351 }
alpar@9	2352 l += aij->val * aij->col->lb;
alpar@9	2353 }
alpar@9	2354 else /* aij->val < 0.0 */
alpar@9	2355 { if (aij->col->ub == +DBL_MAX)
alpar@9	2356 { l = -DBL_MAX;
alpar@9	2357 break;
alpar@9	2358 }
alpar@9	2359 l += aij->val * aij->col->ub;
alpar@9	2360 }
alpar@9	2361 }
alpar@9	2362 /* compute implied upper bound U'[p]; see (4) */
alpar@9	2363 u = 0.0;
alpar@9	2364 for (aij = p->ptr; aij != NULL; aij = aij->r_next)
alpar@9	2365 { if (aij->val > 0.0)
alpar@9	2366 { if (aij->col->ub == +DBL_MAX)
alpar@9	2367 { u = +DBL_MAX;
alpar@9	2368 break;
alpar@9	2369 }
alpar@9	2370 u += aij->val * aij->col->ub;
alpar@9	2371 }
alpar@9	2372 else /* aij->val < 0.0 */
alpar@9	2373 { if (aij->col->lb == -DBL_MAX)
alpar@9	2374 { u = +DBL_MAX;
alpar@9	2375 break;
alpar@9	2376 }
alpar@9	2377 u += aij->val * aij->col->lb;
alpar@9	2378 }
alpar@9	2379 }
alpar@9	2380 /* column bounds are assumed correct, so L'[p] <= U'[p] */
alpar@9	2381 /* check if row lower bound is consistent */
alpar@9	2382 if (p->lb != -DBL_MAX)
alpar@9	2383 { eps = 1e-3 + 1e-6 * fabs(p->lb);
alpar@9	2384 if (p->lb - eps > u)
alpar@9	2385 { ret = 0x33;
alpar@9	2386 goto done;
alpar@9	2387 }
alpar@9	2388 }
alpar@9	2389 /* check if row upper bound is consistent */
alpar@9	2390 if (p->ub != +DBL_MAX)
alpar@9	2391 { eps = 1e-3 + 1e-6 * fabs(p->ub);
alpar@9	2392 if (p->ub + eps < l)
alpar@9	2393 { ret = 0x33;
alpar@9	2394 goto done;
alpar@9	2395 }
alpar@9	2396 }
alpar@9	2397 /* check if row lower bound can be active/forcing */
alpar@9	2398 if (p->lb != -DBL_MAX)
alpar@9	2399 { eps = 1e-9 + 1e-12 * fabs(p->lb);
alpar@9	2400 if (p->lb - eps > l)
alpar@9	2401 { if (p->lb + eps <= u)
alpar@9	2402 ret \|= 0x01;
alpar@9	2403 else
alpar@9	2404 ret \|= 0x02;
alpar@9	2405 }
alpar@9	2406 }
alpar@9	2407 /* check if row upper bound can be active/forcing */
alpar@9	2408 if (p->ub != +DBL_MAX)
alpar@9	2409 { eps = 1e-9 + 1e-12 * fabs(p->ub);
alpar@9	2410 if (p->ub + eps < u)
alpar@9	2411 { /* check if the upper bound is forcing */
alpar@9	2412 if (p->ub - eps >= l)
alpar@9	2413 ret \|= 0x10;
alpar@9	2414 else
alpar@9	2415 ret \|= 0x20;
alpar@9	2416 }
alpar@9	2417 }
alpar@9	2418 done: return ret;
alpar@9	2419 }
alpar@9	2420
alpar@9	2421 /***********************************************************************
alpar@9	2422 * NAME
alpar@9	2423 *
alpar@9	2424 * npp_inactive_bound - remove row lower/upper inactive bound
alpar@9	2425 *
alpar@9	2426 * SYNOPSIS
alpar@9	2427 *
alpar@9	2428 * #include "glpnpp.h"
alpar@9	2429 * void npp_inactive_bound(NPP npp, NPPROW p, int which);
alpar@9	2430 *
alpar@9	2431 * DESCRIPTION
alpar@9	2432 *
alpar@9	2433 * The routine npp_inactive_bound removes lower (if which = 0) or upper
alpar@9	2434 * (if which = 1) bound of row p:
alpar@9	2435 *
alpar@9	2436 * L[p] <= sum a[p,j] x[j] <= U[p],
alpar@9	2437 *
alpar@9	2438 * which (bound) is assumed to be redundant.
alpar@9	2439 *
alpar@9	2440 * PROBLEM TRANSFORMATION
alpar@9	2441 *
alpar@9	2442 * If which = 0, current lower bound L[p] of row p is assigned -oo.
alpar@9	2443 * If which = 1, current upper bound U[p] of row p is assigned +oo.
alpar@9	2444 *
alpar@9	2445 * RECOVERING BASIC SOLUTION
alpar@9	2446 *
alpar@9	2447 * If in solution to the transformed problem row p is inactive
alpar@9	2448 * constraint (GLP_BS), its status is not changed in solution to the
alpar@9	2449 * original problem. Otherwise, status of row p in solution to the
alpar@9	2450 * original problem is defined by its type before transformation and
alpar@9	2451 * its status in solution to the transformed problem as follows:
alpar@9	2452 *
alpar@9	2453 * +---------------------+-------+---------------+---------------+
alpar@9	2454 * \| Row \| Flag \| Row status in \| Row status in \|
alpar@9	2455 * \| type \| which \| transfmd soln \| original soln \|
alpar@9	2456 * +---------------------+-------+---------------+---------------+
alpar@9	2457 * \| sum >= L[p] \| 0 \| GLP_NF \| GLP_NL \|
alpar@9	2458 * \| sum <= U[p] \| 1 \| GLP_NF \| GLP_NU \|
alpar@9	2459 * \| L[p] <= sum <= U[p] \| 0 \| GLP_NU \| GLP_NU \|
alpar@9	2460 * \| L[p] <= sum <= U[p] \| 1 \| GLP_NL \| GLP_NL \|
alpar@9	2461 * \| sum = L[p] = U[p] \| 0 \| GLP_NU \| GLP_NS \|
alpar@9	2462 * \| sum = L[p] = U[p] \| 1 \| GLP_NL \| GLP_NS \|
alpar@9	2463 * +---------------------+-------+---------------+---------------+
alpar@9	2464 *
alpar@9	2465 * RECOVERING INTERIOR-POINT SOLUTION
alpar@9	2466 *
alpar@9	2467 * None needed.
alpar@9	2468 *
alpar@9	2469 * RECOVERING MIP SOLUTION
alpar@9	2470 *
alpar@9	2471 * None needed. */
alpar@9	2472
alpar@9	2473 struct inactive_bound
alpar@9	2474 { /* row inactive bound */
alpar@9	2475 int p;
alpar@9	2476 /* row reference number */
alpar@9	2477 char stat;
alpar@9	2478 /* row status (if active constraint) */
alpar@9	2479 };
alpar@9	2480
alpar@9	2481 static int rcv_inactive_bound(NPP npp, void info);
alpar@9	2482
alpar@9	2483 void npp_inactive_bound(NPP npp, NPPROW p, int which)
alpar@9	2484 { /* remove row lower/upper inactive bound */
alpar@9	2485 struct inactive_bound *info;
alpar@9	2486 if (npp->sol == GLP_SOL)
alpar@9	2487 { /* create transformation stack entry */
alpar@9	2488 info = npp_push_tse(npp,
alpar@9	2489 rcv_inactive_bound, sizeof(struct inactive_bound));
alpar@9	2490 info->p = p->i;
alpar@9	2491 if (p->ub == +DBL_MAX)
alpar@9	2492 info->stat = GLP_NL;
alpar@9	2493 else if (p->lb == -DBL_MAX)
alpar@9	2494 info->stat = GLP_NU;
alpar@9	2495 else if (p->lb != p->ub)
alpar@9	2496 info->stat = (char)(which == 0 ? GLP_NU : GLP_NL);
alpar@9	2497 else
alpar@9	2498 info->stat = GLP_NS;
alpar@9	2499 }
alpar@9	2500 /* remove row inactive bound */
alpar@9	2501 if (which == 0)
alpar@9	2502 { xassert(p->lb != -DBL_MAX);
alpar@9	2503 p->lb = -DBL_MAX;
alpar@9	2504 }
alpar@9	2505 else if (which == 1)
alpar@9	2506 { xassert(p->ub != +DBL_MAX);
alpar@9	2507 p->ub = +DBL_MAX;
alpar@9	2508 }
alpar@9	2509 else
alpar@9	2510 xassert(which != which);
alpar@9	2511 return;
alpar@9	2512 }
alpar@9	2513
alpar@9	2514 static int rcv_inactive_bound(NPP npp, void _info)
alpar@9	2515 { /* recover row status */
alpar@9	2516 struct inactive_bound *info = _info;
alpar@9	2517 if (npp->sol != GLP_SOL)
alpar@9	2518 { npp_error();
alpar@9	2519 return 1;
alpar@9	2520 }
alpar@9	2521 if (npp->r_stat[info->p] == GLP_BS)
alpar@9	2522 npp->r_stat[info->p] = GLP_BS;
alpar@9	2523 else
alpar@9	2524 npp->r_stat[info->p] = info->stat;
alpar@9	2525 return 0;
alpar@9	2526 }
alpar@9	2527
alpar@9	2528 /***********************************************************************
alpar@9	2529 * NAME
alpar@9	2530 *
alpar@9	2531 * npp_implied_bounds - determine implied column bounds
alpar@9	2532 *
alpar@9	2533 * SYNOPSIS
alpar@9	2534 *
alpar@9	2535 * #include "glpnpp.h"
alpar@9	2536 * void npp_implied_bounds(NPP npp, NPPROW p);
alpar@9	2537 *
alpar@9	2538 * DESCRIPTION
alpar@9	2539 *
alpar@9	2540 * The routine npp_implied_bounds inspects general row (constraint) p:
alpar@9	2541 *
alpar@9	2542 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
alpar@9	2543 *
alpar@9	2544 * l[j] <= x[j] <= u[j], (2)
alpar@9	2545 *
alpar@9	2546 * where L[p] <= U[p] and l[j] <= u[j] for all a[p,j] != 0, to compute
alpar@9	2547 * implied bounds of columns (variables x[j]) in this row.
alpar@9	2548 *
alpar@9	2549 * The routine stores implied column bounds l'[j] and u'[j] in column
alpar@9	2550 * descriptors (NPPCOL); it does not change current column bounds l[j]
alpar@9	2551 * and u[j]. (Implied column bounds can be then used to strengthen the
alpar@9	2552 * current column bounds; see the routines npp_implied_lower and
alpar@9	2553 * npp_implied_upper).
alpar@9	2554 *
alpar@9	2555 * ALGORITHM
alpar@9	2556 *
alpar@9	2557 * Current column bounds (2) define implied lower and upper bounds of
alpar@9	2558 * row (1) as follows:
alpar@9	2559 *
alpar@9	2560 * L'[p] = inf sum a[p,j] x[j] =
alpar@9	2561 * j
alpar@9	2562 * (3)
alpar@9	2563 * = sum a[p,j] l[j] + sum a[p,j] u[j],
alpar@9	2564 * j in Jp j in Jn
alpar@9	2565 *
alpar@9	2566 * U'[p] = sup sum a[p,j] x[j] =
alpar@9	2567 * (4)
alpar@9	2568 * = sum a[p,j] u[j] + sum a[p,j] l[j],
alpar@9	2569 * j in Jp j in Jn
alpar@9	2570 *
alpar@9	2571 * Jp = {j: a[p,j] > 0}, Jn = {j: a[p,j] < 0}. (5)
alpar@9	2572 *
alpar@9	2573 * (Note that bounds of all columns in row p are assumed to be correct,
alpar@9	2574 * so L'[p] <= U'[p].)
alpar@9	2575 *
alpar@9	2576 * If L[p] > L'[p] and/or U[p] < U'[p], the lower and/or upper bound of
alpar@9	2577 * row (1) can be active, in which case such row defines implied bounds
alpar@9	2578 * of its variables.
alpar@9	2579 *
alpar@9	2580 * Let x[k] be some variable having in row (1) coefficient a[p,k] != 0.
alpar@9	2581 * Consider a case when row lower bound can be active (L[p] > L'[p]):
alpar@9	2582 *
alpar@9	2583 * sum a[p,j] x[j] >= L[p] ==>
alpar@9	2584 * j
alpar@9	2585 *
alpar@9	2586 * sum a[p,j] x[j] + a[p,k] x[k] >= L[p] ==>
alpar@9	2587 * j!=k
alpar@9	2588 * (6)
alpar@9	2589 * a[p,k] x[k] >= L[p] - sum a[p,j] x[j] ==>
alpar@9	2590 * j!=k
alpar@9	2591 *
alpar@9	2592 * a[p,k] x[k] >= L[p,k],
alpar@9	2593 *
alpar@9	2594 * where
alpar@9	2595 *
alpar@9	2596 * L[p,k] = inf(L[p] - sum a[p,j] x[j]) =
alpar@9	2597 * j!=k
alpar@9	2598 *
alpar@9	2599 * = L[p] - sup sum a[p,j] x[j] = (7)
alpar@9	2600 * j!=k
alpar@9	2601 *
alpar@9	2602 * = L[p] - sum a[p,j] u[j] - sum a[p,j] l[j].
alpar@9	2603 * j in Jp\{k} j in Jn\{k}
alpar@9	2604 *
alpar@9	2605 * Thus:
alpar@9	2606 *
alpar@9	2607 * x[k] >= l'[k] = L[p,k] / a[p,k], if a[p,k] > 0, (8)
alpar@9	2608 *
alpar@9	2609 * x[k] <= u'[k] = L[p,k] / a[p,k], if a[p,k] < 0. (9)
alpar@9	2610 *
alpar@9	2611 * where l'[k] and u'[k] are implied lower and upper bounds of variable
alpar@9	2612 * x[k], resp.
alpar@9	2613 *
alpar@9	2614 * Now consider a similar case when row upper bound can be active
alpar@9	2615 * (U[p] < U'[p]):
alpar@9	2616 *
alpar@9	2617 * sum a[p,j] x[j] <= U[p] ==>
alpar@9	2618 * j
alpar@9	2619 *
alpar@9	2620 * sum a[p,j] x[j] + a[p,k] x[k] <= U[p] ==>
alpar@9	2621 * j!=k
alpar@9	2622 * (10)
alpar@9	2623 * a[p,k] x[k] <= U[p] - sum a[p,j] x[j] ==>
alpar@9	2624 * j!=k
alpar@9	2625 *
alpar@9	2626 * a[p,k] x[k] <= U[p,k],
alpar@9	2627 *
alpar@9	2628 * where:
alpar@9	2629 *
alpar@9	2630 * U[p,k] = sup(U[p] - sum a[p,j] x[j]) =
alpar@9	2631 * j!=k
alpar@9	2632 *
alpar@9	2633 * = U[p] - inf sum a[p,j] x[j] = (11)
alpar@9	2634 * j!=k
alpar@9	2635 *
alpar@9	2636 * = U[p] - sum a[p,j] l[j] - sum a[p,j] u[j].
alpar@9	2637 * j in Jp\{k} j in Jn\{k}
alpar@9	2638 *
alpar@9	2639 * Thus:
alpar@9	2640 *
alpar@9	2641 * x[k] <= u'[k] = U[p,k] / a[p,k], if a[p,k] > 0, (12)
alpar@9	2642 *
alpar@9	2643 * x[k] >= l'[k] = U[p,k] / a[p,k], if a[p,k] < 0. (13)
alpar@9	2644 *
alpar@9	2645 * Note that in formulae (8), (9), (12), and (13) coefficient a[p,k]
alpar@9	2646 * must not be too small in magnitude relatively to other non-zero
alpar@9	2647 * coefficients in row (1), i.e. the following condition must hold:
alpar@9	2648 *
alpar@9	2649 * \|a[p,k]\| >= eps * max(1, \|a[p,j]\|), (14)
alpar@9	2650 * j
alpar@9	2651 *
alpar@9	2652 * where eps is a relative tolerance for constraint coefficients.
alpar@9	2653 * Otherwise the implied column bounds can be numerical inreliable. For
alpar@9	2654 * example, using formula (8) for the following inequality constraint:
alpar@9	2655 *
alpar@9	2656 * 1e-12 x1 - x2 - x3 >= 0,
alpar@9	2657 *
alpar@9	2658 * where x1 >= -1, x2, x3, >= 0, may lead to numerically unreliable
alpar@9	2659 * conclusion that x1 >= 0.
alpar@9	2660 *
alpar@9	2661 * Using formulae (8), (9), (12), and (13) to compute implied bounds
alpar@9	2662 * for one variable requires \|J\| operations, where J = {j: a[p,j] != 0},
alpar@9	2663 * because this needs computing L[p,k] and U[p,k]. Thus, computing
alpar@9	2664 * implied bounds for all variables in row (1) would require \|J\|^2
alpar@9	2665 * operations, that is not a good technique. However, the total number
alpar@9	2666 * of operations can be reduced to \|J\| as follows.
alpar@9	2667 *
alpar@9	2668 * Let a[p,k] > 0. Then from (7) and (11) we have:
alpar@9	2669 *
alpar@9	2670 * L[p,k] = L[p] - (U'[p] - a[p,k] u[k]) =
alpar@9	2671 *
alpar@9	2672 * = L[p] - U'[p] + a[p,k] u[k],
alpar@9	2673 *
alpar@9	2674 * U[p,k] = U[p] - (L'[p] - a[p,k] l[k]) =
alpar@9	2675 *
alpar@9	2676 * = U[p] - L'[p] + a[p,k] l[k],
alpar@9	2677 *
alpar@9	2678 * where L'[p] and U'[p] are implied row lower and upper bounds defined
alpar@9	2679 * by formulae (3) and (4). Substituting these expressions into (8) and
alpar@9	2680 * (12) gives:
alpar@9	2681 *
alpar@9	2682 * l'[k] = L[p,k] / a[p,k] = u[k] + (L[p] - U'[p]) / a[p,k], (15)
alpar@9	2683 *
alpar@9	2684 * u'[k] = U[p,k] / a[p,k] = l[k] + (U[p] - L'[p]) / a[p,k]. (16)
alpar@9	2685 *
alpar@9	2686 * Similarly, if a[p,k] < 0, according to (7) and (11) we have:
alpar@9	2687 *
alpar@9	2688 * L[p,k] = L[p] - (U'[p] - a[p,k] l[k]) =
alpar@9	2689 *
alpar@9	2690 * = L[p] - U'[p] + a[p,k] l[k],
alpar@9	2691 *
alpar@9	2692 * U[p,k] = U[p] - (L'[p] - a[p,k] u[k]) =
alpar@9	2693 *
alpar@9	2694 * = U[p] - L'[p] + a[p,k] u[k],
alpar@9	2695 *
alpar@9	2696 * and substituting these expressions into (8) and (12) gives:
alpar@9	2697 *
alpar@9	2698 * l'[k] = U[p,k] / a[p,k] = u[k] + (U[p] - L'[p]) / a[p,k], (17)
alpar@9	2699 *
alpar@9	2700 * u'[k] = L[p,k] / a[p,k] = l[k] + (L[p] - U'[p]) / a[p,k]. (18)
alpar@9	2701 *
alpar@9	2702 * Note that formulae (15)-(18) can be used only if L'[p] and U'[p]
alpar@9	2703 * exist. However, if for some variable x[j] it happens that l[j] = -oo
alpar@9	2704 * and/or u[j] = +oo, values of L'[p] (if a[p,j] > 0) and/or U'[p] (if
alpar@9	2705 * a[p,j] < 0) are undefined. Consider, therefore, the most general
alpar@9	2706 * situation, when some column bounds (2) may not exist.
alpar@9	2707 *
alpar@9	2708 * Let:
alpar@9	2709 *
alpar@9	2710 * J' = {j : (a[p,j] > 0 and l[j] = -oo) or
alpar@9	2711 * (19)
alpar@9	2712 * (a[p,j] < 0 and u[j] = +oo)}.
alpar@9	2713 *
alpar@9	2714 * Then (assuming that row upper bound U[p] can be active) the following
alpar@9	2715 * three cases are possible:
alpar@9	2716 *
alpar@9	2717 * 1) \|J'\| = 0. In this case L'[p] exists, thus, for all variables x[j]
alpar@9	2718 * in row (1) we can use formulae (16) and (17);
alpar@9	2719 *
alpar@9	2720 * 2) J' = {k}. In this case L'[p] = -oo, however, U[p,k] (11) exists,
alpar@9	2721 * so for variable x[k] we can use formulae (12) and (13). Note that
alpar@9	2722 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] < 0)
alpar@9	2723 * or u'[j] = +oo (if a[p,j] > 0);
alpar@9	2724 *
alpar@9	2725 * 3) \|J'\| > 1. In this case for all variables x[j] in row [1] we have
alpar@9	2726 * l'[j] = -oo (if a[p,j] < 0) or u'[j] = +oo (if a[p,j] > 0).
alpar@9	2727 *
alpar@9	2728 * Similarly, let:
alpar@9	2729 *
alpar@9	2730 * J'' = {j : (a[p,j] > 0 and u[j] = +oo) or
alpar@9	2731 * (20)
alpar@9	2732 * (a[p,j] < 0 and l[j] = -oo)}.
alpar@9	2733 *
alpar@9	2734 * Then (assuming that row lower bound L[p] can be active) the following
alpar@9	2735 * three cases are possible:
alpar@9	2736 *
alpar@9	2737 * 1) \|J''\| = 0. In this case U'[p] exists, thus, for all variables x[j]
alpar@9	2738 * in row (1) we can use formulae (15) and (18);
alpar@9	2739 *
alpar@9	2740 * 2) J'' = {k}. In this case U'[p] = +oo, however, L[p,k] (7) exists,
alpar@9	2741 * so for variable x[k] we can use formulae (8) and (9). Note that
alpar@9	2742 * for all other variables x[j] (j != k) l'[j] = -oo (if a[p,j] > 0)
alpar@9	2743 * or u'[j] = +oo (if a[p,j] < 0);
alpar@9	2744 *
alpar@9	2745 * 3) \|J''\| > 1. In this case for all variables x[j] in row (1) we have
alpar@9	2746 * l'[j] = -oo (if a[p,j] > 0) or u'[j] = +oo (if a[p,j] < 0). */
alpar@9	2747
alpar@9	2748 void npp_implied_bounds(NPP npp, NPPROW p)
alpar@9	2749 { NPPAIJ apj, apk;
alpar@9	2750 double big, eps, temp;
alpar@9	2751 xassert(npp == npp);
alpar@9	2752 /* initialize implied bounds for all variables and determine
alpar@9	2753 maximal magnitude of row coefficients a[p,j] */
alpar@9	2754 big = 1.0;
alpar@9	2755 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2756 { apj->col->ll.ll = -DBL_MAX, apj->col->uu.uu = +DBL_MAX;
alpar@9	2757 if (big < fabs(apj->val)) big = fabs(apj->val);
alpar@9	2758 }
alpar@9	2759 eps = 1e-6 * big;
alpar@9	2760 /* process row lower bound (assuming that it can be active) */
alpar@9	2761 if (p->lb != -DBL_MAX)
alpar@9	2762 { apk = NULL;
alpar@9	2763 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2764 { if (apj->val > 0.0 && apj->col->ub == +DBL_MAX \|\|
alpar@9	2765 apj->val < 0.0 && apj->col->lb == -DBL_MAX)
alpar@9	2766 { if (apk == NULL)
alpar@9	2767 apk = apj;
alpar@9	2768 else
alpar@9	2769 goto skip1;
alpar@9	2770 }
alpar@9	2771 }
alpar@9	2772 /* if a[p,k] = NULL then \|J'\| = 0 else J' = { k } */
alpar@9	2773 temp = p->lb;
alpar@9	2774 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2775 { if (apj == apk)
alpar@9	2776 /* skip a[p,k] */;
alpar@9	2777 else if (apj->val > 0.0)
alpar@9	2778 temp -= apj->val * apj->col->ub;
alpar@9	2779 else /* apj->val < 0.0 */
alpar@9	2780 temp -= apj->val * apj->col->lb;
alpar@9	2781 }
alpar@9	2782 /* compute column implied bounds */
alpar@9	2783 if (apk == NULL)
alpar@9	2784 { /* temp = L[p] - U'[p] */
alpar@9	2785 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2786 { if (apj->val >= +eps)
alpar@9	2787 { /* l'[j] := u[j] + (L[p] - U'[p]) / a[p,j] */
alpar@9	2788 apj->col->ll.ll = apj->col->ub + temp / apj->val;
alpar@9	2789 }
alpar@9	2790 else if (apj->val <= -eps)
alpar@9	2791 { /* u'[j] := l[j] + (L[p] - U'[p]) / a[p,j] */
alpar@9	2792 apj->col->uu.uu = apj->col->lb + temp / apj->val;
alpar@9	2793 }
alpar@9	2794 }
alpar@9	2795 }
alpar@9	2796 else
alpar@9	2797 { /* temp = L[p,k] */
alpar@9	2798 if (apk->val >= +eps)
alpar@9	2799 { /* l'[k] := L[p,k] / a[p,k] */
alpar@9	2800 apk->col->ll.ll = temp / apk->val;
alpar@9	2801 }
alpar@9	2802 else if (apk->val <= -eps)
alpar@9	2803 { /* u'[k] := L[p,k] / a[p,k] */
alpar@9	2804 apk->col->uu.uu = temp / apk->val;
alpar@9	2805 }
alpar@9	2806 }
alpar@9	2807 skip1: ;
alpar@9	2808 }
alpar@9	2809 /* process row upper bound (assuming that it can be active) */
alpar@9	2810 if (p->ub != +DBL_MAX)
alpar@9	2811 { apk = NULL;
alpar@9	2812 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2813 { if (apj->val > 0.0 && apj->col->lb == -DBL_MAX \|\|
alpar@9	2814 apj->val < 0.0 && apj->col->ub == +DBL_MAX)
alpar@9	2815 { if (apk == NULL)
alpar@9	2816 apk = apj;
alpar@9	2817 else
alpar@9	2818 goto skip2;
alpar@9	2819 }
alpar@9	2820 }
alpar@9	2821 /* if a[p,k] = NULL then \|J''\| = 0 else J'' = { k } */
alpar@9	2822 temp = p->ub;
alpar@9	2823 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2824 { if (apj == apk)
alpar@9	2825 /* skip a[p,k] */;
alpar@9	2826 else if (apj->val > 0.0)
alpar@9	2827 temp -= apj->val * apj->col->lb;
alpar@9	2828 else /* apj->val < 0.0 */
alpar@9	2829 temp -= apj->val * apj->col->ub;
alpar@9	2830 }
alpar@9	2831 /* compute column implied bounds */
alpar@9	2832 if (apk == NULL)
alpar@9	2833 { /* temp = U[p] - L'[p] */
alpar@9	2834 for (apj = p->ptr; apj != NULL; apj = apj->r_next)
alpar@9	2835 { if (apj->val >= +eps)
alpar@9	2836 { /* u'[j] := l[j] + (U[p] - L'[p]) / a[p,j] */
alpar@9	2837 apj->col->uu.uu = apj->col->lb + temp / apj->val;
alpar@9	2838 }
alpar@9	2839 else if (apj->val <= -eps)
alpar@9	2840 { /* l'[j] := u[j] + (U[p] - L'[p]) / a[p,j] */
alpar@9	2841 apj->col->ll.ll = apj->col->ub + temp / apj->val;
alpar@9	2842 }
alpar@9	2843 }
alpar@9	2844 }
alpar@9	2845 else
alpar@9	2846 { /* temp = U[p,k] */
alpar@9	2847 if (apk->val >= +eps)
alpar@9	2848 { /* u'[k] := U[p,k] / a[p,k] */
alpar@9	2849 apk->col->uu.uu = temp / apk->val;
alpar@9	2850 }
alpar@9	2851 else if (apk->val <= -eps)
alpar@9	2852 { /* l'[k] := U[p,k] / a[p,k] */
alpar@9	2853 apk->col->ll.ll = temp / apk->val;
alpar@9	2854 }
alpar@9	2855 }
alpar@9	2856 skip2: ;
alpar@9	2857 }
alpar@9	2858 return;
alpar@9	2859 }
alpar@9	2860
alpar@9	2861 /* eof */