lemon-project-template-glpk
view 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 |
line source
1 /* glpscl.c (problem scaling routines) */
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 ***********************************************************************/
25 #include "glpapi.h"
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
476 /* eof */