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