COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpapi12.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 78.6 KB
Line 
1/* glpapi12.c (basis factorization and simplex tableau routines) */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpapi.h"
26
27/***********************************************************************
28*  NAME
29*
30*  glp_bf_exists - check if the basis factorization exists
31*
32*  SYNOPSIS
33*
34*  int glp_bf_exists(glp_prob *lp);
35*
36*  RETURNS
37*
38*  If the basis factorization for the current basis associated with
39*  the specified problem object exists and therefore is available for
40*  computations, the routine glp_bf_exists returns non-zero. Otherwise
41*  the routine returns zero. */
42
43int glp_bf_exists(glp_prob *lp)
44{     int ret;
45      ret = (lp->m == 0 || lp->valid);
46      return ret;
47}
48
49/***********************************************************************
50*  NAME
51*
52*  glp_factorize - compute the basis factorization
53*
54*  SYNOPSIS
55*
56*  int glp_factorize(glp_prob *lp);
57*
58*  DESCRIPTION
59*
60*  The routine glp_factorize computes the basis factorization for the
61*  current basis associated with the specified problem object.
62*
63*  RETURNS
64*
65*  0  The basis factorization has been successfully computed.
66*
67*  GLP_EBADB
68*     The basis matrix is invalid, i.e. the number of basic (auxiliary
69*     and structural) variables differs from the number of rows in the
70*     problem object.
71*
72*  GLP_ESING
73*     The basis matrix is singular within the working precision.
74*
75*  GLP_ECOND
76*     The basis matrix is ill-conditioned. */
77
78static int b_col(void *info, int j, int ind[], double val[])
79{     glp_prob *lp = info;
80      int m = lp->m;
81      GLPAIJ *aij;
82      int k, len;
83      xassert(1 <= j && j <= m);
84      /* determine the ordinal number of basic auxiliary or structural
85         variable x[k] corresponding to basic variable xB[j] */
86      k = lp->head[j];
87      /* build j-th column of the basic matrix, which is k-th column of
88         the scaled augmented matrix (I | -R*A*S) */
89      if (k <= m)
90      {  /* x[k] is auxiliary variable */
91         len = 1;
92         ind[1] = k;
93         val[1] = 1.0;
94      }
95      else
96      {  /* x[k] is structural variable */
97         len = 0;
98         for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
99         {  len++;
100            ind[len] = aij->row->i;
101            val[len] = - aij->row->rii * aij->val * aij->col->sjj;
102         }
103      }
104      return len;
105}
106
107static void copy_bfcp(glp_prob *lp);
108
109int glp_factorize(glp_prob *lp)
110{     int m = lp->m;
111      int n = lp->n;
112      GLPROW **row = lp->row;
113      GLPCOL **col = lp->col;
114      int *head = lp->head;
115      int j, k, stat, ret;
116      /* invalidate the basis factorization */
117      lp->valid = 0;
118      /* build the basis header */
119      j = 0;
120      for (k = 1; k <= m+n; k++)
121      {  if (k <= m)
122         {  stat = row[k]->stat;
123            row[k]->bind = 0;
124         }
125         else
126         {  stat = col[k-m]->stat;
127            col[k-m]->bind = 0;
128         }
129         if (stat == GLP_BS)
130         {  j++;
131            if (j > m)
132            {  /* too many basic variables */
133               ret = GLP_EBADB;
134               goto fini;
135            }
136            head[j] = k;
137            if (k <= m)
138               row[k]->bind = j;
139            else
140               col[k-m]->bind = j;
141         }
142      }
143      if (j < m)
144      {  /* too few basic variables */
145         ret = GLP_EBADB;
146         goto fini;
147      }
148      /* try to factorize the basis matrix */
149      if (m > 0)
150      {  if (lp->bfd == NULL)
151         {  lp->bfd = bfd_create_it();
152            copy_bfcp(lp);
153         }
154         switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
155         {  case 0:
156               /* ok */
157               break;
158            case BFD_ESING:
159               /* singular matrix */
160               ret = GLP_ESING;
161               goto fini;
162            case BFD_ECOND:
163               /* ill-conditioned matrix */
164               ret = GLP_ECOND;
165               goto fini;
166            default:
167               xassert(lp != lp);
168         }
169         lp->valid = 1;
170      }
171      /* factorization successful */
172      ret = 0;
173fini: /* bring the return code to the calling program */
174      return ret;
175}
176
177/***********************************************************************
178*  NAME
179*
180*  glp_bf_updated - check if the basis factorization has been updated
181*
182*  SYNOPSIS
183*
184*  int glp_bf_updated(glp_prob *lp);
185*
186*  RETURNS
187*
188*  If the basis factorization has been just computed from scratch, the
189*  routine glp_bf_updated returns zero. Otherwise, if the factorization
190*  has been updated one or more times, the routine returns non-zero. */
191
192int glp_bf_updated(glp_prob *lp)
193{     int cnt;
194      if (!(lp->m == 0 || lp->valid))
195         xerror("glp_bf_update: basis factorization does not exist\n");
196#if 0 /* 15/XI-2009 */
197      cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
198#else
199      cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
200#endif
201      return cnt;
202}
203
204/***********************************************************************
205*  NAME
206*
207*  glp_get_bfcp - retrieve basis factorization control parameters
208*
209*  SYNOPSIS
210*
211*  void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
212*
213*  DESCRIPTION
214*
215*  The routine glp_get_bfcp retrieves control parameters, which are
216*  used on computing and updating the basis factorization associated
217*  with the specified problem object.
218*
219*  Current values of control parameters are stored by the routine in
220*  a glp_bfcp structure, which the parameter parm points to. */
221
222void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm)
223{     glp_bfcp *bfcp = lp->bfcp;
224      if (bfcp == NULL)
225      {  parm->type = GLP_BF_FT;
226         parm->lu_size = 0;
227         parm->piv_tol = 0.10;
228         parm->piv_lim = 4;
229         parm->suhl = GLP_ON;
230         parm->eps_tol = 1e-15;
231         parm->max_gro = 1e+10;
232         parm->nfs_max = 100;
233         parm->upd_tol = 1e-6;
234         parm->nrs_max = 100;
235         parm->rs_size = 0;
236      }
237      else
238         memcpy(parm, bfcp, sizeof(glp_bfcp));
239      return;
240}
241
242/***********************************************************************
243*  NAME
244*
245*  glp_set_bfcp - change basis factorization control parameters
246*
247*  SYNOPSIS
248*
249*  void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
250*
251*  DESCRIPTION
252*
253*  The routine glp_set_bfcp changes control parameters, which are used
254*  by internal GLPK routines in computing and updating the basis
255*  factorization associated with the specified problem object.
256*
257*  New values of the control parameters should be passed in a structure
258*  glp_bfcp, which the parameter parm points to.
259*
260*  The parameter parm can be specified as NULL, in which case all
261*  control parameters are reset to their default values. */
262
263#if 0 /* 15/XI-2009 */
264static void copy_bfcp(glp_prob *lp)
265{     glp_bfcp _parm, *parm = &_parm;
266      BFD *bfd = lp->bfd;
267      glp_get_bfcp(lp, parm);
268      xassert(bfd != NULL);
269      bfd->type = parm->type;
270      bfd->lu_size = parm->lu_size;
271      bfd->piv_tol = parm->piv_tol;
272      bfd->piv_lim = parm->piv_lim;
273      bfd->suhl = parm->suhl;
274      bfd->eps_tol = parm->eps_tol;
275      bfd->max_gro = parm->max_gro;
276      bfd->nfs_max = parm->nfs_max;
277      bfd->upd_tol = parm->upd_tol;
278      bfd->nrs_max = parm->nrs_max;
279      bfd->rs_size = parm->rs_size;
280      return;
281}
282#else
283static void copy_bfcp(glp_prob *lp)
284{     glp_bfcp _parm, *parm = &_parm;
285      glp_get_bfcp(lp, parm);
286      bfd_set_parm(lp->bfd, parm);
287      return;
288}
289#endif
290
291void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
292{     glp_bfcp *bfcp = lp->bfcp;
293      if (parm == NULL)
294      {  /* reset to default values */
295         if (bfcp != NULL)
296            xfree(bfcp), lp->bfcp = NULL;
297      }
298      else
299      {  /* set to specified values */
300         if (bfcp == NULL)
301            bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp));
302         memcpy(bfcp, parm, sizeof(glp_bfcp));
303         if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG ||
304               bfcp->type == GLP_BF_GR))
305            xerror("glp_set_bfcp: type = %d; invalid parameter\n",
306               bfcp->type);
307         if (bfcp->lu_size < 0)
308            xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
309               bfcp->lu_size);
310         if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
311            xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
312               bfcp->piv_tol);
313         if (bfcp->piv_lim < 1)
314            xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
315               bfcp->piv_lim);
316         if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
317            xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
318               bfcp->suhl);
319         if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
320            xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
321               bfcp->eps_tol);
322         if (bfcp->max_gro < 1.0)
323            xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
324               bfcp->max_gro);
325         if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
326            xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
327               bfcp->nfs_max);
328         if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
329            xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
330               bfcp->upd_tol);
331         if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
332            xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
333               bfcp->nrs_max);
334         if (bfcp->rs_size < 0)
335            xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
336               bfcp->nrs_max);
337         if (bfcp->rs_size == 0)
338            bfcp->rs_size = 20 * bfcp->nrs_max;
339      }
340      if (lp->bfd != NULL) copy_bfcp(lp);
341      return;
342}
343
344/***********************************************************************
345*  NAME
346*
347*  glp_get_bhead - retrieve the basis header information
348*
349*  SYNOPSIS
350*
351*  int glp_get_bhead(glp_prob *lp, int k);
352*
353*  DESCRIPTION
354*
355*  The routine glp_get_bhead returns the basis header information for
356*  the current basis associated with the specified problem object.
357*
358*  RETURNS
359*
360*  If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
361*  routine returns i. Otherwise, if xB[k] is j-th structural variable
362*  (1 <= j <= n), the routine returns m+j. Here m is the number of rows
363*  and n is the number of columns in the problem object. */
364
365int glp_get_bhead(glp_prob *lp, int k)
366{     if (!(lp->m == 0 || lp->valid))
367         xerror("glp_get_bhead: basis factorization does not exist\n");
368      if (!(1 <= k && k <= lp->m))
369         xerror("glp_get_bhead: k = %d; index out of range\n", k);
370      return lp->head[k];
371}
372
373/***********************************************************************
374*  NAME
375*
376*  glp_get_row_bind - retrieve row index in the basis header
377*
378*  SYNOPSIS
379*
380*  int glp_get_row_bind(glp_prob *lp, int i);
381*
382*  RETURNS
383*
384*  The routine glp_get_row_bind returns the index k of basic variable
385*  xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
386*  in the current basis associated with the specified problem object,
387*  where m is the number of rows. However, if i-th auxiliary variable
388*  is non-basic, the routine returns zero. */
389
390int glp_get_row_bind(glp_prob *lp, int i)
391{     if (!(lp->m == 0 || lp->valid))
392         xerror("glp_get_row_bind: basis factorization does not exist\n"
393            );
394      if (!(1 <= i && i <= lp->m))
395         xerror("glp_get_row_bind: i = %d; row number out of range\n",
396            i);
397      return lp->row[i]->bind;
398}
399
400/***********************************************************************
401*  NAME
402*
403*  glp_get_col_bind - retrieve column index in the basis header
404*
405*  SYNOPSIS
406*
407*  int glp_get_col_bind(glp_prob *lp, int j);
408*
409*  RETURNS
410*
411*  The routine glp_get_col_bind returns the index k of basic variable
412*  xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
413*  in the current basis associated with the specified problem object,
414*  where m is the number of rows, n is the number of columns. However,
415*  if j-th structural variable is non-basic, the routine returns zero.*/
416
417int glp_get_col_bind(glp_prob *lp, int j)
418{     if (!(lp->m == 0 || lp->valid))
419         xerror("glp_get_col_bind: basis factorization does not exist\n"
420            );
421      if (!(1 <= j && j <= lp->n))
422         xerror("glp_get_col_bind: j = %d; column number out of range\n"
423            , j);
424      return lp->col[j]->bind;
425}
426
427/***********************************************************************
428*  NAME
429*
430*  glp_ftran - perform forward transformation (solve system B*x = b)
431*
432*  SYNOPSIS
433*
434*  void glp_ftran(glp_prob *lp, double x[]);
435*
436*  DESCRIPTION
437*
438*  The routine glp_ftran performs forward transformation, i.e. solves
439*  the system B*x = b, where B is the basis matrix corresponding to the
440*  current basis for the specified problem object, x is the vector of
441*  unknowns to be computed, b is the vector of right-hand sides.
442*
443*  On entry elements of the vector b should be stored in dense format
444*  in locations x[1], ..., x[m], where m is the number of rows. On exit
445*  the routine stores elements of the vector x in the same locations.
446*
447*  SCALING/UNSCALING
448*
449*  Let A~ = (I | -A) is the augmented constraint matrix of the original
450*  (unscaled) problem. In the scaled LP problem instead the matrix A the
451*  scaled matrix A" = R*A*S is actually used, so
452*
453*     A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
454*                                                                    (1)
455*         = R*(I | A)*S~ = R*A~*S~,
456*
457*  is the scaled augmented constraint matrix, where R and S are diagonal
458*  scaling matrices used to scale rows and columns of the matrix A, and
459*
460*     S~ = diag(inv(R) | S)                                          (2)
461*
462*  is an augmented diagonal scaling matrix.
463*
464*  By definition:
465*
466*     A~ = (B | N),                                                  (3)
467*
468*  where B is the basic matrix, which consists of basic columns of the
469*  augmented constraint matrix A~, and N is a matrix, which consists of
470*  non-basic columns of A~. From (1) it follows that:
471*
472*     A~" = (B" | N") = (R*B*SB | R*N*SN),                           (4)
473*
474*  where SB and SN are parts of the augmented scaling matrix S~, which
475*  correspond to basic and non-basic variables, respectively. Therefore
476*
477*     B" = R*B*SB,                                                   (5)
478*
479*  which is the scaled basis matrix. */
480
481void glp_ftran(glp_prob *lp, double x[])
482{     int m = lp->m;
483      GLPROW **row = lp->row;
484      GLPCOL **col = lp->col;
485      int i, k;
486      /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
487         B"*x" = b", where b" = R*b, x = SB*x" */
488      if (!(m == 0 || lp->valid))
489         xerror("glp_ftran: basis factorization does not exist\n");
490      /* b" := R*b */
491      for (i = 1; i <= m; i++)
492         x[i] *= row[i]->rii;
493      /* x" := inv(B")*b" */
494      if (m > 0) bfd_ftran(lp->bfd, x);
495      /* x := SB*x" */
496      for (i = 1; i <= m; i++)
497      {  k = lp->head[i];
498         if (k <= m)
499            x[i] /= row[k]->rii;
500         else
501            x[i] *= col[k-m]->sjj;
502      }
503      return;
504}
505
506/***********************************************************************
507*  NAME
508*
509*  glp_btran - perform backward transformation (solve system B'*x = b)
510*
511*  SYNOPSIS
512*
513*  void glp_btran(glp_prob *lp, double x[]);
514*
515*  DESCRIPTION
516*
517*  The routine glp_btran performs backward transformation, i.e. solves
518*  the system B'*x = b, where B' is a matrix transposed to the basis
519*  matrix corresponding to the current basis for the specified problem
520*  problem object, x is the vector of unknowns to be computed, b is the
521*  vector of right-hand sides.
522*
523*  On entry elements of the vector b should be stored in dense format
524*  in locations x[1], ..., x[m], where m is the number of rows. On exit
525*  the routine stores elements of the vector x in the same locations.
526*
527*  SCALING/UNSCALING
528*
529*  See comments to the routine glp_ftran. */
530
531void glp_btran(glp_prob *lp, double x[])
532{     int m = lp->m;
533      GLPROW **row = lp->row;
534      GLPCOL **col = lp->col;
535      int i, k;
536      /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
537         (B")'*x" = b", where b" = SB*b, x = R*x" */
538      if (!(m == 0 || lp->valid))
539         xerror("glp_btran: basis factorization does not exist\n");
540      /* b" := SB*b */
541      for (i = 1; i <= m; i++)
542      {  k = lp->head[i];
543         if (k <= m)
544            x[i] /= row[k]->rii;
545         else
546            x[i] *= col[k-m]->sjj;
547      }
548      /* x" := inv[(B")']*b" */
549      if (m > 0) bfd_btran(lp->bfd, x);
550      /* x := R*x" */
551      for (i = 1; i <= m; i++)
552         x[i] *= row[i]->rii;
553      return;
554}
555
556/***********************************************************************
557*  NAME
558*
559*  glp_warm_up - "warm up" LP basis
560*
561*  SYNOPSIS
562*
563*  int glp_warm_up(glp_prob *P);
564*
565*  DESCRIPTION
566*
567*  The routine glp_warm_up "warms up" the LP basis for the specified
568*  problem object using current statuses assigned to rows and columns
569*  (that is, to auxiliary and structural variables).
570*
571*  This operation includes computing factorization of the basis matrix
572*  (if it does not exist), computing primal and dual components of basic
573*  solution, and determining the solution status.
574*
575*  RETURNS
576*
577*  0  The operation has been successfully performed.
578*
579*  GLP_EBADB
580*     The basis matrix is invalid, i.e. the number of basic (auxiliary
581*     and structural) variables differs from the number of rows in the
582*     problem object.
583*
584*  GLP_ESING
585*     The basis matrix is singular within the working precision.
586*
587*  GLP_ECOND
588*     The basis matrix is ill-conditioned. */
589
590int glp_warm_up(glp_prob *P)
591{     GLPROW *row;
592      GLPCOL *col;
593      GLPAIJ *aij;
594      int i, j, type, ret;
595      double eps, temp, *work;
596      /* invalidate basic solution */
597      P->pbs_stat = P->dbs_stat = GLP_UNDEF;
598      P->obj_val = 0.0;
599      P->some = 0;
600      for (i = 1; i <= P->m; i++)
601      {  row = P->row[i];
602         row->prim = row->dual = 0.0;
603      }
604      for (j = 1; j <= P->n; j++)
605      {  col = P->col[j];
606         col->prim = col->dual = 0.0;
607      }
608      /* compute the basis factorization, if necessary */
609      if (!glp_bf_exists(P))
610      {  ret = glp_factorize(P);
611         if (ret != 0) goto done;
612      }
613      /* allocate working array */
614      work = xcalloc(1+P->m, sizeof(double));
615      /* determine and store values of non-basic variables, compute
616         vector (- N * xN) */
617      for (i = 1; i <= P->m; i++)
618         work[i] = 0.0;
619      for (i = 1; i <= P->m; i++)
620      {  row = P->row[i];
621         if (row->stat == GLP_BS)
622            continue;
623         else if (row->stat == GLP_NL)
624            row->prim = row->lb;
625         else if (row->stat == GLP_NU)
626            row->prim = row->ub;
627         else if (row->stat == GLP_NF)
628            row->prim = 0.0;
629         else if (row->stat == GLP_NS)
630            row->prim = row->lb;
631         else
632            xassert(row != row);
633         /* N[j] is i-th column of matrix (I|-A) */
634         work[i] -= row->prim;
635      }
636      for (j = 1; j <= P->n; j++)
637      {  col = P->col[j];
638         if (col->stat == GLP_BS)
639            continue;
640         else if (col->stat == GLP_NL)
641            col->prim = col->lb;
642         else if (col->stat == GLP_NU)
643            col->prim = col->ub;
644         else if (col->stat == GLP_NF)
645            col->prim = 0.0;
646         else if (col->stat == GLP_NS)
647            col->prim = col->lb;
648         else
649            xassert(col != col);
650         /* N[j] is (m+j)-th column of matrix (I|-A) */
651         if (col->prim != 0.0)
652         {  for (aij = col->ptr; aij != NULL; aij = aij->c_next)
653               work[aij->row->i] += aij->val * col->prim;
654         }
655      }
656      /* compute vector of basic variables xB = - inv(B) * N * xN */
657      glp_ftran(P, work);
658      /* store values of basic variables, check primal feasibility */
659      P->pbs_stat = GLP_FEAS;
660      for (i = 1; i <= P->m; i++)
661      {  row = P->row[i];
662         if (row->stat != GLP_BS)
663            continue;
664         row->prim = work[row->bind];
665         type = row->type;
666         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
667         {  eps = 1e-6 + 1e-9 * fabs(row->lb);
668            if (row->prim < row->lb - eps)
669               P->pbs_stat = GLP_INFEAS;
670         }
671         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
672         {  eps = 1e-6 + 1e-9 * fabs(row->ub);
673            if (row->prim > row->ub + eps)
674               P->pbs_stat = GLP_INFEAS;
675         }
676      }
677      for (j = 1; j <= P->n; j++)
678      {  col = P->col[j];
679         if (col->stat != GLP_BS)
680            continue;
681         col->prim = work[col->bind];
682         type = col->type;
683         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
684         {  eps = 1e-6 + 1e-9 * fabs(col->lb);
685            if (col->prim < col->lb - eps)
686               P->pbs_stat = GLP_INFEAS;
687         }
688         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
689         {  eps = 1e-6 + 1e-9 * fabs(col->ub);
690            if (col->prim > col->ub + eps)
691               P->pbs_stat = GLP_INFEAS;
692         }
693      }
694      /* compute value of the objective function */
695      P->obj_val = P->c0;
696      for (j = 1; j <= P->n; j++)
697      {  col = P->col[j];
698         P->obj_val += col->coef * col->prim;
699      }
700      /* build vector cB of objective coefficients at basic variables */
701      for (i = 1; i <= P->m; i++)
702         work[i] = 0.0;
703      for (j = 1; j <= P->n; j++)
704      {  col = P->col[j];
705         if (col->stat == GLP_BS)
706            work[col->bind] = col->coef;
707      }
708      /* compute vector of simplex multipliers pi = inv(B') * cB */
709      glp_btran(P, work);
710      /* compute and store reduced costs of non-basic variables d[j] =
711         c[j] - N'[j] * pi, check dual feasibility */
712      P->dbs_stat = GLP_FEAS;
713      for (i = 1; i <= P->m; i++)
714      {  row = P->row[i];
715         if (row->stat == GLP_BS)
716         {  row->dual = 0.0;
717            continue;
718         }
719         /* N[j] is i-th column of matrix (I|-A) */
720         row->dual = - work[i];
721         type = row->type;
722         temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
723         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
724             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
725            P->dbs_stat = GLP_INFEAS;
726      }
727      for (j = 1; j <= P->n; j++)
728      {  col = P->col[j];
729         if (col->stat == GLP_BS)
730         {  col->dual = 0.0;
731            continue;
732         }
733         /* N[j] is (m+j)-th column of matrix (I|-A) */
734         col->dual = col->coef;
735         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
736            col->dual += aij->val * work[aij->row->i];
737         type = col->type;
738         temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
739         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
740             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
741            P->dbs_stat = GLP_INFEAS;
742      }
743      /* free working array */
744      xfree(work);
745      ret = 0;
746done: return ret;
747}
748
749/***********************************************************************
750*  NAME
751*
752*  glp_eval_tab_row - compute row of the simplex tableau
753*
754*  SYNOPSIS
755*
756*  int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
757*
758*  DESCRIPTION
759*
760*  The routine glp_eval_tab_row computes a row of the current simplex
761*  tableau for the basic variable, which is specified by the number k:
762*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
763*  x[k] is (k-m)-th structural variable, where m is number of rows, and
764*  n is number of columns. The current basis must be available.
765*
766*  The routine stores column indices and numerical values of non-zero
767*  elements of the computed row using sparse format to the locations
768*  ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
769*  0 <= len <= n is number of non-zeros returned on exit.
770*
771*  Element indices stored in the array ind have the same sense as the
772*  index k, i.e. indices 1 to m denote auxiliary variables and indices
773*  m+1 to m+n denote structural ones (all these variables are obviously
774*  non-basic by definition).
775*
776*  The computed row shows how the specified basic variable x[k] = xB[i]
777*  depends on non-basic variables:
778*
779*     xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
780*
781*  where alfa[i,j] are elements of the simplex table row, xN[j] are
782*  non-basic (auxiliary and structural) variables.
783*
784*  RETURNS
785*
786*  The routine returns number of non-zero elements in the simplex table
787*  row stored in the arrays ind and val.
788*
789*  BACKGROUND
790*
791*  The system of equality constraints of the LP problem is:
792*
793*     xR = A * xS,                                                   (1)
794*
795*  where xR is the vector of auxliary variables, xS is the vector of
796*  structural variables, A is the matrix of constraint coefficients.
797*
798*  The system (1) can be written in homogenous form as follows:
799*
800*     A~ * x = 0,                                                    (2)
801*
802*  where A~ = (I | -A) is the augmented constraint matrix (has m rows
803*  and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
804*  structural) variables.
805*
806*  By definition for the current basis we have:
807*
808*     A~ = (B | N),                                                  (3)
809*
810*  where B is the basis matrix. Thus, the system (2) can be written as:
811*
812*     B * xB + N * xN = 0.                                           (4)
813*
814*  From (4) it follows that:
815*
816*     xB = A^ * xN,                                                  (5)
817*
818*  where the matrix
819*
820*     A^ = - inv(B) * N                                              (6)
821*
822*  is called the simplex table.
823*
824*  It is understood that i-th row of the simplex table is:
825*
826*     e * A^ = - e * inv(B) * N,                                     (7)
827*
828*  where e is a unity vector with e[i] = 1.
829*
830*  To compute i-th row of the simplex table the routine first computes
831*  i-th row of the inverse:
832*
833*     rho = inv(B') * e,                                             (8)
834*
835*  where B' is a matrix transposed to B, and then computes elements of
836*  i-th row of the simplex table as scalar products:
837*
838*     alfa[i,j] = - rho * N[j]   for all j,                          (9)
839*
840*  where N[j] is a column of the augmented constraint matrix A~, which
841*  corresponds to some non-basic auxiliary or structural variable. */
842
843int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
844{     int m = lp->m;
845      int n = lp->n;
846      int i, t, len, lll, *iii;
847      double alfa, *rho, *vvv;
848      if (!(m == 0 || lp->valid))
849         xerror("glp_eval_tab_row: basis factorization does not exist\n"
850            );
851      if (!(1 <= k && k <= m+n))
852         xerror("glp_eval_tab_row: k = %d; variable number out of range"
853            , k);
854      /* determine xB[i] which corresponds to x[k] */
855      if (k <= m)
856         i = glp_get_row_bind(lp, k);
857      else
858         i = glp_get_col_bind(lp, k-m);
859      if (i == 0)
860         xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
861      xassert(1 <= i && i <= m);
862      /* allocate working arrays */
863      rho = xcalloc(1+m, sizeof(double));
864      iii = xcalloc(1+m, sizeof(int));
865      vvv = xcalloc(1+m, sizeof(double));
866      /* compute i-th row of the inverse; see (8) */
867      for (t = 1; t <= m; t++) rho[t] = 0.0;
868      rho[i] = 1.0;
869      glp_btran(lp, rho);
870      /* compute i-th row of the simplex table */
871      len = 0;
872      for (k = 1; k <= m+n; k++)
873      {  if (k <= m)
874         {  /* x[k] is auxiliary variable, so N[k] is a unity column */
875            if (glp_get_row_stat(lp, k) == GLP_BS) continue;
876            /* compute alfa[i,j]; see (9) */
877            alfa = - rho[k];
878         }
879         else
880         {  /* x[k] is structural variable, so N[k] is a column of the
881               original constraint matrix A with negative sign */
882            if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
883            /* compute alfa[i,j]; see (9) */
884            lll = glp_get_mat_col(lp, k-m, iii, vvv);
885            alfa = 0.0;
886            for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
887         }
888         /* store alfa[i,j] */
889         if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
890      }
891      xassert(len <= n);
892      /* free working arrays */
893      xfree(rho);
894      xfree(iii);
895      xfree(vvv);
896      /* return to the calling program */
897      return len;
898}
899
900/***********************************************************************
901*  NAME
902*
903*  glp_eval_tab_col - compute column of the simplex tableau
904*
905*  SYNOPSIS
906*
907*  int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
908*
909*  DESCRIPTION
910*
911*  The routine glp_eval_tab_col computes a column of the current simplex
912*  table for the non-basic variable, which is specified by the number k:
913*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
914*  x[k] is (k-m)-th structural variable, where m is number of rows, and
915*  n is number of columns. The current basis must be available.
916*
917*  The routine stores row indices and numerical values of non-zero
918*  elements of the computed column using sparse format to the locations
919*  ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
920*  0 <= len <= m is number of non-zeros returned on exit.
921*
922*  Element indices stored in the array ind have the same sense as the
923*  index k, i.e. indices 1 to m denote auxiliary variables and indices
924*  m+1 to m+n denote structural ones (all these variables are obviously
925*  basic by the definition).
926*
927*  The computed column shows how basic variables depend on the specified
928*  non-basic variable x[k] = xN[j]:
929*
930*     xB[1] = ... + alfa[1,j]*xN[j] + ...
931*     xB[2] = ... + alfa[2,j]*xN[j] + ...
932*              . . . . . .
933*     xB[m] = ... + alfa[m,j]*xN[j] + ...
934*
935*  where alfa[i,j] are elements of the simplex table column, xB[i] are
936*  basic (auxiliary and structural) variables.
937*
938*  RETURNS
939*
940*  The routine returns number of non-zero elements in the simplex table
941*  column stored in the arrays ind and val.
942*
943*  BACKGROUND
944*
945*  As it was explained in comments to the routine glp_eval_tab_row (see
946*  above) the simplex table is the following matrix:
947*
948*     A^ = - inv(B) * N.                                             (1)
949*
950*  Therefore j-th column of the simplex table is:
951*
952*     A^ * e = - inv(B) * N * e = - inv(B) * N[j],                   (2)
953*
954*  where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
955*  is a column of the augmented constraint matrix A~, which corresponds
956*  to the given non-basic auxiliary or structural variable. */
957
958int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
959{     int m = lp->m;
960      int n = lp->n;
961      int t, len, stat;
962      double *col;
963      if (!(m == 0 || lp->valid))
964         xerror("glp_eval_tab_col: basis factorization does not exist\n"
965            );
966      if (!(1 <= k && k <= m+n))
967         xerror("glp_eval_tab_col: k = %d; variable number out of range"
968            , k);
969      if (k <= m)
970         stat = glp_get_row_stat(lp, k);
971      else
972         stat = glp_get_col_stat(lp, k-m);
973      if (stat == GLP_BS)
974         xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
975            k);
976      /* obtain column N[k] with negative sign */
977      col = xcalloc(1+m, sizeof(double));
978      for (t = 1; t <= m; t++) col[t] = 0.0;
979      if (k <= m)
980      {  /* x[k] is auxiliary variable, so N[k] is a unity column */
981         col[k] = -1.0;
982      }
983      else
984      {  /* x[k] is structural variable, so N[k] is a column of the
985            original constraint matrix A with negative sign */
986         len = glp_get_mat_col(lp, k-m, ind, val);
987         for (t = 1; t <= len; t++) col[ind[t]] = val[t];
988      }
989      /* compute column of the simplex table, which corresponds to the
990         specified non-basic variable x[k] */
991      glp_ftran(lp, col);
992      len = 0;
993      for (t = 1; t <= m; t++)
994      {  if (col[t] != 0.0)
995         {  len++;
996            ind[len] = glp_get_bhead(lp, t);
997            val[len] = col[t];
998         }
999      }
1000      xfree(col);
1001      /* return to the calling program */
1002      return len;
1003}
1004
1005/***********************************************************************
1006*  NAME
1007*
1008*  glp_transform_row - transform explicitly specified row
1009*
1010*  SYNOPSIS
1011*
1012*  int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
1013*
1014*  DESCRIPTION
1015*
1016*  The routine glp_transform_row performs the same operation as the
1017*  routine glp_eval_tab_row with exception that the row to be
1018*  transformed is specified explicitly as a sparse vector.
1019*
1020*  The explicitly specified row may be thought as a linear form:
1021*
1022*     x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],             (1)
1023*
1024*  where x is an auxiliary variable for this row, a[j] are coefficients
1025*  of the linear form, x[m+j] are structural variables.
1026*
1027*  On entry column indices and numerical values of non-zero elements of
1028*  the row should be stored in locations ind[1], ..., ind[len] and
1029*  val[1], ..., val[len], where len is the number of non-zero elements.
1030*
1031*  This routine uses the system of equality constraints and the current
1032*  basis in order to express the auxiliary variable x in (1) through the
1033*  current non-basic variables (as if the transformed row were added to
1034*  the problem object and its auxiliary variable were basic), i.e. the
1035*  resultant row has the form:
1036*
1037*     x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n],       (2)
1038*
1039*  where xN[j] are non-basic (auxiliary or structural) variables, n is
1040*  the number of columns in the LP problem object.
1041*
1042*  On exit the routine stores indices and numerical values of non-zero
1043*  elements of the resultant row (2) in locations ind[1], ..., ind[len']
1044*  and val[1], ..., val[len'], where 0 <= len' <= n is the number of
1045*  non-zero elements in the resultant row returned by the routine. Note
1046*  that indices (numbers) of non-basic variables stored in the array ind
1047*  correspond to original ordinal numbers of variables: indices 1 to m
1048*  mean auxiliary variables and indices m+1 to m+n mean structural ones.
1049*
1050*  RETURNS
1051*
1052*  The routine returns len', which is the number of non-zero elements in
1053*  the resultant row stored in the arrays ind and val.
1054*
1055*  BACKGROUND
1056*
1057*  The explicitly specified row (1) is transformed in the same way as it
1058*  were the objective function row.
1059*
1060*  From (1) it follows that:
1061*
1062*     x = aB * xB + aN * xN,                                         (3)
1063*
1064*  where xB is the vector of basic variables, xN is the vector of
1065*  non-basic variables.
1066*
1067*  The simplex table, which corresponds to the current basis, is:
1068*
1069*     xB = [-inv(B) * N] * xN.                                       (4)
1070*
1071*  Therefore substituting xB from (4) to (3) we have:
1072*
1073*     x = aB * [-inv(B) * N] * xN + aN * xN =
1074*                                                                    (5)
1075*       = rho * (-N) * xN + aN * xN = alfa * xN,
1076*
1077*  where:
1078*
1079*     rho = inv(B') * aB,                                            (6)
1080*
1081*  and
1082*
1083*     alfa = aN + rho * (-N)                                         (7)
1084*
1085*  is the resultant row computed by the routine. */
1086
1087int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
1088{     int i, j, k, m, n, t, lll, *iii;
1089      double alfa, *a, *aB, *rho, *vvv;
1090      if (!glp_bf_exists(P))
1091         xerror("glp_transform_row: basis factorization does not exist "
1092            "\n");
1093      m = glp_get_num_rows(P);
1094      n = glp_get_num_cols(P);
1095      /* unpack the row to be transformed to the array a */
1096      a = xcalloc(1+n, sizeof(double));
1097      for (j = 1; j <= n; j++) a[j] = 0.0;
1098      if (!(0 <= len && len <= n))
1099         xerror("glp_transform_row: len = %d; invalid row length\n",
1100            len);
1101      for (t = 1; t <= len; t++)
1102      {  j = ind[t];
1103         if (!(1 <= j && j <= n))
1104            xerror("glp_transform_row: ind[%d] = %d; column index out o"
1105               "f range\n", t, j);
1106         if (val[t] == 0.0)
1107            xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
1108               "t allowed\n", t);
1109         if (a[j] != 0.0)
1110            xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
1111               "ndices not allowed\n", t, j);
1112         a[j] = val[t];
1113      }
1114      /* construct the vector aB */
1115      aB = xcalloc(1+m, sizeof(double));
1116      for (i = 1; i <= m; i++)
1117      {  k = glp_get_bhead(P, i);
1118         /* xB[i] is k-th original variable */
1119         xassert(1 <= k && k <= m+n);
1120         aB[i] = (k <= m ? 0.0 : a[k-m]);
1121      }
1122      /* solve the system B'*rho = aB to compute the vector rho */
1123      rho = aB, glp_btran(P, rho);
1124      /* compute coefficients at non-basic auxiliary variables */
1125      len = 0;
1126      for (i = 1; i <= m; i++)
1127      {  if (glp_get_row_stat(P, i) != GLP_BS)
1128         {  alfa = - rho[i];
1129            if (alfa != 0.0)
1130            {  len++;
1131               ind[len] = i;
1132               val[len] = alfa;
1133            }
1134         }
1135      }
1136      /* compute coefficients at non-basic structural variables */
1137      iii = xcalloc(1+m, sizeof(int));
1138      vvv = xcalloc(1+m, sizeof(double));
1139      for (j = 1; j <= n; j++)
1140      {  if (glp_get_col_stat(P, j) != GLP_BS)
1141         {  alfa = a[j];
1142            lll = glp_get_mat_col(P, j, iii, vvv);
1143            for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
1144            if (alfa != 0.0)
1145            {  len++;
1146               ind[len] = m+j;
1147               val[len] = alfa;
1148            }
1149         }
1150      }
1151      xassert(len <= n);
1152      xfree(iii);
1153      xfree(vvv);
1154      xfree(aB);
1155      xfree(a);
1156      return len;
1157}
1158
1159/***********************************************************************
1160*  NAME
1161*
1162*  glp_transform_col - transform explicitly specified column
1163*
1164*  SYNOPSIS
1165*
1166*  int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
1167*
1168*  DESCRIPTION
1169*
1170*  The routine glp_transform_col performs the same operation as the
1171*  routine glp_eval_tab_col with exception that the column to be
1172*  transformed is specified explicitly as a sparse vector.
1173*
1174*  The explicitly specified column may be thought as if it were added
1175*  to the original system of equality constraints:
1176*
1177*     x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
1178*     x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x            (1)
1179*        .  .  .  .  .  .  .  .  .  .  .  .  .  .  .
1180*     x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
1181*
1182*  where x[i] are auxiliary variables, x[m+j] are structural variables,
1183*  x is a structural variable for the explicitly specified column, a[i]
1184*  are constraint coefficients for x.
1185*
1186*  On entry row indices and numerical values of non-zero elements of
1187*  the column should be stored in locations ind[1], ..., ind[len] and
1188*  val[1], ..., val[len], where len is the number of non-zero elements.
1189*
1190*  This routine uses the system of equality constraints and the current
1191*  basis in order to express the current basic variables through the
1192*  structural variable x in (1) (as if the transformed column were added
1193*  to the problem object and the variable x were non-basic), i.e. the
1194*  resultant column has the form:
1195*
1196*     xB[1] = ... + alfa[1]*x
1197*     xB[2] = ... + alfa[2]*x                                        (2)
1198*        .  .  .  .  .  .
1199*     xB[m] = ... + alfa[m]*x
1200*
1201*  where xB are basic (auxiliary and structural) variables, m is the
1202*  number of rows in the problem object.
1203*
1204*  On exit the routine stores indices and numerical values of non-zero
1205*  elements of the resultant column (2) in locations ind[1], ...,
1206*  ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
1207*  number of non-zero element in the resultant column returned by the
1208*  routine. Note that indices (numbers) of basic variables stored in
1209*  the array ind correspond to original ordinal numbers of variables:
1210*  indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
1211*  structural ones.
1212*
1213*  RETURNS
1214*
1215*  The routine returns len', which is the number of non-zero elements
1216*  in the resultant column stored in the arrays ind and val.
1217*
1218*  BACKGROUND
1219*
1220*  The explicitly specified column (1) is transformed in the same way
1221*  as any other column of the constraint matrix using the formula:
1222*
1223*     alfa = inv(B) * a,                                             (3)
1224*
1225*  where alfa is the resultant column computed by the routine. */
1226
1227int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
1228{     int i, m, t;
1229      double *a, *alfa;
1230      if (!glp_bf_exists(P))
1231         xerror("glp_transform_col: basis factorization does not exist "
1232            "\n");
1233      m = glp_get_num_rows(P);
1234      /* unpack the column to be transformed to the array a */
1235      a = xcalloc(1+m, sizeof(double));
1236      for (i = 1; i <= m; i++) a[i] = 0.0;
1237      if (!(0 <= len && len <= m))
1238         xerror("glp_transform_col: len = %d; invalid column length\n",
1239            len);
1240      for (t = 1; t <= len; t++)
1241      {  i = ind[t];
1242         if (!(1 <= i && i <= m))
1243            xerror("glp_transform_col: ind[%d] = %d; row index out of r"
1244               "ange\n", t, i);
1245         if (val[t] == 0.0)
1246            xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
1247               "t allowed\n", t);
1248         if (a[i] != 0.0)
1249            xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
1250               "ces not allowed\n", t, i);
1251         a[i] = val[t];
1252      }
1253      /* solve the system B*a = alfa to compute the vector alfa */
1254      alfa = a, glp_ftran(P, alfa);
1255      /* store resultant coefficients */
1256      len = 0;
1257      for (i = 1; i <= m; i++)
1258      {  if (alfa[i] != 0.0)
1259         {  len++;
1260            ind[len] = glp_get_bhead(P, i);
1261            val[len] = alfa[i];
1262         }
1263      }
1264      xfree(a);
1265      return len;
1266}
1267
1268/***********************************************************************
1269*  NAME
1270*
1271*  glp_prim_rtest - perform primal ratio test
1272*
1273*  SYNOPSIS
1274*
1275*  int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1276*     const double val[], int dir, double eps);
1277*
1278*  DESCRIPTION
1279*
1280*  The routine glp_prim_rtest performs the primal ratio test using an
1281*  explicitly specified column of the simplex table.
1282*
1283*  The current basic solution associated with the LP problem object
1284*  must be primal feasible.
1285*
1286*  The explicitly specified column of the simplex table shows how the
1287*  basic variables xB depend on some non-basic variable x (which is not
1288*  necessarily presented in the problem object):
1289*
1290*     xB[1] = ... + alfa[1] * x + ...
1291*     xB[2] = ... + alfa[2] * x + ...                                (*)
1292*         .  .  .  .  .  .  .  .
1293*     xB[m] = ... + alfa[m] * x + ...
1294*
1295*  The column (*) is specifed on entry to the routine using the sparse
1296*  format. Ordinal numbers of basic variables xB[i] should be placed in
1297*  locations ind[1], ..., ind[len], where ordinal number 1 to m denote
1298*  auxiliary variables, and ordinal numbers m+1 to m+n denote structural
1299*  variables. The corresponding non-zero coefficients alfa[i] should be
1300*  placed in locations val[1], ..., val[len]. The arrays ind and val are
1301*  not changed on exit.
1302*
1303*  The parameter dir specifies direction in which the variable x changes
1304*  on entering the basis: +1 means increasing, -1 means decreasing.
1305*
1306*  The parameter eps is an absolute tolerance (small positive number)
1307*  used by the routine to skip small alfa[j] of the row (*).
1308*
1309*  The routine determines which basic variable (among specified in
1310*  ind[1], ..., ind[len]) should leave the basis in order to keep primal
1311*  feasibility.
1312*
1313*  RETURNS
1314*
1315*  The routine glp_prim_rtest returns the index piv in the arrays ind
1316*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
1317*  If the adjacent basic solution is primal unbounded and therefore the
1318*  choice cannot be made, the routine returns zero.
1319*
1320*  COMMENTS
1321*
1322*  If the non-basic variable x is presented in the LP problem object,
1323*  the column (*) can be computed with the routine glp_eval_tab_col;
1324*  otherwise it can be computed with the routine glp_transform_col. */
1325
1326int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1327      const double val[], int dir, double eps)
1328{     int k, m, n, piv, t, type, stat;
1329      double alfa, big, beta, lb, ub, temp, teta;
1330      if (glp_get_prim_stat(P) != GLP_FEAS)
1331         xerror("glp_prim_rtest: basic solution is not primal feasible "
1332            "\n");
1333      if (!(dir == +1 || dir == -1))
1334         xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
1335      if (!(0.0 < eps && eps < 1.0))
1336         xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
1337      m = glp_get_num_rows(P);
1338      n = glp_get_num_cols(P);
1339      /* initial settings */
1340      piv = 0, teta = DBL_MAX, big = 0.0;
1341      /* walk through the entries of the specified column */
1342      for (t = 1; t <= len; t++)
1343      {  /* get the ordinal number of basic variable */
1344         k = ind[t];
1345         if (!(1 <= k && k <= m+n))
1346            xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
1347               "f range\n", t, k);
1348         /* determine type, bounds, status and primal value of basic
1349            variable xB[i] = x[k] in the current basic solution */
1350         if (k <= m)
1351         {  type = glp_get_row_type(P, k);
1352            lb = glp_get_row_lb(P, k);
1353            ub = glp_get_row_ub(P, k);
1354            stat = glp_get_row_stat(P, k);
1355            beta = glp_get_row_prim(P, k);
1356         }
1357         else
1358         {  type = glp_get_col_type(P, k-m);
1359            lb = glp_get_col_lb(P, k-m);
1360            ub = glp_get_col_ub(P, k-m);
1361            stat = glp_get_col_stat(P, k-m);
1362            beta = glp_get_col_prim(P, k-m);
1363         }
1364         if (stat != GLP_BS)
1365            xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
1366               "t allowed\n", t, k);
1367         /* determine influence coefficient at basic variable xB[i]
1368            in the explicitly specified column and turn to the case of
1369            increasing the variable x in order to simplify the program
1370            logic */
1371         alfa = (dir > 0 ? + val[t] : - val[t]);
1372         /* analyze main cases */
1373         if (type == GLP_FR)
1374         {  /* xB[i] is free variable */
1375            continue;
1376         }
1377         else if (type == GLP_LO)
1378lo:      {  /* xB[i] has an lower bound */
1379            if (alfa > - eps) continue;
1380            temp = (lb - beta) / alfa;
1381         }
1382         else if (type == GLP_UP)
1383up:      {  /* xB[i] has an upper bound */
1384            if (alfa < + eps) continue;
1385            temp = (ub - beta) / alfa;
1386         }
1387         else if (type == GLP_DB)
1388         {  /* xB[i] has both lower and upper bounds */
1389            if (alfa < 0.0) goto lo; else goto up;
1390         }
1391         else if (type == GLP_FX)
1392         {  /* xB[i] is fixed variable */
1393            if (- eps < alfa && alfa < + eps) continue;
1394            temp = 0.0;
1395         }
1396         else
1397            xassert(type != type);
1398         /* if the value of the variable xB[i] violates its lower or
1399            upper bound (slightly, because the current basis is assumed
1400            to be primal feasible), temp is negative; we can think this
1401            happens due to round-off errors and the value is exactly on
1402            the bound; this allows replacing temp by zero */
1403         if (temp < 0.0) temp = 0.0;
1404         /* apply the minimal ratio test */
1405         if (teta > temp || teta == temp && big < fabs(alfa))
1406            piv = t, teta = temp, big = fabs(alfa);
1407      }
1408      /* return index of the pivot element chosen */
1409      return piv;
1410}
1411
1412/***********************************************************************
1413*  NAME
1414*
1415*  glp_dual_rtest - perform dual ratio test
1416*
1417*  SYNOPSIS
1418*
1419*  int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1420*     const double val[], int dir, double eps);
1421*
1422*  DESCRIPTION
1423*
1424*  The routine glp_dual_rtest performs the dual ratio test using an
1425*  explicitly specified row of the simplex table.
1426*
1427*  The current basic solution associated with the LP problem object
1428*  must be dual feasible.
1429*
1430*  The explicitly specified row of the simplex table is a linear form
1431*  that shows how some basic variable x (which is not necessarily
1432*  presented in the problem object) depends on non-basic variables xN:
1433*
1434*     x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
1435*
1436*  The row (*) is specified on entry to the routine using the sparse
1437*  format. Ordinal numbers of non-basic variables xN[j] should be placed
1438*  in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
1439*  denote auxiliary variables, and ordinal numbers m+1 to m+n denote
1440*  structural variables. The corresponding non-zero coefficients alfa[j]
1441*  should be placed in locations val[1], ..., val[len]. The arrays ind
1442*  and val are not changed on exit.
1443*
1444*  The parameter dir specifies direction in which the variable x changes
1445*  on leaving the basis: +1 means that x goes to its lower bound, and -1
1446*  means that x goes to its upper bound.
1447*
1448*  The parameter eps is an absolute tolerance (small positive number)
1449*  used by the routine to skip small alfa[j] of the row (*).
1450*
1451*  The routine determines which non-basic variable (among specified in
1452*  ind[1], ..., ind[len]) should enter the basis in order to keep dual
1453*  feasibility.
1454*
1455*  RETURNS
1456*
1457*  The routine glp_dual_rtest returns the index piv in the arrays ind
1458*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
1459*  If the adjacent basic solution is dual unbounded and therefore the
1460*  choice cannot be made, the routine returns zero.
1461*
1462*  COMMENTS
1463*
1464*  If the basic variable x is presented in the LP problem object, the
1465*  row (*) can be computed with the routine glp_eval_tab_row; otherwise
1466*  it can be computed with the routine glp_transform_row. */
1467
1468int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1469      const double val[], int dir, double eps)
1470{     int k, m, n, piv, t, stat;
1471      double alfa, big, cost, obj, temp, teta;
1472      if (glp_get_dual_stat(P) != GLP_FEAS)
1473         xerror("glp_dual_rtest: basic solution is not dual feasible\n")
1474            ;
1475      if (!(dir == +1 || dir == -1))
1476         xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
1477      if (!(0.0 < eps && eps < 1.0))
1478         xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
1479      m = glp_get_num_rows(P);
1480      n = glp_get_num_cols(P);
1481      /* take into account optimization direction */
1482      obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
1483      /* initial settings */
1484      piv = 0, teta = DBL_MAX, big = 0.0;
1485      /* walk through the entries of the specified row */
1486      for (t = 1; t <= len; t++)
1487      {  /* get ordinal number of non-basic variable */
1488         k = ind[t];
1489         if (!(1 <= k && k <= m+n))
1490            xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
1491               "f range\n", t, k);
1492         /* determine status and reduced cost of non-basic variable
1493            x[k] = xN[j] in the current basic solution */
1494         if (k <= m)
1495         {  stat = glp_get_row_stat(P, k);
1496            cost = glp_get_row_dual(P, k);
1497         }
1498         else
1499         {  stat = glp_get_col_stat(P, k-m);
1500            cost = glp_get_col_dual(P, k-m);
1501         }
1502         if (stat == GLP_BS)
1503            xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
1504               "lowed\n", t, k);
1505         /* determine influence coefficient at non-basic variable xN[j]
1506            in the explicitly specified row and turn to the case of
1507            increasing the variable x in order to simplify the program
1508            logic */
1509         alfa = (dir > 0 ? + val[t] : - val[t]);
1510         /* analyze main cases */
1511         if (stat == GLP_NL)
1512         {  /* xN[j] is on its lower bound */
1513            if (alfa < + eps) continue;
1514            temp = (obj * cost) / alfa;
1515         }
1516         else if (stat == GLP_NU)
1517         {  /* xN[j] is on its upper bound */
1518            if (alfa > - eps) continue;
1519            temp = (obj * cost) / alfa;
1520         }
1521         else if (stat == GLP_NF)
1522         {  /* xN[j] is non-basic free variable */
1523            if (- eps < alfa && alfa < + eps) continue;
1524            temp = 0.0;
1525         }
1526         else if (stat == GLP_NS)
1527         {  /* xN[j] is non-basic fixed variable */
1528            continue;
1529         }
1530         else
1531            xassert(stat != stat);
1532         /* if the reduced cost of the variable xN[j] violates its zero
1533            bound (slightly, because the current basis is assumed to be
1534            dual feasible), temp is negative; we can think this happens
1535            due to round-off errors and the reduced cost is exact zero;
1536            this allows replacing temp by zero */
1537         if (temp < 0.0) temp = 0.0;
1538         /* apply the minimal ratio test */
1539         if (teta > temp || teta == temp && big < fabs(alfa))
1540            piv = t, teta = temp, big = fabs(alfa);
1541      }
1542      /* return index of the pivot element chosen */
1543      return piv;
1544}
1545
1546/***********************************************************************
1547*  NAME
1548*
1549*  glp_analyze_row - simulate one iteration of dual simplex method
1550*
1551*  SYNOPSIS
1552*
1553*  int glp_analyze_row(glp_prob *P, int len, const int ind[],
1554*     const double val[], int type, double rhs, double eps, int *piv,
1555*     double *x, double *dx, double *y, double *dy, double *dz);
1556*
1557*  DESCRIPTION
1558*
1559*  Let the current basis be optimal or dual feasible, and there be
1560*  specified a row (constraint), which is violated by the current basic
1561*  solution. The routine glp_analyze_row simulates one iteration of the
1562*  dual simplex method to determine some information on the adjacent
1563*  basis (see below), where the specified row becomes active constraint
1564*  (i.e. its auxiliary variable becomes non-basic).
1565*
1566*  The current basic solution associated with the problem object passed
1567*  to the routine must be dual feasible, and its primal components must
1568*  be defined.
1569*
1570*  The row to be analyzed must be previously transformed either with
1571*  the routine glp_eval_tab_row (if the row is in the problem object)
1572*  or with the routine glp_transform_row (if the row is external, i.e.
1573*  not in the problem object). This is needed to express the row only
1574*  through (auxiliary and structural) variables, which are non-basic in
1575*  the current basis:
1576*
1577*     y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
1578*
1579*  where y is an auxiliary variable of the row, alfa[j] is an influence
1580*  coefficient, xN[j] is a non-basic variable.
1581*
1582*  The row is passed to the routine in sparse format. Ordinal numbers
1583*  of non-basic variables are stored in locations ind[1], ..., ind[len],
1584*  where numbers 1 to m denote auxiliary variables while numbers m+1 to
1585*  m+n denote structural variables. Corresponding non-zero coefficients
1586*  alfa[j] are stored in locations val[1], ..., val[len]. The arrays
1587*  ind and val are ot changed on exit.
1588*
1589*  The parameters type and rhs specify the row type and its right-hand
1590*  side as follows:
1591*
1592*     type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
1593*
1594*     type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
1595*
1596*  The parameter eps is an absolute tolerance (small positive number)
1597*  used by the routine to skip small coefficients alfa[j] on performing
1598*  the dual ratio test.
1599*
1600*  If the operation was successful, the routine stores the following
1601*  information to corresponding location (if some parameter is NULL,
1602*  its value is not stored):
1603*
1604*  piv   index in the array ind and val, 1 <= piv <= len, determining
1605*        the non-basic variable, which would enter the adjacent basis;
1606*
1607*  x     value of the non-basic variable in the current basis;
1608*
1609*  dx    difference between values of the non-basic variable in the
1610*        adjacent and current bases, dx = x.new - x.old;
1611*
1612*  y     value of the row (i.e. of its auxiliary variable) in the
1613*        current basis;
1614*
1615*  dy    difference between values of the row in the adjacent and
1616*        current bases, dy = y.new - y.old;
1617*
1618*  dz    difference between values of the objective function in the
1619*        adjacent and current bases, dz = z.new - z.old. Note that in
1620*        case of minimization dz >= 0, and in case of maximization
1621*        dz <= 0, i.e. in the adjacent basis the objective function
1622*        always gets worse (degrades). */
1623
1624int _glp_analyze_row(glp_prob *P, int len, const int ind[],
1625      const double val[], int type, double rhs, double eps, int *_piv,
1626      double *_x, double *_dx, double *_y, double *_dy, double *_dz)
1627{     int t, k, dir, piv, ret = 0;
1628      double x, dx, y, dy, dz;
1629      if (P->pbs_stat == GLP_UNDEF)
1630         xerror("glp_analyze_row: primal basic solution components are "
1631            "undefined\n");
1632      if (P->dbs_stat != GLP_FEAS)
1633         xerror("glp_analyze_row: basic solution is not dual feasible\n"
1634            );
1635      /* compute the row value y = sum alfa[j] * xN[j] in the current
1636         basis */
1637      if (!(0 <= len && len <= P->n))
1638         xerror("glp_analyze_row: len = %d; invalid row length\n", len);
1639      y = 0.0;
1640      for (t = 1; t <= len; t++)
1641      {  /* determine value of x[k] = xN[j] in the current basis */
1642         k = ind[t];
1643         if (!(1 <= k && k <= P->m+P->n))
1644            xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
1645               " of range\n", t, k);
1646         if (k <= P->m)
1647         {  /* x[k] is auxiliary variable */
1648            if (P->row[k]->stat == GLP_BS)
1649               xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
1650                  "ariable is not allowed\n", t, k);
1651            x = P->row[k]->prim;
1652         }
1653         else
1654         {  /* x[k] is structural variable */
1655            if (P->col[k-P->m]->stat == GLP_BS)
1656               xerror("glp_analyze_row: ind[%d] = %d; basic structural "
1657                  "variable is not allowed\n", t, k);
1658            x = P->col[k-P->m]->prim;
1659         }
1660         y += val[t] * x;
1661      }
1662      /* check if the row is primal infeasible in the current basis,
1663         i.e. the constraint is violated at the current point */
1664      if (type == GLP_LO)
1665      {  if (y >= rhs)
1666         {  /* the constraint is not violated */
1667            ret = 1;
1668            goto done;
1669         }
1670         /* in the adjacent basis y goes to its lower bound */
1671         dir = +1;
1672      }
1673      else if (type == GLP_UP)
1674      {  if (y <= rhs)
1675         {  /* the constraint is not violated */
1676            ret = 1;
1677            goto done;
1678         }
1679         /* in the adjacent basis y goes to its upper bound */
1680         dir = -1;
1681      }
1682      else
1683         xerror("glp_analyze_row: type = %d; invalid parameter\n",
1684            type);
1685      /* compute dy = y.new - y.old */
1686      dy = rhs - y;
1687      /* perform dual ratio test to determine which non-basic variable
1688         should enter the adjacent basis to keep it dual feasible */
1689      piv = glp_dual_rtest(P, len, ind, val, dir, eps);
1690      if (piv == 0)
1691      {  /* no dual feasible adjacent basis exists */
1692         ret = 2;
1693         goto done;
1694      }
1695      /* non-basic variable x[k] = xN[j] should enter the basis */
1696      k = ind[piv];
1697      xassert(1 <= k && k <= P->m+P->n);
1698      /* determine its value in the current basis */
1699      if (k <= P->m)
1700         x = P->row[k]->prim;
1701      else
1702         x = P->col[k-P->m]->prim;
1703      /* compute dx = x.new - x.old = dy / alfa[j] */
1704      xassert(val[piv] != 0.0);
1705      dx = dy / val[piv];
1706      /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
1707         cost of xN[j] in the current basis */
1708      if (k <= P->m)
1709         dz = P->row[k]->dual * dx;
1710      else
1711         dz = P->col[k-P->m]->dual * dx;
1712      /* store the analysis results */
1713      if (_piv != NULL) *_piv = piv;
1714      if (_x   != NULL) *_x   = x;
1715      if (_dx  != NULL) *_dx  = dx;
1716      if (_y   != NULL) *_y   = y;
1717      if (_dy  != NULL) *_dy  = dy;
1718      if (_dz  != NULL) *_dz  = dz;
1719done: return ret;
1720}
1721
1722#if 0
1723int main(void)
1724{     /* example program for the routine glp_analyze_row */
1725      glp_prob *P;
1726      glp_smcp parm;
1727      int i, k, len, piv, ret, ind[1+100];
1728      double rhs, x, dx, y, dy, dz, val[1+100];
1729      P = glp_create_prob();
1730      /* read plan.mps (see glpk/examples) */
1731      ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
1732      glp_assert(ret == 0);
1733      /* and solve it to optimality */
1734      ret = glp_simplex(P, NULL);
1735      glp_assert(ret == 0);
1736      glp_assert(glp_get_status(P) == GLP_OPT);
1737      /* the optimal objective value is 296.217 */
1738      /* we would like to know what happens if we would add a new row
1739         (constraint) to plan.mps:
1740         .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
1741      /* first, we specify this new row */
1742      glp_create_index(P);
1743      len = 0;
1744      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1745      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1746      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1747      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1748      rhs = 12;
1749      /* then we can compute value of the row (i.e. of its auxiliary
1750         variable) in the current basis to see if the constraint is
1751         violated */
1752      y = 0.0;
1753      for (k = 1; k <= len; k++)
1754         y += val[k] * glp_get_col_prim(P, ind[k]);
1755      glp_printf("y = %g\n", y);
1756      /* this prints y = 15.1372, so the constraint is violated, since
1757         we require that y <= rhs = 12 */
1758      /* now we transform the row to express it only through non-basic
1759         (auxiliary and artificial) variables */
1760      len = glp_transform_row(P, len, ind, val);
1761      /* finally, we simulate one step of the dual simplex method to
1762         obtain necessary information for the adjacent basis */
1763      ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
1764         &x, &dx, &y, &dy, &dz);
1765      glp_assert(ret == 0);
1766      glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
1767         ind[piv], x, dx, y, dy, dz);
1768      /* this prints dz = 5.64418 and means that in the adjacent basis
1769         the objective function would be 296.217 + 5.64418 = 301.861 */
1770      /* now we actually include the row into the problem object; note
1771         that the arrays ind and val are clobbered, so we need to build
1772         them once again */
1773      len = 0;
1774      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1775      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1776      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1777      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1778      rhs = 12;
1779      i = glp_add_rows(P, 1);
1780      glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
1781      glp_set_mat_row(P, i, len, ind, val);
1782      /* and perform one dual simplex iteration */
1783      glp_init_smcp(&parm);
1784      parm.meth = GLP_DUAL;
1785      parm.it_lim = 1;
1786      glp_simplex(P, &parm);
1787      /* the current objective value is 301.861 */
1788      return 0;
1789}
1790#endif
1791
1792/***********************************************************************
1793*  NAME
1794*
1795*  glp_analyze_bound - analyze active bound of non-basic variable
1796*
1797*  SYNOPSIS
1798*
1799*  void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
1800*     double *limit2, int *var2);
1801*
1802*  DESCRIPTION
1803*
1804*  The routine glp_analyze_bound analyzes the effect of varying the
1805*  active bound of specified non-basic variable.
1806*
1807*  The non-basic variable is specified by the parameter k, where
1808*  1 <= k <= m means auxiliary variable of corresponding row while
1809*  m+1 <= k <= m+n means structural variable (column).
1810*
1811*  Note that the current basic solution must be optimal, and the basis
1812*  factorization must exist.
1813*
1814*  Results of the analysis have the following meaning.
1815*
1816*  value1 is the minimal value of the active bound, at which the basis
1817*  still remains primal feasible and thus optimal. -DBL_MAX means that
1818*  the active bound has no lower limit.
1819*
1820*  var1 is the ordinal number of an auxiliary (1 to m) or structural
1821*  (m+1 to n) basic variable, which reaches its bound first and thereby
1822*  limits further decreasing the active bound being analyzed.
1823*  if value1 = -DBL_MAX, var1 is set to 0.
1824*
1825*  value2 is the maximal value of the active bound, at which the basis
1826*  still remains primal feasible and thus optimal. +DBL_MAX means that
1827*  the active bound has no upper limit.
1828*
1829*  var2 is the ordinal number of an auxiliary (1 to m) or structural
1830*  (m+1 to n) basic variable, which reaches its bound first and thereby
1831*  limits further increasing the active bound being analyzed.
1832*  if value2 = +DBL_MAX, var2 is set to 0. */
1833
1834void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
1835      double *value2, int *var2)
1836{     GLPROW *row;
1837      GLPCOL *col;
1838      int m, n, stat, kase, p, len, piv, *ind;
1839      double x, new_x, ll, uu, xx, delta, *val;
1840      /* sanity checks */
1841      if (P == NULL || P->magic != GLP_PROB_MAGIC)
1842         xerror("glp_analyze_bound: P = %p; invalid problem object\n",
1843            P);
1844      m = P->m, n = P->n;
1845      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
1846         xerror("glp_analyze_bound: optimal basic solution required\n");
1847      if (!(m == 0 || P->valid))
1848         xerror("glp_analyze_bound: basis factorization required\n");
1849      if (!(1 <= k && k <= m+n))
1850         xerror("glp_analyze_bound: k = %d; variable number out of rang"
1851            "e\n", k);
1852      /* retrieve information about the specified non-basic variable
1853         x[k] whose active bound is to be analyzed */
1854      if (k <= m)
1855      {  row = P->row[k];
1856         stat = row->stat;
1857         x = row->prim;
1858      }
1859      else
1860      {  col = P->col[k-m];
1861         stat = col->stat;
1862         x = col->prim;
1863      }
1864      if (stat == GLP_BS)
1865         xerror("glp_analyze_bound: k = %d; basic variable not allowed "
1866            "\n", k);
1867      /* allocate working arrays */
1868      ind = xcalloc(1+m, sizeof(int));
1869      val = xcalloc(1+m, sizeof(double));
1870      /* compute column of the simplex table corresponding to the
1871         non-basic variable x[k] */
1872      len = glp_eval_tab_col(P, k, ind, val);
1873      xassert(0 <= len && len <= m);
1874      /* perform analysis */
1875      for (kase = -1; kase <= +1; kase += 2)
1876      {  /* kase < 0 means active bound of x[k] is decreasing;
1877            kase > 0 means active bound of x[k] is increasing */
1878         /* use the primal ratio test to determine some basic variable
1879            x[p] which reaches its bound first */
1880         piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
1881         if (piv == 0)
1882         {  /* nothing limits changing the active bound of x[k] */
1883            p = 0;
1884            new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
1885            goto store;
1886         }
1887         /* basic variable x[p] limits changing the active bound of
1888            x[k]; determine its value in the current basis */
1889         xassert(1 <= piv && piv <= len);
1890         p = ind[piv];
1891         if (p <= m)
1892         {  row = P->row[p];
1893            ll = glp_get_row_lb(P, row->i);
1894            uu = glp_get_row_ub(P, row->i);
1895            stat = row->stat;
1896            xx = row->prim;
1897         }
1898         else
1899         {  col = P->col[p-m];
1900            ll = glp_get_col_lb(P, col->j);
1901            uu = glp_get_col_ub(P, col->j);
1902            stat = col->stat;
1903            xx = col->prim;
1904         }
1905         xassert(stat == GLP_BS);
1906         /* determine delta x[p] = bound of x[p] - value of x[p] */
1907         if (kase < 0 && val[piv] > 0.0 ||
1908             kase > 0 && val[piv] < 0.0)
1909         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
1910            xassert(ll != -DBL_MAX);
1911            delta = ll - xx;
1912         }
1913         else
1914         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
1915            xassert(uu != +DBL_MAX);
1916            delta = uu - xx;
1917         }
1918         /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
1919            delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
1920            x[k] in the adjacent basis */
1921         xassert(val[piv] != 0.0);
1922         new_x = x + delta / val[piv];
1923store:   /* store analysis results */
1924         if (kase < 0)
1925         {  if (value1 != NULL) *value1 = new_x;
1926            if (var1 != NULL) *var1 = p;
1927         }
1928         else
1929         {  if (value2 != NULL) *value2 = new_x;
1930            if (var2 != NULL) *var2 = p;
1931         }
1932      }
1933      /* free working arrays */
1934      xfree(ind);
1935      xfree(val);
1936      return;
1937}
1938
1939/***********************************************************************
1940*  NAME
1941*
1942*  glp_analyze_coef - analyze objective coefficient at basic variable
1943*
1944*  SYNOPSIS
1945*
1946*  void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1947*     double *value1, double *coef2, int *var2, double *value2);
1948*
1949*  DESCRIPTION
1950*
1951*  The routine glp_analyze_coef analyzes the effect of varying the
1952*  objective coefficient at specified basic variable.
1953*
1954*  The basic variable is specified by the parameter k, where
1955*  1 <= k <= m means auxiliary variable of corresponding row while
1956*  m+1 <= k <= m+n means structural variable (column).
1957*
1958*  Note that the current basic solution must be optimal, and the basis
1959*  factorization must exist.
1960*
1961*  Results of the analysis have the following meaning.
1962*
1963*  coef1 is the minimal value of the objective coefficient, at which
1964*  the basis still remains dual feasible and thus optimal. -DBL_MAX
1965*  means that the objective coefficient has no lower limit.
1966*
1967*  var1 is the ordinal number of an auxiliary (1 to m) or structural
1968*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1969*  bound first and thereby limits further decreasing the objective
1970*  coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
1971*
1972*  value1 is value of the basic variable being analyzed in an adjacent
1973*  basis, which is defined as follows. Let the objective coefficient
1974*  reaches its minimal value (coef1) and continues decreasing. Then the
1975*  reduced cost of the limiting non-basic variable (var1) becomes dual
1976*  infeasible and the current basis becomes non-optimal that forces the
1977*  limiting non-basic variable to enter the basis replacing there some
1978*  basic variable that leaves the basis to keep primal feasibility.
1979*  Should note that on determining the adjacent basis current bounds
1980*  of the basic variable being analyzed are ignored as if it were free
1981*  (unbounded) variable, so it cannot leave the basis. It may happen
1982*  that no dual feasible adjacent basis exists, in which case value1 is
1983*  set to -DBL_MAX or +DBL_MAX.
1984*
1985*  coef2 is the maximal value of the objective coefficient, at which
1986*  the basis still remains dual feasible and thus optimal. +DBL_MAX
1987*  means that the objective coefficient has no upper limit.
1988*
1989*  var2 is the ordinal number of an auxiliary (1 to m) or structural
1990*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1991*  bound first and thereby limits further increasing the objective
1992*  coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
1993*
1994*  value2 is value of the basic variable being analyzed in an adjacent
1995*  basis, which is defined exactly in the same way as value1 above with
1996*  exception that now the objective coefficient is increasing. */
1997
1998void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1999      double *value1, double *coef2, int *var2, double *value2)
2000{     GLPROW *row; GLPCOL *col;
2001      int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
2002         *cind, *rind;
2003      double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
2004         *rval, *cval;
2005      /* sanity checks */
2006      if (P == NULL || P->magic != GLP_PROB_MAGIC)
2007         xerror("glp_analyze_coef: P = %p; invalid problem object\n",
2008            P);
2009      m = P->m, n = P->n;
2010      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
2011         xerror("glp_analyze_coef: optimal basic solution required\n");
2012      if (!(m == 0 || P->valid))
2013         xerror("glp_analyze_coef: basis factorization required\n");
2014      if (!(1 <= k && k <= m+n))
2015         xerror("glp_analyze_coef: k = %d; variable number out of range"
2016            "\n", k);
2017      /* retrieve information about the specified basic variable x[k]
2018         whose objective coefficient c[k] is to be analyzed */
2019      if (k <= m)
2020      {  row = P->row[k];
2021         type = row->type;
2022         lb = row->lb;
2023         ub = row->ub;
2024         coef = 0.0;
2025         stat = row->stat;
2026         x = row->prim;
2027      }
2028      else
2029      {  col = P->col[k-m];
2030         type = col->type;
2031         lb = col->lb;
2032         ub = col->ub;
2033         coef = col->coef;
2034         stat = col->stat;
2035         x = col->prim;
2036      }
2037      if (stat != GLP_BS)
2038         xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
2039            "ed\n", k);
2040      /* allocate working arrays */
2041      cind = xcalloc(1+m, sizeof(int));
2042      cval = xcalloc(1+m, sizeof(double));
2043      rind = xcalloc(1+n, sizeof(int));
2044      rval = xcalloc(1+n, sizeof(double));
2045      /* compute row of the simplex table corresponding to the basic
2046         variable x[k] */
2047      rlen = glp_eval_tab_row(P, k, rind, rval);
2048      xassert(0 <= rlen && rlen <= n);
2049      /* perform analysis */
2050      for (kase = -1; kase <= +1; kase += 2)
2051      {  /* kase < 0 means objective coefficient c[k] is decreasing;
2052            kase > 0 means objective coefficient c[k] is increasing */
2053         /* note that decreasing c[k] is equivalent to increasing dual
2054            variable lambda[k] and vice versa; we need to correctly set
2055            the dir flag as required by the routine glp_dual_rtest */
2056         if (P->dir == GLP_MIN)
2057            dir = - kase;
2058         else if (P->dir == GLP_MAX)
2059            dir = + kase;
2060         else
2061            xassert(P != P);
2062         /* use the dual ratio test to determine non-basic variable
2063            x[q] whose reduced cost d[q] reaches zero bound first */
2064         rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
2065         if (rpiv == 0)
2066         {  /* nothing limits changing c[k] */
2067            lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
2068            q = 0;
2069            /* x[k] keeps its current value */
2070            new_x = x;
2071            goto store;
2072         }
2073         /* non-basic variable x[q] limits changing coefficient c[k];
2074            determine its status and reduced cost d[k] in the current
2075            basis */
2076         xassert(1 <= rpiv && rpiv <= rlen);
2077         q = rind[rpiv];
2078         xassert(1 <= q && q <= m+n);
2079         if (q <= m)
2080         {  row = P->row[q];
2081            stat = row->stat;
2082            d = row->dual;
2083         }
2084         else
2085         {  col = P->col[q-m];
2086            stat = col->stat;
2087            d = col->dual;
2088         }
2089         /* note that delta d[q] = new d[q] - d[q] = - d[q], because
2090            new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
2091            delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
2092         xassert(rval[rpiv] != 0.0);
2093         delta = - d / rval[rpiv];
2094         /* compute new c[k] = c[k] + delta c[k], which is the limiting
2095            value of the objective coefficient c[k] */
2096         lim_coef = coef + delta;
2097         /* let c[k] continue decreasing/increasing that makes d[q]
2098            dual infeasible and forces x[q] to enter the basis;
2099            to perform the primal ratio test we need to know in which
2100            direction x[q] changes on entering the basis; we determine
2101            that analyzing the sign of delta d[q] (see above), since
2102            d[q] may be close to zero having wrong sign */
2103         /* let, for simplicity, the problem is minimization */
2104         if (kase < 0 && rval[rpiv] > 0.0 ||
2105             kase > 0 && rval[rpiv] < 0.0)
2106         {  /* delta d[q] < 0, so d[q] being non-negative will become
2107               negative, so x[q] will increase */
2108            dir = +1;
2109         }
2110         else
2111         {  /* delta d[q] > 0, so d[q] being non-positive will become
2112               positive, so x[q] will decrease */
2113            dir = -1;
2114         }
2115         /* if the problem is maximization, correct the direction */
2116         if (P->dir == GLP_MAX) dir = - dir;
2117         /* check that we didn't make a silly mistake */
2118         if (dir > 0)
2119            xassert(stat == GLP_NL || stat == GLP_NF);
2120         else
2121            xassert(stat == GLP_NU || stat == GLP_NF);
2122         /* compute column of the simplex table corresponding to the
2123            non-basic variable x[q] */
2124         clen = glp_eval_tab_col(P, q, cind, cval);
2125         /* make x[k] temporarily free (unbounded) */
2126         if (k <= m)
2127         {  row = P->row[k];
2128            row->type = GLP_FR;
2129            row->lb = row->ub = 0.0;
2130         }
2131         else
2132         {  col = P->col[k-m];
2133            col->type = GLP_FR;
2134            col->lb = col->ub = 0.0;
2135         }
2136         /* use the primal ratio test to determine some basic variable
2137            which leaves the basis */
2138         cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
2139         /* restore original bounds of the basic variable x[k] */
2140         if (k <= m)
2141         {  row = P->row[k];
2142            row->type = type;
2143            row->lb = lb, row->ub = ub;
2144         }
2145         else
2146         {  col = P->col[k-m];
2147            col->type = type;
2148            col->lb = lb, col->ub = ub;
2149         }
2150         if (cpiv == 0)
2151         {  /* non-basic variable x[q] can change unlimitedly */
2152            if (dir < 0 && rval[rpiv] > 0.0 ||
2153                dir > 0 && rval[rpiv] < 0.0)
2154            {  /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
2155               new_x = -DBL_MAX;
2156            }
2157            else
2158            {  /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
2159               new_x = +DBL_MAX;
2160            }
2161            goto store;
2162         }
2163         /* some basic variable x[p] limits changing non-basic variable
2164            x[q] in the adjacent basis */
2165         xassert(1 <= cpiv && cpiv <= clen);
2166         p = cind[cpiv];
2167         xassert(1 <= p && p <= m+n);
2168         xassert(p != k);
2169         if (p <= m)
2170         {  row = P->row[p];
2171            xassert(row->stat == GLP_BS);
2172            ll = glp_get_row_lb(P, row->i);
2173            uu = glp_get_row_ub(P, row->i);
2174            xx = row->prim;
2175         }
2176         else
2177         {  col = P->col[p-m];
2178            xassert(col->stat == GLP_BS);
2179            ll = glp_get_col_lb(P, col->j);
2180            uu = glp_get_col_ub(P, col->j);
2181            xx = col->prim;
2182         }
2183         /* determine delta x[p] = new x[p] - x[p] */
2184         if (dir < 0 && cval[cpiv] > 0.0 ||
2185             dir > 0 && cval[cpiv] < 0.0)
2186         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
2187            xassert(ll != -DBL_MAX);
2188            delta = ll - xx;
2189         }
2190         else
2191         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
2192            xassert(uu != +DBL_MAX);
2193            delta = uu - xx;
2194         }
2195         /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
2196            delta x[q] = delta x[p] / alfa[p,q] */
2197         xassert(cval[cpiv] != 0.0);
2198         new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
2199store:   /* store analysis results */
2200         if (kase < 0)
2201         {  if (coef1 != NULL) *coef1 = lim_coef;
2202            if (var1 != NULL) *var1 = q;
2203            if (value1 != NULL) *value1 = new_x;
2204         }
2205         else
2206         {  if (coef2 != NULL) *coef2 = lim_coef;
2207            if (var2 != NULL) *var2 = q;
2208            if (value2 != NULL) *value2 = new_x;
2209         }
2210      }
2211      /* free working arrays */
2212      xfree(cind);
2213      xfree(cval);
2214      xfree(rind);
2215      xfree(rval);
2216      return;
2217}
2218
2219/* eof */
Note: See TracBrowser for help on using the repository browser.