diff -r d59bea55db9b -r c445c931472f src/glpscl.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glpscl.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,476 @@ +/* glpscl.c */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics, +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. +* E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "glpapi.h" + +/*********************************************************************** +* min_row_aij - determine minimal |a[i,j]| in i-th row +* +* This routine returns minimal magnitude of (non-zero) constraint +* coefficients in i-th row of the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If i-th row of the matrix is empty, the routine returns 1. */ + +static double min_row_aij(glp_prob *lp, int i, int scaled) +{ GLPAIJ *aij; + double min_aij, temp; + xassert(1 <= i && i <= lp->m); + min_aij = 1.0; + for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next) + { temp = fabs(aij->val); + if (scaled) temp *= (aij->row->rii * aij->col->sjj); + if (aij->r_prev == NULL || min_aij > temp) + min_aij = temp; + } + return min_aij; +} + +/*********************************************************************** +* max_row_aij - determine maximal |a[i,j]| in i-th row +* +* This routine returns maximal magnitude of (non-zero) constraint +* coefficients in i-th row of the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If i-th row of the matrix is empty, the routine returns 1. */ + +static double max_row_aij(glp_prob *lp, int i, int scaled) +{ GLPAIJ *aij; + double max_aij, temp; + xassert(1 <= i && i <= lp->m); + max_aij = 1.0; + for (aij = lp->row[i]->ptr; aij != NULL; aij = aij->r_next) + { temp = fabs(aij->val); + if (scaled) temp *= (aij->row->rii * aij->col->sjj); + if (aij->r_prev == NULL || max_aij < temp) + max_aij = temp; + } + return max_aij; +} + +/*********************************************************************** +* min_col_aij - determine minimal |a[i,j]| in j-th column +* +* This routine returns minimal magnitude of (non-zero) constraint +* coefficients in j-th column of the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If j-th column of the matrix is empty, the routine returns 1. */ + +static double min_col_aij(glp_prob *lp, int j, int scaled) +{ GLPAIJ *aij; + double min_aij, temp; + xassert(1 <= j && j <= lp->n); + min_aij = 1.0; + for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next) + { temp = fabs(aij->val); + if (scaled) temp *= (aij->row->rii * aij->col->sjj); + if (aij->c_prev == NULL || min_aij > temp) + min_aij = temp; + } + return min_aij; +} + +/*********************************************************************** +* max_col_aij - determine maximal |a[i,j]| in j-th column +* +* This routine returns maximal magnitude of (non-zero) constraint +* coefficients in j-th column of the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If j-th column of the matrix is empty, the routine returns 1. */ + +static double max_col_aij(glp_prob *lp, int j, int scaled) +{ GLPAIJ *aij; + double max_aij, temp; + xassert(1 <= j && j <= lp->n); + max_aij = 1.0; + for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next) + { temp = fabs(aij->val); + if (scaled) temp *= (aij->row->rii * aij->col->sjj); + if (aij->c_prev == NULL || max_aij < temp) + max_aij = temp; + } + return max_aij; +} + +/*********************************************************************** +* min_mat_aij - determine minimal |a[i,j]| in constraint matrix +* +* This routine returns minimal magnitude of (non-zero) constraint +* coefficients in the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If the matrix is empty, the routine returns 1. */ + +static double min_mat_aij(glp_prob *lp, int scaled) +{ int i; + double min_aij, temp; + min_aij = 1.0; + for (i = 1; i <= lp->m; i++) + { temp = min_row_aij(lp, i, scaled); + if (i == 1 || min_aij > temp) + min_aij = temp; + } + return min_aij; +} + +/*********************************************************************** +* max_mat_aij - determine maximal |a[i,j]| in constraint matrix +* +* This routine returns maximal magnitude of (non-zero) constraint +* coefficients in the constraint matrix. +* +* If the parameter scaled is zero, the original constraint matrix A is +* assumed. Otherwise, the scaled constraint matrix R*A*S is assumed. +* +* If the matrix is empty, the routine returns 1. */ + +static double max_mat_aij(glp_prob *lp, int scaled) +{ int i; + double max_aij, temp; + max_aij = 1.0; + for (i = 1; i <= lp->m; i++) + { temp = max_row_aij(lp, i, scaled); + if (i == 1 || max_aij < temp) + max_aij = temp; + } + return max_aij; +} + +/*********************************************************************** +* eq_scaling - perform equilibration scaling +* +* This routine performs equilibration scaling of rows and columns of +* the constraint matrix. +* +* If the parameter flag is zero, the routine scales rows at first and +* then columns. Otherwise, the routine scales columns and then rows. +* +* Rows are scaled as follows: +* +* n +* a'[i,j] = a[i,j] / max |a[i,j]|, i = 1,...,m. +* j=1 +* +* This makes the infinity (maximum) norm of each row of the matrix +* equal to 1. +* +* Columns are scaled as follows: +* +* n +* a'[i,j] = a[i,j] / max |a[i,j]|, j = 1,...,n. +* i=1 +* +* This makes the infinity (maximum) norm of each column of the matrix +* equal to 1. */ + +static void eq_scaling(glp_prob *lp, int flag) +{ int i, j, pass; + double temp; + xassert(flag == 0 || flag == 1); + for (pass = 0; pass <= 1; pass++) + { if (pass == flag) + { /* scale rows */ + for (i = 1; i <= lp->m; i++) + { temp = max_row_aij(lp, i, 1); + glp_set_rii(lp, i, glp_get_rii(lp, i) / temp); + } + } + else + { /* scale columns */ + for (j = 1; j <= lp->n; j++) + { temp = max_col_aij(lp, j, 1); + glp_set_sjj(lp, j, glp_get_sjj(lp, j) / temp); + } + } + } + return; +} + +/*********************************************************************** +* gm_scaling - perform geometric mean scaling +* +* This routine performs geometric mean scaling of rows and columns of +* the constraint matrix. +* +* If the parameter flag is zero, the routine scales rows at first and +* then columns. Otherwise, the routine scales columns and then rows. +* +* Rows are scaled as follows: +* +* a'[i,j] = a[i,j] / sqrt(alfa[i] * beta[i]), i = 1,...,m, +* +* where: +* n n +* alfa[i] = min |a[i,j]|, beta[i] = max |a[i,j]|. +* j=1 j=1 +* +* This allows decreasing the ratio beta[i] / alfa[i] for each row of +* the matrix. +* +* Columns are scaled as follows: +* +* a'[i,j] = a[i,j] / sqrt(alfa[j] * beta[j]), j = 1,...,n, +* +* where: +* m m +* alfa[j] = min |a[i,j]|, beta[j] = max |a[i,j]|. +* i=1 i=1 +* +* This allows decreasing the ratio beta[j] / alfa[j] for each column +* of the matrix. */ + +static void gm_scaling(glp_prob *lp, int flag) +{ int i, j, pass; + double temp; + xassert(flag == 0 || flag == 1); + for (pass = 0; pass <= 1; pass++) + { if (pass == flag) + { /* scale rows */ + for (i = 1; i <= lp->m; i++) + { temp = min_row_aij(lp, i, 1) * max_row_aij(lp, i, 1); + glp_set_rii(lp, i, glp_get_rii(lp, i) / sqrt(temp)); + } + } + else + { /* scale columns */ + for (j = 1; j <= lp->n; j++) + { temp = min_col_aij(lp, j, 1) * max_col_aij(lp, j, 1); + glp_set_sjj(lp, j, glp_get_sjj(lp, j) / sqrt(temp)); + } + } + } + return; +} + +/*********************************************************************** +* max_row_ratio - determine worst scaling "quality" for rows +* +* This routine returns the worst scaling "quality" for rows of the +* currently scaled constraint matrix: +* +* m +* ratio = max ratio[i], +* i=1 +* where: +* n n +* ratio[i] = max |a[i,j]| / min |a[i,j]|, 1 <= i <= m, +* j=1 j=1 +* +* is the scaling "quality" of i-th row. */ + +static double max_row_ratio(glp_prob *lp) +{ int i; + double ratio, temp; + ratio = 1.0; + for (i = 1; i <= lp->m; i++) + { temp = max_row_aij(lp, i, 1) / min_row_aij(lp, i, 1); + if (i == 1 || ratio < temp) ratio = temp; + } + return ratio; +} + +/*********************************************************************** +* max_col_ratio - determine worst scaling "quality" for columns +* +* This routine returns the worst scaling "quality" for columns of the +* currently scaled constraint matrix: +* +* n +* ratio = max ratio[j], +* j=1 +* where: +* m m +* ratio[j] = max |a[i,j]| / min |a[i,j]|, 1 <= j <= n, +* i=1 i=1 +* +* is the scaling "quality" of j-th column. */ + +static double max_col_ratio(glp_prob *lp) +{ int j; + double ratio, temp; + ratio = 1.0; + for (j = 1; j <= lp->n; j++) + { temp = max_col_aij(lp, j, 1) / min_col_aij(lp, j, 1); + if (j == 1 || ratio < temp) ratio = temp; + } + return ratio; +} + +/*********************************************************************** +* gm_iterate - perform iterative geometric mean scaling +* +* This routine performs iterative geometric mean scaling of rows and +* columns of the constraint matrix. +* +* The parameter it_max specifies the maximal number of iterations. +* Recommended value of it_max is 15. +* +* The parameter tau specifies a minimal improvement of the scaling +* "quality" on each iteration, 0 < tau < 1. It means than the scaling +* process continues while the following condition is satisfied: +* +* ratio[k] <= tau * ratio[k-1], +* +* where ratio = max |a[i,j]| / min |a[i,j]| is the scaling "quality" +* to be minimized, k is the iteration number. Recommended value of tau +* is 0.90. */ + +static void gm_iterate(glp_prob *lp, int it_max, double tau) +{ int k, flag; + double ratio = 0.0, r_old; + /* if the scaling "quality" for rows is better than for columns, + the rows are scaled first; otherwise, the columns are scaled + first */ + flag = (max_row_ratio(lp) > max_col_ratio(lp)); + for (k = 1; k <= it_max; k++) + { /* save the scaling "quality" from previous iteration */ + r_old = ratio; + /* determine the current scaling "quality" */ + ratio = max_mat_aij(lp, 1) / min_mat_aij(lp, 1); +#if 0 + xprintf("k = %d; ratio = %g\n", k, ratio); +#endif + /* if improvement is not enough, terminate scaling */ + if (k > 1 && ratio > tau * r_old) break; + /* otherwise, perform another iteration */ + gm_scaling(lp, flag); + } + return; +} + +/*********************************************************************** +* NAME +* +* scale_prob - scale problem data +* +* SYNOPSIS +* +* #include "glpscl.h" +* void scale_prob(glp_prob *lp, int flags); +* +* DESCRIPTION +* +* The routine scale_prob performs automatic scaling of problem data +* for the specified problem object. */ + +static void scale_prob(glp_prob *lp, int flags) +{ static const char *fmt = + "%s: min|aij| = %10.3e max|aij| = %10.3e ratio = %10.3e\n"; + double min_aij, max_aij, ratio; + xprintf("Scaling...\n"); + /* cancel the current scaling effect */ + glp_unscale_prob(lp); + /* report original scaling "quality" */ + min_aij = min_mat_aij(lp, 1); + max_aij = max_mat_aij(lp, 1); + ratio = max_aij / min_aij; + xprintf(fmt, " A", min_aij, max_aij, ratio); + /* check if the problem is well scaled */ + if (min_aij >= 0.10 && max_aij <= 10.0) + { xprintf("Problem data seem to be well scaled\n"); + /* skip scaling, if required */ + if (flags & GLP_SF_SKIP) goto done; + } + /* perform iterative geometric mean scaling, if required */ + if (flags & GLP_SF_GM) + { gm_iterate(lp, 15, 0.90); + min_aij = min_mat_aij(lp, 1); + max_aij = max_mat_aij(lp, 1); + ratio = max_aij / min_aij; + xprintf(fmt, "GM", min_aij, max_aij, ratio); + } + /* perform equilibration scaling, if required */ + if (flags & GLP_SF_EQ) + { eq_scaling(lp, max_row_ratio(lp) > max_col_ratio(lp)); + min_aij = min_mat_aij(lp, 1); + max_aij = max_mat_aij(lp, 1); + ratio = max_aij / min_aij; + xprintf(fmt, "EQ", min_aij, max_aij, ratio); + } + /* round scale factors to nearest power of two, if required */ + if (flags & GLP_SF_2N) + { int i, j; + for (i = 1; i <= lp->m; i++) + glp_set_rii(lp, i, round2n(glp_get_rii(lp, i))); + for (j = 1; j <= lp->n; j++) + glp_set_sjj(lp, j, round2n(glp_get_sjj(lp, j))); + min_aij = min_mat_aij(lp, 1); + max_aij = max_mat_aij(lp, 1); + ratio = max_aij / min_aij; + xprintf(fmt, "2N", min_aij, max_aij, ratio); + } +done: return; +} + +/*********************************************************************** +* NAME +* +* glp_scale_prob - scale problem data +* +* SYNOPSIS +* +* void glp_scale_prob(glp_prob *lp, int flags); +* +* DESCRIPTION +* +* The routine glp_scale_prob performs automatic scaling of problem +* data for the specified problem object. +* +* The parameter flags specifies scaling options used by the routine. +* Options can be combined with the bitwise OR operator and may be the +* following: +* +* GLP_SF_GM perform geometric mean scaling; +* GLP_SF_EQ perform equilibration scaling; +* GLP_SF_2N round scale factors to nearest power of two; +* GLP_SF_SKIP skip scaling, if the problem is well scaled. +* +* The parameter flags may be specified as GLP_SF_AUTO, in which case +* the routine chooses scaling options automatically. */ + +void glp_scale_prob(glp_prob *lp, int flags) +{ if (flags & ~(GLP_SF_GM | GLP_SF_EQ | GLP_SF_2N | GLP_SF_SKIP | + GLP_SF_AUTO)) + xerror("glp_scale_prob: flags = 0x%02X; invalid scaling option" + "s\n", flags); + if (flags & GLP_SF_AUTO) + flags = (GLP_SF_GM | GLP_SF_EQ | GLP_SF_SKIP); + scale_prob(lp, flags); + return; +} + +/* eof */