1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpscl.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,476 @@
1.4 +/* glpscl.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +
1.30 +/***********************************************************************
1.31 +* min_row_aij - determine minimal |a[i,j]| in i-th row
1.32 +*
1.33 +* This routine returns minimal magnitude of (non-zero) constraint
1.34 +* coefficients in i-th row of the constraint matrix.
1.35 +*
1.36 +* If the parameter scaled is zero, the original constraint matrix A is
1.37 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.38 +*
1.39 +* If i-th row of the matrix is empty, the routine returns 1. */
1.40 +
1.41 +static double min_row_aij(glp_prob *lp, int i, int scaled)
1.42 +{ GLPAIJ *aij;
1.43 + double min_aij, temp;
1.44 + xassert(1 <= i && i <= lp->m);
1.45 + min_aij = 1.0;
1.46 + for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
1.47 + { temp = fabs(aij->val);
1.48 + if (scaled) temp *= (aij->row->rii * aij->col->sjj);
1.49 + if (aij->r_prev == NULL || min_aij > temp)
1.50 + min_aij = temp;
1.51 + }
1.52 + return min_aij;
1.53 +}
1.54 +
1.55 +/***********************************************************************
1.56 +* max_row_aij - determine maximal |a[i,j]| in i-th row
1.57 +*
1.58 +* This routine returns maximal magnitude of (non-zero) constraint
1.59 +* coefficients in i-th row of the constraint matrix.
1.60 +*
1.61 +* If the parameter scaled is zero, the original constraint matrix A is
1.62 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.63 +*
1.64 +* If i-th row of the matrix is empty, the routine returns 1. */
1.65 +
1.66 +static double max_row_aij(glp_prob *lp, int i, int scaled)
1.67 +{ GLPAIJ *aij;
1.68 + double max_aij, temp;
1.69 + xassert(1 <= i && i <= lp->m);
1.70 + max_aij = 1.0;
1.71 + for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next)
1.72 + { temp = fabs(aij->val);
1.73 + if (scaled) temp *= (aij->row->rii * aij->col->sjj);
1.74 + if (aij->r_prev == NULL || max_aij < temp)
1.75 + max_aij = temp;
1.76 + }
1.77 + return max_aij;
1.78 +}
1.79 +
1.80 +/***********************************************************************
1.81 +* min_col_aij - determine minimal |a[i,j]| in j-th column
1.82 +*
1.83 +* This routine returns minimal magnitude of (non-zero) constraint
1.84 +* coefficients in j-th column of the constraint matrix.
1.85 +*
1.86 +* If the parameter scaled is zero, the original constraint matrix A is
1.87 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.88 +*
1.89 +* If j-th column of the matrix is empty, the routine returns 1. */
1.90 +
1.91 +static double min_col_aij(glp_prob *lp, int j, int scaled)
1.92 +{ GLPAIJ *aij;
1.93 + double min_aij, temp;
1.94 + xassert(1 <= j && j <= lp->n);
1.95 + min_aij = 1.0;
1.96 + for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
1.97 + { temp = fabs(aij->val);
1.98 + if (scaled) temp *= (aij->row->rii * aij->col->sjj);
1.99 + if (aij->c_prev == NULL || min_aij > temp)
1.100 + min_aij = temp;
1.101 + }
1.102 + return min_aij;
1.103 +}
1.104 +
1.105 +/***********************************************************************
1.106 +* max_col_aij - determine maximal |a[i,j]| in j-th column
1.107 +*
1.108 +* This routine returns maximal magnitude of (non-zero) constraint
1.109 +* coefficients in j-th column of the constraint matrix.
1.110 +*
1.111 +* If the parameter scaled is zero, the original constraint matrix A is
1.112 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.113 +*
1.114 +* If j-th column of the matrix is empty, the routine returns 1. */
1.115 +
1.116 +static double max_col_aij(glp_prob *lp, int j, int scaled)
1.117 +{ GLPAIJ *aij;
1.118 + double max_aij, temp;
1.119 + xassert(1 <= j && j <= lp->n);
1.120 + max_aij = 1.0;
1.121 + for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
1.122 + { temp = fabs(aij->val);
1.123 + if (scaled) temp *= (aij->row->rii * aij->col->sjj);
1.124 + if (aij->c_prev == NULL || max_aij < temp)
1.125 + max_aij = temp;
1.126 + }
1.127 + return max_aij;
1.128 +}
1.129 +
1.130 +/***********************************************************************
1.131 +* min_mat_aij - determine minimal |a[i,j]| in constraint matrix
1.132 +*
1.133 +* This routine returns minimal magnitude of (non-zero) constraint
1.134 +* coefficients in the constraint matrix.
1.135 +*
1.136 +* If the parameter scaled is zero, the original constraint matrix A is
1.137 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.138 +*
1.139 +* If the matrix is empty, the routine returns 1. */
1.140 +
1.141 +static double min_mat_aij(glp_prob *lp, int scaled)
1.142 +{ int i;
1.143 + double min_aij, temp;
1.144 + min_aij = 1.0;
1.145 + for (i = 1; i <= lp->m; i++)
1.146 + { temp = min_row_aij(lp, i, scaled);
1.147 + if (i == 1 || min_aij > temp)
1.148 + min_aij = temp;
1.149 + }
1.150 + return min_aij;
1.151 +}
1.152 +
1.153 +/***********************************************************************
1.154 +* max_mat_aij - determine maximal |a[i,j]| in constraint matrix
1.155 +*
1.156 +* This routine returns maximal magnitude of (non-zero) constraint
1.157 +* coefficients in the constraint matrix.
1.158 +*
1.159 +* If the parameter scaled is zero, the original constraint matrix A is
1.160 +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed.
1.161 +*
1.162 +* If the matrix is empty, the routine returns 1. */
1.163 +
1.164 +static double max_mat_aij(glp_prob *lp, int scaled)
1.165 +{ int i;
1.166 + double max_aij, temp;
1.167 + max_aij = 1.0;
1.168 + for (i = 1; i <= lp->m; i++)
1.169 + { temp = max_row_aij(lp, i, scaled);
1.170 + if (i == 1 || max_aij < temp)
1.171 + max_aij = temp;
1.172 + }
1.173 + return max_aij;
1.174 +}
1.175 +
1.176 +/***********************************************************************
1.177 +* eq_scaling - perform equilibration scaling
1.178 +*
1.179 +* This routine performs equilibration scaling of rows and columns of
1.180 +* the constraint matrix.
1.181 +*
1.182 +* If the parameter flag is zero, the routine scales rows at first and
1.183 +* then columns. Otherwise, the routine scales columns and then rows.
1.184 +*
1.185 +* Rows are scaled as follows:
1.186 +*
1.187 +* n
1.188 +* a'[i,j] = a[i,j] / max |a[i,j]|, i = 1,...,m.
1.189 +* j=1
1.190 +*
1.191 +* This makes the infinity (maximum) norm of each row of the matrix
1.192 +* equal to 1.
1.193 +*
1.194 +* Columns are scaled as follows:
1.195 +*
1.196 +* n
1.197 +* a'[i,j] = a[i,j] / max |a[i,j]|, j = 1,...,n.
1.198 +* i=1
1.199 +*
1.200 +* This makes the infinity (maximum) norm of each column of the matrix
1.201 +* equal to 1. */
1.202 +
1.203 +static void eq_scaling(glp_prob *lp, int flag)
1.204 +{ int i, j, pass;
1.205 + double temp;
1.206 + xassert(flag == 0 || flag == 1);
1.207 + for (pass = 0; pass <= 1; pass++)
1.208 + { if (pass == flag)
1.209 + { /* scale rows */
1.210 + for (i = 1; i <= lp->m; i++)
1.211 + { temp = max_row_aij(lp, i, 1);
1.212 + glp_set_rii(lp, i, glp_get_rii(lp, i) / temp);
1.213 + }
1.214 + }
1.215 + else
1.216 + { /* scale columns */
1.217 + for (j = 1; j <= lp->n; j++)
1.218 + { temp = max_col_aij(lp, j, 1);
1.219 + glp_set_sjj(lp, j, glp_get_sjj(lp, j) / temp);
1.220 + }
1.221 + }
1.222 + }
1.223 + return;
1.224 +}
1.225 +
1.226 +/***********************************************************************
1.227 +* gm_scaling - perform geometric mean scaling
1.228 +*
1.229 +* This routine performs geometric mean scaling of rows and columns of
1.230 +* the constraint matrix.
1.231 +*
1.232 +* If the parameter flag is zero, the routine scales rows at first and
1.233 +* then columns. Otherwise, the routine scales columns and then rows.
1.234 +*
1.235 +* Rows are scaled as follows:
1.236 +*
1.237 +* a'[i,j] = a[i,j] / sqrt(alfa[i] * beta[i]), i = 1,...,m,
1.238 +*
1.239 +* where:
1.240 +* n n
1.241 +* alfa[i] = min |a[i,j]|, beta[i] = max |a[i,j]|.
1.242 +* j=1 j=1
1.243 +*
1.244 +* This allows decreasing the ratio beta[i] / alfa[i] for each row of
1.245 +* the matrix.
1.246 +*
1.247 +* Columns are scaled as follows:
1.248 +*
1.249 +* a'[i,j] = a[i,j] / sqrt(alfa[j] * beta[j]), j = 1,...,n,
1.250 +*
1.251 +* where:
1.252 +* m m
1.253 +* alfa[j] = min |a[i,j]|, beta[j] = max |a[i,j]|.
1.254 +* i=1 i=1
1.255 +*
1.256 +* This allows decreasing the ratio beta[j] / alfa[j] for each column
1.257 +* of the matrix. */
1.258 +
1.259 +static void gm_scaling(glp_prob *lp, int flag)
1.260 +{ int i, j, pass;
1.261 + double temp;
1.262 + xassert(flag == 0 || flag == 1);
1.263 + for (pass = 0; pass <= 1; pass++)
1.264 + { if (pass == flag)
1.265 + { /* scale rows */
1.266 + for (i = 1; i <= lp->m; i++)
1.267 + { temp = min_row_aij(lp, i, 1) * max_row_aij(lp, i, 1);
1.268 + glp_set_rii(lp, i, glp_get_rii(lp, i) / sqrt(temp));
1.269 + }
1.270 + }
1.271 + else
1.272 + { /* scale columns */
1.273 + for (j = 1; j <= lp->n; j++)
1.274 + { temp = min_col_aij(lp, j, 1) * max_col_aij(lp, j, 1);
1.275 + glp_set_sjj(lp, j, glp_get_sjj(lp, j) / sqrt(temp));
1.276 + }
1.277 + }
1.278 + }
1.279 + return;
1.280 +}
1.281 +
1.282 +/***********************************************************************
1.283 +* max_row_ratio - determine worst scaling "quality" for rows
1.284 +*
1.285 +* This routine returns the worst scaling "quality" for rows of the
1.286 +* currently scaled constraint matrix:
1.287 +*
1.288 +* m
1.289 +* ratio = max ratio[i],
1.290 +* i=1
1.291 +* where:
1.292 +* n n
1.293 +* ratio[i] = max |a[i,j]| / min |a[i,j]|, 1 <= i <= m,
1.294 +* j=1 j=1
1.295 +*
1.296 +* is the scaling "quality" of i-th row. */
1.297 +
1.298 +static double max_row_ratio(glp_prob *lp)
1.299 +{ int i;
1.300 + double ratio, temp;
1.301 + ratio = 1.0;
1.302 + for (i = 1; i <= lp->m; i++)
1.303 + { temp = max_row_aij(lp, i, 1) / min_row_aij(lp, i, 1);
1.304 + if (i == 1 || ratio < temp) ratio = temp;
1.305 + }
1.306 + return ratio;
1.307 +}
1.308 +
1.309 +/***********************************************************************
1.310 +* max_col_ratio - determine worst scaling "quality" for columns
1.311 +*
1.312 +* This routine returns the worst scaling "quality" for columns of the
1.313 +* currently scaled constraint matrix:
1.314 +*
1.315 +* n
1.316 +* ratio = max ratio[j],
1.317 +* j=1
1.318 +* where:
1.319 +* m m
1.320 +* ratio[j] = max |a[i,j]| / min |a[i,j]|, 1 <= j <= n,
1.321 +* i=1 i=1
1.322 +*
1.323 +* is the scaling "quality" of j-th column. */
1.324 +
1.325 +static double max_col_ratio(glp_prob *lp)
1.326 +{ int j;
1.327 + double ratio, temp;
1.328 + ratio = 1.0;
1.329 + for (j = 1; j <= lp->n; j++)
1.330 + { temp = max_col_aij(lp, j, 1) / min_col_aij(lp, j, 1);
1.331 + if (j == 1 || ratio < temp) ratio = temp;
1.332 + }
1.333 + return ratio;
1.334 +}
1.335 +
1.336 +/***********************************************************************
1.337 +* gm_iterate - perform iterative geometric mean scaling
1.338 +*
1.339 +* This routine performs iterative geometric mean scaling of rows and
1.340 +* columns of the constraint matrix.
1.341 +*
1.342 +* The parameter it_max specifies the maximal number of iterations.
1.343 +* Recommended value of it_max is 15.
1.344 +*
1.345 +* The parameter tau specifies a minimal improvement of the scaling
1.346 +* "quality" on each iteration, 0 < tau < 1. It means than the scaling
1.347 +* process continues while the following condition is satisfied:
1.348 +*
1.349 +* ratio[k] <= tau * ratio[k-1],
1.350 +*
1.351 +* where ratio = max |a[i,j]| / min |a[i,j]| is the scaling "quality"
1.352 +* to be minimized, k is the iteration number. Recommended value of tau
1.353 +* is 0.90. */
1.354 +
1.355 +static void gm_iterate(glp_prob *lp, int it_max, double tau)
1.356 +{ int k, flag;
1.357 + double ratio = 0.0, r_old;
1.358 + /* if the scaling "quality" for rows is better than for columns,
1.359 + the rows are scaled first; otherwise, the columns are scaled
1.360 + first */
1.361 + flag = (max_row_ratio(lp) > max_col_ratio(lp));
1.362 + for (k = 1; k <= it_max; k++)
1.363 + { /* save the scaling "quality" from previous iteration */
1.364 + r_old = ratio;
1.365 + /* determine the current scaling "quality" */
1.366 + ratio = max_mat_aij(lp, 1) / min_mat_aij(lp, 1);
1.367 +#if 0
1.368 + xprintf("k = %d; ratio = %g\n", k, ratio);
1.369 +#endif
1.370 + /* if improvement is not enough, terminate scaling */
1.371 + if (k > 1 && ratio > tau * r_old) break;
1.372 + /* otherwise, perform another iteration */
1.373 + gm_scaling(lp, flag);
1.374 + }
1.375 + return;
1.376 +}
1.377 +
1.378 +/***********************************************************************
1.379 +* NAME
1.380 +*
1.381 +* scale_prob - scale problem data
1.382 +*
1.383 +* SYNOPSIS
1.384 +*
1.385 +* #include "glpscl.h"
1.386 +* void scale_prob(glp_prob *lp, int flags);
1.387 +*
1.388 +* DESCRIPTION
1.389 +*
1.390 +* The routine scale_prob performs automatic scaling of problem data
1.391 +* for the specified problem object. */
1.392 +
1.393 +static void scale_prob(glp_prob *lp, int flags)
1.394 +{ static const char *fmt =
1.395 + "%s: min|aij| = %10.3e max|aij| = %10.3e ratio = %10.3e\n";
1.396 + double min_aij, max_aij, ratio;
1.397 + xprintf("Scaling...\n");
1.398 + /* cancel the current scaling effect */
1.399 + glp_unscale_prob(lp);
1.400 + /* report original scaling "quality" */
1.401 + min_aij = min_mat_aij(lp, 1);
1.402 + max_aij = max_mat_aij(lp, 1);
1.403 + ratio = max_aij / min_aij;
1.404 + xprintf(fmt, " A", min_aij, max_aij, ratio);
1.405 + /* check if the problem is well scaled */
1.406 + if (min_aij >= 0.10 && max_aij <= 10.0)
1.407 + { xprintf("Problem data seem to be well scaled\n");
1.408 + /* skip scaling, if required */
1.409 + if (flags & GLP_SF_SKIP) goto done;
1.410 + }
1.411 + /* perform iterative geometric mean scaling, if required */
1.412 + if (flags & GLP_SF_GM)
1.413 + { gm_iterate(lp, 15, 0.90);
1.414 + min_aij = min_mat_aij(lp, 1);
1.415 + max_aij = max_mat_aij(lp, 1);
1.416 + ratio = max_aij / min_aij;
1.417 + xprintf(fmt, "GM", min_aij, max_aij, ratio);
1.418 + }
1.419 + /* perform equilibration scaling, if required */
1.420 + if (flags & GLP_SF_EQ)
1.421 + { eq_scaling(lp, max_row_ratio(lp) > max_col_ratio(lp));
1.422 + min_aij = min_mat_aij(lp, 1);
1.423 + max_aij = max_mat_aij(lp, 1);
1.424 + ratio = max_aij / min_aij;
1.425 + xprintf(fmt, "EQ", min_aij, max_aij, ratio);
1.426 + }
1.427 + /* round scale factors to nearest power of two, if required */
1.428 + if (flags & GLP_SF_2N)
1.429 + { int i, j;
1.430 + for (i = 1; i <= lp->m; i++)
1.431 + glp_set_rii(lp, i, round2n(glp_get_rii(lp, i)));
1.432 + for (j = 1; j <= lp->n; j++)
1.433 + glp_set_sjj(lp, j, round2n(glp_get_sjj(lp, j)));
1.434 + min_aij = min_mat_aij(lp, 1);
1.435 + max_aij = max_mat_aij(lp, 1);
1.436 + ratio = max_aij / min_aij;
1.437 + xprintf(fmt, "2N", min_aij, max_aij, ratio);
1.438 + }
1.439 +done: return;
1.440 +}
1.441 +
1.442 +/***********************************************************************
1.443 +* NAME
1.444 +*
1.445 +* glp_scale_prob - scale problem data
1.446 +*
1.447 +* SYNOPSIS
1.448 +*
1.449 +* void glp_scale_prob(glp_prob *lp, int flags);
1.450 +*
1.451 +* DESCRIPTION
1.452 +*
1.453 +* The routine glp_scale_prob performs automatic scaling of problem
1.454 +* data for the specified problem object.
1.455 +*
1.456 +* The parameter flags specifies scaling options used by the routine.
1.457 +* Options can be combined with the bitwise OR operator and may be the
1.458 +* following:
1.459 +*
1.460 +* GLP_SF_GM perform geometric mean scaling;
1.461 +* GLP_SF_EQ perform equilibration scaling;
1.462 +* GLP_SF_2N round scale factors to nearest power of two;
1.463 +* GLP_SF_SKIP skip scaling, if the problem is well scaled.
1.464 +*
1.465 +* The parameter flags may be specified as GLP_SF_AUTO, in which case
1.466 +* the routine chooses scaling options automatically. */
1.467 +
1.468 +void glp_scale_prob(glp_prob *lp, int flags)
1.469 +{ if (flags & ~(GLP_SF_GM | GLP_SF_EQ | GLP_SF_2N | GLP_SF_SKIP |
1.470 + GLP_SF_AUTO))
1.471 + xerror("glp_scale_prob: flags = 0x%02X; invalid scaling option"
1.472 + "s\n", flags);
1.473 + if (flags & GLP_SF_AUTO)
1.474 + flags = (GLP_SF_GM | GLP_SF_EQ | GLP_SF_SKIP);
1.475 + scale_prob(lp, flags);
1.476 + return;
1.477 +}
1.478 +
1.479 +/* eof */