ldbl_price.c

Go to the documentation of this file.

/****************************************************************************/
/*                                                                          */
/*  This file is part of QSopt_ex.                                          */
/*                                                                          */
/*  (c) Copyright 2006 by David Applegate, William Cook, Sanjeeb Dash,      */
/*  and Daniel Espinoza                                                     */
/*                                                                          */
/*  Sanjeeb Dash ownership of copyright in QSopt_ex is derived from his     */
/*  copyright in QSopt.                                                     */
/*                                                                          */
/*  This code may be used under the terms of the GNU General Public License */
/*  (Version 2.1 or later) as published by the Free Software Foundation.    */
/*                                                                          */
/*  Alternatively, use is granted for research purposes only.               */
/*                                                                          */
/*  It is your choice of which of these two licenses you are operating      */
/*  under.                                                                  */
/*                                                                          */
/*  We make no guarantees about the correctness or usefulness of this code. */
/*                                                                          */
/****************************************************************************/

/* RCS_INFO = "$RCSfile: ldbl_price.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
//static int TRACE = 0;
#include "qs_config.h"
#include "stddefs.h"
#include "ldbl_qsopt.h"
#include "ldbl_lpdefs.h"
#include "ldbl_fct.h"
#include "ldbl_price.h"
#include "ldbl_basis.h"
#include "ldbl_iqsutil.h"
#include "ldbl_dstruct.h"
#ifdef USEDMALLOC
#include "dmalloc.h"
#endif

#define  ldbl_MULTIP 1
#define  ldbl_PRICE_DEBUG 0

static void ldbl_update_d_scaleinf (
  ldbl_price_info * const p,
  ldbl_heap * const h,
  int const j,
  long double inf,
  int const prule),
  ldbl_update_p_scaleinf (
  ldbl_price_info * const p,
  ldbl_heap * const h,
  int const i,
  long double inf,
  int const prule);

static void ldbl_compute_dualI_inf (
  ldbl_lpinfo * const lp,
  int const j,
  long double * const inf),
  ldbl_compute_dualII_inf (
  ldbl_lpinfo * const lp,
  int const j,
  long double * const inf),
  ldbl_compute_primalI_inf (
  ldbl_lpinfo * const lp,
  int const i,
  long double * const inf),
  ldbl_compute_primalII_inf (
  ldbl_lpinfo * const lp,
  int const i,
  long double * const inf);

void ldbl_ILLprice_free_heap (
  ldbl_price_info * const pinf)
{
  ldbl_ILLheap_free (&(pinf->h));
}

int ldbl_ILLprice_build_heap (
  ldbl_price_info * const pinf,
  int const nkeys,
  long double * keylist)
{
  ldbl_ILLheap_init (&(pinf->h));
  ldbl_EGlpNumSet (pinf->htrigger,
              1.0 +
              (double) nkeys / (PARAM_HEAP_RATIO * ldbl_ILLutil_our_log2 (nkeys)));
  return ldbl_ILLheap_build (&(pinf->h), nkeys, keylist);
}

int ldbl_ILLprice_test_for_heap (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const nkeys,
  long double * keylist,
  int const algo,
  int const upd)
{
  ldbl_heap *const h = &(pinf->h);
  int rval = 0;
  long double ravg;

  if (upd != 0)
  {
    ldbl_EGlpNumInitVar (ravg);
    if (algo == PRIMAL_SIMPLEX)
      ldbl_EGlpNumCopy (ravg, lp->cnts->za_ravg);
    else
      ldbl_EGlpNumCopy (ravg, lp->cnts->y_ravg);
    if (ldbl_EGlpNumIsLeq (ravg, pinf->htrigger))
      pinf->hineff--;
    else
    {
      ldbl_EGlpNumDivUiTo (ravg, 2U);
      if (ldbl_EGlpNumIsLess (pinf->htrigger, ravg))
        pinf->hineff++;
    }
    ldbl_EGlpNumClearVar (ravg);
  }
  if (h->hexist == 0 && pinf->hineff <= 0)
  {
    rval = ldbl_ILLprice_build_heap (pinf, nkeys, keylist);
    CHECKRVALG (rval, CLEANUP);
  }
  else if (h->hexist != 0 && pinf->hineff >= PARAM_HEAP_UTRIGGER)
  {
    ldbl_ILLprice_free_heap (pinf);
    /*
     * printf ("freeing ldbl_heap ..\n");
     * printf ("iter = %d, ravg = %.2f, trigger = %.2f\n",
     * lp->cnts->tot_iter, ravg, pinf->htrigger);
     */
  }

CLEANUP:
  if (rval)
    ldbl_ILLprice_free_heap (pinf);
  return rval;
}

void ldbl_ILLprice_init_pricing_info (
  ldbl_price_info * const pinf)
{
  pinf->p_strategy = -1;
  pinf->d_strategy = -1;
  pinf->pI_price = -1;
  pinf->pII_price = -1;
  pinf->dI_price = -1;
  pinf->dII_price = -1;
  pinf->cur_price = -1;
  pinf->p_scaleinf = 0;
  pinf->d_scaleinf = 0;
  pinf->pdinfo.norms = 0;
  pinf->pdinfo.refframe = 0;
  pinf->psinfo.norms = 0;
  pinf->ddinfo.norms = 0;
  pinf->ddinfo.refframe = 0;
  pinf->dsinfo.norms = 0;
  pinf->dmpinfo.gstart = pinf->pmpinfo.gstart = 0;
  pinf->dmpinfo.gshift = pinf->pmpinfo.gshift = 0;
  pinf->dmpinfo.gsize = pinf->pmpinfo.gsize = 0;
  pinf->dmpinfo.bucket = pinf->pmpinfo.bucket = 0;
  pinf->dmpinfo.perm = pinf->pmpinfo.perm = 0;
  pinf->dmpinfo.infeas = pinf->pmpinfo.infeas = 0;
  ldbl_ILLheap_init (&(pinf->h));
  ldbl_EGlpNumZero (pinf->htrigger);
  pinf->hineff = 0;
}

void ldbl_ILLprice_free_pricing_info (
  ldbl_price_info * const pinf)
{
  ldbl_EGlpNumFreeArray (pinf->p_scaleinf);
  ldbl_EGlpNumFreeArray (pinf->d_scaleinf);
  ldbl_EGlpNumFreeArray (pinf->pdinfo.norms);
  ILL_IFFREE (pinf->pdinfo.refframe, int);
  ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
  ldbl_EGlpNumFreeArray (pinf->ddinfo.norms);
  ILL_IFFREE (pinf->ddinfo.refframe, int);
  ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);

  ldbl_ILLprice_free_mpartial_info (&(pinf->pmpinfo));
  ldbl_ILLprice_free_mpartial_info (&(pinf->dmpinfo));
  ldbl_ILLprice_free_heap (pinf);
}

int ldbl_ILLprice_build_pricing_info (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const phase)
{
  int rval = 0;
  int p_price = -1;
  int d_price = -1;

  switch (phase)
  {
  case PRIMAL_PHASEI:
    p_price = pinf->pI_price;
    break;
  case PRIMAL_PHASEII:
    p_price = pinf->pII_price;
    break;
  case DUAL_PHASEI:
    d_price = pinf->dI_price;
    break;
  case DUAL_PHASEII:
    d_price = pinf->dII_price;
    break;
  }

  if (p_price != -1)
  {
    pinf->cur_price = p_price;

    if (p_price == QS_PRICE_PDANTZIG || p_price == QS_PRICE_PDEVEX ||
        p_price == QS_PRICE_PSTEEP)
    {
      pinf->p_strategy = COMPLETE_PRICING;
      ldbl_EGlpNumFreeArray (pinf->d_scaleinf);
      pinf->d_scaleinf = ldbl_EGlpNumAllocArray (lp->nnbasic);
    }
    else if (p_price == QS_PRICE_PMULTPARTIAL)
      pinf->p_strategy = MULTI_PART_PRICING;

    switch (p_price)
    {
    case QS_PRICE_PDEVEX:
      if (pinf->pdinfo.norms)
        return rval;
      rval = ldbl_ILLprice_build_pdevex_norms (lp, &(pinf->pdinfo), 0);
      CHECKRVALG(rval,CLEANUP);
      break;
    case QS_PRICE_PSTEEP:
      if (pinf->psinfo.norms)
        return rval;
      rval = ldbl_ILLprice_build_psteep_norms (lp, &(pinf->psinfo));
      CHECKRVALG(rval,CLEANUP);
      break;
    case QS_PRICE_PMULTPARTIAL:
      rval = ldbl_ILLprice_build_mpartial_info (lp, pinf, COL_PRICING);
      CHECKRVALG(rval,CLEANUP);
      break;
    }
  }
  else if (d_price != -1)
  {
    pinf->cur_price = d_price;

    if (d_price == QS_PRICE_DDANTZIG || d_price == QS_PRICE_DSTEEP ||
        d_price == QS_PRICE_DDEVEX)
    {
      pinf->d_strategy = COMPLETE_PRICING;
      ldbl_EGlpNumFreeArray (pinf->p_scaleinf);
      pinf->p_scaleinf = ldbl_EGlpNumAllocArray (lp->nrows);
    }
    else if (d_price == QS_PRICE_DMULTPARTIAL)
      pinf->d_strategy = MULTI_PART_PRICING;

    switch (d_price)
    {
    case QS_PRICE_DSTEEP:
      if (pinf->dsinfo.norms)
        return rval;
      rval = ldbl_ILLprice_build_dsteep_norms (lp, &(pinf->dsinfo));
      CHECKRVALG(rval,CLEANUP);
      break;
    case QS_PRICE_DMULTPARTIAL:
      rval = ldbl_ILLprice_build_mpartial_info (lp, pinf, ROW_PRICING);
      CHECKRVALG(rval,CLEANUP);
      break;
    case QS_PRICE_DDEVEX:
      if (pinf->ddinfo.norms)
        return rval;
      rval = ldbl_ILLprice_build_ddevex_norms (lp, &(pinf->ddinfo), 0);
      CHECKRVALG(rval,CLEANUP);
      break;
    }
  }

CLEANUP:
  if (rval)
    ldbl_ILLprice_free_pricing_info (pinf);
  EG_RETURN(rval);
}

int ldbl_ILLprice_update_pricing_info (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const phase,
  ldbl_svector * const wz,
  int const eindex,
  int const lindex,
  long double y)
{
  int rval = 0;
  int p_price = -1;
  int d_price = -1;

  switch (phase)
  {
  case PRIMAL_PHASEI:
    p_price = pinf->pI_price;
    break;
  case PRIMAL_PHASEII:
    p_price = pinf->pII_price;
    break;
  case DUAL_PHASEI:
    d_price = pinf->dI_price;
    break;
  case DUAL_PHASEII:
    d_price = pinf->dII_price;
    break;
  }

  if (p_price != -1)
  {
    if (p_price == QS_PRICE_PDEVEX)
    {
      rval = ldbl_ILLprice_update_pdevex_norms (lp, &(pinf->pdinfo), eindex, y);
      CHECKRVALG(rval,CLEANUP);
    }
    else if (p_price == QS_PRICE_PSTEEP)
      ldbl_ILLprice_update_psteep_norms (lp, &(pinf->psinfo), wz, eindex, y);
  }
  else if (d_price != -1)
  {
    if (d_price == QS_PRICE_DSTEEP)
      ldbl_ILLprice_update_dsteep_norms (lp, &(pinf->dsinfo), wz, lindex, y);
    else if (d_price == QS_PRICE_DDEVEX)
    {
      rval = ldbl_ILLprice_update_ddevex_norms (lp, &(pinf->ddinfo), lindex, y);
      CHECKRVALG(rval,CLEANUP);
    }
  }
CLEANUP:
  EG_RETURN(rval);
}

int ldbl_ILLprice_get_price (
  ldbl_price_info * const p,
  int const phase)
{
  int pri = -1;

  switch (phase)
  {
  case PRIMAL_PHASEI:
    return p->pI_price;
  case PRIMAL_PHASEII:
    return p->pII_price;
  case DUAL_PHASEI:
    return p->dI_price;
  case DUAL_PHASEII:
    return p->dII_price;
  }
  return pri;
}

void ldbl_ILLprice_free_mpartial_info (
  ldbl_mpart_info * p)
{
  ILL_IFFREE (p->gstart, int);
  ILL_IFFREE (p->gshift, int);
  ILL_IFFREE (p->gsize, int);
  ILL_IFFREE (p->bucket, int);
  ldbl_EGlpNumFreeArray (p->infeas);
  ILL_IFFREE (p->perm, int);
}

int ldbl_ILLprice_build_mpartial_info (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const pricetype)
{
  int i = 0;
  int rval = 0;
  int extra = 0;
  int nelems;
  ldbl_mpart_info *p;

  p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
  p->k = 50;
  p->cgroup = 0;
  nelems = (pricetype == COL_PRICING) ? lp->nnbasic : lp->nrows;

  if (nelems % p->k)
    extra = nelems - p->k * (nelems / p->k);
  p->ngroups = nelems / p->k;
  if (extra != 0)
    p->ngroups++;

  ILL_SAFE_MALLOC (p->gstart, p->ngroups, int);
  ILL_SAFE_MALLOC (p->gshift, p->ngroups, int);
  ILL_SAFE_MALLOC (p->gsize, p->ngroups, int);
  ILL_SAFE_MALLOC (p->bucket, 2 * p->k, int);
  p->infeas = ldbl_EGlpNumAllocArray (2 * p->k);
  ILL_SAFE_MALLOC (p->perm, 2 * p->k, int);

  p->bsize = 0;

  if (extra != 0)
  {
    p->gstart[0] = 0;
    p->gshift[0] = 1;
    p->gsize[0] = extra;
    for (i = 1; i < p->ngroups; i++)
    {
      p->gstart[i] = extra + i - 1;
      p->gshift[i] = p->ngroups - 1;
      p->gsize[i] = p->k;
    }
  }
  else
  {
    for (i = 0; i < p->ngroups; i++)
    {
      p->gstart[i] = i;
      p->gshift[i] = p->ngroups;
      p->gsize[i] = p->k;
    }
  }

CLEANUP:
  if (rval)
    ldbl_ILLprice_free_mpartial_info (p);
  EG_RETURN(rval);
}

void ldbl_ILLprice_init_mpartial_price (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const phase,
  int const pricetype)
{
  int i;
  ldbl_mpart_info *p;

  p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
  p->bsize = 0;
  i = p->cgroup;
  do
  {
    ldbl_ILLprice_mpartial_group (lp, p, phase, i, pricetype);
    i = (i + 1) % p->ngroups;
  } while (i != p->cgroup && p->bsize <= p->k);
  p->cgroup = i;
}

void ldbl_ILLprice_update_mpartial_price (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const phase,
  int const pricetype)
{
  int i = 0;
  int csize = 0;
  long double infeas;
  ldbl_mpart_info *p;
  ldbl_price_res pr;

  p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
  ldbl_EGlpNumInitVar (pr.dinfeas);
  ldbl_EGlpNumInitVar (pr.pinfeas);
  ldbl_EGlpNumInitVar (infeas);

#ifdef ldbl_MULTIP
  i = 0;
  while (i < p->bsize)
  {
    if (pricetype == COL_PRICING)
    {
      ldbl_ILLprice_column (lp, p->bucket[i], phase, &pr);
      ldbl_EGlpNumCopy (infeas, pr.dinfeas);
    }
    else
    {
      ldbl_ILLprice_row (lp, p->bucket[i], phase, &pr);
      ldbl_EGlpNumCopy (infeas, pr.pinfeas);
    }
    if (!ldbl_EGlpNumIsNeqqZero (infeas))
    {
      p->bucket[i] = p->bucket[p->bsize - 1];
      p->bsize--;
    }
    else
    {
      ldbl_EGlpNumCopy (p->infeas[i], infeas);
      i++;
    }
  }
  if (p->bsize > 0)
  {
    for (i = 0; i < p->bsize; i++)
      p->perm[i] = i;
    ldbl_EGutilPermSort ((size_t) (p->bsize), p->perm,
                    (const long double * const) p->infeas);

    csize = QSMIN (p->bsize, p->k);
    for (i = csize - 1; i >= 0; i--)
      lp->iwork[p->bucket[p->perm[i]]] = 1;

    for (i = 0, csize = 0; i < p->bsize; i++)
      if (lp->iwork[p->bucket[i]] == 1)
      {
        ldbl_EGlpNumCopy (p->infeas[csize], p->infeas[i]);
        p->bucket[csize] = p->bucket[i];
        csize++;
      }
    p->bsize = csize;
  }
#else
  p->bsize = 0;
#endif

  i = p->cgroup;
  do
  {
    ldbl_ILLprice_mpartial_group (lp, p, phase, i, pricetype);
    i = (i + 1) % p->ngroups;
  } while (i != p->cgroup && p->bsize <= p->k);
  p->cgroup = i;

#ifdef ldbl_MULTIP
  for (i = 0; i < csize; i++)
    lp->iwork[p->bucket[i]] = 0;
#endif
  ldbl_EGlpNumClearVar (infeas);
  ldbl_EGlpNumClearVar (pr.pinfeas);
  ldbl_EGlpNumClearVar (pr.dinfeas);
}

void ldbl_ILLprice_delete_onempart_price (
  /*ldbl_lpinfo * const lp,*/
  ldbl_price_info * const pinf,
  int const indx,
  int const pricetype)
{
  int i = 0;
  ldbl_mpart_info *p;

  p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);

  for (i = 0; i < p->bsize; i++)
    if (p->bucket[i] == indx)
    {
      p->bucket[i] = p->bucket[p->bsize - 1];
      ldbl_EGlpNumCopy (p->infeas[i], p->infeas[p->bsize - 1]);
      p->bsize--;
      break;
    }
}

void ldbl_ILLprice_mpartial_group (
  ldbl_lpinfo * const lp,
  ldbl_mpart_info * const p,
  int const phase,
  int const g,
  int const pricetype)
{
  int i, ix;
  int gstart = p->gstart[g];
  int gsize = p->gsize[g];
  int gshift = p->gshift[g];
  long double infeas;
  ldbl_price_res pr;

  ldbl_EGlpNumInitVar (pr.dinfeas);
  ldbl_EGlpNumInitVar (pr.pinfeas);
  ldbl_EGlpNumInitVar (infeas);

  for (i = 0, ix = gstart; i < gsize; i++, ix += gshift)
  {
#ifdef ldbl_MULTIP
    if (lp->iwork[ix])
      continue;
#endif
    if (pricetype == COL_PRICING)
    {
      ldbl_ILLprice_column (lp, ix, phase, &pr);
      ldbl_EGlpNumCopy (infeas, pr.dinfeas);
    }
    else
    {
      ldbl_ILLprice_row (lp, ix, phase, &pr);
      ldbl_EGlpNumCopy (infeas, pr.pinfeas);
    }
    if (ldbl_EGlpNumIsNeqqZero (infeas))
    {
      ldbl_EGlpNumCopy (p->infeas[p->bsize], infeas);
      p->bucket[p->bsize] = ix;
      p->bsize++;
    }
  }
  ldbl_EGlpNumClearVar (infeas);
  ldbl_EGlpNumClearVar (pr.dinfeas);
  ldbl_EGlpNumClearVar (pr.pinfeas);
}

void ldbl_ILLprice_column (
  ldbl_lpinfo * const lp,
  int const ix,
  int const phase,
  ldbl_price_res * const pr)
{
  int i;
  int col;
  int mcnt;
  int mbeg;
  long double sum;

  ldbl_EGlpNumZero (pr->dinfeas);
  col = lp->nbaz[ix];
  if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFIXED)
    return;
  ldbl_EGlpNumInitVar (sum);
  ldbl_EGlpNumZero (sum);
  mcnt = lp->matcnt[col];
  mbeg = lp->matbeg[col];

  if (phase == PRIMAL_PHASEII)
  {
    for (i = 0; i < mcnt; i++)
      ldbl_EGlpNumAddInnProdTo (sum, lp->piz[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
    ldbl_EGlpNumCopyDiff (lp->dz[ix], lp->cz[col], sum);
    ldbl_compute_dualII_inf (lp, ix, &(pr->dinfeas));
  }
  else
  {
    for (i = 0; i < mcnt; i++)
      ldbl_EGlpNumAddInnProdTo (sum, lp->pIpiz[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
    ldbl_EGlpNumCopyNeg (lp->pIdz[ix], sum);
    ldbl_compute_dualI_inf (lp, ix, &(pr->dinfeas));
  }
  ldbl_EGlpNumClearVar (sum);
}

void ldbl_ILLprice_row (
  ldbl_lpinfo * const lp,
  int const ix,
  int const phase,
  ldbl_price_res * const pr)
{
  if (phase == DUAL_PHASEII)
    ldbl_compute_primalII_inf (lp, ix, &(pr->pinfeas));
  else
    ldbl_compute_primalI_inf (lp, ix, &(pr->pinfeas));
}

int ldbl_ILLprice_build_pdevex_norms (
  ldbl_lpinfo * const lp,
  ldbl_p_devex_info * const pdinfo,
  int const reinit)
{
  int j;
  int rval = 0;

  if (reinit == 0)
  {
    pdinfo->ninit = 0;
    pdinfo->norms = ldbl_EGlpNumAllocArray (lp->nnbasic);
    ILL_SAFE_MALLOC (pdinfo->refframe, lp->ncols, int);
  }

  if (reinit != 0)
    pdinfo->ninit++;

  for (j = 0; j < lp->ncols; j++)
  {
    if (lp->vstat[j] == STAT_BASIC)
      pdinfo->refframe[j] = 0;
    else
    {
      ldbl_EGlpNumOne (pdinfo->norms[lp->vindex[j]]);
      pdinfo->refframe[j] = 1;
    }
  }

CLEANUP:
  if (rval)
  {
    ldbl_EGlpNumFreeArray (pdinfo->norms);
    ILL_IFFREE (pdinfo->refframe, int);
  }
  EG_RETURN(rval);
}

int ldbl_ILLprice_update_pdevex_norms (
  ldbl_lpinfo * const lp,
  ldbl_p_devex_info * const pdinfo,
  int const eindex,
  long double yl)
{
  int i, j;
  long double normj;
  long double zAj;
  long double ntmp, ntmp2;

  ldbl_EGlpNumInitVar (normj);
  ldbl_EGlpNumInitVar (zAj);
  ldbl_EGlpNumInitVar (ntmp);
  ldbl_EGlpNumInitVar (ntmp2);
  ldbl_EGlpNumZero (normj);

  for (i = 0; i < lp->yjz.nzcnt; i++)
    if (pdinfo->refframe[lp->baz[lp->yjz.indx[i]]])
      ldbl_EGlpNumAddInnProdTo (normj, lp->yjz.coef[i], lp->yjz.coef[i]);

  if (pdinfo->refframe[lp->nbaz[eindex]])
    ldbl_EGlpNumAddTo (normj, ldbl_oneLpNum);

  ldbl_EGlpNumSet(ntmp,1000.0);
  ldbl_EGlpNumSet(ntmp2,0.001);
  ldbl_EGlpNumMultTo(ntmp,pdinfo->norms[eindex]);
  ldbl_EGlpNumMultTo(ntmp2,pdinfo->norms[eindex]);
  if (ldbl_EGlpNumIsLess (normj, ntmp2) || ldbl_EGlpNumIsLess (ntmp, normj))
  {
    ldbl_EGlpNumClearVar (zAj);
    ldbl_EGlpNumClearVar (normj);
    ldbl_EGlpNumClearVar (ntmp);
    ldbl_EGlpNumClearVar (ntmp2);
    return ldbl_ILLprice_build_pdevex_norms (lp, pdinfo, 1);
  }

  for (i = 0; i < lp->zA.nzcnt; i++)
  {
    j = lp->zA.indx[i];
    ldbl_EGlpNumCopyFrac (zAj, lp->zA.coef[i], yl);
    ldbl_EGlpNumMultTo (zAj, zAj);
    ldbl_EGlpNumMultTo (zAj, normj);
    if (ldbl_EGlpNumIsLess (pdinfo->norms[j], zAj))
      ldbl_EGlpNumCopy (pdinfo->norms[j], zAj);
  }
  ldbl_EGlpNumDivTo (normj, yl);
  ldbl_EGlpNumDivTo (normj, yl);
  if (ldbl_EGlpNumIsLess (normj, ldbl_oneLpNum))
    ldbl_EGlpNumCopy (pdinfo->norms[eindex], ldbl_oneLpNum);
  else
    ldbl_EGlpNumCopy (pdinfo->norms[eindex], normj);
  ldbl_EGlpNumClearVar (zAj);
  ldbl_EGlpNumClearVar (normj);
  ldbl_EGlpNumClearVar (ntmp);
  ldbl_EGlpNumClearVar (ntmp2);
  return 0;
}

int ldbl_ILLprice_build_psteep_norms (
  ldbl_lpinfo * const lp,
  ldbl_p_steep_info * const psinfo)
{
  int j;
  int rval = 0;
  ldbl_svector yz;

  ldbl_ILLsvector_init (&yz);
  rval = ldbl_ILLsvector_alloc (&yz, lp->nrows);
  CHECKRVALG(rval,CLEANUP);
  psinfo->norms = ldbl_EGlpNumAllocArray (lp->nnbasic);

  for (j = 0; j < lp->nnbasic; j++)
  {
    rval = ILLstring_report (NULL, &lp->O->reporter);
    CHECKRVALG(rval,CLEANUP);
    ldbl_ILLfct_compute_yz (lp, &yz, 0, lp->nbaz[j]);
    ldbl_EGlpNumInnProd (psinfo->norms[j], yz.coef, yz.coef, (size_t) yz.nzcnt);
    ldbl_EGlpNumAddTo (psinfo->norms[j], ldbl_oneLpNum);
  }

CLEANUP:
  ldbl_ILLsvector_free (&yz);
  if (rval)
    ldbl_EGlpNumFreeArray (psinfo->norms);

  EG_RETURN(rval);
}

void ldbl_ILLprice_update_psteep_norms (
  ldbl_lpinfo * const lp,
  ldbl_p_steep_info * const psinfo,
  ldbl_svector * const wz,
  int const eindex,
  long double yl)
{
  int i, j, k;
  int mcnt, mbeg;
  long double normj,ntmp;
  long double zAj, wAj;
  long double *v = 0;

  ldbl_EGlpNumInitVar (normj);
  ldbl_EGlpNumInitVar (zAj);
  ldbl_EGlpNumInitVar (ntmp);
  ldbl_EGlpNumInitVar (wAj);
  ldbl_EGlpNumInnProd (normj, lp->yjz.coef, lp->yjz.coef, (size_t) (lp->yjz.nzcnt));
  ldbl_EGlpNumAddTo (normj, ldbl_oneLpNum);

#if 0
  Bico - remove warnings for dist
    if (fabs ((normj - psinfo->norms[eindex]) / normj) > 1000.0 /* 0.01 */ )
    {
      printf ("warning: incorrect norm values\n");
      printf ("anorm = %.6f, pnorm = %.6f\n", normj, psinfo->norms[eindex]);
      fflush (stdout);
    }
#endif

  ldbl_ILLfct_load_workvector (lp, wz);
  v = lp->work.coef;

  for (k = 0; k < lp->zA.nzcnt; k++)
  {
    j = lp->zA.indx[k];
    ldbl_EGlpNumCopy (zAj, lp->zA.coef[k]);
    ldbl_EGlpNumZero (wAj);
    mcnt = lp->matcnt[lp->nbaz[j]];
    mbeg = lp->matbeg[lp->nbaz[j]];
    for (i = 0; i < mcnt; i++)
      ldbl_EGlpNumAddInnProdTo (wAj, lp->matval[mbeg + i], v[lp->matind[mbeg + i]]);

    /* compute ntmp = (zAj * ((zAj * normj / yl) - (2.0 * wAj))) / yl; */ 
    ldbl_EGlpNumCopy(ntmp,zAj);
    ldbl_EGlpNumMultTo(ntmp,normj);
    ldbl_EGlpNumDivTo(ntmp,yl);
    ldbl_EGlpNumSubTo(ntmp,wAj);
    ldbl_EGlpNumSubTo(ntmp,wAj);
    ldbl_EGlpNumMultTo(ntmp,zAj);
    ldbl_EGlpNumDivTo(ntmp,yl);
    /* set psinfo->norms[j] += (zAj * ((zAj * normj / yl) - (2.0 * wAj))) / yl; */
    ldbl_EGlpNumAddTo(psinfo->norms[j],ntmp);
    if (ldbl_EGlpNumIsLess (psinfo->norms[j], ldbl_oneLpNum))
      ldbl_EGlpNumOne (psinfo->norms[j]);
  }

  ldbl_EGlpNumCopyFrac (psinfo->norms[eindex], normj, yl);
  ldbl_EGlpNumDivTo (psinfo->norms[eindex], yl);
  if (ldbl_EGlpNumIsLess (psinfo->norms[eindex], ldbl_oneLpNum))
    ldbl_EGlpNumOne (psinfo->norms[eindex]);

  ldbl_ILLfct_zero_workvector (lp);
  ldbl_EGlpNumClearVar (wAj);
  ldbl_EGlpNumClearVar (zAj);
  ldbl_EGlpNumClearVar (normj);
  ldbl_EGlpNumClearVar (ntmp);
}

int ldbl_ILLprice_build_ddevex_norms (
  ldbl_lpinfo * const lp,
  ldbl_d_devex_info * const ddinfo,
  int const reinit)
{
  int i;
  int rval = 0;

  if (reinit == 0)
  {
    ddinfo->ninit = 0;
    ddinfo->norms = ldbl_EGlpNumAllocArray (lp->nrows);
    ILL_SAFE_MALLOC (ddinfo->refframe, lp->ncols, int);
  }
  if (reinit != 0)
    ddinfo->ninit++;

  for (i = 0; i < lp->ncols; i++)
    ddinfo->refframe[i] = (lp->vstat[i] == STAT_BASIC) ? 1 : 0;

  for (i = 0; i < lp->nrows; i++)
    ldbl_EGlpNumOne (ddinfo->norms[i]);

CLEANUP:
  if (rval)
  {
    ldbl_EGlpNumFreeArray (ddinfo->norms);
    ILL_IFFREE (ddinfo->refframe, int);
  }
  EG_RETURN(rval);
}

int ldbl_ILLprice_update_ddevex_norms (
  ldbl_lpinfo * const lp,
  ldbl_d_devex_info * const ddinfo,
  int const lindex,
  long double yl)
{
  int i, r;
  long double normi;
  long double yr;
  long double ntmp,ntmp2;

  ldbl_EGlpNumInitVar (ntmp);
  ldbl_EGlpNumInitVar (ntmp2);
  ldbl_EGlpNumInitVar (normi);
  ldbl_EGlpNumInitVar (yr);
  ldbl_EGlpNumZero (normi);

  for (i = 0; i < lp->zA.nzcnt; i++)
    if (ddinfo->refframe[lp->nbaz[lp->zA.indx[i]]])
      ldbl_EGlpNumAddInnProdTo (normi, lp->zA.coef[i], lp->zA.coef[i]);

  if (ddinfo->refframe[lp->baz[lindex]])
    ldbl_EGlpNumAddTo (normi, ldbl_oneLpNum);

  ldbl_EGlpNumSet(ntmp,1000.0);
  ldbl_EGlpNumSet(ntmp2,0.001);
  ldbl_EGlpNumMultTo(ntmp,ddinfo->norms[lindex]);
  ldbl_EGlpNumMultTo(ntmp2,ddinfo->norms[lindex]);
  if (ldbl_EGlpNumIsLess(normi, ntmp2) || ldbl_EGlpNumIsLess(ntmp, normi))
  {
    ldbl_EGlpNumClearVar (normi);
    ldbl_EGlpNumClearVar (yr);
    ldbl_EGlpNumClearVar (ntmp);
    ldbl_EGlpNumClearVar (ntmp2);
    return ldbl_ILLprice_build_ddevex_norms (lp, ddinfo, 1);
  }

  for (i = 0; i < lp->yjz.nzcnt; i++)
  {
    r = lp->yjz.indx[i];
    ldbl_EGlpNumCopy(yr, lp->yjz.coef[i]);
    ldbl_EGlpNumCopy(ntmp,yr);
    ldbl_EGlpNumMultTo(ntmp,yr);
    ldbl_EGlpNumMultTo(ntmp,normi);
    ldbl_EGlpNumDivTo(ntmp,yl);
    ldbl_EGlpNumDivTo(ntmp,yl);
    if (ldbl_EGlpNumIsLess (ddinfo->norms[r], ntmp))
      ldbl_EGlpNumCopy (ddinfo->norms[r], ntmp);
  }
  ldbl_EGlpNumCopy (ddinfo->norms[lindex], normi);
  ldbl_EGlpNumDivTo(ddinfo->norms[lindex], yl);
  ldbl_EGlpNumDivTo(ddinfo->norms[lindex], yl);
  if (ldbl_EGlpNumIsLess (ddinfo->norms[lindex], ldbl_oneLpNum))
    ldbl_EGlpNumOne (ddinfo->norms[lindex]);
  ldbl_EGlpNumClearVar (normi);
  ldbl_EGlpNumClearVar (yr);
  ldbl_EGlpNumClearVar (ntmp);
  ldbl_EGlpNumClearVar (ntmp2);
  return 0;
}

int ldbl_ILLprice_build_dsteep_norms (
  ldbl_lpinfo * const lp,
  ldbl_d_steep_info * const dsinfo)
{
  int i;
  int rval = 0;
  ldbl_svector z;

  ldbl_ILLsvector_init (&z);
  rval = ldbl_ILLsvector_alloc (&z, lp->nrows);
  CHECKRVALG(rval,CLEANUP);
  dsinfo->norms = ldbl_EGlpNumAllocArray (lp->nrows);

  for (i = 0; i < lp->nrows; i++)
  {
    rval = ILLstring_report (NULL, &lp->O->reporter);
    CHECKRVALG(rval,CLEANUP);

    ldbl_ILLfct_compute_zz (lp, &z, i);

    ldbl_EGlpNumInnProd (dsinfo->norms[i], z.coef, z.coef, (size_t) z.nzcnt);
    if (ldbl_EGlpNumIsLess (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM))
      ldbl_EGlpNumCopy (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM);
  }

CLEANUP:
  ldbl_ILLsvector_free (&z);
  if (rval)
    ldbl_EGlpNumFreeArray (dsinfo->norms);

  EG_RETURN(rval);
}

int ldbl_ILLprice_get_dsteep_norms (
  ldbl_lpinfo * const lp,
  int const count,
  int *const rowind,
  long double * const norms)
{
  int i;
  int rval = 0;
  ldbl_svector z;

  ldbl_ILLsvector_init (&z);
  rval = ldbl_ILLsvector_alloc (&z, lp->nrows);
  CHECKRVALG(rval,CLEANUP);

  for (i = 0; i < count; i++)
  {
    ldbl_ILLfct_compute_zz (lp, &z, rowind[i]);
    ldbl_EGlpNumInnProd (norms[i], z.coef, z.coef, (size_t) z.nzcnt);
  }

CLEANUP:
  ldbl_ILLsvector_free (&z);
  EG_RETURN(rval);
}

void ldbl_ILLprice_update_dsteep_norms (
  ldbl_lpinfo * const lp,
  ldbl_d_steep_info * const dsinfo,
  ldbl_svector * const wz,
  int const lindex,
  long double yl)
{
  int i, k;
  long double yij;
  long double norml;
  long double *v = 0;
  long double ntmp;

  ldbl_EGlpNumInitVar (ntmp);
  ldbl_EGlpNumInitVar (norml);
  ldbl_EGlpNumInitVar (yij);
  ldbl_EGlpNumInnProd (norml, lp->zz.coef, lp->zz.coef, (size_t) (lp->zz.nzcnt));

#if 0
  Bico - remove warnings for dist
    if (fabs ((norml - dsinfo->norms[lindex]) / norml) > 1000.0 /*0.01 */ )
    {
      printf ("warning: incorrect dnorm values\n");
      printf ("anorm = %.6f, pnorm = %.6f\n", norml, dsinfo->norms[lindex]);
      fflush (stdout);
    }
#endif

  ldbl_ILLfct_load_workvector (lp, wz);
  v = lp->work.coef;

  for (k = 0; k < lp->yjz.nzcnt; k++)
  {
    i = lp->yjz.indx[k];
    ldbl_EGlpNumCopy (yij, lp->yjz.coef[k]);
    /* compute in ntmp (yij * ((yij * norml / yl) - (2.0 * v[i]))) / yl; */
    ldbl_EGlpNumCopy(ntmp,yij);
    ldbl_EGlpNumMultTo(ntmp,norml);
    ldbl_EGlpNumDivTo(ntmp,yl);
    ldbl_EGlpNumSubTo(ntmp,v[i]);
    ldbl_EGlpNumSubTo(ntmp,v[i]);
    ldbl_EGlpNumMultTo (ntmp, yij);
    ldbl_EGlpNumDivTo (ntmp, yl);
    /* set dsinfo->norms[i] += (yij * ((yij * norml / yl) - (2.0 * v[i]))) / yl;*/
    ldbl_EGlpNumAddTo(dsinfo->norms[i], ntmp);
    if (ldbl_EGlpNumIsLess (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM))
      ldbl_EGlpNumCopy (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM);
  }
  ldbl_EGlpNumCopyFrac (dsinfo->norms[lindex], norml, yl);
  ldbl_EGlpNumDivTo (dsinfo->norms[lindex], yl);
  if (ldbl_EGlpNumIsLess (dsinfo->norms[lindex], ldbl_PARAM_MIN_DNORM))
    ldbl_EGlpNumCopy (dsinfo->norms[lindex], ldbl_PARAM_MIN_DNORM);

  ldbl_ILLfct_zero_workvector (lp);
  ldbl_EGlpNumClearVar (norml);
  ldbl_EGlpNumClearVar (ntmp);
  ldbl_EGlpNumClearVar (yij);
}

static void ldbl_update_d_scaleinf (
  ldbl_price_info * const p,
  ldbl_heap * const h,
  int const j,
  long double inf,
  int const prule)
{
  if (!ldbl_EGlpNumIsNeqqZero (inf))
  {
    ldbl_EGlpNumZero (p->d_scaleinf[j]);
    if (h->hexist != 0 && h->loc[j] != -1)
      ldbl_ILLheap_delete (h, j);
  }
  else
  {
    if (prule == QS_PRICE_PDANTZIG)
      ldbl_EGlpNumCopy (p->d_scaleinf[j], inf);
    else if (prule == QS_PRICE_PDEVEX)
      ldbl_EGlpNumCopySqrOver (p->d_scaleinf[j], inf, p->pdinfo.norms[j]);
    else if (prule == QS_PRICE_PSTEEP)
      ldbl_EGlpNumCopySqrOver (p->d_scaleinf[j], inf, p->psinfo.norms[j]);

    if (h->hexist != 0)
    {
      if (h->loc[j] == -1)
        ldbl_ILLheap_insert (h, j);
      else
        ldbl_ILLheap_modify (h, j);
    }
  }
}

static void ldbl_compute_dualI_inf (
  ldbl_lpinfo * const lp,
  const int j,
  long double * const inf)
{
  int col = lp->nbaz[j];
  int vt = lp->vtype[col];
  int vs = lp->vstat[col];
  long double*dj = &(lp->pIdz[j]);
  long double*ftol = &(lp->tol->id_tol);
  ldbl_EGlpNumZero (*inf);
  if (vt != VARTIFICIAL && vt != VFIXED)
  {
    if( ldbl_EGlpNumIsSumLess(*dj,*ftol,ldbl_zeroLpNum) && (vs == STAT_LOWER || vs == STAT_ZERO))
      ldbl_EGlpNumCopyNeg(*inf,*dj);
    else if (ldbl_EGlpNumIsLess(*ftol, *dj) && (vs == STAT_UPPER || vs == STAT_ZERO))
      ldbl_EGlpNumCopy (*inf, *dj);
  }
}

static void ldbl_compute_dualII_inf (
  ldbl_lpinfo * const lp,
  int const j,
  long double * const inf)
{
  int col = lp->nbaz[j];
  int vt = lp->vtype[col];
  int vs = lp->vstat[col];
  long double*dj = &(lp->dz[j]);
  long double*ftol = &(lp->tol->dfeas_tol);
  ldbl_EGlpNumZero (*inf);
  if (vt != VARTIFICIAL && vt != VFIXED)
  {
    if( ldbl_EGlpNumIsSumLess(*dj,*ftol,ldbl_zeroLpNum) && (vs == STAT_LOWER || vs == STAT_ZERO))
      ldbl_EGlpNumCopyNeg(*inf,*dj);
    else if (ldbl_EGlpNumIsLess(*ftol,*dj) && (vs == STAT_UPPER || vs == STAT_ZERO))
      ldbl_EGlpNumCopy (*inf, *dj);
  }
}

void ldbl_ILLprice_compute_dual_inf (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const p,
  int *const ix,
  int const icnt,
  int const phase)
{
  int i;
  int price;
  long double inf;
  ldbl_heap *h = &(p->h);

  price = (phase == PRIMAL_PHASEI) ? p->pI_price : p->pII_price;
  ldbl_EGlpNumInitVar (inf);
  ldbl_EGlpNumZero (inf);

  if (phase == PRIMAL_PHASEI)
  {
    if (ix == NULL)
      for (i = 0; i < lp->nnbasic; i++)
      {
        ldbl_compute_dualI_inf (lp, i, &(inf));
        ldbl_update_d_scaleinf (p, h, i, inf, price);
      }
    else
      for (i = 0; i < icnt; i++)
      {
        ldbl_compute_dualI_inf (lp, ix[i], &(inf));
        ldbl_update_d_scaleinf (p, h, ix[i], inf, price);
      }
  }
  else if (phase == PRIMAL_PHASEII)
  {
    if (ix == NULL)
      for (i = 0; i < lp->nnbasic; i++)
      {
        ldbl_compute_dualII_inf (lp, i, &inf);
        ldbl_update_d_scaleinf (p, h, i, inf, price);
      }
    else
      for (i = 0; i < icnt; i++)
      {
        ldbl_compute_dualII_inf (lp, ix[i], &inf);
        ldbl_update_d_scaleinf (p, h, ix[i], inf, price);
      }
  }
  ldbl_EGlpNumClearVar (inf);
}

void ldbl_ILLprice_primal (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  ldbl_price_res * const pr,
  int const phase)
{
  int j, vs;
  long double d_e, d_max;
  long double *ftol = &(lp->tol->dfeas_tol);
  ldbl_heap *const h = &(pinf->h);

  ldbl_EGlpNumInitVar (d_e);
  ldbl_EGlpNumInitVar (d_max);
  pr->eindex = -1;
  ldbl_EGlpNumZero(d_max);

#if USEHEAP > 0
  ldbl_ILLprice_test_for_heap (lp, pinf, lp->nnbasic, pinf->d_scaleinf,
                          PRIMAL_SIMPLEX, 1);
#endif

  if (pinf->p_strategy == COMPLETE_PRICING)
  {
    if (h->hexist)
    {
      pr->eindex = ldbl_ILLheap_findmin (h);
      if (pr->eindex != -1)
        ldbl_ILLheap_delete (h, pr->eindex);
    }
    else
    {
      for (j = 0; j < lp->nnbasic; j++)
      {
        if (ldbl_EGlpNumIsLess (d_max, pinf->d_scaleinf[j]))
        {
          ldbl_EGlpNumCopy (d_max, pinf->d_scaleinf[j]);
          pr->eindex = j;
        }
      }
    }
  }
  else if (pinf->p_strategy == MULTI_PART_PRICING)
  {
    for (j = 0; j < pinf->pmpinfo.bsize; j++)
    {
      if (ldbl_EGlpNumIsLess (d_max, pinf->pmpinfo.infeas[j]))
      {
        ldbl_EGlpNumCopy (d_max, pinf->pmpinfo.infeas[j]);
        pr->eindex = pinf->pmpinfo.bucket[j];
      }
    }
  }

  if (pr->eindex < 0)
    pr->price_stat = PRICE_OPTIMAL;
  else
  {
    if (phase == PRIMAL_PHASEI)
      ldbl_EGlpNumCopy (d_e, lp->pIdz[pr->eindex]);
    else
      ldbl_EGlpNumCopy (d_e, lp->dz[pr->eindex]);
    vs = lp->vstat[lp->nbaz[pr->eindex]];

    pr->price_stat = PRICE_NONOPTIMAL;
    if (vs == STAT_UPPER || (vs == STAT_ZERO && ldbl_EGlpNumIsLess (*ftol, d_e)))
      pr->dir = VDECREASE;
    else
      pr->dir = VINCREASE;
  }
  ldbl_EGlpNumClearVar (d_e);
  ldbl_EGlpNumClearVar (d_max);
}

static void ldbl_update_p_scaleinf (
  ldbl_price_info * const p,
  ldbl_heap * const h,
  int const i,
  long double inf,
  int const prule)
{
  if (!ldbl_EGlpNumIsNeqqZero (inf))
  {
    ldbl_EGlpNumZero (p->p_scaleinf[i]);
    if (h->hexist != 0 && h->loc[i] != -1)
      ldbl_ILLheap_delete (h, i);
  }
  else
  {
    if (prule == QS_PRICE_DDANTZIG)
      ldbl_EGlpNumCopy (p->p_scaleinf[i], inf);
    else if (prule == QS_PRICE_DSTEEP)
      ldbl_EGlpNumCopySqrOver (p->p_scaleinf[i], inf, p->dsinfo.norms[i]);
    else if (prule == QS_PRICE_DDEVEX)
      ldbl_EGlpNumCopySqrOver (p->p_scaleinf[i], inf, p->ddinfo.norms[i]);

    if (h->hexist != 0)
    {
      if (h->loc[i] == -1)
        ldbl_ILLheap_insert (h, i);
      else
        ldbl_ILLheap_modify (h, i);
    }
  }
}

static void ldbl_compute_primalI_inf (
  ldbl_lpinfo * const lp,
  int const i,
  long double * const inf)
{
  int const col = lp->baz[i];
  long double*x = &(lp->xbz[i]);
  long double*l = &(lp->lz[col]);
  long double*u = &(lp->uz[col]);
  long double*ftol = &(lp->tol->ip_tol);
  ldbl_EGlpNumZero (*inf);

  if (ldbl_EGlpNumIsLess (*ftol, *x) && ldbl_EGlpNumIsNeqq (*u, ldbl_INFTY))
    ldbl_EGlpNumCopy (*inf, *x);
  else if (ldbl_EGlpNumIsNeqq (*l, ldbl_NINFTY) && ldbl_EGlpNumIsSumLess (*x, *ftol,ldbl_zeroLpNum))
    ldbl_EGlpNumCopy (*inf, *x);
}

static void ldbl_compute_primalII_inf (
  ldbl_lpinfo * const lp,
  int const i,
  long double * const inf)
{
  int const col = lp->baz[i];
  long double*x = &(lp->xbz[i]);
  long double*l = &(lp->lz[col]);
  long double*u = &(lp->uz[col]);
  long double*ftol = &(lp->tol->pfeas_tol);
  ldbl_EGlpNumZero (*inf);

  if (ldbl_EGlpNumIsNeqq (*u, ldbl_INFTY) && ldbl_EGlpNumIsSumLess (*u, *ftol, *x))
    ldbl_EGlpNumCopyDiff (*inf, *x, *u);
  else if (ldbl_EGlpNumIsNeqq (*l, ldbl_NINFTY) && ldbl_EGlpNumIsSumLess (*x, *ftol, *l))
    ldbl_EGlpNumCopyDiff (*inf, *l, *x);
}

void ldbl_ILLprice_compute_primal_inf (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const p,
  int *const ix,
  int const icnt,
  int const phase)
{
  int i;
  int price;
  long double inf;
  ldbl_heap *h = &(p->h);

  price = (phase == DUAL_PHASEI) ? p->dI_price : p->dII_price;
  ldbl_EGlpNumInitVar (inf);
  ldbl_EGlpNumZero (inf);

  if (phase == DUAL_PHASEI)
  {
    if (ix == NULL)
      for (i = 0; i < lp->nrows; i++)
      {
        ldbl_compute_primalI_inf (lp, i, &inf);
        ldbl_update_p_scaleinf (p, h, i, inf, price);
      }
    else
      for (i = 0; i < icnt; i++)
      {
        ldbl_compute_primalI_inf (lp, ix[i], &inf);
        ldbl_update_p_scaleinf (p, h, ix[i], inf, price);
      }
  }
  else if (phase == DUAL_PHASEII)
  {
    if (ix == NULL)
      for (i = 0; i < lp->nrows; i++)
      {
        ldbl_compute_primalII_inf (lp, i, &inf);
        ldbl_update_p_scaleinf (p, h, i, inf, price);
      }
    else
      for (i = 0; i < icnt; i++)
      {
        ldbl_compute_primalII_inf (lp, ix[i], &inf);
        ldbl_update_p_scaleinf (p, h, ix[i], inf, price);
      }
  }
  ldbl_EGlpNumClearVar (inf);
}

void ldbl_ILLprice_dual (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  int const phase,
  ldbl_price_res * const pr)
{
  int i;
  long double p_max;
  long double ubound;
  long double*ftol = &(lp->tol->pfeas_tol);
  ldbl_heap *const h = &(pinf->h);

  ldbl_EGlpNumInitVar (p_max);
  ldbl_EGlpNumInitVar (ubound);
  pr->lindex = -1;
  ldbl_EGlpNumZero(p_max);

#if USEHEAP > 0
  ldbl_ILLprice_test_for_heap (lp, pinf, lp->nrows, pinf->p_scaleinf, DUAL_SIMPLEX,
                          1);
#endif

  if (pinf->d_strategy == COMPLETE_PRICING)
  {
    if (h->hexist)
    {
      pr->lindex = ldbl_ILLheap_findmin (h);
      if (pr->lindex != -1)
        ldbl_ILLheap_delete (h, pr->lindex);
    }
    else
    {
      for (i = 0; i < lp->nrows; i++)
      {
        if (ldbl_EGlpNumIsLess (p_max, pinf->p_scaleinf[i]))
        {
          ldbl_EGlpNumCopy (p_max, pinf->p_scaleinf[i]);
          pr->lindex = i;
        }
      }
    }
  }
  else if (pinf->d_strategy == MULTI_PART_PRICING)
  {
    for (i = 0; i < pinf->dmpinfo.bsize; i++)
    {
      if (ldbl_EGlpNumIsLess (p_max, pinf->dmpinfo.infeas[i]))
      {
        ldbl_EGlpNumCopy (p_max, pinf->dmpinfo.infeas[i]);
        pr->lindex = pinf->dmpinfo.bucket[i];
      }
    }
  }

  if (pr->lindex < 0)
    pr->price_stat = PRICE_OPTIMAL;
  else
  {
    pr->price_stat = NONOPTIMAL;

    if (ldbl_EGlpNumIsNeqq (lp->uz[lp->baz[pr->lindex]], ldbl_INFTY))
    {
      if (phase == DUAL_PHASEI)
        ldbl_EGlpNumZero(ubound);
      else
        ldbl_EGlpNumCopy(ubound,lp->uz[lp->baz[pr->lindex]]);
      if (ldbl_EGlpNumIsSumLess (*ftol, ubound, lp->xbz[pr->lindex]))
        pr->lvstat = STAT_UPPER;
      else
        pr->lvstat = STAT_LOWER;
    }
    else
      pr->lvstat = STAT_LOWER;
  }
  ldbl_EGlpNumClearVar (p_max);
  ldbl_EGlpNumClearVar (ubound);
}

int ldbl_ILLprice_get_rownorms (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  long double * const rnorms)
{
  int rval = 0;
  int i;

  if (pinf->dsinfo.norms == NULL)
  {
    rval = ldbl_ILLprice_build_dsteep_norms (lp, &(pinf->dsinfo));
    CHECKRVALG(rval,CLEANUP);
  }
  for (i = 0; i < lp->nrows; i++)
    ldbl_EGlpNumCopy (rnorms[i], pinf->dsinfo.norms[i]);

CLEANUP:
  if (rval)
    ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);

  return rval;
}

int ldbl_ILLprice_get_colnorms (
  ldbl_lpinfo * const lp,
  ldbl_price_info * const pinf,
  long double * const cnorms)
{
  int rval = 0;
  int i, j;

  if (pinf->psinfo.norms == NULL)
  {
    rval = ldbl_ILLprice_build_psteep_norms (lp, &(pinf->psinfo));
    CHECKRVALG(rval,CLEANUP);
  }
  for (i = 0; i < lp->nrows; i++)
    ldbl_EGlpNumZero (cnorms[lp->baz[i]]);
  for (j = 0; j < lp->nnbasic; j++)
    ldbl_EGlpNumCopy (cnorms[lp->nbaz[j]], pinf->psinfo.norms[j]);

CLEANUP:
  if (rval)
    ldbl_EGlpNumFreeArray (pinf->psinfo.norms);

  return rval;
}

int ldbl_ILLprice_get_newnorms (
  ldbl_lpinfo * const lp,
  int const nelems,
  long double * const norms,
  int *const matcnt,
  int *const matbeg,
  int *const matind,
  long double * const matval,
  int const option)
{
  int i, j;
  int rval = 0;
  ldbl_svector a;
  ldbl_svector y;

  ldbl_ILLsvector_init (&y);
  rval = ldbl_ILLsvector_alloc (&y, lp->nrows);
  CHECKRVALG(rval,CLEANUP);

  for (j = 0; j < nelems; j++)
  {
    a.nzcnt = matcnt[j];
    a.indx = &(matind[matbeg[j]]);
    a.coef = &(matval[matbeg[j]]);

    if (option == COLUMN_SOLVE)
      ldbl_ILLbasis_column_solve (lp, &a, &y);
    else
      ldbl_ILLbasis_row_solve (lp, &a, &y);

    ldbl_EGlpNumOne (norms[j]);
    for (i = 0; i < y.nzcnt; i++)
      ldbl_EGlpNumAddInnProdTo (norms[j], y.coef[i], y.coef[i]);
  }

CLEANUP:
  ldbl_ILLsvector_free (&y);
  EG_RETURN(rval);
}

int ldbl_ILLprice_get_new_rownorms (
  ldbl_lpinfo * const lp,
  int const newrows,
  long double * const rnorms,
  int *const rmatcnt,
  int *const rmatbeg,
  int *const rmatind,
  long double * const rmatval)
{
  return ldbl_ILLprice_get_newnorms (lp, newrows, rnorms, rmatcnt, rmatbeg, rmatind,
                                rmatval, ROW_SOLVE);
}

int ldbl_ILLprice_get_new_colnorms (
  ldbl_lpinfo * const lp,
  int const newrows,
  long double * const rnorms,
  int *const matcnt,
  int *const matbeg,
  int *const matind,
  long double * const matval)
{
  return ldbl_ILLprice_get_newnorms (lp, newrows, rnorms, matcnt, matbeg, matind,
                                matval, COLUMN_SOLVE);
}

int ldbl_ILLprice_load_rownorms (
  ldbl_lpinfo * const lp,
  long double * const rnorms,
  ldbl_price_info * const pinf)
{
  int i;
  int rval = 0;

  ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);
  pinf->dsinfo.norms = ldbl_EGlpNumAllocArray (lp->nrows);

  for (i = 0; i < lp->nrows; i++)
  {
    ldbl_EGlpNumCopy (pinf->dsinfo.norms[i], rnorms[i]);
    if (ldbl_EGlpNumIsLess (pinf->dsinfo.norms[i], ldbl_PARAM_MIN_DNORM))
      ldbl_EGlpNumCopy (pinf->dsinfo.norms[i], ldbl_PARAM_MIN_DNORM);
  }

  EG_RETURN(rval);
}

int ldbl_ILLprice_load_colnorms (
  ldbl_lpinfo * const lp,
  long double * const cnorms,
  ldbl_price_info * const pinf)
{
  int j;
  int rval = 0;

  ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
  pinf->psinfo.norms = ldbl_EGlpNumAllocArray (lp->nnbasic);

  for (j = 0; j < lp->nnbasic; j++)
  {
    ldbl_EGlpNumCopy (pinf->psinfo.norms[j], cnorms[lp->nbaz[j]]);
    if (ldbl_EGlpNumIsLess (pinf->psinfo.norms[j], ldbl_oneLpNum))
      ldbl_EGlpNumOne (pinf->psinfo.norms[j]);
  }

  EG_RETURN(rval);
}

#if ldbl_PRICE_DEBUG > 0
void ldbl_test_dsteep_norms (
  ldbl_lpinfo * lp,
  ldbl_price_info * p)
{
  int i, errn = 0;
  long double *pn = ldbl_EGlpNumAllocArray(lp->nrows);
  long double err, diff;
  ldbl_EGlpNumZero (err);

  ldbl_EGlpNumInitVar (err);
  ldbl_EGlpNumInitVar (diff);

  ldbl_ILLprice_get_dsteep_norms (lp, lp->yjz.nzcnt, lp->yjz.indx, pn);
  for (i = 0; i < lp->yjz.nzcnt; i++)
  {
    ldbl_EGlpNumCopyDiff (diff, pn[i], p->dsinfo.norms[lp->yjz.indx[i]]);
    ldbl_EGlpNumCopyAbs(diff,diff);
    if (ldbl_EGlpNumIsLess (ldbl_PFEAS_TOLER, diff))
    {
      errn++;
      ldbl_EGlpNumAddTo (err, diff);
      ldbl_EGlpNumCopy (p->dsinfo.norms[lp->yjz.indx[i]], pn[i]);
    }
  }
  if (errn)
    printf ("%d: dnorm errn = %d, err = %.6f\n", lp->cnts->tot_iter, errn,
            ldbl_EGlpNumToLf (err));
  ldbl_EGlpNumFreeArray (pn);
  ldbl_EGlpNumClearVar (diff);
  ldbl_EGlpNumClearVar (err);
}
#endif

---