/* 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: <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 "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 */