#include <string.h>
#include "commonlib.h"
#include "lp_lib.h"
#include "lp_report.h"
#include "lp_scale.h"

#ifdef FORTIFY
# include "lp_fortify.h"
#endif


/*
    Scaling routines for lp_solve v5.0+
   ----------------------------------------------------------------------------------
    Author:        Kjell Eikland
    Contact:       kjell.eikland@broadpark.no
    License terms: LGPL.

    Requires:      lp_lib.h, lp_scale.h

    Release notes:
    v5.0.0  1 January 2004      Significantly expanded and repackaged scaling
                                routines.
    v5.0.1  20 February 2004    Modified rounding behaviour in several areas.
    v5.1.0  20 July 2004        Reworked with flexible matrix storage model.
    v5.2.0  20 February 2005    Converted to matrix storage model without the OF.

   ----------------------------------------------------------------------------------
*/

/* First define scaling and unscaling primitives */

REAL scaled_value(lprec *lp, REAL value, int index)
{
  if(fabs(value) < lp->infinite) {
    if(lp->scaling_used) {
      if(index > lp->rows)
        value /= lp->scalars[index];
      else
        value *= lp->scalars[index];
    }
  }
  else
    value = my_sign(value)*lp->infinite;
  return(value);
}

REAL unscaled_value(lprec *lp, REAL value, int index)
{
  if(fabs(value) < lp->infinite) {
    if(lp->scaling_used) {
      if(index > lp->rows)
        value *= lp->scalars[index];
      else
        value /= lp->scalars[index];
    }
  }
  else
    value = my_sign(value)*lp->infinite;
  return(value);
}

STATIC REAL scaled_mat(lprec *lp, REAL value, int rownr, int colnr)
{
  if(lp->scaling_used)
    value *= lp->scalars[rownr] * lp->scalars[lp->rows + colnr];
  return( value );
}

STATIC REAL unscaled_mat(lprec *lp, REAL value, int rownr, int colnr)
{
  if(lp->scaling_used)
    value /= lp->scalars[rownr] * lp->scalars[lp->rows + colnr];
  return( value );
}

/* Compute the scale factor by the formulae:
      FALSE: SUM (log |Aij|) ^ 2
      TRUE:  SUM (log |Aij| - RowScale[i] - ColScale[j]) ^ 2 */
REAL CurtisReidMeasure(lprec *lp, MYBOOL _Advanced, REAL *FRowScale, REAL *FColScale)
{
  int      i, nz;
  REAL     absvalue, logvalue;
  register REAL result;
  MATrec   *mat = lp->matA;
  REAL     *value;
  int      *rownr, *colnr;

  /* Do OF part */
  result = 0;
  for(i = 1; i <= lp->columns; i++) {
    absvalue = fabs(lp->orig_obj[i]);
    if(absvalue > 0) {
      logvalue = log(absvalue);
      if(_Advanced)
        logvalue -= FRowScale[0] + FColScale[i];
      result += logvalue*logvalue;
    }
  }

  /* Do constraint matrix part */
  mat_validate(mat);
  value = &(COL_MAT_VALUE(0));
  rownr = &(COL_MAT_ROWNR(0));
  colnr = &(COL_MAT_COLNR(0));
  nz = get_nonzeros(lp);
  for(i = 0; i < nz;
      i++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
    absvalue = fabs(*value);
    if(absvalue > 0) {
      logvalue = log(absvalue);
      if(_Advanced)
        logvalue -= FRowScale[*rownr] + FColScale[*colnr];
      result += logvalue*logvalue;
    }
  }
  return( result );
}

/* Implementation of the Curtis-Reid scaling based on the paper
   "On the Automatic Scaling of Matrices for Gaussian
    Elimination," Journal of the Institute of Mathematics and
    Its Applications (1972) 10, 118-124.

    Solve the system | M   E | (r)   (sigma)
                     |       | ( ) = (     )
                     | E^T N | (c)   ( tau )

    by the conjugate gradient method (clever recurrences).

    E is the matrix A with all elements = 1

    M is diagonal matrix of row    counts (RowCount)
    N is diagonal matrix of column counts (ColCount)

    sigma is the vector of row    logarithm sums (RowSum)
    tau   is the vector of column logarithm sums (ColSum)

    r, c are resulting row and column scalings (RowScale, ColScale) */

int CurtisReidScales(lprec *lp, MYBOOL _Advanced, REAL *FRowScale, REAL *FColScale)
{
  int    i, row, col, ent, nz;
  REAL   *RowScalem2, *ColScalem2,
         *RowSum, *ColSum,
         *residual_even, *residual_odd;
  REAL   sk,   qk,     ek,
         skm1, qkm1,   ekm1,
         qkm2, qkqkm1, ekm2, ekekm1,
         absvalue, logvalue,
         StopTolerance;
  int    *RowCount, *ColCount, colMax;
  int    Result;
  MATrec *mat = lp->matA;
  REAL   *value;
  int    *rownr, *colnr;

  if(CurtisReidMeasure(lp, _Advanced, FRowScale, FColScale)<0.1*get_nonzeros(lp))
  return(0);

  /* Allocate temporary memory and find RowSum and ColSum measures */
  nz = get_nonzeros(lp);
  colMax = lp->columns;

  allocREAL(lp, &RowSum, lp->rows+1, TRUE);
  allocINT(lp,  &RowCount, lp->rows+1, TRUE);
  allocREAL(lp, &residual_odd, lp->rows+1, TRUE);

  allocREAL(lp, &ColSum, colMax+1, TRUE);
  allocINT(lp,  &ColCount, colMax+1, TRUE);
  allocREAL(lp, &residual_even, colMax+1, TRUE);

  allocREAL(lp, &RowScalem2, lp->rows+1, FALSE);
  allocREAL(lp, &ColScalem2, colMax+1, FALSE);

  /* Set origin for row scaling */
  for(i = 1; i <= colMax; i++) {
    absvalue=fabs(lp->orig_obj[i]);
    if(absvalue>0) {
      logvalue = log(absvalue);
      ColSum[i] += logvalue;
      RowSum[0] += logvalue;
      ColCount[i]++;
      RowCount[0]++;
    }
  }

  value = &(COL_MAT_VALUE(0));
  rownr = &(COL_MAT_ROWNR(0));
  colnr = &(COL_MAT_COLNR(0));
  for(i = 0; i < nz;
      i++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
    absvalue=fabs(*value);
    if(absvalue>0) {
      logvalue = log(absvalue);
      ColSum[*colnr] += logvalue;
      RowSum[*rownr] += logvalue;
      ColCount[*colnr]++;
      RowCount[*rownr]++;
    }
  }

  /* Make sure we dont't have division by zero errors */
  for(row = 0; row <= lp->rows; row++)
    if(RowCount[row] == 0)
      RowCount[row] = 1;
  for(col = 1; col <= colMax; col++)
    if(ColCount[col] == 0)
      ColCount[col] = 1;

  /* Initialize to RowScale = RowCount-1 RowSum
                   ColScale = 0.0
                   residual = ColSum - ET RowCount-1 RowSum */

  StopTolerance= MAX(lp->scalelimit-floor(lp->scalelimit), DEF_SCALINGEPS);
  StopTolerance *= (REAL) nz;
  for(row = 0; row <= lp->rows; row++) {
    FRowScale[row] = RowSum[row] / (REAL) RowCount[row];
    RowScalem2[row] = FRowScale[row];
  }

  /* Compute initial residual */
  for(col = 1; col <= colMax; col++) {
    FColScale[col] = 0;
    ColScalem2[col] = 0;
    residual_even[col] = ColSum[col];

    if(lp->orig_obj[col] != 0)
      residual_even[col] -= RowSum[0] / (REAL) RowCount[0];

    i = mat->col_end[col-1];
    rownr = &(COL_MAT_ROWNR(i));
    ent = mat->col_end[col];
    for(; i < ent;
        i++, rownr += matRowColStep) {
      residual_even[col] -= RowSum[*rownr] / (REAL) RowCount[*rownr];
    }
  }

  /* Compute sk */
  sk = 0;
  skm1 = 0;
  for(col = 1; col <= colMax; col++)
    sk += (residual_even[col]*residual_even[col]) / (REAL) ColCount[col];

  Result = 0;
  qk=1; qkm1=0; qkm2=0;
  ek=0; ekm1=0; ekm2=0;

  while(sk>StopTolerance) {
  /* Given the values of residual and sk, construct
     ColScale (when pass is even)
     RowScale (when pass is odd)  */

    qkqkm1 = qk * qkm1;
    ekekm1 = ek * ekm1;
    if((Result % 2) == 0) { /* pass is even; construct RowScale[pass+1] */
      if(Result != 0) {
        for(row = 0; row <= lp->rows; row++)
          RowScalem2[row] = FRowScale[row];
        if(qkqkm1 != 0) {
          for(row = 0; row <= lp->rows; row++)
            FRowScale[row]*=(1 + ekekm1 / qkqkm1);
          for(row = 0; row<=lp->rows; row++)
            FRowScale[row]+=(residual_odd[row] / (qkqkm1 * (REAL) RowCount[row]) -
                             RowScalem2[row] * ekekm1 / qkqkm1);
        }
      }
    }
    else { /* pass is odd; construct ColScale[pass+1] */
      for(col = 1; col <= colMax; col++)
        ColScalem2[col] = FColScale[col];
      if(qkqkm1 != 0) {
        for(col = 1; col <= colMax; col++)
          FColScale[col] *= (1 + ekekm1 / qkqkm1);
        for(col = 1; col <= colMax; col++)
          FColScale[col] += (residual_even[col] / ((REAL) ColCount[col] * qkqkm1) -
                             ColScalem2[col] * ekekm1 / qkqkm1);
      }
    }

    /* update residual and sk (pass + 1) */
    if((Result % 2) == 0) { /* even */
       /* residual */
      for(row = 0; row <= lp->rows; row++)
        residual_odd[row] *= ek;

      for(i = 1; i <= colMax; i++)
        if(lp->orig_obj[i] != 0)
          residual_odd[0] += (residual_even[i] / (REAL) ColCount[i]);

      rownr = &(COL_MAT_ROWNR(0));
      colnr = &(COL_MAT_COLNR(0));
      for(i = 0; i < nz;
          i++, rownr += matRowColStep, colnr += matRowColStep) {
        residual_odd[*rownr] += (residual_even[*colnr] / (REAL) ColCount[*colnr]);
      }
      for(row = 0; row <= lp->rows; row++)
        residual_odd[row] *= (-1 / qk);

      /* sk */
      skm1 = sk;
      sk = 0;
      for(row = 0; row <= lp->rows; row++)
        sk += (residual_odd[row]*residual_odd[row]) / (REAL) RowCount[row];
    }
    else { /* odd */
      /* residual */
      for(col = 1; col <= colMax; col++)
        residual_even[col] *= ek;

      for(i = 1; i <= colMax; i++)
        if(lp->orig_obj[i] != 0)
          residual_even[i] += (residual_odd[0] / (REAL) RowCount[0]);

      rownr = &(COL_MAT_ROWNR(0));
      colnr = &(COL_MAT_COLNR(0));
      for(i = 0; i < nz;
          i++, rownr += matRowColStep, colnr += matRowColStep) {
        residual_even[*colnr] += (residual_odd[*rownr] / (REAL) RowCount[*rownr]);
      }
      for(col = 1; col <= colMax; col++)
        residual_even[col] *= (-1 / qk);

      /* sk */
      skm1 = sk;
      sk = 0;
      for(col = 1; col <= colMax; col++)
        sk += (residual_even[col]*residual_even[col]) / (REAL) ColCount[col];
    }

    /* Compute ek and qk */
    ekm2=ekm1;
    ekm1=ek;
    ek=qk * sk / skm1;

    qkm2=qkm1;
    qkm1=qk;
    qk=1-ek;

    Result++;
  }

  /* Synchronize the RowScale and ColScale vectors */
  ekekm1 = ek * ekm1;
  if(qkm1 != 0) {
  if((Result % 2) == 0) { /* pass is even, compute RowScale */
    for(row = 0; row<=lp->rows; row++)
      FRowScale[row]*=(1.0 + ekekm1 / qkm1);
    for(row = 0; row<=lp->rows; row++)
      FRowScale[row]+=(residual_odd[row] / (qkm1 * (REAL) RowCount[row]) -
                      RowScalem2[row] * ekekm1 / qkm1);
  }
  else { /* pass is odd, compute ColScale */
    for(col=1; col<=colMax; col++)
      FColScale[col]*=(1 + ekekm1 / qkm1);
    for(col=1; col<=colMax; col++)
      FColScale[col]+=(residual_even[col] / ((REAL) ColCount[col] * qkm1) -
                       ColScalem2[col] * ekekm1 / qkm1);
  }
  }

  /* Do validation, if indicated */
  if(FALSE && mat_validate(mat)){
    double check, error;

    /* CHECK: M RowScale + E ColScale = RowSum */
    error = 0;
    for(row = 0; row <= lp->rows; row++) {
      check = (REAL) RowCount[row] * FRowScale[row];
      if(row == 0) {
        for(i = 1; i <= colMax; i++) {
          if(lp->orig_obj[i] != 0)
            check += FColScale[i];
        }
      }
      else {
        i = mat->row_end[row-1];
        ent = mat->row_end[row];
        for(; i < ent; i++) {
          col = ROW_MAT_COLNR(i);
          check += FColScale[col];
        }
      }
      check -= RowSum[row];
      error += check*check;
    }

    /* CHECK: E^T RowScale + N ColScale = ColSum */
    error = 0;
    for(col = 1; col <= colMax; col++) {
      check = (REAL) ColCount[col] * FColScale[col];

      if(lp->orig_obj[col] != 0)
        check += FRowScale[0];

      i = mat->col_end[col-1];
      ent = mat->col_end[col];
      rownr = &(COL_MAT_ROWNR(i));
      for(; i < ent;
          i++, rownr += matRowColStep) {
        check += FRowScale[*rownr];
      }
      check -= ColSum[col];
      error += check*check;
    }
  }

  /* Convert to scaling factors (rounding to nearest power
     of 2 can optionally be done as a separate step later) */
  for(col = 1; col <= colMax; col++) {
    absvalue = exp(-FColScale[col]);
    if(absvalue < MIN_SCALAR) absvalue = MIN_SCALAR;
    if(absvalue > MAX_SCALAR) absvalue = MAX_SCALAR;
    if(!is_int(lp,col) || is_integerscaling(lp))
        FColScale[col] = absvalue;
    else
        FColScale[col] = 1;
  }
  for(row = 0; row <= lp->rows; row++) {
    absvalue = exp(-FRowScale[row]);
    if(absvalue < MIN_SCALAR) absvalue = MIN_SCALAR;
    if(absvalue > MAX_SCALAR) absvalue = MAX_SCALAR;
    FRowScale[row] = absvalue;
  }

 /* free temporary memory */
  FREE(RowSum);
  FREE(ColSum);
  FREE(RowCount);
  FREE(ColCount);
  FREE(residual_even);
  FREE(residual_odd);
  FREE(RowScalem2);
  FREE(ColScalem2);

  return(Result);

}

STATIC MYBOOL scaleCR(lprec *lp, REAL *scaledelta)
{
  REAL *scalechange = NULL;
  int  Result;

  if(!lp->scaling_used) {
    allocREAL(lp, &lp->scalars, lp->sum_alloc + 1, FALSE);
    for(Result = 0; Result <= lp->sum; Result++)
      lp->scalars[Result] = 1;
    lp->scaling_used = TRUE;
  }

  if(scaledelta == NULL)
    allocREAL(lp, &scalechange, lp->sum + 1, FALSE);
  else
    scalechange = scaledelta;

  Result=CurtisReidScales(lp, FALSE, scalechange, &scalechange[lp->rows]);
  if(Result>0) {

    /* Do the scaling*/
    if(scale_updaterows(lp, scalechange, TRUE) ||
       scale_updatecolumns(lp, &scalechange[lp->rows], TRUE))
      lp->scalemode |= SCALE_CURTISREID;

    set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
  }

  if(scaledelta == NULL)
    FREE(scalechange);

  return((MYBOOL) (Result > 0));
}

STATIC MYBOOL transform_for_scale(lprec *lp, REAL *value)
{
  MYBOOL Accept = TRUE;
  *value = fabs(*value);
#ifdef Paranoia
  if(*value < lp->epsmachine) {
    Accept = FALSE;
    report(lp, SEVERE, "transform_for_scale: A zero-valued entry was passed\n");
  }
  else
#endif
  if(is_scalemode(lp, SCALE_LOGARITHMIC))
    *value = log(*value);
  else if(is_scalemode(lp, SCALE_QUADRATIC))
    (*value) *= (*value);
  return( Accept );
}

STATIC void accumulate_for_scale(lprec *lp, REAL *min, REAL *max, REAL value)
{
  if(transform_for_scale(lp, &value)) {
    if(is_scaletype(lp, SCALE_MEAN)) {
      *max += value;
      *min += 1;
    }
    else {
      SETMAX(*max, value);
      SETMIN(*min, value);
    }
  }
}

STATIC REAL minmax_to_scale(lprec *lp, REAL min, REAL max, int itemcount)
{
  REAL scale;

  /* Initialize according to transformation / weighting model */
  if(is_scalemode(lp, SCALE_LOGARITHMIC))
    scale = 0;
  else
    scale = 1;
  if(itemcount <= 0)
    return(scale);

  /* Compute base scalar according to chosen scaling type */
  if(is_scaletype(lp, SCALE_MEAN)) {
    if(min > 0)
      scale = max / min;
  }
  else if(is_scaletype(lp, SCALE_RANGE))
    scale = (max + min) / 2;
  else if(is_scaletype(lp, SCALE_GEOMETRIC))
    scale = sqrt(min*max);
  else if(is_scaletype(lp, SCALE_EXTREME))
    scale = max;

  /* Compute final scalar according to transformation / weighting model */
  if(is_scalemode(lp, SCALE_LOGARITHMIC))
    scale = exp(-scale);
  else if(is_scalemode(lp, SCALE_QUADRATIC)) {
    if(scale == 0)
      scale = 1;
    else
      scale = 1 / sqrt(scale);
  }
  else {
    if(scale == 0)
      scale = 1;
    else
      scale = 1 / scale;
  }

  /* Make sure we are within acceptable scaling ranges */
  SETMAX(scale, MIN_SCALAR);
  SETMIN(scale, MAX_SCALAR);

  return(scale);
}

STATIC REAL roundPower2(REAL scale)
/* Purpose is to round a number to it nearest power of 2; in a system
   with binary number representation, this avoids rounding errors when
   scale is used to normalize another value */
{
  long int power2;
  MYBOOL   isSmall = FALSE;

  if(scale == 1)
    return( scale );

  /* Obtain the fractional power of 2 */
  if(scale < 2) {
    scale = 2 / scale;
    isSmall = TRUE;
  }
  else
    scale /= 2;
  scale = log(scale)/log(2.0);

  /* Find the desired nearest power of two and compute the associated scalar */
  power2 = (long) ceil(scale-0.5);
  scale = 1 << power2;
  if(isSmall)
    scale = 1.0 / scale;

  return( scale );

}

STATIC MYBOOL scale_updatecolumns(lprec *lp, REAL *scalechange, MYBOOL updateonly)
{
  int i, j;

  /* Verify that the scale change is significant (different from the unit) */
  for(i = lp->columns; i > 0; i--)
    if(fabs(scalechange[i]-1) > lp->epsprimal)
      break;
  if(i <= 0)
    return( FALSE );

 /* Update the pre-existing column scalar */
  if(updateonly)
    for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++)
      lp->scalars[j] *= scalechange[i];
  else
    for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++)
      lp->scalars[j] = scalechange[i];

  return( TRUE );
}

STATIC MYBOOL scale_updaterows(lprec *lp, REAL *scalechange, MYBOOL updateonly)
{
  int i;

  /* Verify that the scale change is significant (different from the unit) */
  for(i = lp->rows; i >= 0; i--) {
    if(fabs(scalechange[i]-1) > lp->epsprimal)
      break;
  }
  if(i < 0)
    return( FALSE );

 /* Update the pre-existing row scalar */
  if(updateonly)
    for(i = 0; i <= lp->rows; i++)
      lp->scalars[i] *= scalechange[i];
  else
    for(i = 0; i <= lp->rows; i++)
      lp->scalars[i] = scalechange[i];

  return( TRUE );
}

STATIC MYBOOL scale_columns(lprec *lp, REAL *scaledelta)
{
  int     i,j, colMax, nz;
  REAL    *scalechange;
  REAL    *value;
  int     *colnr;
  MATrec  *mat = lp->matA;

  /* Check that columns are in fact targeted */
  if((lp->scalemode & SCALE_ROWSONLY) != 0)
    return( TRUE );

  if(scaledelta == NULL)
    scalechange = &lp->scalars[lp->rows];
  else
    scalechange = &scaledelta[lp->rows];

  colMax = lp->columns;

  /* Scale matrix entries (including any Lagrangean constraints) */
  for(i = 1; i <= lp->columns; i++) {
    lp->orig_obj[i] *= scalechange[i];
  }

  mat_validate(lp->matA);
  nz = get_nonzeros(lp);
  value = &(COL_MAT_VALUE(0));
  colnr = &(COL_MAT_COLNR(0));
  for(i = 0; i < nz;
      i++, value += matValueStep, colnr += matRowColStep) {
    (*value) *= scalechange[*colnr];
  }

  /* Scale variable bounds as well */
  for(i = 1, j = lp->rows + 1; j <= lp->sum; i++, j++) {
    if(lp->orig_lowbo[j] > -lp->infinite)
      lp->orig_lowbo[j] /= scalechange[i];
    if(lp->orig_upbo[j] < lp->infinite)
      lp->orig_upbo[j] /= scalechange[i];
    if(lp->sc_lobound[i] != 0)
      lp->sc_lobound[i] /= scalechange[i];
  }

  lp->columns_scaled = TRUE;
  set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);

  return( TRUE );
}

STATIC MYBOOL scale_rows(lprec *lp, REAL *scaledelta)
{
  int     i, j, nz, colMax;
  REAL    *scalechange;
  REAL    *value;
  int     *rownr;
  MATrec  *mat = lp->matA;


  /* Check that rows are in fact targeted */
  if((lp->scalemode & SCALE_COLSONLY) != 0)
    return( TRUE );

  if(scaledelta == NULL)
    scalechange = lp->scalars;
  else
    scalechange = scaledelta;

  colMax = lp->columns;

  /* First row-scale the matrix (including the objective function) */
  for(i = 1; i <= colMax; i++) {
    lp->orig_obj[i] *= scalechange[0];
  }

  nz = get_nonzeros(lp);
  value = &(COL_MAT_VALUE(0));
  rownr = &(COL_MAT_ROWNR(0));
  for(i = 0; i < nz;
      i++, value += matValueStep, rownr += matRowColStep) {
    (*value) *= scalechange[*rownr];
  }

  /* ...and scale the rhs and the row bounds (RANGES in MPS!!) */
  for(i = 0; i <= lp->rows; i++) {
    if(fabs(lp->orig_rhs[i]) < lp->infinite)
      lp->orig_rhs[i] *= scalechange[i];

    j = lp->presolve_undo->var_to_orig[i];
    if(j != 0)
      lp->presolve_undo->fixed_rhs[j] *= scalechange[i];

    if(lp->orig_upbo[i] < lp->infinite)     /* This is the range */
      lp->orig_upbo[i] *= scalechange[i];

    if((lp->orig_lowbo[i] != 0) && (fabs(lp->orig_lowbo[i]) < lp->infinite))
      lp->orig_lowbo[i] *= scalechange[i];
  }

  set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);

  return( TRUE );
}

STATIC REAL scale(lprec *lp, REAL *scaledelta)
{
  int     i, j, nz, row_count, nzOF = 0;
  REAL    *row_max, *row_min, *scalechange = NULL, absval;
  REAL    col_max, col_min;
  MYBOOL  rowscaled, colscaled;
  MATrec  *mat = lp->matA;
  REAL    *value;
  int     *rownr;

  if(is_scaletype(lp, SCALE_NONE))
    return(0.0);

  if(!lp->scaling_used) {
    allocREAL(lp, &lp->scalars, lp->sum_alloc + 1, FALSE);
    for(i = 0; i <= lp->sum; i++) {
      lp->scalars[i] = 1;
    }
    lp->scaling_used = TRUE;
  }
#ifdef Paranoia
  else
    for(i = 0; i <= lp->sum; i++) {
      if(lp->scalars[i] == 0)
        report(lp, SEVERE, "scale: Zero-valued scalar found at index %d\n", i);
    }
#endif
  if(scaledelta == NULL)
    allocREAL(lp, &scalechange, lp->sum + 1, FALSE);
  else
    scalechange = scaledelta;

 /* Must initialize due to computation of scaling statistic - KE */
  for(i = 0; i <= lp->sum; i++)
    scalechange[i] = 1;

  row_count = lp->rows;
  allocREAL(lp, &row_max, row_count + 1, TRUE);
  allocREAL(lp, &row_min, row_count + 1, FALSE);

  /* Initialise min and max values of rows */
  for(i = 0; i <= row_count; i++) {
    if(is_scaletype(lp, SCALE_MEAN))
      row_min[i] = 0;             /* Carries the count of elements */
    else
      row_min[i] = lp->infinite;  /* Carries the minimum element */
  }

  /* Calculate row scaling data */
  for(j = 1; j <= lp->columns; j++) {

    absval = lp->orig_obj[j];
    if(absval != 0) {
      absval = scaled_mat(lp, absval, 0, j);
      accumulate_for_scale(lp, &row_min[0], &row_max[0], absval);
      nzOF++;
    }

    i = mat->col_end[j - 1];
    value = &(COL_MAT_VALUE(i));
    rownr = &(COL_MAT_ROWNR(i));
    nz = mat->col_end[j];
    for(; i < nz;
        i++, value += matValueStep, rownr += matRowColStep) {
      absval = scaled_mat(lp, *value, *rownr, j);
      accumulate_for_scale(lp, &row_min[*rownr], &row_max[*rownr], absval);
    }
  }

  /* Calculate scale factors for rows */
  i = 0;
  for(; i <= lp->rows; i++) {
    if(i == 0)
      nz = nzOF;
    else
      nz = mat_rowlength(lp->matA, i);
    absval = minmax_to_scale(lp, row_min[i], row_max[i], nz); /* nz instead of nzOF KJEI 20/05/2010 */
    if(absval == 0)
      absval = 1;
    scalechange[i] = absval;
  }

  FREE(row_max);
  FREE(row_min);

  /* Row-scale the matrix (including the objective function and Lagrangean constraints) */
  rowscaled = scale_updaterows(lp, scalechange, TRUE);

  /* Calculate column scales */
  i = 1;
  for(j = 1; j <= lp->columns; j++) {
    if(is_int(lp,j) && !is_integerscaling(lp)) { /* do not scale integer columns */
      scalechange[lp->rows + j] = 1;
    }
    else {
      col_max = 0;
      if(is_scaletype(lp, SCALE_MEAN))
        col_min = 0;
      else
        col_min = lp->infinite;

      absval = lp->orig_obj[j];
      if(absval != 0) {
        absval = scaled_mat(lp, absval, 0, j);
        accumulate_for_scale(lp, &col_min, &col_max, absval);
      }

      i = mat->col_end[j - 1];
      value = &(COL_MAT_VALUE(i));
      rownr = &(COL_MAT_ROWNR(i));
      nz = mat->col_end[j];
      for(; i < nz;
          i++, value += matValueStep, rownr += matRowColStep) {
        absval = scaled_mat(lp, *value, *rownr, j);
        accumulate_for_scale(lp, &col_min, &col_max, absval);
      }
      nz = mat_collength(lp->matA, j);
      if(fabs(lp->orig_obj[j]) > 0)
        nz++;
      scalechange[lp->rows + j] = minmax_to_scale(lp, col_min, col_max, nz);
    }
  }

  /* ... and then column-scale the already row-scaled matrix */
  colscaled = scale_updatecolumns(lp, &scalechange[lp->rows], TRUE);

  /* Create a geometric mean-type measure of the extent of scaling performed; */
  /* ideally, upon successive calls to scale() the value should converge to 0 */
  if(rowscaled || colscaled) {
    col_max = 0;
    for(j = 1; j <= lp->columns; j++)
      col_max += log(scalechange[lp->rows + j]);
    col_max = exp(col_max/lp->columns);

    i = 0;
    col_min = 0;
    for(; i <= lp->rows; i++)
      col_min += log(scalechange[i]);
    col_min = exp(col_min/row_count);
  }
  else {
    col_max = 1;
    col_min = 1;
  }

  if(scaledelta == NULL)
    FREE(scalechange);

  return(1 - sqrt(col_max*col_min));
}

STATIC MYBOOL finalize_scaling(lprec *lp, REAL *scaledelta)
{
  int i;

  /* Check if we should equilibrate */
  if(is_scalemode(lp, SCALE_EQUILIBRATE) && !is_scaletype(lp, SCALE_CURTISREID)) {
    int oldmode;

    oldmode = lp->scalemode;
    lp->scalemode = SCALE_LINEAR + SCALE_EXTREME;
    scale(lp, scaledelta);
    lp->scalemode = oldmode;
  }

  /* Check if we should prevent rounding errors */
  if(is_scalemode(lp, SCALE_POWER2)) {
    REAL *scalars;
    if(scaledelta == NULL)
      scalars = lp->scalars;
    else
      scalars = scaledelta;

    for(i = 0; i <= lp->sum; i++)
      scalars[i] = roundPower2(scalars[i]);
  }

  /* Then transfer the scalars to the model's data */
  return( scale_rows(lp, scaledelta) && scale_columns(lp, scaledelta) );

}

STATIC REAL auto_scale(lprec *lp)
{
  int    n = 1;
  REAL   scalingmetric = 0, *scalenew = NULL;

  if(lp->scaling_used &&
     ((((lp->scalemode & SCALE_DYNUPDATE) == 0)) || (lp->bb_level > 0)))
    return( scalingmetric);

  if(lp->scalemode != SCALE_NONE) {

    /* Allocate array for incremental scaling if appropriate */
    if((lp->solvecount > 1) && (lp->bb_level < 1) &&
       ((lp->scalemode & SCALE_DYNUPDATE) != 0))
      allocREAL(lp, &scalenew, lp->sum + 1, FALSE);

    if(is_scaletype(lp, SCALE_CURTISREID)) {
      scalingmetric = scaleCR(lp, scalenew);
    }
    else {
      REAL scalinglimit, scalingdelta;
      int  count;

      /* Integer value of scalelimit holds the maximum number of iterations; default to 1 */
      count = (int) floor(lp->scalelimit);
      scalinglimit = lp->scalelimit;
      if((count == 0) || (scalinglimit == 0)) {
        if(scalinglimit > 0)
          count = DEF_SCALINGLIMIT;  /* A non-zero convergence has been given, default to max 5 iterations */
        else
          count = 1;
      }
      else
        scalinglimit -= count;

      /* Scale to desired relative convergence or iteration limit */
      n = 0;
      scalingdelta = 1.0;
      scalingmetric = 1.0;
      while((n < count) && (fabs(scalingdelta) > scalinglimit)) {
        n++;
        scalingdelta = scale(lp, scalenew);
        scalingmetric = scalingmetric*(1+scalingdelta);
      }
      scalingmetric -= 1;
    }
  }

  /* Update the inf norm of the elements of the matrix (excluding the OF) */
  mat_computemax(lp->matA);

  /* Check if we really have to do scaling */
  if(lp->scaling_used && (fabs(scalingmetric) >= lp->epsprimal))
    /* Ok, do it */
    finalize_scaling(lp, scalenew);

  else {

    /* Otherwise reset scaling variables */
    if(lp->scalars != NULL) {
      FREE(lp->scalars);
    }
    lp->scaling_used = FALSE;
    lp->columns_scaled = FALSE;
  }
  if(scalenew != NULL)
    FREE(scalenew);

  return(scalingmetric);
}

STATIC void unscale_columns(lprec *lp)
{
  int     i, j, nz;
  MATrec  *mat = lp->matA;
  REAL    *value;
  int     *rownr, *colnr;

  if(!lp->columns_scaled)
    return;

  /* Unscale OF */
  for(j = 1; j <= lp->columns; j++) {
    lp->orig_obj[j] = unscaled_mat(lp, lp->orig_obj[j], 0, j);
  }

  /* Unscale mat */
  mat_validate(mat);
  nz = get_nonzeros(lp);
  value = &(COL_MAT_VALUE(0));
  rownr = &(COL_MAT_ROWNR(0));
  colnr = &(COL_MAT_COLNR(0));
  for(j = 0; j < nz;
      j++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
    *value = unscaled_mat(lp, *value, *rownr, *colnr);
  }

  /* Unscale bounds as well */
  for(i = lp->rows + 1, j = 1; i <= lp->sum; i++, j++) {
    lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
    lp->orig_upbo[i]  = unscaled_value(lp, lp->orig_upbo[i], i);
    lp->sc_lobound[j]  = unscaled_value(lp, lp->sc_lobound[j], i);
  }

  for(i = lp->rows + 1; i<= lp->sum; i++)
    lp->scalars[i] = 1;

  lp->columns_scaled = FALSE;
  set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
}

void undoscale(lprec *lp)
{
  int     i, j, nz;
  MATrec  *mat = lp->matA;
  REAL    *value;
  int     *rownr, *colnr;

  if(lp->scaling_used) {

    /* Unscale the OF */
    for(j = 1; j <= lp->columns; j++) {
      lp->orig_obj[j] = unscaled_mat(lp, lp->orig_obj[j], 0, j);
    }

    /* Unscale the matrix */
    mat_validate(mat);
    nz = get_nonzeros(lp);
    value = &(COL_MAT_VALUE(0));
    rownr = &(COL_MAT_ROWNR(0));
    colnr = &(COL_MAT_COLNR(0));
    for(j = 0; j < nz;
        j++, value += matValueStep, rownr += matRowColStep, colnr += matRowColStep) {
      *value = unscaled_mat(lp, *value, *rownr, *colnr);
    }

    /* Unscale variable bounds */
    for(i = lp->rows + 1, j = 1; i <= lp->sum; i++, j++) {
      lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
      lp->orig_upbo[i]  = unscaled_value(lp, lp->orig_upbo[i], i);
      lp->sc_lobound[j]  = unscaled_value(lp, lp->sc_lobound[j], i);
    }

    /* Unscale the rhs, upper and lower bounds... */
    for(i = 0; i <= lp->rows; i++) {
      lp->orig_rhs[i] = unscaled_value(lp, lp->orig_rhs[i], i);
      j = lp->presolve_undo->var_to_orig[i];
      if(j != 0)
        lp->presolve_undo->fixed_rhs[j] = unscaled_value(lp, lp->presolve_undo->fixed_rhs[j], i);
      lp->orig_lowbo[i] = unscaled_value(lp, lp->orig_lowbo[i], i);
      lp->orig_upbo[i] = unscaled_value(lp, lp->orig_upbo[i], i);
    }

    FREE(lp->scalars);
    lp->scaling_used = FALSE;
    lp->columns_scaled = FALSE;

    set_action(&lp->spx_action, ACTION_REBASE | ACTION_REINVERT | ACTION_RECOMPUTE);
  }
}