COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpini02.c @ 9:33de93886c88

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

Import GLPK 4.47

File size: 8.8 KB
Line 
1/* glpini02.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, 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
27struct var
28{     /* structural variable */
29      int j;
30      /* ordinal number */
31      double q;
32      /* penalty value */
33};
34
35static int fcmp(const void *ptr1, const void *ptr2)
36{     /* this routine is passed to the qsort() function */
37      struct var *col1 = (void *)ptr1, *col2 = (void *)ptr2;
38      if (col1->q < col2->q) return -1;
39      if (col1->q > col2->q) return +1;
40      return 0;
41}
42
43static int get_column(glp_prob *lp, int j, int ind[], double val[])
44{     /* Bixby's algorithm assumes that the constraint matrix is scaled
45         such that the maximum absolute value in every non-zero row and
46         column is 1 */
47      int k, len;
48      double big;
49      len = glp_get_mat_col(lp, j, ind, val);
50      big = 0.0;
51      for (k = 1; k <= len; k++)
52         if (big < fabs(val[k])) big = fabs(val[k]);
53      if (big == 0.0) big = 1.0;
54      for (k = 1; k <= len; k++) val[k] /= big;
55      return len;
56}
57
58static void cpx_basis(glp_prob *lp)
59{     /* main routine */
60      struct var *C, *C2, *C3, *C4;
61      int m, n, i, j, jk, k, l, ll, t, n2, n3, n4, type, len, *I, *r,
62         *ind;
63      double alpha, gamma, cmax, temp, *v, *val;
64      xprintf("Constructing initial basis...\n");
65      /* determine the number of rows and columns */
66      m = glp_get_num_rows(lp);
67      n = glp_get_num_cols(lp);
68      /* allocate working arrays */
69      C = xcalloc(1+n, sizeof(struct var));
70      I = xcalloc(1+m, sizeof(int));
71      r = xcalloc(1+m, sizeof(int));
72      v = xcalloc(1+m, sizeof(double));
73      ind = xcalloc(1+m, sizeof(int));
74      val = xcalloc(1+m, sizeof(double));
75      /* make all auxiliary variables non-basic */
76      for (i = 1; i <= m; i++)
77      {  if (glp_get_row_type(lp, i) != GLP_DB)
78            glp_set_row_stat(lp, i, GLP_NS);
79         else if (fabs(glp_get_row_lb(lp, i)) <=
80                  fabs(glp_get_row_ub(lp, i)))
81            glp_set_row_stat(lp, i, GLP_NL);
82         else
83            glp_set_row_stat(lp, i, GLP_NU);
84      }
85      /* make all structural variables non-basic */
86      for (j = 1; j <= n; j++)
87      {  if (glp_get_col_type(lp, j) != GLP_DB)
88            glp_set_col_stat(lp, j, GLP_NS);
89         else if (fabs(glp_get_col_lb(lp, j)) <=
90                  fabs(glp_get_col_ub(lp, j)))
91            glp_set_col_stat(lp, j, GLP_NL);
92         else
93            glp_set_col_stat(lp, j, GLP_NU);
94      }
95      /* C2 is a set of free structural variables */
96      n2 = 0, C2 = C + 0;
97      for (j = 1; j <= n; j++)
98      {  type = glp_get_col_type(lp, j);
99         if (type == GLP_FR)
100         {  n2++;
101            C2[n2].j = j;
102            C2[n2].q = 0.0;
103         }
104      }
105      /* C3 is a set of structural variables having excatly one (lower
106         or upper) bound */
107      n3 = 0, C3 = C2 + n2;
108      for (j = 1; j <= n; j++)
109      {  type = glp_get_col_type(lp, j);
110         if (type == GLP_LO)
111         {  n3++;
112            C3[n3].j = j;
113            C3[n3].q = + glp_get_col_lb(lp, j);
114         }
115         else if (type == GLP_UP)
116         {  n3++;
117            C3[n3].j = j;
118            C3[n3].q = - glp_get_col_ub(lp, j);
119         }
120      }
121      /* C4 is a set of structural variables having both (lower and
122         upper) bounds */
123      n4 = 0, C4 = C3 + n3;
124      for (j = 1; j <= n; j++)
125      {  type = glp_get_col_type(lp, j);
126         if (type == GLP_DB)
127         {  n4++;
128            C4[n4].j = j;
129            C4[n4].q = glp_get_col_lb(lp, j) - glp_get_col_ub(lp, j);
130         }
131      }
132      /* compute gamma = max{|c[j]|: 1 <= j <= n} */
133      gamma = 0.0;
134      for (j = 1; j <= n; j++)
135      {  temp = fabs(glp_get_obj_coef(lp, j));
136         if (gamma < temp) gamma = temp;
137      }
138      /* compute cmax */
139      cmax = (gamma == 0.0 ? 1.0 : 1000.0 * gamma);
140      /* compute final penalty for all structural variables within sets
141         C2, C3, and C4 */
142      switch (glp_get_obj_dir(lp))
143      {  case GLP_MIN: temp = +1.0; break;
144         case GLP_MAX: temp = -1.0; break;
145         default: xassert(lp != lp);
146      }
147      for (k = 1; k <= n2+n3+n4; k++)
148      {  j = C[k].j;
149         C[k].q += (temp * glp_get_obj_coef(lp, j)) / cmax;
150      }
151      /* sort structural variables within C2, C3, and C4 in ascending
152         order of penalty value */
153      qsort(C2+1, n2, sizeof(struct var), fcmp);
154      for (k = 1; k < n2; k++) xassert(C2[k].q <= C2[k+1].q);
155      qsort(C3+1, n3, sizeof(struct var), fcmp);
156      for (k = 1; k < n3; k++) xassert(C3[k].q <= C3[k+1].q);
157      qsort(C4+1, n4, sizeof(struct var), fcmp);
158      for (k = 1; k < n4; k++) xassert(C4[k].q <= C4[k+1].q);
159      /*** STEP 1 ***/
160      for (i = 1; i <= m; i++)
161      {  type = glp_get_row_type(lp, i);
162         if (type != GLP_FX)
163         {  /* row i is either free or inequality constraint */
164            glp_set_row_stat(lp, i, GLP_BS);
165            I[i] = 1;
166            r[i] = 1;
167         }
168         else
169         {  /* row i is equality constraint */
170            I[i] = 0;
171            r[i] = 0;
172         }
173         v[i] = +DBL_MAX;
174      }
175      /*** STEP 2 ***/
176      for (k = 1; k <= n2+n3+n4; k++)
177      {  jk = C[k].j;
178         len = get_column(lp, jk, ind, val);
179         /* let alpha = max{|A[l,jk]|: r[l] = 0} and let l' be such
180            that alpha = |A[l',jk]| */
181         alpha = 0.0, ll = 0;
182         for (t = 1; t <= len; t++)
183         {  l = ind[t];
184            if (r[l] == 0 && alpha < fabs(val[t]))
185               alpha = fabs(val[t]), ll = l;
186         }
187         if (alpha >= 0.99)
188         {  /* B := B union {jk} */
189            glp_set_col_stat(lp, jk, GLP_BS);
190            I[ll] = 1;
191            v[ll] = alpha;
192            /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */
193            for (t = 1; t <= len; t++)
194            {  l = ind[t];
195               if (val[t] != 0.0) r[l]++;
196            }
197            /* continue to the next k */
198            continue;
199         }
200         /* if |A[l,jk]| > 0.01 * v[l] for some l, continue to the
201            next k */
202         for (t = 1; t <= len; t++)
203         {  l = ind[t];
204            if (fabs(val[t]) > 0.01 * v[l]) break;
205         }
206         if (t <= len) continue;
207         /* otherwise, let alpha = max{|A[l,jk]|: I[l] = 0} and let l'
208            be such that alpha = |A[l',jk]| */
209         alpha = 0.0, ll = 0;
210         for (t = 1; t <= len; t++)
211         {  l = ind[t];
212            if (I[l] == 0 && alpha < fabs(val[t]))
213               alpha = fabs(val[t]), ll = l;
214         }
215         /* if alpha = 0, continue to the next k */
216         if (alpha == 0.0) continue;
217         /* B := B union {jk} */
218         glp_set_col_stat(lp, jk, GLP_BS);
219         I[ll] = 1;
220         v[ll] = alpha;
221         /* r[l] := r[l] + 1 for all l such that |A[l,jk]| != 0 */
222         for (t = 1; t <= len; t++)
223         {  l = ind[t];
224            if (val[t] != 0.0) r[l]++;
225         }
226      }
227      /*** STEP 3 ***/
228      /* add an artificial variable (auxiliary variable for equality
229         constraint) to cover each remaining uncovered row */
230      for (i = 1; i <= m; i++)
231         if (I[i] == 0) glp_set_row_stat(lp, i, GLP_BS);
232      /* free working arrays */
233      xfree(C);
234      xfree(I);
235      xfree(r);
236      xfree(v);
237      xfree(ind);
238      xfree(val);
239      return;
240}
241
242/***********************************************************************
243*  NAME
244*
245*  glp_cpx_basis - construct Bixby's initial LP basis
246*
247*  SYNOPSIS
248*
249*  void glp_cpx_basis(glp_prob *lp);
250*
251*  DESCRIPTION
252*
253*  The routine glp_cpx_basis constructs an advanced initial basis for
254*  the specified problem object.
255*
256*  The routine is based on Bixby's algorithm described in the paper:
257*
258*  Robert E. Bixby. Implementing the Simplex Method: The Initial Basis.
259*  ORSA Journal on Computing, Vol. 4, No. 3, 1992, pp. 267-84. */
260
261void glp_cpx_basis(glp_prob *lp)
262{     if (lp->m == 0 || lp->n == 0)
263         glp_std_basis(lp);
264      else
265         cpx_basis(lp);
266      return;
267}
268
269/* eof */
Note: See TracBrowser for help on using the repository browser.