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