|
1 /* glpscl.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 /*********************************************************************** |
|
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 * n |
|
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 */ |