ldbl_price.c

Go to the documentation of this file.
00001 /****************************************************************************/
00002 /*                                                                          */
00003 /*  This file is part of QSopt_ex.                                          */
00004 /*                                                                          */
00005 /*  (c) Copyright 2006 by David Applegate, William Cook, Sanjeeb Dash,      */
00006 /*  and Daniel Espinoza                                                     */
00007 /*                                                                          */
00008 /*  Sanjeeb Dash ownership of copyright in QSopt_ex is derived from his     */
00009 /*  copyright in QSopt.                                                     */
00010 /*                                                                          */
00011 /*  This code may be used under the terms of the GNU General Public License */
00012 /*  (Version 2.1 or later) as published by the Free Software Foundation.    */
00013 /*                                                                          */
00014 /*  Alternatively, use is granted for research purposes only.               */
00015 /*                                                                          */
00016 /*  It is your choice of which of these two licenses you are operating      */
00017 /*  under.                                                                  */
00018 /*                                                                          */
00019 /*  We make no guarantees about the correctness or usefulness of this code. */
00020 /*                                                                          */
00021 /****************************************************************************/
00022 
00023 /* RCS_INFO = "$RCSfile: ldbl_price.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 //static int TRACE = 0;
00025 #include "qs_config.h"
00026 #include "stddefs.h"
00027 #include "ldbl_qsopt.h"
00028 #include "ldbl_lpdefs.h"
00029 #include "ldbl_fct.h"
00030 #include "ldbl_price.h"
00031 #include "ldbl_basis.h"
00032 #include "ldbl_iqsutil.h"
00033 #include "ldbl_dstruct.h"
00034 #ifdef USEDMALLOC
00035 #include "dmalloc.h"
00036 #endif
00037 
00038 #define  ldbl_MULTIP 1
00039 #define  ldbl_PRICE_DEBUG 0
00040 
00041 static void ldbl_update_d_scaleinf (
00042   ldbl_price_info * const p,
00043   ldbl_heap * const h,
00044   int const j,
00045   long double inf,
00046   int const prule),
00047   ldbl_update_p_scaleinf (
00048   ldbl_price_info * const p,
00049   ldbl_heap * const h,
00050   int const i,
00051   long double inf,
00052   int const prule);
00053 
00054 static void ldbl_compute_dualI_inf (
00055   ldbl_lpinfo * const lp,
00056   int const j,
00057   long double * const inf),
00058   ldbl_compute_dualII_inf (
00059   ldbl_lpinfo * const lp,
00060   int const j,
00061   long double * const inf),
00062   ldbl_compute_primalI_inf (
00063   ldbl_lpinfo * const lp,
00064   int const i,
00065   long double * const inf),
00066   ldbl_compute_primalII_inf (
00067   ldbl_lpinfo * const lp,
00068   int const i,
00069   long double * const inf);
00070 
00071 void ldbl_ILLprice_free_heap (
00072   ldbl_price_info * const pinf)
00073 {
00074   ldbl_ILLheap_free (&(pinf->h));
00075 }
00076 
00077 int ldbl_ILLprice_build_heap (
00078   ldbl_price_info * const pinf,
00079   int const nkeys,
00080   long double * keylist)
00081 {
00082   ldbl_ILLheap_init (&(pinf->h));
00083   ldbl_EGlpNumSet (pinf->htrigger,
00084               1.0 +
00085               (double) nkeys / (PARAM_HEAP_RATIO * ldbl_ILLutil_our_log2 (nkeys)));
00086   return ldbl_ILLheap_build (&(pinf->h), nkeys, keylist);
00087 }
00088 
00089 int ldbl_ILLprice_test_for_heap (
00090   ldbl_lpinfo * const lp,
00091   ldbl_price_info * const pinf,
00092   int const nkeys,
00093   long double * keylist,
00094   int const algo,
00095   int const upd)
00096 {
00097   ldbl_heap *const h = &(pinf->h);
00098   int rval = 0;
00099   long double ravg;
00100 
00101   if (upd != 0)
00102   {
00103     ldbl_EGlpNumInitVar (ravg);
00104     if (algo == PRIMAL_SIMPLEX)
00105       ldbl_EGlpNumCopy (ravg, lp->cnts->za_ravg);
00106     else
00107       ldbl_EGlpNumCopy (ravg, lp->cnts->y_ravg);
00108     if (ldbl_EGlpNumIsLeq (ravg, pinf->htrigger))
00109       pinf->hineff--;
00110     else
00111     {
00112       ldbl_EGlpNumDivUiTo (ravg, 2U);
00113       if (ldbl_EGlpNumIsLess (pinf->htrigger, ravg))
00114         pinf->hineff++;
00115     }
00116     ldbl_EGlpNumClearVar (ravg);
00117   }
00118   if (h->hexist == 0 && pinf->hineff <= 0)
00119   {
00120     rval = ldbl_ILLprice_build_heap (pinf, nkeys, keylist);
00121     CHECKRVALG (rval, CLEANUP);
00122   }
00123   else if (h->hexist != 0 && pinf->hineff >= PARAM_HEAP_UTRIGGER)
00124   {
00125     ldbl_ILLprice_free_heap (pinf);
00126     /*
00127      * printf ("freeing ldbl_heap ..\n");
00128      * printf ("iter = %d, ravg = %.2f, trigger = %.2f\n",
00129      * lp->cnts->tot_iter, ravg, pinf->htrigger);
00130      */
00131   }
00132 
00133 CLEANUP:
00134   if (rval)
00135     ldbl_ILLprice_free_heap (pinf);
00136   return rval;
00137 }
00138 
00139 void ldbl_ILLprice_init_pricing_info (
00140   ldbl_price_info * const pinf)
00141 {
00142   pinf->p_strategy = -1;
00143   pinf->d_strategy = -1;
00144   pinf->pI_price = -1;
00145   pinf->pII_price = -1;
00146   pinf->dI_price = -1;
00147   pinf->dII_price = -1;
00148   pinf->cur_price = -1;
00149   pinf->p_scaleinf = 0;
00150   pinf->d_scaleinf = 0;
00151   pinf->pdinfo.norms = 0;
00152   pinf->pdinfo.refframe = 0;
00153   pinf->psinfo.norms = 0;
00154   pinf->ddinfo.norms = 0;
00155   pinf->ddinfo.refframe = 0;
00156   pinf->dsinfo.norms = 0;
00157   pinf->dmpinfo.gstart = pinf->pmpinfo.gstart = 0;
00158   pinf->dmpinfo.gshift = pinf->pmpinfo.gshift = 0;
00159   pinf->dmpinfo.gsize = pinf->pmpinfo.gsize = 0;
00160   pinf->dmpinfo.bucket = pinf->pmpinfo.bucket = 0;
00161   pinf->dmpinfo.perm = pinf->pmpinfo.perm = 0;
00162   pinf->dmpinfo.infeas = pinf->pmpinfo.infeas = 0;
00163   ldbl_ILLheap_init (&(pinf->h));
00164   ldbl_EGlpNumZero (pinf->htrigger);
00165   pinf->hineff = 0;
00166 }
00167 
00168 void ldbl_ILLprice_free_pricing_info (
00169   ldbl_price_info * const pinf)
00170 {
00171   ldbl_EGlpNumFreeArray (pinf->p_scaleinf);
00172   ldbl_EGlpNumFreeArray (pinf->d_scaleinf);
00173   ldbl_EGlpNumFreeArray (pinf->pdinfo.norms);
00174   ILL_IFFREE (pinf->pdinfo.refframe, int);
00175   ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
00176   ldbl_EGlpNumFreeArray (pinf->ddinfo.norms);
00177   ILL_IFFREE (pinf->ddinfo.refframe, int);
00178   ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);
00179 
00180   ldbl_ILLprice_free_mpartial_info (&(pinf->pmpinfo));
00181   ldbl_ILLprice_free_mpartial_info (&(pinf->dmpinfo));
00182   ldbl_ILLprice_free_heap (pinf);
00183 }
00184 
00185 int ldbl_ILLprice_build_pricing_info (
00186   ldbl_lpinfo * const lp,
00187   ldbl_price_info * const pinf,
00188   int const phase)
00189 {
00190   int rval = 0;
00191   int p_price = -1;
00192   int d_price = -1;
00193 
00194   switch (phase)
00195   {
00196   case PRIMAL_PHASEI:
00197     p_price = pinf->pI_price;
00198     break;
00199   case PRIMAL_PHASEII:
00200     p_price = pinf->pII_price;
00201     break;
00202   case DUAL_PHASEI:
00203     d_price = pinf->dI_price;
00204     break;
00205   case DUAL_PHASEII:
00206     d_price = pinf->dII_price;
00207     break;
00208   }
00209 
00210   if (p_price != -1)
00211   {
00212     pinf->cur_price = p_price;
00213 
00214     if (p_price == QS_PRICE_PDANTZIG || p_price == QS_PRICE_PDEVEX ||
00215         p_price == QS_PRICE_PSTEEP)
00216     {
00217       pinf->p_strategy = COMPLETE_PRICING;
00218       ldbl_EGlpNumFreeArray (pinf->d_scaleinf);
00219       pinf->d_scaleinf = ldbl_EGlpNumAllocArray (lp->nnbasic);
00220     }
00221     else if (p_price == QS_PRICE_PMULTPARTIAL)
00222       pinf->p_strategy = MULTI_PART_PRICING;
00223 
00224     switch (p_price)
00225     {
00226     case QS_PRICE_PDEVEX:
00227       if (pinf->pdinfo.norms)
00228         return rval;
00229       rval = ldbl_ILLprice_build_pdevex_norms (lp, &(pinf->pdinfo), 0);
00230       CHECKRVALG(rval,CLEANUP);
00231       break;
00232     case QS_PRICE_PSTEEP:
00233       if (pinf->psinfo.norms)
00234         return rval;
00235       rval = ldbl_ILLprice_build_psteep_norms (lp, &(pinf->psinfo));
00236       CHECKRVALG(rval,CLEANUP);
00237       break;
00238     case QS_PRICE_PMULTPARTIAL:
00239       rval = ldbl_ILLprice_build_mpartial_info (lp, pinf, COL_PRICING);
00240       CHECKRVALG(rval,CLEANUP);
00241       break;
00242     }
00243   }
00244   else if (d_price != -1)
00245   {
00246     pinf->cur_price = d_price;
00247 
00248     if (d_price == QS_PRICE_DDANTZIG || d_price == QS_PRICE_DSTEEP ||
00249         d_price == QS_PRICE_DDEVEX)
00250     {
00251       pinf->d_strategy = COMPLETE_PRICING;
00252       ldbl_EGlpNumFreeArray (pinf->p_scaleinf);
00253       pinf->p_scaleinf = ldbl_EGlpNumAllocArray (lp->nrows);
00254     }
00255     else if (d_price == QS_PRICE_DMULTPARTIAL)
00256       pinf->d_strategy = MULTI_PART_PRICING;
00257 
00258     switch (d_price)
00259     {
00260     case QS_PRICE_DSTEEP:
00261       if (pinf->dsinfo.norms)
00262         return rval;
00263       rval = ldbl_ILLprice_build_dsteep_norms (lp, &(pinf->dsinfo));
00264       CHECKRVALG(rval,CLEANUP);
00265       break;
00266     case QS_PRICE_DMULTPARTIAL:
00267       rval = ldbl_ILLprice_build_mpartial_info (lp, pinf, ROW_PRICING);
00268       CHECKRVALG(rval,CLEANUP);
00269       break;
00270     case QS_PRICE_DDEVEX:
00271       if (pinf->ddinfo.norms)
00272         return rval;
00273       rval = ldbl_ILLprice_build_ddevex_norms (lp, &(pinf->ddinfo), 0);
00274       CHECKRVALG(rval,CLEANUP);
00275       break;
00276     }
00277   }
00278 
00279 CLEANUP:
00280   if (rval)
00281     ldbl_ILLprice_free_pricing_info (pinf);
00282   EG_RETURN(rval);
00283 }
00284 
00285 int ldbl_ILLprice_update_pricing_info (
00286   ldbl_lpinfo * const lp,
00287   ldbl_price_info * const pinf,
00288   int const phase,
00289   ldbl_svector * const wz,
00290   int const eindex,
00291   int const lindex,
00292   long double y)
00293 {
00294   int rval = 0;
00295   int p_price = -1;
00296   int d_price = -1;
00297 
00298   switch (phase)
00299   {
00300   case PRIMAL_PHASEI:
00301     p_price = pinf->pI_price;
00302     break;
00303   case PRIMAL_PHASEII:
00304     p_price = pinf->pII_price;
00305     break;
00306   case DUAL_PHASEI:
00307     d_price = pinf->dI_price;
00308     break;
00309   case DUAL_PHASEII:
00310     d_price = pinf->dII_price;
00311     break;
00312   }
00313 
00314   if (p_price != -1)
00315   {
00316     if (p_price == QS_PRICE_PDEVEX)
00317     {
00318       rval = ldbl_ILLprice_update_pdevex_norms (lp, &(pinf->pdinfo), eindex, y);
00319       CHECKRVALG(rval,CLEANUP);
00320     }
00321     else if (p_price == QS_PRICE_PSTEEP)
00322       ldbl_ILLprice_update_psteep_norms (lp, &(pinf->psinfo), wz, eindex, y);
00323   }
00324   else if (d_price != -1)
00325   {
00326     if (d_price == QS_PRICE_DSTEEP)
00327       ldbl_ILLprice_update_dsteep_norms (lp, &(pinf->dsinfo), wz, lindex, y);
00328     else if (d_price == QS_PRICE_DDEVEX)
00329     {
00330       rval = ldbl_ILLprice_update_ddevex_norms (lp, &(pinf->ddinfo), lindex, y);
00331       CHECKRVALG(rval,CLEANUP);
00332     }
00333   }
00334 CLEANUP:
00335   EG_RETURN(rval);
00336 }
00337 
00338 int ldbl_ILLprice_get_price (
00339   ldbl_price_info * const p,
00340   int const phase)
00341 {
00342   int pri = -1;
00343 
00344   switch (phase)
00345   {
00346   case PRIMAL_PHASEI:
00347     return p->pI_price;
00348   case PRIMAL_PHASEII:
00349     return p->pII_price;
00350   case DUAL_PHASEI:
00351     return p->dI_price;
00352   case DUAL_PHASEII:
00353     return p->dII_price;
00354   }
00355   return pri;
00356 }
00357 
00358 void ldbl_ILLprice_free_mpartial_info (
00359   ldbl_mpart_info * p)
00360 {
00361   ILL_IFFREE (p->gstart, int);
00362   ILL_IFFREE (p->gshift, int);
00363   ILL_IFFREE (p->gsize, int);
00364   ILL_IFFREE (p->bucket, int);
00365   ldbl_EGlpNumFreeArray (p->infeas);
00366   ILL_IFFREE (p->perm, int);
00367 }
00368 
00369 int ldbl_ILLprice_build_mpartial_info (
00370   ldbl_lpinfo * const lp,
00371   ldbl_price_info * const pinf,
00372   int const pricetype)
00373 {
00374   int i = 0;
00375   int rval = 0;
00376   int extra = 0;
00377   int nelems;
00378   ldbl_mpart_info *p;
00379 
00380   p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
00381   p->k = 50;
00382   p->cgroup = 0;
00383   nelems = (pricetype == COL_PRICING) ? lp->nnbasic : lp->nrows;
00384 
00385   if (nelems % p->k)
00386     extra = nelems - p->k * (nelems / p->k);
00387   p->ngroups = nelems / p->k;
00388   if (extra != 0)
00389     p->ngroups++;
00390 
00391   ILL_SAFE_MALLOC (p->gstart, p->ngroups, int);
00392   ILL_SAFE_MALLOC (p->gshift, p->ngroups, int);
00393   ILL_SAFE_MALLOC (p->gsize, p->ngroups, int);
00394   ILL_SAFE_MALLOC (p->bucket, 2 * p->k, int);
00395   p->infeas = ldbl_EGlpNumAllocArray (2 * p->k);
00396   ILL_SAFE_MALLOC (p->perm, 2 * p->k, int);
00397 
00398   p->bsize = 0;
00399 
00400   if (extra != 0)
00401   {
00402     p->gstart[0] = 0;
00403     p->gshift[0] = 1;
00404     p->gsize[0] = extra;
00405     for (i = 1; i < p->ngroups; i++)
00406     {
00407       p->gstart[i] = extra + i - 1;
00408       p->gshift[i] = p->ngroups - 1;
00409       p->gsize[i] = p->k;
00410     }
00411   }
00412   else
00413   {
00414     for (i = 0; i < p->ngroups; i++)
00415     {
00416       p->gstart[i] = i;
00417       p->gshift[i] = p->ngroups;
00418       p->gsize[i] = p->k;
00419     }
00420   }
00421 
00422 CLEANUP:
00423   if (rval)
00424     ldbl_ILLprice_free_mpartial_info (p);
00425   EG_RETURN(rval);
00426 }
00427 
00428 void ldbl_ILLprice_init_mpartial_price (
00429   ldbl_lpinfo * const lp,
00430   ldbl_price_info * const pinf,
00431   int const phase,
00432   int const pricetype)
00433 {
00434   int i;
00435   ldbl_mpart_info *p;
00436 
00437   p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
00438   p->bsize = 0;
00439   i = p->cgroup;
00440   do
00441   {
00442     ldbl_ILLprice_mpartial_group (lp, p, phase, i, pricetype);
00443     i = (i + 1) % p->ngroups;
00444   } while (i != p->cgroup && p->bsize <= p->k);
00445   p->cgroup = i;
00446 }
00447 
00448 void ldbl_ILLprice_update_mpartial_price (
00449   ldbl_lpinfo * const lp,
00450   ldbl_price_info * const pinf,
00451   int const phase,
00452   int const pricetype)
00453 {
00454   int i = 0;
00455   int csize = 0;
00456   long double infeas;
00457   ldbl_mpart_info *p;
00458   ldbl_price_res pr;
00459 
00460   p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
00461   ldbl_EGlpNumInitVar (pr.dinfeas);
00462   ldbl_EGlpNumInitVar (pr.pinfeas);
00463   ldbl_EGlpNumInitVar (infeas);
00464 
00465 #ifdef ldbl_MULTIP
00466   i = 0;
00467   while (i < p->bsize)
00468   {
00469     if (pricetype == COL_PRICING)
00470     {
00471       ldbl_ILLprice_column (lp, p->bucket[i], phase, &pr);
00472       ldbl_EGlpNumCopy (infeas, pr.dinfeas);
00473     }
00474     else
00475     {
00476       ldbl_ILLprice_row (lp, p->bucket[i], phase, &pr);
00477       ldbl_EGlpNumCopy (infeas, pr.pinfeas);
00478     }
00479     if (!ldbl_EGlpNumIsNeqqZero (infeas))
00480     {
00481       p->bucket[i] = p->bucket[p->bsize - 1];
00482       p->bsize--;
00483     }
00484     else
00485     {
00486       ldbl_EGlpNumCopy (p->infeas[i], infeas);
00487       i++;
00488     }
00489   }
00490   if (p->bsize > 0)
00491   {
00492     for (i = 0; i < p->bsize; i++)
00493       p->perm[i] = i;
00494     ldbl_EGutilPermSort ((size_t) (p->bsize), p->perm,
00495                     (const long double * const) p->infeas);
00496 
00497     csize = QSMIN (p->bsize, p->k);
00498     for (i = csize - 1; i >= 0; i--)
00499       lp->iwork[p->bucket[p->perm[i]]] = 1;
00500 
00501     for (i = 0, csize = 0; i < p->bsize; i++)
00502       if (lp->iwork[p->bucket[i]] == 1)
00503       {
00504         ldbl_EGlpNumCopy (p->infeas[csize], p->infeas[i]);
00505         p->bucket[csize] = p->bucket[i];
00506         csize++;
00507       }
00508     p->bsize = csize;
00509   }
00510 #else
00511   p->bsize = 0;
00512 #endif
00513 
00514   i = p->cgroup;
00515   do
00516   {
00517     ldbl_ILLprice_mpartial_group (lp, p, phase, i, pricetype);
00518     i = (i + 1) % p->ngroups;
00519   } while (i != p->cgroup && p->bsize <= p->k);
00520   p->cgroup = i;
00521 
00522 #ifdef ldbl_MULTIP
00523   for (i = 0; i < csize; i++)
00524     lp->iwork[p->bucket[i]] = 0;
00525 #endif
00526   ldbl_EGlpNumClearVar (infeas);
00527   ldbl_EGlpNumClearVar (pr.pinfeas);
00528   ldbl_EGlpNumClearVar (pr.dinfeas);
00529 }
00530 
00531 void ldbl_ILLprice_delete_onempart_price (
00532   /*ldbl_lpinfo * const lp,*/
00533   ldbl_price_info * const pinf,
00534   int const indx,
00535   int const pricetype)
00536 {
00537   int i = 0;
00538   ldbl_mpart_info *p;
00539 
00540   p = (pricetype == COL_PRICING) ? &(pinf->pmpinfo) : &(pinf->dmpinfo);
00541 
00542   for (i = 0; i < p->bsize; i++)
00543     if (p->bucket[i] == indx)
00544     {
00545       p->bucket[i] = p->bucket[p->bsize - 1];
00546       ldbl_EGlpNumCopy (p->infeas[i], p->infeas[p->bsize - 1]);
00547       p->bsize--;
00548       break;
00549     }
00550 }
00551 
00552 void ldbl_ILLprice_mpartial_group (
00553   ldbl_lpinfo * const lp,
00554   ldbl_mpart_info * const p,
00555   int const phase,
00556   int const g,
00557   int const pricetype)
00558 {
00559   int i, ix;
00560   int gstart = p->gstart[g];
00561   int gsize = p->gsize[g];
00562   int gshift = p->gshift[g];
00563   long double infeas;
00564   ldbl_price_res pr;
00565 
00566   ldbl_EGlpNumInitVar (pr.dinfeas);
00567   ldbl_EGlpNumInitVar (pr.pinfeas);
00568   ldbl_EGlpNumInitVar (infeas);
00569 
00570   for (i = 0, ix = gstart; i < gsize; i++, ix += gshift)
00571   {
00572 #ifdef ldbl_MULTIP
00573     if (lp->iwork[ix])
00574       continue;
00575 #endif
00576     if (pricetype == COL_PRICING)
00577     {
00578       ldbl_ILLprice_column (lp, ix, phase, &pr);
00579       ldbl_EGlpNumCopy (infeas, pr.dinfeas);
00580     }
00581     else
00582     {
00583       ldbl_ILLprice_row (lp, ix, phase, &pr);
00584       ldbl_EGlpNumCopy (infeas, pr.pinfeas);
00585     }
00586     if (ldbl_EGlpNumIsNeqqZero (infeas))
00587     {
00588       ldbl_EGlpNumCopy (p->infeas[p->bsize], infeas);
00589       p->bucket[p->bsize] = ix;
00590       p->bsize++;
00591     }
00592   }
00593   ldbl_EGlpNumClearVar (infeas);
00594   ldbl_EGlpNumClearVar (pr.dinfeas);
00595   ldbl_EGlpNumClearVar (pr.pinfeas);
00596 }
00597 
00598 void ldbl_ILLprice_column (
00599   ldbl_lpinfo * const lp,
00600   int const ix,
00601   int const phase,
00602   ldbl_price_res * const pr)
00603 {
00604   int i;
00605   int col;
00606   int mcnt;
00607   int mbeg;
00608   long double sum;
00609 
00610   ldbl_EGlpNumZero (pr->dinfeas);
00611   col = lp->nbaz[ix];
00612   if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFIXED)
00613     return;
00614   ldbl_EGlpNumInitVar (sum);
00615   ldbl_EGlpNumZero (sum);
00616   mcnt = lp->matcnt[col];
00617   mbeg = lp->matbeg[col];
00618 
00619   if (phase == PRIMAL_PHASEII)
00620   {
00621     for (i = 0; i < mcnt; i++)
00622       ldbl_EGlpNumAddInnProdTo (sum, lp->piz[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00623     ldbl_EGlpNumCopyDiff (lp->dz[ix], lp->cz[col], sum);
00624     ldbl_compute_dualII_inf (lp, ix, &(pr->dinfeas));
00625   }
00626   else
00627   {
00628     for (i = 0; i < mcnt; i++)
00629       ldbl_EGlpNumAddInnProdTo (sum, lp->pIpiz[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00630     ldbl_EGlpNumCopyNeg (lp->pIdz[ix], sum);
00631     ldbl_compute_dualI_inf (lp, ix, &(pr->dinfeas));
00632   }
00633   ldbl_EGlpNumClearVar (sum);
00634 }
00635 
00636 void ldbl_ILLprice_row (
00637   ldbl_lpinfo * const lp,
00638   int const ix,
00639   int const phase,
00640   ldbl_price_res * const pr)
00641 {
00642   if (phase == DUAL_PHASEII)
00643     ldbl_compute_primalII_inf (lp, ix, &(pr->pinfeas));
00644   else
00645     ldbl_compute_primalI_inf (lp, ix, &(pr->pinfeas));
00646 }
00647 
00648 int ldbl_ILLprice_build_pdevex_norms (
00649   ldbl_lpinfo * const lp,
00650   ldbl_p_devex_info * const pdinfo,
00651   int const reinit)
00652 {
00653   int j;
00654   int rval = 0;
00655 
00656   if (reinit == 0)
00657   {
00658     pdinfo->ninit = 0;
00659     pdinfo->norms = ldbl_EGlpNumAllocArray (lp->nnbasic);
00660     ILL_SAFE_MALLOC (pdinfo->refframe, lp->ncols, int);
00661   }
00662 
00663   if (reinit != 0)
00664     pdinfo->ninit++;
00665 
00666   for (j = 0; j < lp->ncols; j++)
00667   {
00668     if (lp->vstat[j] == STAT_BASIC)
00669       pdinfo->refframe[j] = 0;
00670     else
00671     {
00672       ldbl_EGlpNumOne (pdinfo->norms[lp->vindex[j]]);
00673       pdinfo->refframe[j] = 1;
00674     }
00675   }
00676 
00677 CLEANUP:
00678   if (rval)
00679   {
00680     ldbl_EGlpNumFreeArray (pdinfo->norms);
00681     ILL_IFFREE (pdinfo->refframe, int);
00682   }
00683   EG_RETURN(rval);
00684 }
00685 
00686 int ldbl_ILLprice_update_pdevex_norms (
00687   ldbl_lpinfo * const lp,
00688   ldbl_p_devex_info * const pdinfo,
00689   int const eindex,
00690   long double yl)
00691 {
00692   int i, j;
00693   long double normj;
00694   long double zAj;
00695   long double ntmp, ntmp2;
00696 
00697   ldbl_EGlpNumInitVar (normj);
00698   ldbl_EGlpNumInitVar (zAj);
00699   ldbl_EGlpNumInitVar (ntmp);
00700   ldbl_EGlpNumInitVar (ntmp2);
00701   ldbl_EGlpNumZero (normj);
00702 
00703   for (i = 0; i < lp->yjz.nzcnt; i++)
00704     if (pdinfo->refframe[lp->baz[lp->yjz.indx[i]]])
00705       ldbl_EGlpNumAddInnProdTo (normj, lp->yjz.coef[i], lp->yjz.coef[i]);
00706 
00707   if (pdinfo->refframe[lp->nbaz[eindex]])
00708     ldbl_EGlpNumAddTo (normj, ldbl_oneLpNum);
00709 
00710   ldbl_EGlpNumSet(ntmp,1000.0);
00711   ldbl_EGlpNumSet(ntmp2,0.001);
00712   ldbl_EGlpNumMultTo(ntmp,pdinfo->norms[eindex]);
00713   ldbl_EGlpNumMultTo(ntmp2,pdinfo->norms[eindex]);
00714   if (ldbl_EGlpNumIsLess (normj, ntmp2) || ldbl_EGlpNumIsLess (ntmp, normj))
00715   {
00716     ldbl_EGlpNumClearVar (zAj);
00717     ldbl_EGlpNumClearVar (normj);
00718     ldbl_EGlpNumClearVar (ntmp);
00719     ldbl_EGlpNumClearVar (ntmp2);
00720     return ldbl_ILLprice_build_pdevex_norms (lp, pdinfo, 1);
00721   }
00722 
00723   for (i = 0; i < lp->zA.nzcnt; i++)
00724   {
00725     j = lp->zA.indx[i];
00726     ldbl_EGlpNumCopyFrac (zAj, lp->zA.coef[i], yl);
00727     ldbl_EGlpNumMultTo (zAj, zAj);
00728     ldbl_EGlpNumMultTo (zAj, normj);
00729     if (ldbl_EGlpNumIsLess (pdinfo->norms[j], zAj))
00730       ldbl_EGlpNumCopy (pdinfo->norms[j], zAj);
00731   }
00732   ldbl_EGlpNumDivTo (normj, yl);
00733   ldbl_EGlpNumDivTo (normj, yl);
00734   if (ldbl_EGlpNumIsLess (normj, ldbl_oneLpNum))
00735     ldbl_EGlpNumCopy (pdinfo->norms[eindex], ldbl_oneLpNum);
00736   else
00737     ldbl_EGlpNumCopy (pdinfo->norms[eindex], normj);
00738   ldbl_EGlpNumClearVar (zAj);
00739   ldbl_EGlpNumClearVar (normj);
00740   ldbl_EGlpNumClearVar (ntmp);
00741   ldbl_EGlpNumClearVar (ntmp2);
00742   return 0;
00743 }
00744 
00745 int ldbl_ILLprice_build_psteep_norms (
00746   ldbl_lpinfo * const lp,
00747   ldbl_p_steep_info * const psinfo)
00748 {
00749   int j;
00750   int rval = 0;
00751   ldbl_svector yz;
00752 
00753   ldbl_ILLsvector_init (&yz);
00754   rval = ldbl_ILLsvector_alloc (&yz, lp->nrows);
00755   CHECKRVALG(rval,CLEANUP);
00756   psinfo->norms = ldbl_EGlpNumAllocArray (lp->nnbasic);
00757 
00758   for (j = 0; j < lp->nnbasic; j++)
00759   {
00760     rval = ILLstring_report (NULL, &lp->O->reporter);
00761     CHECKRVALG(rval,CLEANUP);
00762     ldbl_ILLfct_compute_yz (lp, &yz, 0, lp->nbaz[j]);
00763     ldbl_EGlpNumInnProd (psinfo->norms[j], yz.coef, yz.coef, (size_t) yz.nzcnt);
00764     ldbl_EGlpNumAddTo (psinfo->norms[j], ldbl_oneLpNum);
00765   }
00766 
00767 CLEANUP:
00768   ldbl_ILLsvector_free (&yz);
00769   if (rval)
00770     ldbl_EGlpNumFreeArray (psinfo->norms);
00771 
00772   EG_RETURN(rval);
00773 }
00774 
00775 void ldbl_ILLprice_update_psteep_norms (
00776   ldbl_lpinfo * const lp,
00777   ldbl_p_steep_info * const psinfo,
00778   ldbl_svector * const wz,
00779   int const eindex,
00780   long double yl)
00781 {
00782   int i, j, k;
00783   int mcnt, mbeg;
00784   long double normj,ntmp;
00785   long double zAj, wAj;
00786   long double *v = 0;
00787 
00788   ldbl_EGlpNumInitVar (normj);
00789   ldbl_EGlpNumInitVar (zAj);
00790   ldbl_EGlpNumInitVar (ntmp);
00791   ldbl_EGlpNumInitVar (wAj);
00792   ldbl_EGlpNumInnProd (normj, lp->yjz.coef, lp->yjz.coef, (size_t) (lp->yjz.nzcnt));
00793   ldbl_EGlpNumAddTo (normj, ldbl_oneLpNum);
00794 
00795 #if 0
00796   Bico - remove warnings for dist
00797     if (fabs ((normj - psinfo->norms[eindex]) / normj) > 1000.0 /* 0.01 */ )
00798     {
00799       printf ("warning: incorrect norm values\n");
00800       printf ("anorm = %.6f, pnorm = %.6f\n", normj, psinfo->norms[eindex]);
00801       fflush (stdout);
00802     }
00803 #endif
00804 
00805   ldbl_ILLfct_load_workvector (lp, wz);
00806   v = lp->work.coef;
00807 
00808   for (k = 0; k < lp->zA.nzcnt; k++)
00809   {
00810     j = lp->zA.indx[k];
00811     ldbl_EGlpNumCopy (zAj, lp->zA.coef[k]);
00812     ldbl_EGlpNumZero (wAj);
00813     mcnt = lp->matcnt[lp->nbaz[j]];
00814     mbeg = lp->matbeg[lp->nbaz[j]];
00815     for (i = 0; i < mcnt; i++)
00816       ldbl_EGlpNumAddInnProdTo (wAj, lp->matval[mbeg + i], v[lp->matind[mbeg + i]]);
00817 
00818     /* compute ntmp = (zAj * ((zAj * normj / yl) - (2.0 * wAj))) / yl; */ 
00819     ldbl_EGlpNumCopy(ntmp,zAj);
00820     ldbl_EGlpNumMultTo(ntmp,normj);
00821     ldbl_EGlpNumDivTo(ntmp,yl);
00822     ldbl_EGlpNumSubTo(ntmp,wAj);
00823     ldbl_EGlpNumSubTo(ntmp,wAj);
00824     ldbl_EGlpNumMultTo(ntmp,zAj);
00825     ldbl_EGlpNumDivTo(ntmp,yl);
00826     /* set psinfo->norms[j] += (zAj * ((zAj * normj / yl) - (2.0 * wAj))) / yl; */
00827     ldbl_EGlpNumAddTo(psinfo->norms[j],ntmp);
00828     if (ldbl_EGlpNumIsLess (psinfo->norms[j], ldbl_oneLpNum))
00829       ldbl_EGlpNumOne (psinfo->norms[j]);
00830   }
00831 
00832   ldbl_EGlpNumCopyFrac (psinfo->norms[eindex], normj, yl);
00833   ldbl_EGlpNumDivTo (psinfo->norms[eindex], yl);
00834   if (ldbl_EGlpNumIsLess (psinfo->norms[eindex], ldbl_oneLpNum))
00835     ldbl_EGlpNumOne (psinfo->norms[eindex]);
00836 
00837   ldbl_ILLfct_zero_workvector (lp);
00838   ldbl_EGlpNumClearVar (wAj);
00839   ldbl_EGlpNumClearVar (zAj);
00840   ldbl_EGlpNumClearVar (normj);
00841   ldbl_EGlpNumClearVar (ntmp);
00842 }
00843 
00844 int ldbl_ILLprice_build_ddevex_norms (
00845   ldbl_lpinfo * const lp,
00846   ldbl_d_devex_info * const ddinfo,
00847   int const reinit)
00848 {
00849   int i;
00850   int rval = 0;
00851 
00852   if (reinit == 0)
00853   {
00854     ddinfo->ninit = 0;
00855     ddinfo->norms = ldbl_EGlpNumAllocArray (lp->nrows);
00856     ILL_SAFE_MALLOC (ddinfo->refframe, lp->ncols, int);
00857   }
00858   if (reinit != 0)
00859     ddinfo->ninit++;
00860 
00861   for (i = 0; i < lp->ncols; i++)
00862     ddinfo->refframe[i] = (lp->vstat[i] == STAT_BASIC) ? 1 : 0;
00863 
00864   for (i = 0; i < lp->nrows; i++)
00865     ldbl_EGlpNumOne (ddinfo->norms[i]);
00866 
00867 CLEANUP:
00868   if (rval)
00869   {
00870     ldbl_EGlpNumFreeArray (ddinfo->norms);
00871     ILL_IFFREE (ddinfo->refframe, int);
00872   }
00873   EG_RETURN(rval);
00874 }
00875 
00876 int ldbl_ILLprice_update_ddevex_norms (
00877   ldbl_lpinfo * const lp,
00878   ldbl_d_devex_info * const ddinfo,
00879   int const lindex,
00880   long double yl)
00881 {
00882   int i, r;
00883   long double normi;
00884   long double yr;
00885   long double ntmp,ntmp2;
00886 
00887   ldbl_EGlpNumInitVar (ntmp);
00888   ldbl_EGlpNumInitVar (ntmp2);
00889   ldbl_EGlpNumInitVar (normi);
00890   ldbl_EGlpNumInitVar (yr);
00891   ldbl_EGlpNumZero (normi);
00892 
00893   for (i = 0; i < lp->zA.nzcnt; i++)
00894     if (ddinfo->refframe[lp->nbaz[lp->zA.indx[i]]])
00895       ldbl_EGlpNumAddInnProdTo (normi, lp->zA.coef[i], lp->zA.coef[i]);
00896 
00897   if (ddinfo->refframe[lp->baz[lindex]])
00898     ldbl_EGlpNumAddTo (normi, ldbl_oneLpNum);
00899 
00900   ldbl_EGlpNumSet(ntmp,1000.0);
00901   ldbl_EGlpNumSet(ntmp2,0.001);
00902   ldbl_EGlpNumMultTo(ntmp,ddinfo->norms[lindex]);
00903   ldbl_EGlpNumMultTo(ntmp2,ddinfo->norms[lindex]);
00904   if (ldbl_EGlpNumIsLess(normi, ntmp2) || ldbl_EGlpNumIsLess(ntmp, normi))
00905   {
00906     ldbl_EGlpNumClearVar (normi);
00907     ldbl_EGlpNumClearVar (yr);
00908     ldbl_EGlpNumClearVar (ntmp);
00909     ldbl_EGlpNumClearVar (ntmp2);
00910     return ldbl_ILLprice_build_ddevex_norms (lp, ddinfo, 1);
00911   }
00912 
00913   for (i = 0; i < lp->yjz.nzcnt; i++)
00914   {
00915     r = lp->yjz.indx[i];
00916     ldbl_EGlpNumCopy(yr, lp->yjz.coef[i]);
00917     ldbl_EGlpNumCopy(ntmp,yr);
00918     ldbl_EGlpNumMultTo(ntmp,yr);
00919     ldbl_EGlpNumMultTo(ntmp,normi);
00920     ldbl_EGlpNumDivTo(ntmp,yl);
00921     ldbl_EGlpNumDivTo(ntmp,yl);
00922     if (ldbl_EGlpNumIsLess (ddinfo->norms[r], ntmp))
00923       ldbl_EGlpNumCopy (ddinfo->norms[r], ntmp);
00924   }
00925   ldbl_EGlpNumCopy (ddinfo->norms[lindex], normi);
00926   ldbl_EGlpNumDivTo(ddinfo->norms[lindex], yl);
00927   ldbl_EGlpNumDivTo(ddinfo->norms[lindex], yl);
00928   if (ldbl_EGlpNumIsLess (ddinfo->norms[lindex], ldbl_oneLpNum))
00929     ldbl_EGlpNumOne (ddinfo->norms[lindex]);
00930   ldbl_EGlpNumClearVar (normi);
00931   ldbl_EGlpNumClearVar (yr);
00932   ldbl_EGlpNumClearVar (ntmp);
00933   ldbl_EGlpNumClearVar (ntmp2);
00934   return 0;
00935 }
00936 
00937 int ldbl_ILLprice_build_dsteep_norms (
00938   ldbl_lpinfo * const lp,
00939   ldbl_d_steep_info * const dsinfo)
00940 {
00941   int i;
00942   int rval = 0;
00943   ldbl_svector z;
00944 
00945   ldbl_ILLsvector_init (&z);
00946   rval = ldbl_ILLsvector_alloc (&z, lp->nrows);
00947   CHECKRVALG(rval,CLEANUP);
00948   dsinfo->norms = ldbl_EGlpNumAllocArray (lp->nrows);
00949 
00950   for (i = 0; i < lp->nrows; i++)
00951   {
00952     rval = ILLstring_report (NULL, &lp->O->reporter);
00953     CHECKRVALG(rval,CLEANUP);
00954 
00955     ldbl_ILLfct_compute_zz (lp, &z, i);
00956 
00957     ldbl_EGlpNumInnProd (dsinfo->norms[i], z.coef, z.coef, (size_t) z.nzcnt);
00958     if (ldbl_EGlpNumIsLess (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM))
00959       ldbl_EGlpNumCopy (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM);
00960   }
00961 
00962 CLEANUP:
00963   ldbl_ILLsvector_free (&z);
00964   if (rval)
00965     ldbl_EGlpNumFreeArray (dsinfo->norms);
00966 
00967   EG_RETURN(rval);
00968 }
00969 
00970 int ldbl_ILLprice_get_dsteep_norms (
00971   ldbl_lpinfo * const lp,
00972   int const count,
00973   int *const rowind,
00974   long double * const norms)
00975 {
00976   int i;
00977   int rval = 0;
00978   ldbl_svector z;
00979 
00980   ldbl_ILLsvector_init (&z);
00981   rval = ldbl_ILLsvector_alloc (&z, lp->nrows);
00982   CHECKRVALG(rval,CLEANUP);
00983 
00984   for (i = 0; i < count; i++)
00985   {
00986     ldbl_ILLfct_compute_zz (lp, &z, rowind[i]);
00987     ldbl_EGlpNumInnProd (norms[i], z.coef, z.coef, (size_t) z.nzcnt);
00988   }
00989 
00990 CLEANUP:
00991   ldbl_ILLsvector_free (&z);
00992   EG_RETURN(rval);
00993 }
00994 
00995 void ldbl_ILLprice_update_dsteep_norms (
00996   ldbl_lpinfo * const lp,
00997   ldbl_d_steep_info * const dsinfo,
00998   ldbl_svector * const wz,
00999   int const lindex,
01000   long double yl)
01001 {
01002   int i, k;
01003   long double yij;
01004   long double norml;
01005   long double *v = 0;
01006   long double ntmp;
01007 
01008   ldbl_EGlpNumInitVar (ntmp);
01009   ldbl_EGlpNumInitVar (norml);
01010   ldbl_EGlpNumInitVar (yij);
01011   ldbl_EGlpNumInnProd (norml, lp->zz.coef, lp->zz.coef, (size_t) (lp->zz.nzcnt));
01012 
01013 #if 0
01014   Bico - remove warnings for dist
01015     if (fabs ((norml - dsinfo->norms[lindex]) / norml) > 1000.0 /*0.01 */ )
01016     {
01017       printf ("warning: incorrect dnorm values\n");
01018       printf ("anorm = %.6f, pnorm = %.6f\n", norml, dsinfo->norms[lindex]);
01019       fflush (stdout);
01020     }
01021 #endif
01022 
01023   ldbl_ILLfct_load_workvector (lp, wz);
01024   v = lp->work.coef;
01025 
01026   for (k = 0; k < lp->yjz.nzcnt; k++)
01027   {
01028     i = lp->yjz.indx[k];
01029     ldbl_EGlpNumCopy (yij, lp->yjz.coef[k]);
01030     /* compute in ntmp (yij * ((yij * norml / yl) - (2.0 * v[i]))) / yl; */
01031     ldbl_EGlpNumCopy(ntmp,yij);
01032     ldbl_EGlpNumMultTo(ntmp,norml);
01033     ldbl_EGlpNumDivTo(ntmp,yl);
01034     ldbl_EGlpNumSubTo(ntmp,v[i]);
01035     ldbl_EGlpNumSubTo(ntmp,v[i]);
01036     ldbl_EGlpNumMultTo (ntmp, yij);
01037     ldbl_EGlpNumDivTo (ntmp, yl);
01038     /* set dsinfo->norms[i] += (yij * ((yij * norml / yl) - (2.0 * v[i]))) / yl;*/
01039     ldbl_EGlpNumAddTo(dsinfo->norms[i], ntmp);
01040     if (ldbl_EGlpNumIsLess (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM))
01041       ldbl_EGlpNumCopy (dsinfo->norms[i], ldbl_PARAM_MIN_DNORM);
01042   }
01043   ldbl_EGlpNumCopyFrac (dsinfo->norms[lindex], norml, yl);
01044   ldbl_EGlpNumDivTo (dsinfo->norms[lindex], yl);
01045   if (ldbl_EGlpNumIsLess (dsinfo->norms[lindex], ldbl_PARAM_MIN_DNORM))
01046     ldbl_EGlpNumCopy (dsinfo->norms[lindex], ldbl_PARAM_MIN_DNORM);
01047 
01048   ldbl_ILLfct_zero_workvector (lp);
01049   ldbl_EGlpNumClearVar (norml);
01050   ldbl_EGlpNumClearVar (ntmp);
01051   ldbl_EGlpNumClearVar (yij);
01052 }
01053 
01054 static void ldbl_update_d_scaleinf (
01055   ldbl_price_info * const p,
01056   ldbl_heap * const h,
01057   int const j,
01058   long double inf,
01059   int const prule)
01060 {
01061   if (!ldbl_EGlpNumIsNeqqZero (inf))
01062   {
01063     ldbl_EGlpNumZero (p->d_scaleinf[j]);
01064     if (h->hexist != 0 && h->loc[j] != -1)
01065       ldbl_ILLheap_delete (h, j);
01066   }
01067   else
01068   {
01069     if (prule == QS_PRICE_PDANTZIG)
01070       ldbl_EGlpNumCopy (p->d_scaleinf[j], inf);
01071     else if (prule == QS_PRICE_PDEVEX)
01072       ldbl_EGlpNumCopySqrOver (p->d_scaleinf[j], inf, p->pdinfo.norms[j]);
01073     else if (prule == QS_PRICE_PSTEEP)
01074       ldbl_EGlpNumCopySqrOver (p->d_scaleinf[j], inf, p->psinfo.norms[j]);
01075 
01076     if (h->hexist != 0)
01077     {
01078       if (h->loc[j] == -1)
01079         ldbl_ILLheap_insert (h, j);
01080       else
01081         ldbl_ILLheap_modify (h, j);
01082     }
01083   }
01084 }
01085 
01086 static void ldbl_compute_dualI_inf (
01087   ldbl_lpinfo * const lp,
01088   const int j,
01089   long double * const inf)
01090 {
01091   int col = lp->nbaz[j];
01092   int vt = lp->vtype[col];
01093   int vs = lp->vstat[col];
01094   long double*dj = &(lp->pIdz[j]);
01095   long double*ftol = &(lp->tol->id_tol);
01096   ldbl_EGlpNumZero (*inf);
01097   if (vt != VARTIFICIAL && vt != VFIXED)
01098   {
01099     if( ldbl_EGlpNumIsSumLess(*dj,*ftol,ldbl_zeroLpNum) && (vs == STAT_LOWER || vs == STAT_ZERO))
01100       ldbl_EGlpNumCopyNeg(*inf,*dj);
01101     else if (ldbl_EGlpNumIsLess(*ftol, *dj) && (vs == STAT_UPPER || vs == STAT_ZERO))
01102       ldbl_EGlpNumCopy (*inf, *dj);
01103   }
01104 }
01105 
01106 static void ldbl_compute_dualII_inf (
01107   ldbl_lpinfo * const lp,
01108   int const j,
01109   long double * const inf)
01110 {
01111   int col = lp->nbaz[j];
01112   int vt = lp->vtype[col];
01113   int vs = lp->vstat[col];
01114   long double*dj = &(lp->dz[j]);
01115   long double*ftol = &(lp->tol->dfeas_tol);
01116   ldbl_EGlpNumZero (*inf);
01117   if (vt != VARTIFICIAL && vt != VFIXED)
01118   {
01119     if( ldbl_EGlpNumIsSumLess(*dj,*ftol,ldbl_zeroLpNum) && (vs == STAT_LOWER || vs == STAT_ZERO))
01120       ldbl_EGlpNumCopyNeg(*inf,*dj);
01121     else if (ldbl_EGlpNumIsLess(*ftol,*dj) && (vs == STAT_UPPER || vs == STAT_ZERO))
01122       ldbl_EGlpNumCopy (*inf, *dj);
01123   }
01124 }
01125 
01126 void ldbl_ILLprice_compute_dual_inf (
01127   ldbl_lpinfo * const lp,
01128   ldbl_price_info * const p,
01129   int *const ix,
01130   int const icnt,
01131   int const phase)
01132 {
01133   int i;
01134   int price;
01135   long double inf;
01136   ldbl_heap *h = &(p->h);
01137 
01138   price = (phase == PRIMAL_PHASEI) ? p->pI_price : p->pII_price;
01139   ldbl_EGlpNumInitVar (inf);
01140   ldbl_EGlpNumZero (inf);
01141 
01142   if (phase == PRIMAL_PHASEI)
01143   {
01144     if (ix == NULL)
01145       for (i = 0; i < lp->nnbasic; i++)
01146       {
01147         ldbl_compute_dualI_inf (lp, i, &(inf));
01148         ldbl_update_d_scaleinf (p, h, i, inf, price);
01149       }
01150     else
01151       for (i = 0; i < icnt; i++)
01152       {
01153         ldbl_compute_dualI_inf (lp, ix[i], &(inf));
01154         ldbl_update_d_scaleinf (p, h, ix[i], inf, price);
01155       }
01156   }
01157   else if (phase == PRIMAL_PHASEII)
01158   {
01159     if (ix == NULL)
01160       for (i = 0; i < lp->nnbasic; i++)
01161       {
01162         ldbl_compute_dualII_inf (lp, i, &inf);
01163         ldbl_update_d_scaleinf (p, h, i, inf, price);
01164       }
01165     else
01166       for (i = 0; i < icnt; i++)
01167       {
01168         ldbl_compute_dualII_inf (lp, ix[i], &inf);
01169         ldbl_update_d_scaleinf (p, h, ix[i], inf, price);
01170       }
01171   }
01172   ldbl_EGlpNumClearVar (inf);
01173 }
01174 
01175 void ldbl_ILLprice_primal (
01176   ldbl_lpinfo * const lp,
01177   ldbl_price_info * const pinf,
01178   ldbl_price_res * const pr,
01179   int const phase)
01180 {
01181   int j, vs;
01182   long double d_e, d_max;
01183   long double *ftol = &(lp->tol->dfeas_tol);
01184   ldbl_heap *const h = &(pinf->h);
01185 
01186   ldbl_EGlpNumInitVar (d_e);
01187   ldbl_EGlpNumInitVar (d_max);
01188   pr->eindex = -1;
01189   ldbl_EGlpNumZero(d_max);
01190 
01191 #if USEHEAP > 0
01192   ldbl_ILLprice_test_for_heap (lp, pinf, lp->nnbasic, pinf->d_scaleinf,
01193                           PRIMAL_SIMPLEX, 1);
01194 #endif
01195 
01196   if (pinf->p_strategy == COMPLETE_PRICING)
01197   {
01198     if (h->hexist)
01199     {
01200       pr->eindex = ldbl_ILLheap_findmin (h);
01201       if (pr->eindex != -1)
01202         ldbl_ILLheap_delete (h, pr->eindex);
01203     }
01204     else
01205     {
01206       for (j = 0; j < lp->nnbasic; j++)
01207       {
01208         if (ldbl_EGlpNumIsLess (d_max, pinf->d_scaleinf[j]))
01209         {
01210           ldbl_EGlpNumCopy (d_max, pinf->d_scaleinf[j]);
01211           pr->eindex = j;
01212         }
01213       }
01214     }
01215   }
01216   else if (pinf->p_strategy == MULTI_PART_PRICING)
01217   {
01218     for (j = 0; j < pinf->pmpinfo.bsize; j++)
01219     {
01220       if (ldbl_EGlpNumIsLess (d_max, pinf->pmpinfo.infeas[j]))
01221       {
01222         ldbl_EGlpNumCopy (d_max, pinf->pmpinfo.infeas[j]);
01223         pr->eindex = pinf->pmpinfo.bucket[j];
01224       }
01225     }
01226   }
01227 
01228   if (pr->eindex < 0)
01229     pr->price_stat = PRICE_OPTIMAL;
01230   else
01231   {
01232     if (phase == PRIMAL_PHASEI)
01233       ldbl_EGlpNumCopy (d_e, lp->pIdz[pr->eindex]);
01234     else
01235       ldbl_EGlpNumCopy (d_e, lp->dz[pr->eindex]);
01236     vs = lp->vstat[lp->nbaz[pr->eindex]];
01237 
01238     pr->price_stat = PRICE_NONOPTIMAL;
01239     if (vs == STAT_UPPER || (vs == STAT_ZERO && ldbl_EGlpNumIsLess (*ftol, d_e)))
01240       pr->dir = VDECREASE;
01241     else
01242       pr->dir = VINCREASE;
01243   }
01244   ldbl_EGlpNumClearVar (d_e);
01245   ldbl_EGlpNumClearVar (d_max);
01246 }
01247 
01248 static void ldbl_update_p_scaleinf (
01249   ldbl_price_info * const p,
01250   ldbl_heap * const h,
01251   int const i,
01252   long double inf,
01253   int const prule)
01254 {
01255   if (!ldbl_EGlpNumIsNeqqZero (inf))
01256   {
01257     ldbl_EGlpNumZero (p->p_scaleinf[i]);
01258     if (h->hexist != 0 && h->loc[i] != -1)
01259       ldbl_ILLheap_delete (h, i);
01260   }
01261   else
01262   {
01263     if (prule == QS_PRICE_DDANTZIG)
01264       ldbl_EGlpNumCopy (p->p_scaleinf[i], inf);
01265     else if (prule == QS_PRICE_DSTEEP)
01266       ldbl_EGlpNumCopySqrOver (p->p_scaleinf[i], inf, p->dsinfo.norms[i]);
01267     else if (prule == QS_PRICE_DDEVEX)
01268       ldbl_EGlpNumCopySqrOver (p->p_scaleinf[i], inf, p->ddinfo.norms[i]);
01269 
01270     if (h->hexist != 0)
01271     {
01272       if (h->loc[i] == -1)
01273         ldbl_ILLheap_insert (h, i);
01274       else
01275         ldbl_ILLheap_modify (h, i);
01276     }
01277   }
01278 }
01279 
01280 static void ldbl_compute_primalI_inf (
01281   ldbl_lpinfo * const lp,
01282   int const i,
01283   long double * const inf)
01284 {
01285   int const col = lp->baz[i];
01286   long double*x = &(lp->xbz[i]);
01287   long double*l = &(lp->lz[col]);
01288   long double*u = &(lp->uz[col]);
01289   long double*ftol = &(lp->tol->ip_tol);
01290   ldbl_EGlpNumZero (*inf);
01291 
01292   if (ldbl_EGlpNumIsLess (*ftol, *x) && ldbl_EGlpNumIsNeqq (*u, ldbl_INFTY))
01293     ldbl_EGlpNumCopy (*inf, *x);
01294   else if (ldbl_EGlpNumIsNeqq (*l, ldbl_NINFTY) && ldbl_EGlpNumIsSumLess (*x, *ftol,ldbl_zeroLpNum))
01295     ldbl_EGlpNumCopy (*inf, *x);
01296 }
01297 
01298 static void ldbl_compute_primalII_inf (
01299   ldbl_lpinfo * const lp,
01300   int const i,
01301   long double * const inf)
01302 {
01303   int const col = lp->baz[i];
01304   long double*x = &(lp->xbz[i]);
01305   long double*l = &(lp->lz[col]);
01306   long double*u = &(lp->uz[col]);
01307   long double*ftol = &(lp->tol->pfeas_tol);
01308   ldbl_EGlpNumZero (*inf);
01309 
01310   if (ldbl_EGlpNumIsNeqq (*u, ldbl_INFTY) && ldbl_EGlpNumIsSumLess (*u, *ftol, *x))
01311     ldbl_EGlpNumCopyDiff (*inf, *x, *u);
01312   else if (ldbl_EGlpNumIsNeqq (*l, ldbl_NINFTY) && ldbl_EGlpNumIsSumLess (*x, *ftol, *l))
01313     ldbl_EGlpNumCopyDiff (*inf, *l, *x);
01314 }
01315 
01316 void ldbl_ILLprice_compute_primal_inf (
01317   ldbl_lpinfo * const lp,
01318   ldbl_price_info * const p,
01319   int *const ix,
01320   int const icnt,
01321   int const phase)
01322 {
01323   int i;
01324   int price;
01325   long double inf;
01326   ldbl_heap *h = &(p->h);
01327 
01328   price = (phase == DUAL_PHASEI) ? p->dI_price : p->dII_price;
01329   ldbl_EGlpNumInitVar (inf);
01330   ldbl_EGlpNumZero (inf);
01331 
01332   if (phase == DUAL_PHASEI)
01333   {
01334     if (ix == NULL)
01335       for (i = 0; i < lp->nrows; i++)
01336       {
01337         ldbl_compute_primalI_inf (lp, i, &inf);
01338         ldbl_update_p_scaleinf (p, h, i, inf, price);
01339       }
01340     else
01341       for (i = 0; i < icnt; i++)
01342       {
01343         ldbl_compute_primalI_inf (lp, ix[i], &inf);
01344         ldbl_update_p_scaleinf (p, h, ix[i], inf, price);
01345       }
01346   }
01347   else if (phase == DUAL_PHASEII)
01348   {
01349     if (ix == NULL)
01350       for (i = 0; i < lp->nrows; i++)
01351       {
01352         ldbl_compute_primalII_inf (lp, i, &inf);
01353         ldbl_update_p_scaleinf (p, h, i, inf, price);
01354       }
01355     else
01356       for (i = 0; i < icnt; i++)
01357       {
01358         ldbl_compute_primalII_inf (lp, ix[i], &inf);
01359         ldbl_update_p_scaleinf (p, h, ix[i], inf, price);
01360       }
01361   }
01362   ldbl_EGlpNumClearVar (inf);
01363 }
01364 
01365 void ldbl_ILLprice_dual (
01366   ldbl_lpinfo * const lp,
01367   ldbl_price_info * const pinf,
01368   int const phase,
01369   ldbl_price_res * const pr)
01370 {
01371   int i;
01372   long double p_max;
01373   long double ubound;
01374   long double*ftol = &(lp->tol->pfeas_tol);
01375   ldbl_heap *const h = &(pinf->h);
01376 
01377   ldbl_EGlpNumInitVar (p_max);
01378   ldbl_EGlpNumInitVar (ubound);
01379   pr->lindex = -1;
01380   ldbl_EGlpNumZero(p_max);
01381 
01382 #if USEHEAP > 0
01383   ldbl_ILLprice_test_for_heap (lp, pinf, lp->nrows, pinf->p_scaleinf, DUAL_SIMPLEX,
01384                           1);
01385 #endif
01386 
01387   if (pinf->d_strategy == COMPLETE_PRICING)
01388   {
01389     if (h->hexist)
01390     {
01391       pr->lindex = ldbl_ILLheap_findmin (h);
01392       if (pr->lindex != -1)
01393         ldbl_ILLheap_delete (h, pr->lindex);
01394     }
01395     else
01396     {
01397       for (i = 0; i < lp->nrows; i++)
01398       {
01399         if (ldbl_EGlpNumIsLess (p_max, pinf->p_scaleinf[i]))
01400         {
01401           ldbl_EGlpNumCopy (p_max, pinf->p_scaleinf[i]);
01402           pr->lindex = i;
01403         }
01404       }
01405     }
01406   }
01407   else if (pinf->d_strategy == MULTI_PART_PRICING)
01408   {
01409     for (i = 0; i < pinf->dmpinfo.bsize; i++)
01410     {
01411       if (ldbl_EGlpNumIsLess (p_max, pinf->dmpinfo.infeas[i]))
01412       {
01413         ldbl_EGlpNumCopy (p_max, pinf->dmpinfo.infeas[i]);
01414         pr->lindex = pinf->dmpinfo.bucket[i];
01415       }
01416     }
01417   }
01418 
01419   if (pr->lindex < 0)
01420     pr->price_stat = PRICE_OPTIMAL;
01421   else
01422   {
01423     pr->price_stat = NONOPTIMAL;
01424 
01425     if (ldbl_EGlpNumIsNeqq (lp->uz[lp->baz[pr->lindex]], ldbl_INFTY))
01426     {
01427       if (phase == DUAL_PHASEI)
01428         ldbl_EGlpNumZero(ubound);
01429       else
01430         ldbl_EGlpNumCopy(ubound,lp->uz[lp->baz[pr->lindex]]);
01431       if (ldbl_EGlpNumIsSumLess (*ftol, ubound, lp->xbz[pr->lindex]))
01432         pr->lvstat = STAT_UPPER;
01433       else
01434         pr->lvstat = STAT_LOWER;
01435     }
01436     else
01437       pr->lvstat = STAT_LOWER;
01438   }
01439   ldbl_EGlpNumClearVar (p_max);
01440   ldbl_EGlpNumClearVar (ubound);
01441 }
01442 
01443 int ldbl_ILLprice_get_rownorms (
01444   ldbl_lpinfo * const lp,
01445   ldbl_price_info * const pinf,
01446   long double * const rnorms)
01447 {
01448   int rval = 0;
01449   int i;
01450 
01451   if (pinf->dsinfo.norms == NULL)
01452   {
01453     rval = ldbl_ILLprice_build_dsteep_norms (lp, &(pinf->dsinfo));
01454     CHECKRVALG(rval,CLEANUP);
01455   }
01456   for (i = 0; i < lp->nrows; i++)
01457     ldbl_EGlpNumCopy (rnorms[i], pinf->dsinfo.norms[i]);
01458 
01459 CLEANUP:
01460   if (rval)
01461     ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);
01462 
01463   return rval;
01464 }
01465 
01466 int ldbl_ILLprice_get_colnorms (
01467   ldbl_lpinfo * const lp,
01468   ldbl_price_info * const pinf,
01469   long double * const cnorms)
01470 {
01471   int rval = 0;
01472   int i, j;
01473 
01474   if (pinf->psinfo.norms == NULL)
01475   {
01476     rval = ldbl_ILLprice_build_psteep_norms (lp, &(pinf->psinfo));
01477     CHECKRVALG(rval,CLEANUP);
01478   }
01479   for (i = 0; i < lp->nrows; i++)
01480     ldbl_EGlpNumZero (cnorms[lp->baz[i]]);
01481   for (j = 0; j < lp->nnbasic; j++)
01482     ldbl_EGlpNumCopy (cnorms[lp->nbaz[j]], pinf->psinfo.norms[j]);
01483 
01484 CLEANUP:
01485   if (rval)
01486     ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
01487 
01488   return rval;
01489 }
01490 
01491 int ldbl_ILLprice_get_newnorms (
01492   ldbl_lpinfo * const lp,
01493   int const nelems,
01494   long double * const norms,
01495   int *const matcnt,
01496   int *const matbeg,
01497   int *const matind,
01498   long double * const matval,
01499   int const option)
01500 {
01501   int i, j;
01502   int rval = 0;
01503   ldbl_svector a;
01504   ldbl_svector y;
01505 
01506   ldbl_ILLsvector_init (&y);
01507   rval = ldbl_ILLsvector_alloc (&y, lp->nrows);
01508   CHECKRVALG(rval,CLEANUP);
01509 
01510   for (j = 0; j < nelems; j++)
01511   {
01512     a.nzcnt = matcnt[j];
01513     a.indx = &(matind[matbeg[j]]);
01514     a.coef = &(matval[matbeg[j]]);
01515 
01516     if (option == COLUMN_SOLVE)
01517       ldbl_ILLbasis_column_solve (lp, &a, &y);
01518     else
01519       ldbl_ILLbasis_row_solve (lp, &a, &y);
01520 
01521     ldbl_EGlpNumOne (norms[j]);
01522     for (i = 0; i < y.nzcnt; i++)
01523       ldbl_EGlpNumAddInnProdTo (norms[j], y.coef[i], y.coef[i]);
01524   }
01525 
01526 CLEANUP:
01527   ldbl_ILLsvector_free (&y);
01528   EG_RETURN(rval);
01529 }
01530 
01531 int ldbl_ILLprice_get_new_rownorms (
01532   ldbl_lpinfo * const lp,
01533   int const newrows,
01534   long double * const rnorms,
01535   int *const rmatcnt,
01536   int *const rmatbeg,
01537   int *const rmatind,
01538   long double * const rmatval)
01539 {
01540   return ldbl_ILLprice_get_newnorms (lp, newrows, rnorms, rmatcnt, rmatbeg, rmatind,
01541                                 rmatval, ROW_SOLVE);
01542 }
01543 
01544 int ldbl_ILLprice_get_new_colnorms (
01545   ldbl_lpinfo * const lp,
01546   int const newrows,
01547   long double * const rnorms,
01548   int *const matcnt,
01549   int *const matbeg,
01550   int *const matind,
01551   long double * const matval)
01552 {
01553   return ldbl_ILLprice_get_newnorms (lp, newrows, rnorms, matcnt, matbeg, matind,
01554                                 matval, COLUMN_SOLVE);
01555 }
01556 
01557 int ldbl_ILLprice_load_rownorms (
01558   ldbl_lpinfo * const lp,
01559   long double * const rnorms,
01560   ldbl_price_info * const pinf)
01561 {
01562   int i;
01563   int rval = 0;
01564 
01565   ldbl_EGlpNumFreeArray (pinf->dsinfo.norms);
01566   pinf->dsinfo.norms = ldbl_EGlpNumAllocArray (lp->nrows);
01567 
01568   for (i = 0; i < lp->nrows; i++)
01569   {
01570     ldbl_EGlpNumCopy (pinf->dsinfo.norms[i], rnorms[i]);
01571     if (ldbl_EGlpNumIsLess (pinf->dsinfo.norms[i], ldbl_PARAM_MIN_DNORM))
01572       ldbl_EGlpNumCopy (pinf->dsinfo.norms[i], ldbl_PARAM_MIN_DNORM);
01573   }
01574 
01575   EG_RETURN(rval);
01576 }
01577 
01578 int ldbl_ILLprice_load_colnorms (
01579   ldbl_lpinfo * const lp,
01580   long double * const cnorms,
01581   ldbl_price_info * const pinf)
01582 {
01583   int j;
01584   int rval = 0;
01585 
01586   ldbl_EGlpNumFreeArray (pinf->psinfo.norms);
01587   pinf->psinfo.norms = ldbl_EGlpNumAllocArray (lp->nnbasic);
01588 
01589   for (j = 0; j < lp->nnbasic; j++)
01590   {
01591     ldbl_EGlpNumCopy (pinf->psinfo.norms[j], cnorms[lp->nbaz[j]]);
01592     if (ldbl_EGlpNumIsLess (pinf->psinfo.norms[j], ldbl_oneLpNum))
01593       ldbl_EGlpNumOne (pinf->psinfo.norms[j]);
01594   }
01595 
01596   EG_RETURN(rval);
01597 }
01598 
01599 #if ldbl_PRICE_DEBUG > 0
01600 void ldbl_test_dsteep_norms (
01601   ldbl_lpinfo * lp,
01602   ldbl_price_info * p)
01603 {
01604   int i, errn = 0;
01605   long double *pn = ldbl_EGlpNumAllocArray(lp->nrows);
01606   long double err, diff;
01607   ldbl_EGlpNumZero (err);
01608 
01609   ldbl_EGlpNumInitVar (err);
01610   ldbl_EGlpNumInitVar (diff);
01611 
01612   ldbl_ILLprice_get_dsteep_norms (lp, lp->yjz.nzcnt, lp->yjz.indx, pn);
01613   for (i = 0; i < lp->yjz.nzcnt; i++)
01614   {
01615     ldbl_EGlpNumCopyDiff (diff, pn[i], p->dsinfo.norms[lp->yjz.indx[i]]);
01616     ldbl_EGlpNumCopyAbs(diff,diff);
01617     if (ldbl_EGlpNumIsLess (ldbl_PFEAS_TOLER, diff))
01618     {
01619       errn++;
01620       ldbl_EGlpNumAddTo (err, diff);
01621       ldbl_EGlpNumCopy (p->dsinfo.norms[lp->yjz.indx[i]], pn[i]);
01622     }
01623   }
01624   if (errn)
01625     printf ("%d: dnorm errn = %d, err = %.6f\n", lp->cnts->tot_iter, errn,
01626             ldbl_EGlpNumToLf (err));
01627   ldbl_EGlpNumFreeArray (pn);
01628   ldbl_EGlpNumClearVar (diff);
01629   ldbl_EGlpNumClearVar (err);
01630 }
01631 #endif

Generated on Thu Mar 29 09:32:31 2012 for QSopt_ex by  doxygen 1.4.7