You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
478 lines
16 KiB
478 lines
16 KiB
/* glpscl.c (problem scaling routines) */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
|
* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* 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 <http://www.gnu.org/licenses/>.
|
|
***********************************************************************/
|
|
|
|
#include "env.h"
|
|
#include "misc.h"
|
|
#include "prob.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:
|
|
*
|
|
* m
|
|
* 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 */
|