COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpssx01.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 27.1 KB
Line 
1/* glpssx01.c */
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 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 "glpenv.h"
26#include "glpssx.h"
27#define xfault xerror
28
29/*----------------------------------------------------------------------
30// ssx_create - create simplex solver workspace.
31//
32// This routine creates the workspace used by simplex solver routines,
33// and returns a pointer to it.
34//
35// Parameters m, n, and nnz specify, respectively, the number of rows,
36// columns, and non-zero constraint coefficients.
37//
38// This routine only allocates the memory for the workspace components,
39// so the workspace needs to be saturated by data. */
40
41SSX *ssx_create(int m, int n, int nnz)
42{     SSX *ssx;
43      int i, j, k;
44      if (m < 1)
45         xfault("ssx_create: m = %d; invalid number of rows\n", m);
46      if (n < 1)
47         xfault("ssx_create: n = %d; invalid number of columns\n", n);
48      if (nnz < 0)
49         xfault("ssx_create: nnz = %d; invalid number of non-zero const"
50            "raint coefficients\n", nnz);
51      ssx = xmalloc(sizeof(SSX));
52      ssx->m = m;
53      ssx->n = n;
54      ssx->type = xcalloc(1+m+n, sizeof(int));
55      ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));
56      for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);
57      ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));
58      for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);
59      ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));
60      for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);
61      ssx->A_ptr = xcalloc(1+n+1, sizeof(int));
62      ssx->A_ptr[n+1] = nnz+1;
63      ssx->A_ind = xcalloc(1+nnz, sizeof(int));
64      ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));
65      for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);
66      ssx->stat = xcalloc(1+m+n, sizeof(int));
67      ssx->Q_row = xcalloc(1+m+n, sizeof(int));
68      ssx->Q_col = xcalloc(1+m+n, sizeof(int));
69      ssx->binv = bfx_create_binv();
70      ssx->bbar = xcalloc(1+m, sizeof(mpq_t));
71      for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);
72      ssx->pi = xcalloc(1+m, sizeof(mpq_t));
73      for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);
74      ssx->cbar = xcalloc(1+n, sizeof(mpq_t));
75      for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);
76      ssx->rho = xcalloc(1+m, sizeof(mpq_t));
77      for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);
78      ssx->ap = xcalloc(1+n, sizeof(mpq_t));
79      for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);
80      ssx->aq = xcalloc(1+m, sizeof(mpq_t));
81      for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);
82      mpq_init(ssx->delta);
83      return ssx;
84}
85
86/*----------------------------------------------------------------------
87// ssx_factorize - factorize the current basis matrix.
88//
89// This routine computes factorization of the current basis matrix B
90// and returns the singularity flag. If the matrix B is non-singular,
91// the flag is zero, otherwise non-zero. */
92
93static int basis_col(void *info, int j, int ind[], mpq_t val[])
94{     /* this auxiliary routine provides row indices and numeric values
95         of non-zero elements in j-th column of the matrix B */
96      SSX *ssx = info;
97      int m = ssx->m;
98      int n = ssx->n;
99      int *A_ptr = ssx->A_ptr;
100      int *A_ind = ssx->A_ind;
101      mpq_t *A_val = ssx->A_val;
102      int *Q_col = ssx->Q_col;
103      int k, len, ptr;
104      xassert(1 <= j && j <= m);
105      k = Q_col[j]; /* x[k] = xB[j] */
106      xassert(1 <= k && k <= m+n);
107      /* j-th column of the matrix B is k-th column of the augmented
108         constraint matrix (I | -A) */
109      if (k <= m)
110      {  /* it is a column of the unity matrix I */
111         len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1);
112      }
113      else
114      {  /* it is a column of the original constraint matrix -A */
115         len = 0;
116         for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
117         {  len++;
118            ind[len] = A_ind[ptr];
119            mpq_neg(val[len], A_val[ptr]);
120         }
121      }
122      return len;
123}
124
125int ssx_factorize(SSX *ssx)
126{     int ret;
127      ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx);
128      return ret;
129}
130
131/*----------------------------------------------------------------------
132// ssx_get_xNj - determine value of non-basic variable.
133//
134// This routine determines the value of non-basic variable xN[j] in the
135// current basic solution defined as follows:
136//
137//    0,             if xN[j] is free variable
138//    lN[j],         if xN[j] is on its lower bound
139//    uN[j],         if xN[j] is on its upper bound
140//    lN[j] = uN[j], if xN[j] is fixed variable
141//
142// where lN[j] and uN[j] are lower and upper bounds of xN[j]. */
143
144void ssx_get_xNj(SSX *ssx, int j, mpq_t x)
145{     int m = ssx->m;
146      int n = ssx->n;
147      mpq_t *lb = ssx->lb;
148      mpq_t *ub = ssx->ub;
149      int *stat = ssx->stat;
150      int *Q_col = ssx->Q_col;
151      int k;
152      xassert(1 <= j && j <= n);
153      k = Q_col[m+j]; /* x[k] = xN[j] */
154      xassert(1 <= k && k <= m+n);
155      switch (stat[k])
156      {  case SSX_NL:
157            /* xN[j] is on its lower bound */
158            mpq_set(x, lb[k]); break;
159         case SSX_NU:
160            /* xN[j] is on its upper bound */
161            mpq_set(x, ub[k]); break;
162         case SSX_NF:
163            /* xN[j] is free variable */
164            mpq_set_si(x, 0, 1); break;
165         case SSX_NS:
166            /* xN[j] is fixed variable */
167            mpq_set(x, lb[k]); break;
168         default:
169            xassert(stat != stat);
170      }
171      return;
172}
173
174/*----------------------------------------------------------------------
175// ssx_eval_bbar - compute values of basic variables.
176//
177// This routine computes values of basic variables xB in the current
178// basic solution as follows:
179//
180//    beta = - inv(B) * N * xN,
181//
182// where B is the basis matrix, N is the matrix of non-basic columns,
183// xN is a vector of current values of non-basic variables. */
184
185void ssx_eval_bbar(SSX *ssx)
186{     int m = ssx->m;
187      int n = ssx->n;
188      mpq_t *coef = ssx->coef;
189      int *A_ptr = ssx->A_ptr;
190      int *A_ind = ssx->A_ind;
191      mpq_t *A_val = ssx->A_val;
192      int *Q_col = ssx->Q_col;
193      mpq_t *bbar = ssx->bbar;
194      int i, j, k, ptr;
195      mpq_t x, temp;
196      mpq_init(x);
197      mpq_init(temp);
198      /* bbar := 0 */
199      for (i = 1; i <= m; i++)
200         mpq_set_si(bbar[i], 0, 1);
201      /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */
202      for (j = 1; j <= n; j++)
203      {  ssx_get_xNj(ssx, j, x);
204         if (mpq_sgn(x) == 0) continue;
205         k = Q_col[m+j]; /* x[k] = xN[j] */
206         if (k <= m)
207         {  /* N[j] is a column of the unity matrix I */
208            mpq_sub(bbar[k], bbar[k], x);
209         }
210         else
211         {  /* N[j] is a column of the original constraint matrix -A */
212            for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
213            {  mpq_mul(temp, A_val[ptr], x);
214               mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp);
215            }
216         }
217      }
218      /* bbar := inv(B) * bbar */
219      bfx_ftran(ssx->binv, bbar, 0);
220#if 1
221      /* compute value of the objective function */
222      /* bbar[0] := c[0] */
223      mpq_set(bbar[0], coef[0]);
224      /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */
225      for (i = 1; i <= m; i++)
226      {  k = Q_col[i]; /* x[k] = xB[i] */
227         if (mpq_sgn(coef[k]) == 0) continue;
228         mpq_mul(temp, coef[k], bbar[i]);
229         mpq_add(bbar[0], bbar[0], temp);
230      }
231      /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */
232      for (j = 1; j <= n; j++)
233      {  k = Q_col[m+j]; /* x[k] = xN[j] */
234         if (mpq_sgn(coef[k]) == 0) continue;
235         ssx_get_xNj(ssx, j, x);
236         mpq_mul(temp, coef[k], x);
237         mpq_add(bbar[0], bbar[0], temp);
238      }
239#endif
240      mpq_clear(x);
241      mpq_clear(temp);
242      return;
243}
244
245/*----------------------------------------------------------------------
246// ssx_eval_pi - compute values of simplex multipliers.
247//
248// This routine computes values of simplex multipliers (shadow prices)
249// pi in the current basic solution as follows:
250//
251//    pi = inv(B') * cB,
252//
253// where B' is a matrix transposed to the basis matrix B, cB is a vector
254// of objective coefficients at basic variables xB. */
255
256void ssx_eval_pi(SSX *ssx)
257{     int m = ssx->m;
258      mpq_t *coef = ssx->coef;
259      int *Q_col = ssx->Q_col;
260      mpq_t *pi = ssx->pi;
261      int i;
262      /* pi := cB */
263      for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]);
264      /* pi := inv(B') * cB */
265      bfx_btran(ssx->binv, pi);
266      return;
267}
268
269/*----------------------------------------------------------------------
270// ssx_eval_dj - compute reduced cost of non-basic variable.
271//
272// This routine computes reduced cost d[j] of non-basic variable xN[j]
273// in the current basic solution as follows:
274//
275//    d[j] = cN[j] - N[j] * pi,
276//
277// where cN[j] is an objective coefficient at xN[j], N[j] is a column
278// of the augmented constraint matrix (I | -A) corresponding to xN[j],
279// pi is the vector of simplex multipliers (shadow prices). */
280
281void ssx_eval_dj(SSX *ssx, int j, mpq_t dj)
282{     int m = ssx->m;
283      int n = ssx->n;
284      mpq_t *coef = ssx->coef;
285      int *A_ptr = ssx->A_ptr;
286      int *A_ind = ssx->A_ind;
287      mpq_t *A_val = ssx->A_val;
288      int *Q_col = ssx->Q_col;
289      mpq_t *pi = ssx->pi;
290      int k, ptr, end;
291      mpq_t temp;
292      mpq_init(temp);
293      xassert(1 <= j && j <= n);
294      k = Q_col[m+j]; /* x[k] = xN[j] */
295      xassert(1 <= k && k <= m+n);
296      /* j-th column of the matrix N is k-th column of the augmented
297         constraint matrix (I | -A) */
298      if (k <= m)
299      {  /* it is a column of the unity matrix I */
300         mpq_sub(dj, coef[k], pi[k]);
301      }
302      else
303      {  /* it is a column of the original constraint matrix -A */
304         mpq_set(dj, coef[k]);
305         for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++)
306         {  mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]);
307            mpq_add(dj, dj, temp);
308         }
309      }
310      mpq_clear(temp);
311      return;
312}
313
314/*----------------------------------------------------------------------
315// ssx_eval_cbar - compute reduced costs of all non-basic variables.
316//
317// This routine computes the vector of reduced costs pi in the current
318// basic solution for all non-basic variables, including fixed ones. */
319
320void ssx_eval_cbar(SSX *ssx)
321{     int n = ssx->n;
322      mpq_t *cbar = ssx->cbar;
323      int j;
324      for (j = 1; j <= n; j++)
325         ssx_eval_dj(ssx, j, cbar[j]);
326      return;
327}
328
329/*----------------------------------------------------------------------
330// ssx_eval_rho - compute p-th row of the inverse.
331//
332// This routine computes p-th row of the matrix inv(B), where B is the
333// current basis matrix.
334//
335// p-th row of the inverse is computed using the following formula:
336//
337//    rho = inv(B') * e[p],
338//
339// where B' is a matrix transposed to B, e[p] is a unity vector, which
340// contains one in p-th position. */
341
342void ssx_eval_rho(SSX *ssx)
343{     int m = ssx->m;
344      int p = ssx->p;
345      mpq_t *rho = ssx->rho;
346      int i;
347      xassert(1 <= p && p <= m);
348      /* rho := 0 */
349      for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1);
350      /* rho := e[p] */
351      mpq_set_si(rho[p], 1, 1);
352      /* rho := inv(B') * rho */
353      bfx_btran(ssx->binv, rho);
354      return;
355}
356
357/*----------------------------------------------------------------------
358// ssx_eval_row - compute pivot row of the simplex table.
359//
360// This routine computes p-th (pivot) row of the current simplex table
361// A~ = - inv(B) * N using the following formula:
362//
363//    A~[p] = - N' * inv(B') * e[p] = - N' * rho[p],
364//
365// where N' is a matrix transposed to the matrix N, rho[p] is p-th row
366// of the inverse inv(B). */
367
368void ssx_eval_row(SSX *ssx)
369{     int m = ssx->m;
370      int n = ssx->n;
371      int *A_ptr = ssx->A_ptr;
372      int *A_ind = ssx->A_ind;
373      mpq_t *A_val = ssx->A_val;
374      int *Q_col = ssx->Q_col;
375      mpq_t *rho = ssx->rho;
376      mpq_t *ap = ssx->ap;
377      int j, k, ptr;
378      mpq_t temp;
379      mpq_init(temp);
380      for (j = 1; j <= n; j++)
381      {  /* ap[j] := - N'[j] * rho (inner product) */
382         k = Q_col[m+j]; /* x[k] = xN[j] */
383         if (k <= m)
384            mpq_neg(ap[j], rho[k]);
385         else
386         {  mpq_set_si(ap[j], 0, 1);
387            for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
388            {  mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]);
389               mpq_add(ap[j], ap[j], temp);
390            }
391         }
392      }
393      mpq_clear(temp);
394      return;
395}
396
397/*----------------------------------------------------------------------
398// ssx_eval_col - compute pivot column of the simplex table.
399//
400// This routine computes q-th (pivot) column of the current simplex
401// table A~ = - inv(B) * N using the following formula:
402//
403//    A~[q] = - inv(B) * N[q],
404//
405// where N[q] is q-th column of the matrix N corresponding to chosen
406// non-basic variable xN[q]. */
407
408void ssx_eval_col(SSX *ssx)
409{     int m = ssx->m;
410      int n = ssx->n;
411      int *A_ptr = ssx->A_ptr;
412      int *A_ind = ssx->A_ind;
413      mpq_t *A_val = ssx->A_val;
414      int *Q_col = ssx->Q_col;
415      int q = ssx->q;
416      mpq_t *aq = ssx->aq;
417      int i, k, ptr;
418      xassert(1 <= q && q <= n);
419      /* aq := 0 */
420      for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1);
421      /* aq := N[q] */
422      k = Q_col[m+q]; /* x[k] = xN[q] */
423      if (k <= m)
424      {  /* N[q] is a column of the unity matrix I */
425         mpq_set_si(aq[k], 1, 1);
426      }
427      else
428      {  /* N[q] is a column of the original constraint matrix -A */
429         for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
430            mpq_neg(aq[A_ind[ptr]], A_val[ptr]);
431      }
432      /* aq := inv(B) * aq */
433      bfx_ftran(ssx->binv, aq, 1);
434      /* aq := - aq */
435      for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]);
436      return;
437}
438
439/*----------------------------------------------------------------------
440// ssx_chuzc - choose pivot column.
441//
442// This routine chooses non-basic variable xN[q] whose reduced cost
443// indicates possible improving of the objective function to enter it
444// in the basis.
445//
446// Currently the standard (textbook) pricing is used, i.e. that
447// non-basic variable is preferred which has greatest reduced cost (in
448// magnitude).
449//
450// If xN[q] has been chosen, the routine stores its number q and also
451// sets the flag q_dir that indicates direction in which xN[q] has to
452// change (+1 means increasing, -1 means decreasing).
453//
454// If the choice cannot be made, because the current basic solution is
455// dual feasible, the routine sets the number q to 0. */
456
457void ssx_chuzc(SSX *ssx)
458{     int m = ssx->m;
459      int n = ssx->n;
460      int dir = (ssx->dir == SSX_MIN ? +1 : -1);
461      int *Q_col = ssx->Q_col;
462      int *stat = ssx->stat;
463      mpq_t *cbar = ssx->cbar;
464      int j, k, s, q, q_dir;
465      double best, temp;
466      /* nothing is chosen so far */
467      q = 0, q_dir = 0, best = 0.0;
468      /* look through the list of non-basic variables */
469      for (j = 1; j <= n; j++)
470      {  k = Q_col[m+j]; /* x[k] = xN[j] */
471         s = dir * mpq_sgn(cbar[j]);
472         if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 ||
473             (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0)
474         {  /* reduced cost of xN[j] indicates possible improving of
475               the objective function */
476            temp = fabs(mpq_get_d(cbar[j]));
477            xassert(temp != 0.0);
478            if (q == 0 || best < temp)
479               q = j, q_dir = - s, best = temp;
480         }
481      }
482      ssx->q = q, ssx->q_dir = q_dir;
483      return;
484}
485
486/*----------------------------------------------------------------------
487// ssx_chuzr - choose pivot row.
488//
489// This routine looks through elements of q-th column of the simplex
490// table and chooses basic variable xB[p] which should leave the basis.
491//
492// The choice is based on the standard (textbook) ratio test.
493//
494// If xB[p] has been chosen, the routine stores its number p and also
495// sets its non-basic status p_stat which should be assigned to xB[p]
496// when it has left the basis and become xN[q].
497//
498// Special case p < 0 means that xN[q] is double-bounded variable and
499// it reaches its opposite bound before any basic variable does that,
500// so the current basis remains unchanged.
501//
502// If the choice cannot be made, because xN[q] can infinitely change in
503// the feasible direction, the routine sets the number p to 0. */
504
505void ssx_chuzr(SSX *ssx)
506{     int m = ssx->m;
507      int n = ssx->n;
508      int *type = ssx->type;
509      mpq_t *lb = ssx->lb;
510      mpq_t *ub = ssx->ub;
511      int *Q_col = ssx->Q_col;
512      mpq_t *bbar = ssx->bbar;
513      int q = ssx->q;
514      mpq_t *aq = ssx->aq;
515      int q_dir = ssx->q_dir;
516      int i, k, s, t, p, p_stat;
517      mpq_t teta, temp;
518      mpq_init(teta);
519      mpq_init(temp);
520      xassert(1 <= q && q <= n);
521      xassert(q_dir == +1 || q_dir == -1);
522      /* nothing is chosen so far */
523      p = 0, p_stat = 0;
524      /* look through the list of basic variables */
525      for (i = 1; i <= m; i++)
526      {  s = q_dir * mpq_sgn(aq[i]);
527         if (s < 0)
528         {  /* xB[i] decreases */
529            k = Q_col[i]; /* x[k] = xB[i] */
530            t = type[k];
531            if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
532            {  /* xB[i] has finite lower bound */
533               mpq_sub(temp, bbar[i], lb[k]);
534               mpq_div(temp, temp, aq[i]);
535               mpq_abs(temp, temp);
536               if (p == 0 || mpq_cmp(teta, temp) > 0)
537               {  p = i;
538                  p_stat = (t == SSX_FX ? SSX_NS : SSX_NL);
539                  mpq_set(teta, temp);
540               }
541            }
542         }
543         else if (s > 0)
544         {  /* xB[i] increases */
545            k = Q_col[i]; /* x[k] = xB[i] */
546            t = type[k];
547            if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
548            {  /* xB[i] has finite upper bound */
549               mpq_sub(temp, bbar[i], ub[k]);
550               mpq_div(temp, temp, aq[i]);
551               mpq_abs(temp, temp);
552               if (p == 0 || mpq_cmp(teta, temp) > 0)
553               {  p = i;
554                  p_stat = (t == SSX_FX ? SSX_NS : SSX_NU);
555                  mpq_set(teta, temp);
556               }
557            }
558         }
559         /* if something has been chosen and the ratio test indicates
560            exact degeneracy, the search can be finished */
561         if (p != 0 && mpq_sgn(teta) == 0) break;
562      }
563      /* if xN[q] is double-bounded, check if it can reach its opposite
564         bound before any basic variable */
565      k = Q_col[m+q]; /* x[k] = xN[q] */
566      if (type[k] == SSX_DB)
567      {  mpq_sub(temp, ub[k], lb[k]);
568         if (p == 0 || mpq_cmp(teta, temp) > 0)
569         {  p = -1;
570            p_stat = -1;
571            mpq_set(teta, temp);
572         }
573      }
574      ssx->p = p;
575      ssx->p_stat = p_stat;
576      /* if xB[p] has been chosen, determine its actual change in the
577         adjacent basis (it has the same sign as q_dir) */
578      if (p != 0)
579      {  xassert(mpq_sgn(teta) >= 0);
580         if (q_dir > 0)
581            mpq_set(ssx->delta, teta);
582         else
583            mpq_neg(ssx->delta, teta);
584      }
585      mpq_clear(teta);
586      mpq_clear(temp);
587      return;
588}
589
590/*----------------------------------------------------------------------
591// ssx_update_bbar - update values of basic variables.
592//
593// This routine recomputes the current values of basic variables for
594// the adjacent basis.
595//
596// The simplex table for the current basis is the following:
597//
598//    xB[i] = sum{j in 1..n} alfa[i,j] * xN[q],  i = 1,...,m
599//
600// therefore
601//
602//    delta xB[i] = alfa[i,q] * delta xN[q],  i = 1,...,m
603//
604// where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the
605// adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of
606// xB[i]. This gives formulae for recomputing values of xB[i]:
607//
608//    xB.new[p] = xN[q] + delta xN[q]
609//
610// (because xN[q] becomes xB[p] in the adjacent basis), and
611//
612//    xB.new[i] = xB[i] + alfa[i,q] * delta xN[q],  i != p
613//
614// for other basic variables. */
615
616void ssx_update_bbar(SSX *ssx)
617{     int m = ssx->m;
618      int n = ssx->n;
619      mpq_t *bbar = ssx->bbar;
620      mpq_t *cbar = ssx->cbar;
621      int p = ssx->p;
622      int q = ssx->q;
623      mpq_t *aq = ssx->aq;
624      int i;
625      mpq_t temp;
626      mpq_init(temp);
627      xassert(1 <= q && q <= n);
628      if (p < 0)
629      {  /* xN[q] is double-bounded and goes to its opposite bound */
630         /* nop */;
631      }
632      else
633      {  /* xN[q] becomes xB[p] in the adjacent basis */
634         /* xB.new[p] = xN[q] + delta xN[q] */
635         xassert(1 <= p && p <= m);
636         ssx_get_xNj(ssx, q, temp);
637         mpq_add(bbar[p], temp, ssx->delta);
638      }
639      /* update values of other basic variables depending on xN[q] */
640      for (i = 1; i <= m; i++)
641      {  if (i == p) continue;
642         /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */
643         if (mpq_sgn(aq[i]) == 0) continue;
644         mpq_mul(temp, aq[i], ssx->delta);
645         mpq_add(bbar[i], bbar[i], temp);
646      }
647#if 1
648      /* update value of the objective function */
649      /* z.new = z + d[q] * delta xN[q] */
650      mpq_mul(temp, cbar[q], ssx->delta);
651      mpq_add(bbar[0], bbar[0], temp);
652#endif
653      mpq_clear(temp);
654      return;
655}
656
657/*----------------------------------------------------------------------
658-- ssx_update_pi - update simplex multipliers.
659--
660-- This routine recomputes the vector of simplex multipliers for the
661-- adjacent basis. */
662
663void ssx_update_pi(SSX *ssx)
664{     int m = ssx->m;
665      int n = ssx->n;
666      mpq_t *pi = ssx->pi;
667      mpq_t *cbar = ssx->cbar;
668      int p = ssx->p;
669      int q = ssx->q;
670      mpq_t *aq = ssx->aq;
671      mpq_t *rho = ssx->rho;
672      int i;
673      mpq_t new_dq, temp;
674      mpq_init(new_dq);
675      mpq_init(temp);
676      xassert(1 <= p && p <= m);
677      xassert(1 <= q && q <= n);
678      /* compute d[q] in the adjacent basis */
679      mpq_div(new_dq, cbar[q], aq[p]);
680      /* update the vector of simplex multipliers */
681      for (i = 1; i <= m; i++)
682      {  if (mpq_sgn(rho[i]) == 0) continue;
683         mpq_mul(temp, new_dq, rho[i]);
684         mpq_sub(pi[i], pi[i], temp);
685      }
686      mpq_clear(new_dq);
687      mpq_clear(temp);
688      return;
689}
690
691/*----------------------------------------------------------------------
692// ssx_update_cbar - update reduced costs of non-basic variables.
693//
694// This routine recomputes the vector of reduced costs of non-basic
695// variables for the adjacent basis. */
696
697void ssx_update_cbar(SSX *ssx)
698{     int m = ssx->m;
699      int n = ssx->n;
700      mpq_t *cbar = ssx->cbar;
701      int p = ssx->p;
702      int q = ssx->q;
703      mpq_t *ap = ssx->ap;
704      int j;
705      mpq_t temp;
706      mpq_init(temp);
707      xassert(1 <= p && p <= m);
708      xassert(1 <= q && q <= n);
709      /* compute d[q] in the adjacent basis */
710      /* d.new[q] = d[q] / alfa[p,q] */
711      mpq_div(cbar[q], cbar[q], ap[q]);
712      /* update reduced costs of other non-basic variables */
713      for (j = 1; j <= n; j++)
714      {  if (j == q) continue;
715         /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */
716         if (mpq_sgn(ap[j]) == 0) continue;
717         mpq_mul(temp, ap[j], cbar[q]);
718         mpq_sub(cbar[j], cbar[j], temp);
719      }
720      mpq_clear(temp);
721      return;
722}
723
724/*----------------------------------------------------------------------
725// ssx_change_basis - change current basis to adjacent one.
726//
727// This routine changes the current basis to the adjacent one swapping
728// basic variable xB[p] and non-basic variable xN[q]. */
729
730void ssx_change_basis(SSX *ssx)
731{     int m = ssx->m;
732      int n = ssx->n;
733      int *type = ssx->type;
734      int *stat = ssx->stat;
735      int *Q_row = ssx->Q_row;
736      int *Q_col = ssx->Q_col;
737      int p = ssx->p;
738      int q = ssx->q;
739      int p_stat = ssx->p_stat;
740      int k, kp, kq;
741      if (p < 0)
742      {  /* special case: xN[q] goes to its opposite bound */
743         xassert(1 <= q && q <= n);
744         k = Q_col[m+q]; /* x[k] = xN[q] */
745         xassert(type[k] == SSX_DB);
746         switch (stat[k])
747         {  case SSX_NL:
748               stat[k] = SSX_NU;
749               break;
750            case SSX_NU:
751               stat[k] = SSX_NL;
752               break;
753            default:
754               xassert(stat != stat);
755         }
756      }
757      else
758      {  /* xB[p] leaves the basis, xN[q] enters the basis */
759         xassert(1 <= p && p <= m);
760         xassert(1 <= q && q <= n);
761         kp = Q_col[p];   /* x[kp] = xB[p] */
762         kq = Q_col[m+q]; /* x[kq] = xN[q] */
763         /* check non-basic status of xB[p] which becomes xN[q] */
764         switch (type[kp])
765         {  case SSX_FR:
766               xassert(p_stat == SSX_NF);
767               break;
768            case SSX_LO:
769               xassert(p_stat == SSX_NL);
770               break;
771            case SSX_UP:
772               xassert(p_stat == SSX_NU);
773               break;
774            case SSX_DB:
775               xassert(p_stat == SSX_NL || p_stat == SSX_NU);
776               break;
777            case SSX_FX:
778               xassert(p_stat == SSX_NS);
779               break;
780            default:
781               xassert(type != type);
782         }
783         /* swap xB[p] and xN[q] */
784         stat[kp] = (char)p_stat, stat[kq] = SSX_BS;
785         Q_row[kp] = m+q, Q_row[kq] = p;
786         Q_col[p] = kq, Q_col[m+q] = kp;
787         /* update factorization of the basis matrix */
788         if (bfx_update(ssx->binv, p))
789         {  if (ssx_factorize(ssx))
790               xassert(("Internal error: basis matrix is singular", 0));
791         }
792      }
793      return;
794}
795
796/*----------------------------------------------------------------------
797// ssx_delete - delete simplex solver workspace.
798//
799// This routine deletes the simplex solver workspace freeing all the
800// memory allocated to this object. */
801
802void ssx_delete(SSX *ssx)
803{     int m = ssx->m;
804      int n = ssx->n;
805      int nnz = ssx->A_ptr[n+1]-1;
806      int i, j, k;
807      xfree(ssx->type);
808      for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]);
809      xfree(ssx->lb);
810      for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]);
811      xfree(ssx->ub);
812      for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]);
813      xfree(ssx->coef);
814      xfree(ssx->A_ptr);
815      xfree(ssx->A_ind);
816      for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]);
817      xfree(ssx->A_val);
818      xfree(ssx->stat);
819      xfree(ssx->Q_row);
820      xfree(ssx->Q_col);
821      bfx_delete_binv(ssx->binv);
822      for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]);
823      xfree(ssx->bbar);
824      for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]);
825      xfree(ssx->pi);
826      for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]);
827      xfree(ssx->cbar);
828      for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]);
829      xfree(ssx->rho);
830      for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]);
831      xfree(ssx->ap);
832      for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]);
833      xfree(ssx->aq);
834      mpq_clear(ssx->delta);
835      xfree(ssx);
836      return;
837}
838
839/* eof */
Note: See TracBrowser for help on using the repository browser.