lemon-project-template-glpk
comparison deps/glpk/src/glpscl.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:ebc11c9b475b |
---|---|
1 /* glpscl.c (problem scaling 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 * min_row_aij - determine minimal |a[i,j]| in i-th row | |
29 * | |
30 * This routine returns minimal magnitude of (non-zero) constraint | |
31 * coefficients in i-th row of the constraint matrix. | |
32 * | |
33 * If the parameter scaled is zero, the original constraint matrix A is | |
34 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
35 * | |
36 * If i-th row of the matrix is empty, the routine returns 1. */ | |
37 | |
38 static double min_row_aij(glp_prob *lp, int i, int scaled) | |
39 { GLPAIJ *aij; | |
40 double min_aij, temp; | |
41 xassert(1 <= i && i <= lp->m); | |
42 min_aij = 1.0; | |
43 for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next) | |
44 { temp = fabs(aij->val); | |
45 if (scaled) temp *= (aij->row->rii * aij->col->sjj); | |
46 if (aij->r_prev == NULL || min_aij > temp) | |
47 min_aij = temp; | |
48 } | |
49 return min_aij; | |
50 } | |
51 | |
52 /*********************************************************************** | |
53 * max_row_aij - determine maximal |a[i,j]| in i-th row | |
54 * | |
55 * This routine returns maximal magnitude of (non-zero) constraint | |
56 * coefficients in i-th row of the constraint matrix. | |
57 * | |
58 * If the parameter scaled is zero, the original constraint matrix A is | |
59 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
60 * | |
61 * If i-th row of the matrix is empty, the routine returns 1. */ | |
62 | |
63 static double max_row_aij(glp_prob *lp, int i, int scaled) | |
64 { GLPAIJ *aij; | |
65 double max_aij, temp; | |
66 xassert(1 <= i && i <= lp->m); | |
67 max_aij = 1.0; | |
68 for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next) | |
69 { temp = fabs(aij->val); | |
70 if (scaled) temp *= (aij->row->rii * aij->col->sjj); | |
71 if (aij->r_prev == NULL || max_aij < temp) | |
72 max_aij = temp; | |
73 } | |
74 return max_aij; | |
75 } | |
76 | |
77 /*********************************************************************** | |
78 * min_col_aij - determine minimal |a[i,j]| in j-th column | |
79 * | |
80 * This routine returns minimal magnitude of (non-zero) constraint | |
81 * coefficients in j-th column of the constraint matrix. | |
82 * | |
83 * If the parameter scaled is zero, the original constraint matrix A is | |
84 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
85 * | |
86 * If j-th column of the matrix is empty, the routine returns 1. */ | |
87 | |
88 static double min_col_aij(glp_prob *lp, int j, int scaled) | |
89 { GLPAIJ *aij; | |
90 double min_aij, temp; | |
91 xassert(1 <= j && j <= lp->n); | |
92 min_aij = 1.0; | |
93 for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next) | |
94 { temp = fabs(aij->val); | |
95 if (scaled) temp *= (aij->row->rii * aij->col->sjj); | |
96 if (aij->c_prev == NULL || min_aij > temp) | |
97 min_aij = temp; | |
98 } | |
99 return min_aij; | |
100 } | |
101 | |
102 /*********************************************************************** | |
103 * max_col_aij - determine maximal |a[i,j]| in j-th column | |
104 * | |
105 * This routine returns maximal magnitude of (non-zero) constraint | |
106 * coefficients in j-th column of the constraint matrix. | |
107 * | |
108 * If the parameter scaled is zero, the original constraint matrix A is | |
109 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
110 * | |
111 * If j-th column of the matrix is empty, the routine returns 1. */ | |
112 | |
113 static double max_col_aij(glp_prob *lp, int j, int scaled) | |
114 { GLPAIJ *aij; | |
115 double max_aij, temp; | |
116 xassert(1 <= j && j <= lp->n); | |
117 max_aij = 1.0; | |
118 for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next) | |
119 { temp = fabs(aij->val); | |
120 if (scaled) temp *= (aij->row->rii * aij->col->sjj); | |
121 if (aij->c_prev == NULL || max_aij < temp) | |
122 max_aij = temp; | |
123 } | |
124 return max_aij; | |
125 } | |
126 | |
127 /*********************************************************************** | |
128 * min_mat_aij - determine minimal |a[i,j]| in constraint matrix | |
129 * | |
130 * This routine returns minimal magnitude of (non-zero) constraint | |
131 * coefficients in the constraint matrix. | |
132 * | |
133 * If the parameter scaled is zero, the original constraint matrix A is | |
134 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
135 * | |
136 * If the matrix is empty, the routine returns 1. */ | |
137 | |
138 static double min_mat_aij(glp_prob *lp, int scaled) | |
139 { int i; | |
140 double min_aij, temp; | |
141 min_aij = 1.0; | |
142 for (i = 1; i <= lp->m; i++) | |
143 { temp = min_row_aij(lp, i, scaled); | |
144 if (i == 1 || min_aij > temp) | |
145 min_aij = temp; | |
146 } | |
147 return min_aij; | |
148 } | |
149 | |
150 /*********************************************************************** | |
151 * max_mat_aij - determine maximal |a[i,j]| in constraint matrix | |
152 * | |
153 * This routine returns maximal magnitude of (non-zero) constraint | |
154 * coefficients in the constraint matrix. | |
155 * | |
156 * If the parameter scaled is zero, the original constraint matrix A is | |
157 * assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. | |
158 * | |
159 * If the matrix is empty, the routine returns 1. */ | |
160 | |
161 static double max_mat_aij(glp_prob *lp, int scaled) | |
162 { int i; | |
163 double max_aij, temp; | |
164 max_aij = 1.0; | |
165 for (i = 1; i <= lp->m; i++) | |
166 { temp = max_row_aij(lp, i, scaled); | |
167 if (i == 1 || max_aij < temp) | |
168 max_aij = temp; | |
169 } | |
170 return max_aij; | |
171 } | |
172 | |
173 /*********************************************************************** | |
174 * eq_scaling - perform equilibration scaling | |
175 * | |
176 * This routine performs equilibration scaling of rows and columns of | |
177 * the constraint matrix. | |
178 * | |
179 * If the parameter flag is zero, the routine scales rows at first and | |
180 * then columns. Otherwise, the routine scales columns and then rows. | |
181 * | |
182 * Rows are scaled as follows: | |
183 * | |
184 * n | |
185 * a'[i,j] = a[i,j] / max |a[i,j]|, i = 1,...,m. | |
186 * j=1 | |
187 * | |
188 * This makes the infinity (maximum) norm of each row of the matrix | |
189 * equal to 1. | |
190 * | |
191 * Columns are scaled as follows: | |
192 * | |
193 * m | |
194 * a'[i,j] = a[i,j] / max |a[i,j]|, j = 1,...,n. | |
195 * i=1 | |
196 * | |
197 * This makes the infinity (maximum) norm of each column of the matrix | |
198 * equal to 1. */ | |
199 | |
200 static void eq_scaling(glp_prob *lp, int flag) | |
201 { int i, j, pass; | |
202 double temp; | |
203 xassert(flag == 0 || flag == 1); | |
204 for (pass = 0; pass <= 1; pass++) | |
205 { if (pass == flag) | |
206 { /* scale rows */ | |
207 for (i = 1; i <= lp->m; i++) | |
208 { temp = max_row_aij(lp, i, 1); | |
209 glp_set_rii(lp, i, glp_get_rii(lp, i) / temp); | |
210 } | |
211 } | |
212 else | |
213 { /* scale columns */ | |
214 for (j = 1; j <= lp->n; j++) | |
215 { temp = max_col_aij(lp, j, 1); | |
216 glp_set_sjj(lp, j, glp_get_sjj(lp, j) / temp); | |
217 } | |
218 } | |
219 } | |
220 return; | |
221 } | |
222 | |
223 /*********************************************************************** | |
224 * gm_scaling - perform geometric mean scaling | |
225 * | |
226 * This routine performs geometric mean scaling of rows and columns of | |
227 * the constraint matrix. | |
228 * | |
229 * If the parameter flag is zero, the routine scales rows at first and | |
230 * then columns. Otherwise, the routine scales columns and then rows. | |
231 * | |
232 * Rows are scaled as follows: | |
233 * | |
234 * a'[i,j] = a[i,j] / sqrt(alfa[i] * beta[i]), i = 1,...,m, | |
235 * | |
236 * where: | |
237 * n n | |
238 * alfa[i] = min |a[i,j]|, beta[i] = max |a[i,j]|. | |
239 * j=1 j=1 | |
240 * | |
241 * This allows decreasing the ratio beta[i] / alfa[i] for each row of | |
242 * the matrix. | |
243 * | |
244 * Columns are scaled as follows: | |
245 * | |
246 * a'[i,j] = a[i,j] / sqrt(alfa[j] * beta[j]), j = 1,...,n, | |
247 * | |
248 * where: | |
249 * m m | |
250 * alfa[j] = min |a[i,j]|, beta[j] = max |a[i,j]|. | |
251 * i=1 i=1 | |
252 * | |
253 * This allows decreasing the ratio beta[j] / alfa[j] for each column | |
254 * of the matrix. */ | |
255 | |
256 static void gm_scaling(glp_prob *lp, int flag) | |
257 { int i, j, pass; | |
258 double temp; | |
259 xassert(flag == 0 || flag == 1); | |
260 for (pass = 0; pass <= 1; pass++) | |
261 { if (pass == flag) | |
262 { /* scale rows */ | |
263 for (i = 1; i <= lp->m; i++) | |
264 { temp = min_row_aij(lp, i, 1) * max_row_aij(lp, i, 1); | |
265 glp_set_rii(lp, i, glp_get_rii(lp, i) / sqrt(temp)); | |
266 } | |
267 } | |
268 else | |
269 { /* scale columns */ | |
270 for (j = 1; j <= lp->n; j++) | |
271 { temp = min_col_aij(lp, j, 1) * max_col_aij(lp, j, 1); | |
272 glp_set_sjj(lp, j, glp_get_sjj(lp, j) / sqrt(temp)); | |
273 } | |
274 } | |
275 } | |
276 return; | |
277 } | |
278 | |
279 /*********************************************************************** | |
280 * max_row_ratio - determine worst scaling "quality" for rows | |
281 * | |
282 * This routine returns the worst scaling "quality" for rows of the | |
283 * currently scaled constraint matrix: | |
284 * | |
285 * m | |
286 * ratio = max ratio[i], | |
287 * i=1 | |
288 * where: | |
289 * n n | |
290 * ratio[i] = max |a[i,j]| / min |a[i,j]|, 1 <= i <= m, | |
291 * j=1 j=1 | |
292 * | |
293 * is the scaling "quality" of i-th row. */ | |
294 | |
295 static double max_row_ratio(glp_prob *lp) | |
296 { int i; | |
297 double ratio, temp; | |
298 ratio = 1.0; | |
299 for (i = 1; i <= lp->m; i++) | |
300 { temp = max_row_aij(lp, i, 1) / min_row_aij(lp, i, 1); | |
301 if (i == 1 || ratio < temp) ratio = temp; | |
302 } | |
303 return ratio; | |
304 } | |
305 | |
306 /*********************************************************************** | |
307 * max_col_ratio - determine worst scaling "quality" for columns | |
308 * | |
309 * This routine returns the worst scaling "quality" for columns of the | |
310 * currently scaled constraint matrix: | |
311 * | |
312 * n | |
313 * ratio = max ratio[j], | |
314 * j=1 | |
315 * where: | |
316 * m m | |
317 * ratio[j] = max |a[i,j]| / min |a[i,j]|, 1 <= j <= n, | |
318 * i=1 i=1 | |
319 * | |
320 * is the scaling "quality" of j-th column. */ | |
321 | |
322 static double max_col_ratio(glp_prob *lp) | |
323 { int j; | |
324 double ratio, temp; | |
325 ratio = 1.0; | |
326 for (j = 1; j <= lp->n; j++) | |
327 { temp = max_col_aij(lp, j, 1) / min_col_aij(lp, j, 1); | |
328 if (j == 1 || ratio < temp) ratio = temp; | |
329 } | |
330 return ratio; | |
331 } | |
332 | |
333 /*********************************************************************** | |
334 * gm_iterate - perform iterative geometric mean scaling | |
335 * | |
336 * This routine performs iterative geometric mean scaling of rows and | |
337 * columns of the constraint matrix. | |
338 * | |
339 * The parameter it_max specifies the maximal number of iterations. | |
340 * Recommended value of it_max is 15. | |
341 * | |
342 * The parameter tau specifies a minimal improvement of the scaling | |
343 * "quality" on each iteration, 0 < tau < 1. It means than the scaling | |
344 * process continues while the following condition is satisfied: | |
345 * | |
346 * ratio[k] <= tau * ratio[k-1], | |
347 * | |
348 * where ratio = max |a[i,j]| / min |a[i,j]| is the scaling "quality" | |
349 * to be minimized, k is the iteration number. Recommended value of tau | |
350 * is 0.90. */ | |
351 | |
352 static void gm_iterate(glp_prob *lp, int it_max, double tau) | |
353 { int k, flag; | |
354 double ratio = 0.0, r_old; | |
355 /* if the scaling "quality" for rows is better than for columns, | |
356 the rows are scaled first; otherwise, the columns are scaled | |
357 first */ | |
358 flag = (max_row_ratio(lp) > max_col_ratio(lp)); | |
359 for (k = 1; k <= it_max; k++) | |
360 { /* save the scaling "quality" from previous iteration */ | |
361 r_old = ratio; | |
362 /* determine the current scaling "quality" */ | |
363 ratio = max_mat_aij(lp, 1) / min_mat_aij(lp, 1); | |
364 #if 0 | |
365 xprintf("k = %d; ratio = %g\n", k, ratio); | |
366 #endif | |
367 /* if improvement is not enough, terminate scaling */ | |
368 if (k > 1 && ratio > tau * r_old) break; | |
369 /* otherwise, perform another iteration */ | |
370 gm_scaling(lp, flag); | |
371 } | |
372 return; | |
373 } | |
374 | |
375 /*********************************************************************** | |
376 * NAME | |
377 * | |
378 * scale_prob - scale problem data | |
379 * | |
380 * SYNOPSIS | |
381 * | |
382 * #include "glpscl.h" | |
383 * void scale_prob(glp_prob *lp, int flags); | |
384 * | |
385 * DESCRIPTION | |
386 * | |
387 * The routine scale_prob performs automatic scaling of problem data | |
388 * for the specified problem object. */ | |
389 | |
390 static void scale_prob(glp_prob *lp, int flags) | |
391 { static const char *fmt = | |
392 "%s: min|aij| = %10.3e max|aij| = %10.3e ratio = %10.3e\n"; | |
393 double min_aij, max_aij, ratio; | |
394 xprintf("Scaling...\n"); | |
395 /* cancel the current scaling effect */ | |
396 glp_unscale_prob(lp); | |
397 /* report original scaling "quality" */ | |
398 min_aij = min_mat_aij(lp, 1); | |
399 max_aij = max_mat_aij(lp, 1); | |
400 ratio = max_aij / min_aij; | |
401 xprintf(fmt, " A", min_aij, max_aij, ratio); | |
402 /* check if the problem is well scaled */ | |
403 if (min_aij >= 0.10 && max_aij <= 10.0) | |
404 { xprintf("Problem data seem to be well scaled\n"); | |
405 /* skip scaling, if required */ | |
406 if (flags & GLP_SF_SKIP) goto done; | |
407 } | |
408 /* perform iterative geometric mean scaling, if required */ | |
409 if (flags & GLP_SF_GM) | |
410 { gm_iterate(lp, 15, 0.90); | |
411 min_aij = min_mat_aij(lp, 1); | |
412 max_aij = max_mat_aij(lp, 1); | |
413 ratio = max_aij / min_aij; | |
414 xprintf(fmt, "GM", min_aij, max_aij, ratio); | |
415 } | |
416 /* perform equilibration scaling, if required */ | |
417 if (flags & GLP_SF_EQ) | |
418 { eq_scaling(lp, max_row_ratio(lp) > max_col_ratio(lp)); | |
419 min_aij = min_mat_aij(lp, 1); | |
420 max_aij = max_mat_aij(lp, 1); | |
421 ratio = max_aij / min_aij; | |
422 xprintf(fmt, "EQ", min_aij, max_aij, ratio); | |
423 } | |
424 /* round scale factors to nearest power of two, if required */ | |
425 if (flags & GLP_SF_2N) | |
426 { int i, j; | |
427 for (i = 1; i <= lp->m; i++) | |
428 glp_set_rii(lp, i, round2n(glp_get_rii(lp, i))); | |
429 for (j = 1; j <= lp->n; j++) | |
430 glp_set_sjj(lp, j, round2n(glp_get_sjj(lp, j))); | |
431 min_aij = min_mat_aij(lp, 1); | |
432 max_aij = max_mat_aij(lp, 1); | |
433 ratio = max_aij / min_aij; | |
434 xprintf(fmt, "2N", min_aij, max_aij, ratio); | |
435 } | |
436 done: return; | |
437 } | |
438 | |
439 /*********************************************************************** | |
440 * NAME | |
441 * | |
442 * glp_scale_prob - scale problem data | |
443 * | |
444 * SYNOPSIS | |
445 * | |
446 * void glp_scale_prob(glp_prob *lp, int flags); | |
447 * | |
448 * DESCRIPTION | |
449 * | |
450 * The routine glp_scale_prob performs automatic scaling of problem | |
451 * data for the specified problem object. | |
452 * | |
453 * The parameter flags specifies scaling options used by the routine. | |
454 * Options can be combined with the bitwise OR operator and may be the | |
455 * following: | |
456 * | |
457 * GLP_SF_GM perform geometric mean scaling; | |
458 * GLP_SF_EQ perform equilibration scaling; | |
459 * GLP_SF_2N round scale factors to nearest power of two; | |
460 * GLP_SF_SKIP skip scaling, if the problem is well scaled. | |
461 * | |
462 * The parameter flags may be specified as GLP_SF_AUTO, in which case | |
463 * the routine chooses scaling options automatically. */ | |
464 | |
465 void glp_scale_prob(glp_prob *lp, int flags) | |
466 { if (flags & ~(GLP_SF_GM | GLP_SF_EQ | GLP_SF_2N | GLP_SF_SKIP | | |
467 GLP_SF_AUTO)) | |
468 xerror("glp_scale_prob: flags = 0x%02X; invalid scaling option" | |
469 "s\n", flags); | |
470 if (flags & GLP_SF_AUTO) | |
471 flags = (GLP_SF_GM | GLP_SF_EQ | GLP_SF_SKIP); | |
472 scale_prob(lp, flags); | |
473 return; | |
474 } | |
475 | |
476 /* eof */ |