|
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 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 struct var |
|
28 { /* structural variable */ |
|
29 int j; |
|
30 /* ordinal number */ |
|
31 double q; |
|
32 /* penalty value */ |
|
33 }; |
|
34 |
|
35 static 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 |
|
43 static 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 |
|
58 static 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 |
|
261 void 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 */ |