dbl_fct.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: dbl_fct.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 static int TRACE = 0;
00025 
00026 //#define dbl_FCT_DEBUG 10
00027 #define dbl_FCT_DEBUG 0
00028 
00029 #include "qs_config.h"
00030 #include "dbl_iqsutil.h"
00031 #include "dbl_lpdefs.h"
00032 #include "stddefs.h"
00033 #include "dbl_basis.h"
00034 #include "dbl_fct.h"
00035 #include "dbl_price.h"
00036 #include "dbl_ratio.h"
00037 #include "dbl_dstruct.h"
00038 
00039 dbl_bndinfo *dbl_ILLfct_new_bndinfo (
00040   void)
00041 {
00042   dbl_bndinfo *nbnd = (dbl_bndinfo *) malloc (sizeof (dbl_bndinfo));
00043 
00044   if (!nbnd)
00045   {
00046     fprintf (stderr, "not enough memory, in %s\n", __func__);
00047     exit (1);
00048   }
00049   dbl_EGlpNumInitVar ((nbnd->pbound));
00050   dbl_EGlpNumInitVar ((nbnd->cbound));
00051   return nbnd;
00052 }
00053 
00054 void dbl_ILLfct_free_bndinfo (
00055   dbl_bndinfo * binfo)
00056 {
00057   dbl_EGlpNumClearVar ((binfo->pbound));
00058   dbl_EGlpNumClearVar ((binfo->cbound));
00059   ILL_IFFREE (binfo, dbl_bndinfo);
00060   return;
00061 }
00062 
00063 static int dbl_compute_zA1 (
00064   dbl_lpinfo * lp,
00065   dbl_svector * z,
00066   dbl_svector * zA,
00067   double ztoler),
00068 /*
00069   compute_zA2 (dbl_lpinfo * lp,
00070                dbl_svector * z,
00071                dbl_svector * zA,
00072                const double* ztoler), */
00073   dbl_compute_zA3 (
00074   dbl_lpinfo * lp,
00075   dbl_svector * z,
00076   dbl_svector * zA,
00077   double ztoler),
00078   dbl_expand_var_bounds (
00079   dbl_lpinfo * lp,
00080   double ftol,
00081   int *chgb),
00082   dbl_expand_var_coefs (
00083   dbl_lpinfo * lp,
00084   double ftol,
00085   int *chgc);
00086 
00087 static void dbl_update_piv_values (
00088   dbl_count_struct * c,
00089   int phase,
00090   const double piv),
00091 /*  copy_vectors (dbl_svector * a,
00092                 dbl_svector * b),*/
00093   dbl_add_vectors (
00094   dbl_lpinfo * lp,
00095   dbl_svector * a,
00096   dbl_svector * b,
00097   dbl_svector * c,
00098   const double t);
00099 
00100 static double dbl_my_rand (
00101   int bound,
00102   ILLrandstate * r);
00103 
00104 
00105 void dbl_ILLfct_load_workvector (
00106   dbl_lpinfo * lp,
00107   dbl_svector * s)
00108 {
00109   int i;
00110 
00111   for (i = 0; i < s->nzcnt; i++)
00112   {
00113     lp->work.indx[i] = s->indx[i];
00114     dbl_EGlpNumCopy (lp->work.coef[s->indx[i]], s->coef[i]);
00115   }
00116   lp->work.nzcnt = s->nzcnt;
00117 }
00118 
00119 void dbl_ILLfct_zero_workvector (
00120   dbl_lpinfo * lp)
00121 {
00122   int i;
00123 
00124   for (i = 0; i < lp->work.nzcnt; i++)
00125     dbl_EGlpNumZero (lp->work.coef[lp->work.indx[i]]);
00126   lp->work.nzcnt = 0;
00127 }
00128 
00129 void dbl_ILLfct_set_variable_type (
00130   dbl_lpinfo * lp)
00131 {
00132   int j;
00133 
00134   for (j = 0; j < lp->ncols; j++)
00135   {
00136 
00137     if (lp->matcnt[j] == 1 && lp->O->rowmap[lp->matind[lp->matbeg[j]]] == j)
00138       lp->vclass[j] = CLASS_LOGICAL;
00139     else
00140       lp->vclass[j] = CLASS_STRUCT;
00141     switch ((dbl_EGlpNumIsEqqual (lp->uz[j], dbl_INFTY) ? 1U : 0U) |
00142             (dbl_EGlpNumIsEqqual (lp->lz[j], dbl_NINFTY) ? 2U : 0U))
00143     {
00144     case 0:
00145       if (dbl_EGlpNumIsLess (lp->lz[j], lp->uz[j]))
00146         lp->vtype[j] = VBOUNDED;
00147       else if (!dbl_EGlpNumIsNeqqZero (lp->lz[j]) &&
00148                (lp->vclass[j] == CLASS_LOGICAL))
00149         lp->vtype[j] = VARTIFICIAL;
00150       else
00151         lp->vtype[j] = VFIXED;
00152       break;
00153     case 3:
00154       lp->vtype[j] = VFREE;
00155       break;
00156     case 1:
00157       lp->vtype[j] = VLOWER;
00158       break;
00159     case 2:
00160       lp->vtype[j] = VUPPER;
00161       break;
00162     }
00163   }
00164 }
00165 
00166 /* compute various vectors */
00167 
00168 void dbl_ILLfct_compute_pobj (
00169   dbl_lpinfo * lp)
00170 {
00171   int i, j;
00172   int col;
00173   double sum;
00174 
00175   dbl_EGlpNumInitVar (sum);
00176   dbl_EGlpNumZero (sum);
00177 
00178   for (i = 0; i < lp->nrows; i++)
00179     dbl_EGlpNumAddInnProdTo (sum, lp->cz[lp->baz[i]], lp->xbz[i]);
00180 
00181   for (j = 0; j < lp->nnbasic; j++)
00182   {
00183     col = lp->nbaz[j];
00184     if (lp->vstat[col] == STAT_UPPER)
00185       dbl_EGlpNumAddInnProdTo (sum, lp->cz[col], lp->uz[col]);
00186     else if (lp->vstat[col] == STAT_LOWER)
00187       dbl_EGlpNumAddInnProdTo (sum, lp->cz[col], lp->lz[col]);
00188   }
00189   dbl_EGlpNumCopy (lp->pobjval, sum);
00190   dbl_EGlpNumCopy (lp->objval, sum);
00191   dbl_EGlpNumClearVar (sum);
00192 }
00193 
00194 void dbl_ILLfct_compute_dobj (
00195   dbl_lpinfo * lp)
00196 {
00197   int i, j;
00198   int col;
00199   double sum;
00200 
00201   dbl_EGlpNumInitVar (sum);
00202   dbl_EGlpNumZero (sum);
00203 
00204   for (i = 0; i < lp->nrows; i++)
00205     dbl_EGlpNumAddInnProdTo (sum, lp->piz[i], lp->bz[i]);
00206 
00207   for (j = 0; j < lp->nnbasic; j++)
00208   {
00209     col = lp->nbaz[j];
00210     if (lp->vstat[col] == STAT_UPPER)
00211       dbl_EGlpNumAddInnProdTo (sum, lp->dz[j], lp->uz[col]);
00212     else if (lp->vstat[col] == STAT_LOWER)
00213       dbl_EGlpNumAddInnProdTo (sum, lp->dz[j], lp->lz[col]);
00214   }
00215   dbl_EGlpNumCopy (lp->dobjval, sum);
00216   dbl_EGlpNumCopy (lp->objval, sum);
00217   dbl_EGlpNumClearVar (sum);
00218 }
00219 
00220 void dbl_ILLfct_compute_xbz (
00221   dbl_lpinfo * lp)
00222 {
00223   int i, j, r;
00224   int col, mcnt, mbeg;
00225   dbl_svector *srhs = &(lp->srhs);
00226   dbl_svector *ssoln = &(lp->ssoln);
00227   double xval;
00228 
00229   dbl_EGlpNumInitVar (xval);
00230 
00231   for (i = 0; i < lp->nrows; i++)
00232   {
00233     dbl_EGlpNumZero (lp->xbz[i]);
00234     dbl_EGlpNumCopy (srhs->coef[i], lp->bz[i]);
00235   }
00236   for (j = 0; j < lp->nnbasic; j++)
00237   {
00238     col = lp->nbaz[j];
00239     dbl_EGlpNumZero (xval);
00240     if (lp->vstat[col] == STAT_UPPER && dbl_EGlpNumIsNeqqZero (lp->uz[col]))
00241       dbl_EGlpNumCopy (xval, lp->uz[col]);
00242     else if (lp->vstat[col] == STAT_LOWER && dbl_EGlpNumIsNeqqZero (lp->lz[col]))
00243       dbl_EGlpNumCopy (xval, lp->lz[col]);
00244 
00245     if (dbl_EGlpNumIsNeqqZero (xval))
00246     {
00247       mcnt = lp->matcnt[col];
00248       mbeg = lp->matbeg[col];
00249       for (i = 0; i < mcnt; i++)
00250         dbl_EGlpNumSubInnProdTo (srhs->coef[lp->matind[mbeg + i]], xval,
00251                              lp->matval[mbeg + i]);
00252     }
00253   }
00254   for (i = 0, r = 0; i < lp->nrows; i++)
00255     if (dbl_EGlpNumIsNeqqZero (srhs->coef[i]))
00256     {
00257       dbl_EGlpNumCopy (srhs->coef[r], srhs->coef[i]);
00258       srhs->indx[r] = i;
00259       r++;
00260     }
00261   srhs->nzcnt = r;
00262 
00263   dbl_ILLbasis_column_solve (lp, srhs, ssoln);
00264   for (i = 0; i < ssoln->nzcnt; i++)
00265     dbl_EGlpNumCopy (lp->xbz[ssoln->indx[i]], ssoln->coef[i]);
00266   dbl_EGlpNumClearVar (xval);
00267 }
00268 
00269 void dbl_ILLfct_compute_piz (
00270   dbl_lpinfo * lp)
00271 {
00272   int i, r;
00273   dbl_svector *srhs = &(lp->srhs);
00274   dbl_svector *ssoln = &(lp->ssoln);
00275 
00276   for (i = 0, r = 0; i < lp->nrows; i++)
00277   {
00278     dbl_EGlpNumZero (lp->piz[i]);
00279     if (dbl_EGlpNumIsNeqqZero (lp->cz[lp->baz[i]]))
00280     {
00281       srhs->indx[r] = i;
00282       dbl_EGlpNumCopy (srhs->coef[r], lp->cz[lp->baz[i]]);
00283       r++;
00284     }
00285   }
00286   srhs->nzcnt = r;
00287 
00288   dbl_ILLbasis_row_solve (lp, srhs, ssoln);
00289   for (i = 0; i < ssoln->nzcnt; i++)
00290     dbl_EGlpNumCopy (lp->piz[ssoln->indx[i]], ssoln->coef[i]);
00291 }
00292 
00293 void dbl_ILLfct_compute_dz (
00294   dbl_lpinfo * lp)
00295 {
00296   int i, j;
00297   int col;
00298   int mcnt, mbeg;
00299   double sum;
00300 
00301   dbl_EGlpNumInitVar (sum);
00302 
00303   for (j = 0; j < lp->nnbasic; j++)
00304   {
00305     dbl_EGlpNumZero (sum);
00306     col = lp->nbaz[j];
00307     mcnt = lp->matcnt[col];
00308     mbeg = lp->matbeg[col];
00309     for (i = 0; i < mcnt; i++)
00310       dbl_EGlpNumAddInnProdTo (sum, lp->piz[lp->matind[mbeg + i]],
00311                            lp->matval[mbeg + i]);
00312     dbl_EGlpNumCopyDiff (lp->dz[j], lp->cz[col], sum);
00313   }
00314   dbl_EGlpNumClearVar (sum);
00315 }
00316 
00317 void dbl_ILLfct_compute_phaseI_xbz (
00318   dbl_lpinfo * lp)
00319 {
00320   int i, j, r;
00321   int col, mcnt, mbeg;
00322   dbl_svector *srhs = &(lp->srhs);
00323   dbl_svector *ssoln = &(lp->ssoln);
00324 
00325   for (i = 0; i < lp->nrows; i++)
00326   {
00327     dbl_EGlpNumZero (lp->xbz[i]);
00328     dbl_EGlpNumZero (srhs->coef[i]);
00329   }
00330   for (j = 0; j < lp->nnbasic; j++)
00331   {
00332     col = lp->nbaz[j];
00333 
00334     if (lp->dfeas[j])
00335     {
00336       mcnt = lp->matcnt[col];
00337       mbeg = lp->matbeg[col];
00338       if (lp->dfeas[j] == -1)
00339         for (i = 0; i < mcnt; i++)
00340           dbl_EGlpNumSubTo (srhs->coef[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00341       else
00342         for (i = 0; i < mcnt; i++)
00343           dbl_EGlpNumAddTo (srhs->coef[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00344     }
00345   }
00346   for (i = 0, r = 0; i < lp->nrows; i++)
00347     if (dbl_EGlpNumIsNeqqZero (srhs->coef[i]))
00348     {
00349       dbl_EGlpNumCopy (srhs->coef[r], srhs->coef[i]);
00350       srhs->indx[r] = i;
00351       r++;
00352     }
00353   srhs->nzcnt = r;
00354 
00355   dbl_ILLbasis_column_solve (lp, srhs, ssoln);
00356   for (i = 0; i < ssoln->nzcnt; i++)
00357     dbl_EGlpNumCopy (lp->xbz[ssoln->indx[i]], ssoln->coef[i]);
00358 }
00359 
00360 void dbl_ILLfct_compute_phaseI_piz (
00361   dbl_lpinfo * lp)
00362 {
00363   int i, r;
00364   dbl_svector *srhs = &(lp->srhs);
00365   dbl_svector *ssoln = &(lp->ssoln);
00366 
00367   for (i = 0, r = 0; i < lp->nrows; i++)
00368   {
00369     dbl_EGlpNumZero (lp->pIpiz[i]);
00370     if (lp->bfeas[i] != 0)
00371     {
00372       srhs->indx[r] = i;
00373       dbl_EGlpNumSet (srhs->coef[r], (double) lp->bfeas[i]);
00374       r++;
00375     }
00376   }
00377   srhs->nzcnt = r;
00378 
00379   dbl_ILLbasis_row_solve (lp, srhs, ssoln);
00380   for (i = 0; i < ssoln->nzcnt; i++)
00381     dbl_EGlpNumCopy (lp->pIpiz[ssoln->indx[i]], ssoln->coef[i]);
00382   dbl_ILLfct_update_counts (lp, CNT_P1PINZ, ssoln->nzcnt, dbl_zeroLpNum);
00383 }
00384 
00385 void dbl_ILLfct_compute_phaseI_dz (
00386   dbl_lpinfo * lp)
00387 {
00388   int i, j;
00389   int col;
00390   int mcnt, mbeg;
00391   double sum;
00392 
00393   dbl_EGlpNumInitVar (sum);
00394   ILL_IFTRACE ("%s\n", __func__);
00395 
00396   for (j = 0; j < lp->nnbasic; j++)
00397   {
00398     dbl_EGlpNumZero (sum);
00399     col = lp->nbaz[j];
00400     mcnt = lp->matcnt[col];
00401     mbeg = lp->matbeg[col];
00402     for (i = 0; i < mcnt; i++)
00403       dbl_EGlpNumAddInnProdTo (sum, lp->pIpiz[lp->matind[mbeg + i]],
00404                            lp->matval[mbeg + i]);
00405     dbl_EGlpNumCopyNeg (lp->pIdz[j], sum);
00406     ILL_IFTRACE ("%d:%d:%lf:%la\n", j, col, dbl_EGlpNumToLf (sum),
00407                  dbl_EGlpNumToLf (sum));
00408   }
00409   dbl_EGlpNumClearVar (sum);
00410 }
00411 
00412 void dbl_ILLfct_compute_yz (
00413   dbl_lpinfo * lp,
00414   dbl_svector * yz,
00415   dbl_svector * updz,
00416   int col)
00417 {
00418   dbl_svector a;
00419 
00420   a.nzcnt = lp->matcnt[col];
00421   a.indx = &(lp->matind[lp->matbeg[col]]);
00422   a.coef = &(lp->matval[lp->matbeg[col]]);
00423 
00424   dbl_ILLfactor_set_factor_dparam (lp->f, QS_FACTOR_SZERO_TOL, dbl_PIVZ_TOLER);
00425   if (updz)
00426     dbl_ILLbasis_column_solve_update (lp, &a, updz, yz);
00427   else
00428     dbl_ILLbasis_column_solve (lp, &a, yz);
00429   dbl_ILLfactor_set_factor_dparam (lp->f, QS_FACTOR_SZERO_TOL, dbl_SZERO_TOLER);
00430 }
00431 
00432 void dbl_ILLfct_compute_zz (
00433   dbl_lpinfo * lp,
00434   dbl_svector * zz,
00435   int row)
00436 {
00437   dbl_ILLfct_compute_binvrow (lp, zz, row, dbl_PIVZ_TOLER);
00438 }
00439 
00440 void dbl_ILLfct_compute_binvrow (
00441   dbl_lpinfo * lp,
00442   dbl_svector * zz,
00443   int row,
00444   double ztoler)
00445 {
00446   dbl_svector a;
00447   double e;
00448 
00449   dbl_EGlpNumInitVar (e);
00450   dbl_EGlpNumOne (e);
00451 
00452   a.nzcnt = 1;
00453   a.coef = &e;
00454   a.indx = &row;
00455 
00456   if (dbl_EGlpNumIsGreatZero (ztoler))
00457     dbl_ILLfactor_set_factor_dparam (lp->f, QS_FACTOR_SZERO_TOL, ztoler);
00458   dbl_ILLbasis_row_solve (lp, &a, zz);
00459   if (dbl_EGlpNumIsGreatZero (ztoler))
00460     dbl_ILLfactor_set_factor_dparam (lp->f, QS_FACTOR_SZERO_TOL, dbl_SZERO_TOLER);
00461   dbl_EGlpNumClearVar (e);
00462 }
00463 
00464 void dbl_ILLfct_compute_psteep_upv (
00465   dbl_lpinfo * lp,
00466   dbl_svector * swz)
00467 {
00468   dbl_ILLbasis_row_solve (lp, &(lp->yjz), swz);
00469 }
00470 
00471 void dbl_ILLfct_compute_dsteep_upv (
00472   dbl_lpinfo * lp,
00473   dbl_svector * swz)
00474 {
00475   dbl_ILLbasis_column_solve (lp, &(lp->zz), swz);
00476 }
00477 
00478 static int dbl_compute_zA1 (
00479   dbl_lpinfo * lp,
00480   dbl_svector * z,
00481   dbl_svector * zA,
00482   double ztoler)
00483 {
00484   int rval = 0;
00485   int i, j, nz = 0;
00486   int col, mcnt, mbeg;
00487   double sum;
00488   double *v = 0;
00489 
00490   dbl_EGlpNumInitVar (sum);
00491   v = dbl_EGlpNumAllocArray (lp->nrows);
00492 
00493   for (i = 0; i < lp->nrows; i++)
00494     dbl_EGlpNumZero (v[i]);
00495   for (i = 0; i < z->nzcnt; i++)
00496     dbl_EGlpNumCopy (v[z->indx[i]], z->coef[i]);
00497 
00498   for (j = 0; j < lp->nnbasic; j++)
00499   {
00500     dbl_EGlpNumZero (sum);
00501     col = lp->nbaz[j];
00502     mcnt = lp->matcnt[col];
00503     mbeg = lp->matbeg[col];
00504     for (i = 0; i < mcnt; i++)
00505       dbl_EGlpNumAddInnProdTo (sum, v[lp->matind[mbeg + i]], lp->matval[mbeg + i]);
00506 
00507     if (dbl_EGlpNumIsNeqZero (sum, ztoler))
00508     {
00509       dbl_EGlpNumCopy (zA->coef[nz], sum);
00510       zA->indx[nz] = j;
00511       nz++;
00512     }
00513   }
00514   zA->nzcnt = nz;
00515 
00516   dbl_EGlpNumClearVar (sum);
00517   dbl_EGlpNumFreeArray (v);
00518   EG_RETURN (rval);
00519 }
00520 
00521 
00522 static int dbl_compute_zA3 (
00523   dbl_lpinfo * lp,
00524   dbl_svector * z,
00525   dbl_svector * zA,
00526   double ztoler)
00527 {
00528   int rval = 0;
00529   int i, j, k, ix;
00530   int nz = 0;
00531   int row, col;
00532   int rcnt, rbeg;
00533   double val;
00534 
00535   dbl_EGlpNumInitVar (val);
00536   k = 0;
00537   for (i = 0; i < z->nzcnt; i++)
00538   {
00539     row = z->indx[i];
00540     dbl_EGlpNumCopy (val, z->coef[i]);
00541     rcnt = lp->rowcnt[row];
00542     rbeg = lp->rowbeg[row];
00543     for (j = 0; j < rcnt; j++)
00544     {
00545       col = lp->rowind[rbeg + j];
00546       if (lp->vstat[col] != STAT_BASIC)
00547       {
00548         ix = lp->vindex[col];
00549         if (lp->iwork[ix] == 0)
00550         {
00551           lp->iwork[ix] = 1;
00552           lp->work.indx[k++] = ix;
00553         }
00554         dbl_EGlpNumAddInnProdTo (lp->work.coef[ix], val, lp->rowval[rbeg + j]);
00555       }
00556     }
00557   }
00558   for (j = 0; j < k; j++)
00559   {
00560     ix = lp->work.indx[j];
00561     dbl_EGlpNumCopy (val, lp->work.coef[ix]);
00562     dbl_EGlpNumZero (lp->work.coef[ix]);
00563     lp->iwork[ix] = 0;
00564     if (dbl_EGlpNumIsNeqZero (val, ztoler))
00565     {
00566       dbl_EGlpNumCopy (zA->coef[nz], val);
00567       zA->indx[nz] = ix;
00568       nz++;
00569     }
00570   }
00571   zA->nzcnt = nz;
00572   dbl_EGlpNumClearVar (val);
00573   EG_RETURN (rval);
00574 }
00575 
00576 int dbl_ILLfct_compute_zA (
00577   dbl_lpinfo * lp,
00578   dbl_svector * z,
00579   dbl_svector * zA)
00580 {
00581   if (z->nzcnt < lp->nrows / 2)
00582     return dbl_compute_zA3 (lp, z, zA, dbl_PIVZ_TOLER);
00583   else
00584     return dbl_compute_zA1 (lp, z, zA, dbl_PIVZ_TOLER);
00585 }
00586 
00587 /* compute v^T A */
00588 void dbl_ILLfct_compute_vA (
00589   dbl_lpinfo * lp,
00590   dbl_svector * v,
00591   double * vA)
00592 {
00593   int i, j;
00594   int row, col;
00595   int rcnt, rbeg;
00596   double val;
00597 
00598   dbl_EGlpNumInitVar (val);
00599 
00600   for (j = 0; j < lp->ncols; j++)
00601     dbl_EGlpNumZero (vA[j]);
00602 
00603   for (i = 0; i < v->nzcnt; i++)
00604   {
00605     row = v->indx[i];
00606     dbl_EGlpNumCopy (val, v->coef[i]);
00607     rcnt = lp->rowcnt[row];
00608     rbeg = lp->rowbeg[row];
00609     for (j = 0; j < rcnt; j++)
00610     {
00611       col = lp->rowind[rbeg + j];
00612       dbl_EGlpNumAddInnProdTo (vA[col], val, lp->rowval[rbeg + j]);
00613     }
00614   }
00615 
00616   for (j = 0; j < lp->ncols; j++)
00617     if (!dbl_EGlpNumIsNeqZero (vA[j], dbl_SZERO_TOLER))
00618       dbl_EGlpNumZero (vA[j]);
00619 
00620   dbl_EGlpNumClearVar (val);
00621   return;
00622 }
00623 
00624 /* update information */
00625 
00626 /*
00627 1) lvstat - new status of leaving var.
00628 */
00629 void dbl_ILLfct_update_basis_info (
00630   dbl_lpinfo * lp,
00631   int eindex,
00632   int lindex,
00633   int lvstat)
00634 {
00635   int evar;
00636   int lvar;
00637 
00638   evar = lp->nbaz[eindex];
00639 
00640   if (lindex >= 0)
00641   {                             /* variable leaves basis */
00642     lvar = lp->baz[lindex];
00643     lp->vstat[evar] = STAT_BASIC;
00644     lp->vstat[lvar] = lvstat;
00645     lp->vindex[evar] = lindex;
00646     lp->vindex[lvar] = eindex;
00647     lp->baz[lindex] = evar;
00648     lp->nbaz[eindex] = lvar;
00649     (lp->basisid)++;
00650   }
00651   else
00652   {
00653     lp->vstat[evar] = (lp->vstat[evar] == STAT_LOWER) ? STAT_UPPER : STAT_LOWER;
00654   }
00655 }
00656 
00657 void dbl_ILLfct_update_xz (
00658   dbl_lpinfo * lp,
00659   double tz,
00660   int eindex,
00661   int lindex)
00662 {
00663   int i, evar, estat;
00664 
00665   ILL_IFTRACE ("%s:%la:%d:%d:%d\n", __func__, dbl_EGlpNumToLf (tz), eindex,
00666                lindex, lp->yjz.nzcnt);
00667 
00668   if (dbl_EGlpNumIsNeqqZero (tz))
00669     for (i = 0; i < lp->yjz.nzcnt; i++)
00670       dbl_EGlpNumSubInnProdTo (lp->xbz[lp->yjz.indx[i]], tz, lp->yjz.coef[i]);
00671 
00672   if (lindex >= 0)
00673   {                             /* variable leaves basis */
00674     evar = lp->nbaz[eindex];
00675     estat = lp->vstat[evar];
00676     if (estat == STAT_LOWER)
00677       dbl_EGlpNumCopySum (lp->xbz[lindex], lp->lz[evar], tz);
00678     else if (estat == STAT_UPPER)
00679       dbl_EGlpNumCopySum (lp->xbz[lindex], lp->uz[evar], tz);
00680     else if (estat == STAT_ZERO)
00681       dbl_EGlpNumCopy (lp->xbz[lindex], tz);
00682   }
00683 }
00684 
00685 void dbl_ILLfct_update_piz (
00686   dbl_lpinfo * lp,
00687   double alpha)
00688 {
00689   int i;
00690 
00691   for (i = 0; i < lp->zz.nzcnt; i++)
00692     dbl_EGlpNumAddInnProdTo (lp->piz[lp->zz.indx[i]], alpha, lp->zz.coef[i]);
00693 }
00694 
00695 void dbl_ILLfct_update_pIpiz (
00696   dbl_lpinfo * lp,
00697   dbl_svector * z,
00698   const double alpha)
00699 {
00700   int i;
00701 
00702   if (!dbl_EGlpNumIsNeqqZero (alpha))
00703     return;
00704   if (dbl_EGlpNumIsEqqual (alpha, dbl_oneLpNum))
00705   {
00706     for (i = 0; i < z->nzcnt; i++)
00707       dbl_EGlpNumAddTo (lp->pIpiz[z->indx[i]], z->coef[i]);
00708   }
00709   else
00710   {
00711     for (i = 0; i < z->nzcnt; i++)
00712       dbl_EGlpNumAddInnProdTo (lp->pIpiz[z->indx[i]], alpha, z->coef[i]);
00713   }
00714 }
00715 
00716 void dbl_ILLfct_update_dz (
00717   dbl_lpinfo * lp,
00718   int eindex,
00719   double alpha)
00720 {
00721   int i;
00722 
00723   for (i = 0; i < lp->zA.nzcnt; i++)
00724     dbl_EGlpNumSubInnProdTo (lp->dz[lp->zA.indx[i]], alpha, lp->zA.coef[i]);
00725   dbl_EGlpNumCopyNeg (lp->dz[eindex], alpha);
00726 }
00727 
00728 void dbl_ILLfct_update_pIdz (
00729   dbl_lpinfo * lp,
00730   dbl_svector * zA,
00731   int eindex,
00732   const double alpha)
00733 {
00734   int i;
00735 
00736   if (!dbl_EGlpNumIsNeqqZero (alpha))
00737     return;
00738 
00739   if (dbl_EGlpNumIsEqqual (alpha, dbl_oneLpNum))
00740   {
00741     for (i = 0; i < zA->nzcnt; i++)
00742       dbl_EGlpNumSubTo (lp->pIdz[zA->indx[i]], zA->coef[i]);
00743   }
00744   else
00745   {
00746     for (i = 0; i < zA->nzcnt; i++)
00747       dbl_EGlpNumSubInnProdTo (lp->pIdz[zA->indx[i]], alpha, zA->coef[i]);
00748   }
00749   if (eindex > -1)
00750     dbl_EGlpNumCopyNeg (lp->pIdz[eindex], alpha);
00751 }
00752 
00753 /* bound and coef shift routines */
00754 
00755 /* scale bound in dbl_my_rand to get more random digits, unless bound is large */
00756 static double dbl_my_rand (
00757   int bound,
00758   ILLrandstate * r)
00759 {
00760   int k = bound, scale = 1;
00761   double v = 0.0;
00762 
00763   if (bound < 100000)
00764   {
00765     k = 20000 * bound;
00766     scale = 20000;
00767   }
00768   v = 1 + (ILLutil_lprand (r) % (k));
00769   return v / (double) scale;
00770 }
00771 
00772 static int dbl_expand_var_bounds (
00773   dbl_lpinfo * lp,
00774   double ftol,
00775   int *chgb)
00776 {
00777   int rval = 0;
00778   int i, col, nchg = 0;
00779   double newb, cftol;
00780   double *x, *l, *u;
00781   ILLrandstate r;
00782 
00783   dbl_EGlpNumInitVar (newb);
00784   dbl_EGlpNumInitVar (cftol);
00785   dbl_EGlpNumCopyAbs (cftol, ftol);
00786   dbl_EGlpNumDivUiTo (cftol, 10);
00787 
00788   ILLutil_sprand (1, &r);
00789 
00790   for (i = 0; i < lp->nrows; i++)
00791   {
00792     col = lp->baz[i];
00793     if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFREE)
00794       continue;
00795     x = &(lp->xbz[i]);
00796     l = &(lp->lz[col]);
00797     u = &(lp->uz[col]);
00798     /* we use newb as temporal variable outside the if's scope */
00799     dbl_EGlpNumCopyDiff (newb, *x, ftol);
00800     if (dbl_EGlpNumIsNeqq (*l, dbl_NINFTY) && dbl_EGlpNumIsLess (newb, *l))
00801     {
00802       dbl_EGlpNumSet (newb, -1.0 * (dbl_my_rand (50, &(lp->rstate)) + 1.0));
00803       dbl_EGlpNumMultTo (newb, cftol);
00804       if (dbl_EGlpNumIsLess (*x, *l))
00805         dbl_EGlpNumAddTo (newb, *x);
00806       else
00807         dbl_EGlpNumAddTo (newb, *l);
00808       rval = dbl_ILLfct_bound_shift (lp, col, BOUND_LOWER, newb);
00809       CHECKRVALG (rval, CLEANUP);
00810       nchg++;
00811     }
00812     dbl_EGlpNumCopySum (newb, *x, ftol);
00813     if (dbl_EGlpNumIsNeqq (*u, dbl_INFTY) && dbl_EGlpNumIsLess (*u, newb))
00814     {
00815       dbl_EGlpNumSet (newb, dbl_my_rand (50, &(lp->rstate)) + 1.0);
00816       dbl_EGlpNumMultTo (newb, cftol);
00817       if (dbl_EGlpNumIsLess (*x, *u))
00818         dbl_EGlpNumAddTo (newb, *u);
00819       else
00820         dbl_EGlpNumAddTo (newb, *x);
00821       rval = dbl_ILLfct_bound_shift (lp, col, BOUND_UPPER, newb);
00822       CHECKRVALG (rval, CLEANUP);
00823       nchg++;
00824     }
00825   }
00826   *chgb = nchg;
00827 
00828 CLEANUP:
00829   dbl_EGlpNumClearVar (newb);
00830   dbl_EGlpNumClearVar (cftol);
00831   EG_RETURN (rval);
00832 }
00833 
00834 static int dbl_expand_phaseI_bounds (
00835   dbl_lpinfo * lp,
00836   int *chgb)
00837 {
00838   int rval = 0;
00839   int i, col, nchg = 0;
00840   double newb, cftol;
00841   double *u, *l, *x;
00842   ILLrandstate r;
00843 
00844   dbl_EGlpNumInitVar (newb);
00845   dbl_EGlpNumInitVar (cftol);
00846   dbl_EGlpNumCopyAbs (cftol, lp->tol->ip_tol);
00847   dbl_EGlpNumDivUiTo (cftol, 10);
00848   ILLutil_sprand (1, &r);
00849 
00850   for (i = 0; i < lp->nrows; i++)
00851   {
00852     col = lp->baz[i];
00853     if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFREE)
00854       continue;
00855     x = &(lp->xbz[i]);
00856     l = &(lp->lz[col]);
00857     u = &(lp->uz[col]);
00858 
00859     if (dbl_EGlpNumIsNeqq (*l, dbl_NINFTY) && dbl_EGlpNumIsEqual (*x, *l, cftol))
00860     {
00861       dbl_EGlpNumSet (newb, dbl_my_rand (50, &(lp->rstate)) + 1.0);
00862       dbl_EGlpNumMultTo (newb, cftol);
00863       dbl_EGlpNumSign (newb);
00864       dbl_EGlpNumAddTo (newb, *l);
00865       rval = dbl_ILLfct_bound_shift (lp, col, BOUND_LOWER, newb);
00866       CHECKRVALG (rval, CLEANUP);
00867       nchg++;
00868     }
00869     if (dbl_EGlpNumIsNeqq (*u, dbl_INFTY) && dbl_EGlpNumIsEqual (*x, *u, cftol))
00870     {
00871       dbl_EGlpNumSet (newb, dbl_my_rand (50, &(lp->rstate)) + 1.0);
00872       dbl_EGlpNumMultTo (newb, cftol);
00873       dbl_EGlpNumAddTo (newb, *u);
00874       rval = dbl_ILLfct_bound_shift (lp, col, BOUND_UPPER, newb);
00875       CHECKRVALG (rval, CLEANUP);
00876       nchg++;
00877     }
00878   }
00879   *chgb = nchg;
00880 
00881 CLEANUP:
00882   dbl_EGlpNumClearVar (newb);
00883   dbl_EGlpNumClearVar (cftol);
00884   EG_RETURN (rval);
00885 }
00886 
00887 int dbl_ILLfct_adjust_viol_bounds (
00888   dbl_lpinfo * lp)
00889 {
00890   int rval = 0;
00891   int chgb = 0;
00892   double tol;
00893 
00894   dbl_EGlpNumInitVar (tol);
00895   dbl_EGlpNumCopyNeg (tol, lp->tol->pfeas_tol);
00896   rval = dbl_expand_var_bounds (lp, tol, &chgb);
00897 #if dbl_FCT_DEBUG > 0
00898   if (rval == 0)
00899     printf ("adjusting %d bounds\n", chgb);
00900 #endif
00901   dbl_EGlpNumClearVar (tol);
00902   EG_RETURN (rval);
00903 }
00904 
00905 int dbl_ILLfct_perturb_bounds (
00906   dbl_lpinfo * lp)
00907 {
00908   int rval = 0;
00909   int chgb = 0;
00910 
00911   rval = dbl_expand_var_bounds (lp, lp->tol->ip_tol, &chgb);
00912 #if dbl_FCT_DEBUG > 0
00913   if (rval == 0)
00914     printf ("perturbing %d bounds\n", chgb);
00915 #endif
00916   EG_RETURN (rval);
00917 }
00918 
00919 int dbl_ILLfct_perturb_phaseI_bounds (
00920   dbl_lpinfo * lp)
00921 {
00922   int rval = 0;
00923   int chgb = 0;
00924 
00925   rval = dbl_expand_phaseI_bounds (lp, &chgb);
00926 #if dbl_FCT_DEBUG > 0
00927   if (rval == 0)
00928     printf ("perturbing %d phase I bounds\n", chgb);
00929 #endif
00930   EG_RETURN (rval);
00931 }
00932 
00933 int dbl_ILLfct_bound_shift (
00934   dbl_lpinfo * lp,
00935   int col,
00936   int bndtype,
00937   double newbnd)
00938 {
00939   int rval = 0;
00940   dbl_bndinfo *nbnd = 0;
00941 
00942   ILL_IFTRACE ("\n%s:%d:%d:%la", __func__, col, bndtype, dbl_EGlpNumToLf (newbnd));
00943   nbnd = dbl_ILLfct_new_bndinfo ();
00944 
00945   nbnd->varnum = col;
00946   nbnd->btype = bndtype;
00947   if (bndtype == BOUND_LOWER)
00948   {
00949     dbl_EGlpNumCopy (nbnd->pbound, lp->lz[col]);
00950     dbl_EGlpNumCopy (nbnd->cbound, newbnd);
00951     dbl_EGlpNumCopy (lp->lz[col], newbnd);
00952   }
00953   else
00954   {
00955     dbl_EGlpNumCopy (nbnd->pbound, lp->uz[col]);
00956     dbl_EGlpNumCopy (nbnd->cbound, newbnd);
00957     dbl_EGlpNumCopy (lp->uz[col], newbnd);
00958   }
00959   ILL_IFTRACE (":%la", dbl_EGlpNumToLf (nbnd->pbound));
00960   if (lp->vtype[col] == VFIXED || lp->vtype[col] == VARTIFICIAL)
00961   {
00962     /* printf ("changing f/a bound\n"); */
00963     if (dbl_EGlpNumIsLess (lp->lz[col], lp->uz[col]))
00964       lp->vtype[col] = VBOUNDED;
00965   }
00966 
00967   nbnd->next = lp->bchanges;
00968   lp->bchanges = nbnd;
00969   lp->nbchange++;
00970 
00971 //CLEANUP:
00972   if (rval)
00973     dbl_ILLfct_free_bndinfo (nbnd);
00974   ILL_IFTRACE ("\n");
00975   EG_RETURN (rval);
00976 }
00977 
00978 void dbl_ILLfct_unroll_bound_change (
00979   dbl_lpinfo * lp)
00980 {
00981   int col;
00982   int changex = 0;
00983   dbl_bndinfo *bptr = lp->bchanges;
00984   dbl_bndinfo *nptr = 0;
00985 
00986   ILL_IFTRACE ("%s:", __func__);
00987 
00988   while (lp->nbchange != 0)
00989   {
00990     col = bptr->varnum;
00991     ILL_IFTRACE (":%d", col);
00992 
00993     if (bptr->btype == BOUND_UPPER)
00994       dbl_EGlpNumCopy (lp->uz[col], bptr->pbound);
00995     else
00996       dbl_EGlpNumCopy (lp->lz[col], bptr->pbound);
00997 
00998     if (lp->vtype[col] == VBOUNDED)
00999     {
01000       if (dbl_EGlpNumIsEqqual (lp->lz[col], lp->uz[col]))
01001         lp->vtype[col] = (!dbl_EGlpNumIsNeqqZero (lp->lz[col])) ?
01002           VARTIFICIAL : VFIXED;
01003     }
01004 
01005     if (lp->vstat[col] != STAT_BASIC)
01006     {
01007       if ((bptr->btype == BOUND_UPPER && lp->vstat[col] == STAT_UPPER) ||
01008           (bptr->btype == BOUND_LOWER && lp->vstat[col] == STAT_LOWER))
01009         changex++;
01010     }
01011     nptr = bptr->next;
01012     dbl_EGlpNumClearVar ((bptr->cbound));
01013     dbl_EGlpNumClearVar ((bptr->pbound));
01014     ILL_IFFREE (bptr, dbl_bndinfo);
01015     bptr = nptr;
01016     lp->nbchange--;
01017   }
01018   lp->bchanges = bptr;
01019   ILL_IFTRACE ("\n");
01020   if (changex)
01021     dbl_ILLfct_compute_xbz (lp);
01022 }
01023 
01024 static int dbl_expand_var_coefs (
01025   dbl_lpinfo * lp,
01026   double ftol,
01027   int *chgc)
01028 {
01029   int rval = 0;
01030   int i, col, vs, vt;
01031   int nchg = 0;
01032   double newc, cftol, mftol[1];
01033   double *c, *dj;
01034   ILLrandstate r;
01035 
01036   dbl_EGlpNumInitVar (newc);
01037   dbl_EGlpNumInitVar (cftol);
01038   dbl_EGlpNumInitVar (mftol[0]);
01039   dbl_EGlpNumCopyAbs (cftol, ftol);
01040   dbl_EGlpNumDivUiTo (cftol, 10);
01041   dbl_EGlpNumCopyNeg (mftol[0], ftol);
01042   ILLutil_sprand (1, &r);
01043 
01044   for (i = 0; i < lp->nnbasic; i++)
01045   {
01046     dj = &(lp->dz[i]);
01047     col = lp->nbaz[i];
01048     c = &(lp->cz[col]);
01049     vs = lp->vstat[col];
01050     vt = lp->vtype[col];
01051 
01052     if (vt == VARTIFICIAL || vt == VFIXED)
01053       continue;
01054     switch (vs)
01055     {
01056     case STAT_ZERO:
01057       dbl_EGlpNumCopyDiff (newc, *c, *dj);
01058       rval = dbl_ILLfct_coef_shift (lp, col, newc);
01059       CHECKRVALG (rval, CLEANUP);
01060       nchg++;
01061       break;
01062     case STAT_LOWER:
01063       if (dbl_EGlpNumIsLess (*dj, ftol))
01064       {
01065         dbl_EGlpNumSet (newc, dbl_my_rand (50, &(lp->rstate)) + 1.0);
01066         dbl_EGlpNumMultTo (newc, cftol);
01067         dbl_EGlpNumAddTo (newc, *c);
01068         if (dbl_EGlpNumIsLessZero (*dj))
01069           dbl_EGlpNumSubTo (newc, *dj);
01070         rval = dbl_ILLfct_coef_shift (lp, col, newc);
01071         CHECKRVALG (rval, CLEANUP);
01072         nchg++;
01073       }
01074       break;
01075     case STAT_UPPER:
01076       if (dbl_EGlpNumIsLess (mftol[0], *dj))
01077       {
01078         dbl_EGlpNumSet (newc, dbl_my_rand (50, &(lp->rstate)) + 1.0);
01079         dbl_EGlpNumMultTo (newc, cftol);
01080         dbl_EGlpNumSign (newc);
01081         dbl_EGlpNumAddTo (newc, *c);
01082         if (dbl_EGlpNumIsGreatZero (*dj))
01083           dbl_EGlpNumSubTo (newc, *dj);
01084         rval = dbl_ILLfct_coef_shift (lp, col, newc);
01085         CHECKRVALG (rval, CLEANUP);
01086         nchg++;
01087       }
01088       break;
01089     default:
01090       break;
01091     }
01092   }
01093   *chgc = nchg;
01094 
01095 CLEANUP:
01096   dbl_EGlpNumClearVar (mftol[0]);
01097   dbl_EGlpNumClearVar (newc);
01098   dbl_EGlpNumClearVar (cftol);
01099   EG_RETURN (rval);
01100 }
01101 
01102 int dbl_ILLfct_adjust_viol_coefs (
01103   dbl_lpinfo * lp)
01104 {
01105   int rval = 0;
01106   int chgc = 0;
01107   double tol;
01108 
01109   dbl_EGlpNumInitVar (tol);
01110   dbl_EGlpNumCopyNeg (tol, lp->tol->dfeas_tol);
01111 
01112   rval = dbl_expand_var_coefs (lp, tol, &chgc);
01113 #if dbl_FCT_DEBUG > 0
01114   if (rval == 0)
01115     printf ("perturbing %d coefs\n", chgc);
01116 #endif
01117   dbl_EGlpNumClearVar (tol);
01118   EG_RETURN (rval);
01119 }
01120 
01121 int dbl_ILLfct_perturb_coefs (
01122   dbl_lpinfo * lp)
01123 {
01124   int rval = 0;
01125   int chgc = 0;
01126 
01127   rval = dbl_expand_var_coefs (lp, lp->tol->id_tol, &chgc);
01128 #if dbl_FCT_DEBUG > 0
01129   if (rval == 0)
01130     printf ("perturbing %d coefs\n", chgc);
01131 #endif
01132   EG_RETURN (rval);
01133 }
01134 
01135 int dbl_ILLfct_coef_shift (
01136   dbl_lpinfo * lp,
01137   int col,
01138   double newcoef)
01139 {
01140   int rval = 0;
01141   dbl_coefinfo *ncoef = 0;
01142 
01143   ILL_SAFE_MALLOC (ncoef, 1, dbl_coefinfo);
01144   dbl_EGlpNumInitVar ((ncoef->pcoef));
01145   dbl_EGlpNumInitVar ((ncoef->ccoef));
01146 
01147   ncoef->varnum = col;
01148   dbl_EGlpNumCopy (ncoef->pcoef, lp->cz[col]);
01149   dbl_EGlpNumCopy (ncoef->ccoef, newcoef);
01150   dbl_EGlpNumCopy (lp->cz[col], newcoef);
01151   ncoef->next = lp->cchanges;
01152   lp->cchanges = ncoef;
01153   dbl_EGlpNumAddTo (lp->dz[lp->vindex[col]], ncoef->ccoef);
01154   dbl_EGlpNumSubTo (lp->dz[lp->vindex[col]], ncoef->pcoef);
01155   lp->ncchange++;
01156 
01157 CLEANUP:
01158   if (rval)
01159   {
01160     dbl_EGlpNumClearVar ((ncoef->pcoef));
01161     dbl_EGlpNumClearVar ((ncoef->ccoef));
01162     ILL_IFFREE (ncoef, dbl_coefinfo);
01163   }
01164   EG_RETURN (rval);
01165 }
01166 
01167 void dbl_ILLfct_unroll_coef_change (
01168   dbl_lpinfo * lp)
01169 {
01170   int bascoef = 0;
01171   dbl_coefinfo *cptr = (dbl_coefinfo *) lp->cchanges;
01172   dbl_coefinfo *nptr = 0;
01173 
01174   while (lp->ncchange != 0)
01175   {
01176     dbl_EGlpNumCopy (lp->cz[cptr->varnum], cptr->pcoef);
01177     if (lp->vstat[cptr->varnum] != STAT_BASIC)
01178     {
01179       dbl_EGlpNumAddTo (lp->dz[lp->vindex[cptr->varnum]], cptr->pcoef);
01180       dbl_EGlpNumSubTo (lp->dz[lp->vindex[cptr->varnum]], cptr->ccoef);
01181     }
01182     else
01183       bascoef++;
01184 
01185     nptr = cptr->next;
01186     dbl_EGlpNumClearVar ((cptr->pcoef));
01187     dbl_EGlpNumClearVar ((cptr->ccoef));
01188     ILL_IFFREE (cptr, dbl_coefinfo);
01189     cptr = nptr;
01190     lp->ncchange--;
01191   }
01192   lp->cchanges = cptr;
01193   if (bascoef)
01194   {
01195     dbl_ILLfct_compute_piz (lp);
01196     dbl_ILLfct_compute_dz (lp);
01197   }
01198 }
01199 
01200 /* feasibility routines */
01201 void dbl_ILLfct_check_pfeasible (
01202   dbl_lpinfo * lp,
01203   dbl_feas_info * fs,
01204   const double ftol)
01205 {
01206   int i, col;
01207   double infeas, err1, err2;
01208 
01209   dbl_EGlpNumInitVar (infeas);
01210   dbl_EGlpNumInitVar (err1);
01211   dbl_EGlpNumInitVar (err2);
01212   dbl_EGlpNumZero (infeas);
01213   fs->pstatus = PRIMAL_FEASIBLE;
01214   dbl_EGlpNumZero (fs->totinfeas);
01215   ILL_IFTRACE ("%s:tol %la\n", __func__, dbl_EGlpNumToLf (ftol));
01216 
01217   for (i = 0; i < lp->nrows; i++)
01218   {
01219     col = lp->baz[i];
01220     dbl_EGlpNumCopyDiff (err1, lp->xbz[i], lp->uz[col]);
01221     dbl_EGlpNumCopyDiff (err2, lp->lz[col], lp->xbz[i]);
01222     if (dbl_EGlpNumIsLess (ftol, err1)
01223         && dbl_EGlpNumIsNeq (lp->uz[col], dbl_INFTY, dbl_oneLpNum))
01224     {
01225       dbl_EGlpNumAddTo (infeas, err1);
01226       WARNINGL (QSE_WLVL, dbl_EGlpNumIsLess (dbl_INFTY, err1),
01227                "This is imposible lu = %15lg xbz = %15lg" " dbl_INFTY = %15lg",
01228                dbl_EGlpNumToLf (lp->uz[col]), dbl_EGlpNumToLf (lp->xbz[i]),
01229                dbl_EGlpNumToLf (dbl_INFTY));
01230       lp->bfeas[i] = 1;
01231     }
01232     else if (dbl_EGlpNumIsLess (ftol, err2)
01233              && dbl_EGlpNumIsNeq (lp->lz[col], dbl_NINFTY, dbl_oneLpNum))
01234     {
01235       dbl_EGlpNumAddTo (infeas, err2);
01236       WARNINGL (QSE_WLVL, dbl_EGlpNumIsLess (dbl_INFTY, err2),
01237                "This is imposible lz = %15lg xbz = %15lg" " dbl_NINFTY = %15lg",
01238                dbl_EGlpNumToLf (lp->lz[col]), dbl_EGlpNumToLf (lp->xbz[i]),
01239                dbl_EGlpNumToLf (dbl_NINFTY));
01240       lp->bfeas[i] = -1;
01241     }
01242     else
01243       lp->bfeas[i] = 0;
01244   }
01245   if (dbl_EGlpNumIsNeqqZero (infeas))
01246   {
01247     fs->pstatus = PRIMAL_INFEASIBLE;
01248     dbl_EGlpNumCopy (fs->totinfeas, infeas);
01249     ILL_IFTRACE ("%s:inf %la\n", __func__, dbl_EGlpNumToLf (infeas));
01250     if (dbl_EGlpNumIsLessZero (fs->totinfeas))
01251     {
01252       printf ("Negative infeasibility, Imposible! %lf %la\n",
01253               dbl_EGlpNumToLf (infeas), dbl_EGlpNumToLf (infeas));
01254     }
01255   }
01256   dbl_EGlpNumCopy (lp->pinfeas, infeas);
01257   dbl_EGlpNumClearVar (infeas);
01258   dbl_EGlpNumClearVar (err1);
01259   dbl_EGlpNumClearVar (err2);
01260 }
01261 
01262 /* feasibility routines */
01263 void dbl_ILLfct_check_pIpfeasible (
01264   dbl_lpinfo * lp,
01265   dbl_feas_info * fs,
01266   double ftol)
01267 {
01268   int i, col;
01269   int ninf = 0;
01270 
01271   fs->pstatus = PRIMAL_FEASIBLE;
01272   dbl_EGlpNumZero (fs->totinfeas);
01273 
01274   for (i = 0; i < lp->nrows; i++)
01275   {
01276     if (!dbl_EGlpNumIsNeqZero (lp->xbz[i], ftol))
01277       continue;
01278     col = lp->baz[i];
01279     if (dbl_EGlpNumIsGreatZero(lp->xbz[i]) &&
01280         dbl_EGlpNumIsNeqq (lp->uz[col], dbl_INFTY))
01281     {
01282       ninf++;
01283     }
01284     else if (dbl_EGlpNumIsLessZero (lp->xbz[i]) &&
01285              dbl_EGlpNumIsNeqq (lp->lz[col], dbl_NINFTY))
01286     {
01287       ninf++;
01288     }
01289   }
01290   if (ninf != 0)
01291     fs->pstatus = PRIMAL_INFEASIBLE;
01292 }
01293 
01294 void dbl_ILLfct_check_dfeasible (
01295   dbl_lpinfo * lp,
01296   dbl_feas_info * fs,
01297   const double ftol)
01298 {
01299   int j, col;
01300   double infeas;
01301 
01302   dbl_EGlpNumInitVar (infeas);
01303   dbl_EGlpNumZero (infeas);
01304   fs->dstatus = DUAL_FEASIBLE;
01305   dbl_EGlpNumZero (fs->totinfeas);
01306 
01307   for (j = 0; j < lp->nnbasic; j++)
01308   {
01309     lp->dfeas[j] = 0;
01310     if (!dbl_EGlpNumIsNeqZero (lp->dz[j], ftol))
01311       continue;
01312     col = lp->nbaz[j];
01313 
01314     if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFIXED)
01315       continue;
01316 
01317     if (dbl_EGlpNumIsLessZero (lp->dz[j]) &&
01318         (lp->vstat[col] == STAT_LOWER || lp->vstat[col] == STAT_ZERO))
01319     {
01320       dbl_EGlpNumSubTo (infeas, lp->dz[j]);
01321       lp->dfeas[j] = -1;
01322     }
01323     else if (dbl_EGlpNumIsGreatZero (lp->dz[j]) &&
01324              (lp->vstat[col] == STAT_UPPER || lp->vstat[col] == STAT_ZERO))
01325     {
01326       dbl_EGlpNumAddTo (infeas, lp->dz[j]);
01327       lp->dfeas[j] = 1;
01328     }
01329   }
01330 
01331   if (dbl_EGlpNumIsNeqqZero (infeas))
01332   {
01333     dbl_EGlpNumCopy (fs->totinfeas, infeas);
01334     fs->dstatus = DUAL_INFEASIBLE;
01335     ILL_IFTRACE ("%s:inf %la\n", __func__, dbl_EGlpNumToLf (infeas));
01336     if (dbl_EGlpNumIsLessZero (fs->totinfeas))
01337     {
01338       printf ("Negative infeasibility, Imposible! %lf %la\n",
01339               dbl_EGlpNumToLf (infeas), dbl_EGlpNumToLf (infeas));
01340     }
01341   }
01342   dbl_EGlpNumCopy (lp->dinfeas, infeas);
01343   dbl_EGlpNumClearVar (infeas);
01344 }
01345 
01346 void dbl_ILLfct_check_pIdfeasible (
01347   dbl_lpinfo * lp,
01348   dbl_feas_info * fs,
01349   double ftol)
01350 {
01351   int j, col;
01352   int ninf = 0;
01353   double *dz = lp->pIdz;
01354 
01355   fs->dstatus = DUAL_FEASIBLE;
01356 
01357   for (j = 0; j < lp->nnbasic; j++)
01358   {
01359     if (!dbl_EGlpNumIsNeqZero (dz[j], ftol))
01360       continue;
01361     col = lp->nbaz[j];
01362     if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFIXED)
01363       continue;
01364 
01365     if (dbl_EGlpNumIsLessZero (dz[j]) &&
01366         (lp->vstat[col] == STAT_LOWER || lp->vstat[col] == STAT_ZERO))
01367       ninf++;
01368     else if (dbl_EGlpNumIsGreatZero (dz[j]) &&
01369              (lp->vstat[col] == STAT_UPPER || lp->vstat[col] == STAT_ZERO))
01370       ninf++;
01371   }
01372 
01373   if (ninf != 0)
01374     fs->dstatus = DUAL_INFEASIBLE;
01375 }
01376 
01377 void dbl_ILLfct_dual_adjust (
01378   dbl_lpinfo * lp,
01379   const double ftol)
01380 {
01381   int j, col;
01382 
01383   for (j = 0; j < lp->nnbasic; j++)
01384   {
01385     if (!dbl_EGlpNumIsNeqZero (lp->dz[j], ftol))
01386       continue;
01387     col = lp->nbaz[j];
01388     if (dbl_EGlpNumIsLessZero (lp->dz[j]) &&
01389         dbl_EGlpNumIsNeqq (lp->uz[col], dbl_INFTY))
01390       lp->vstat[col] = STAT_UPPER;
01391     else if (dbl_EGlpNumIsGreatZero (lp->dz[j]) &&
01392              dbl_EGlpNumIsNeqq (lp->lz[col], dbl_NINFTY))
01393       lp->vstat[col] = STAT_LOWER;
01394   }
01395 }
01396 
01397 void dbl_ILLfct_dphaseI_simple_update (
01398   dbl_lpinfo * lp,
01399   double ftol)
01400 {
01401   int j, col;
01402 
01403   for (j = 0; j < lp->nnbasic; j++)
01404   {
01405     if (!dbl_EGlpNumIsNeqZero (lp->dz[j], ftol))
01406       continue;
01407     col = lp->nbaz[j];
01408     if (dbl_EGlpNumIsLessZero (lp->dz[j] ) && lp->vtype[col] == VBOUNDED)
01409       lp->vstat[col] = STAT_UPPER;
01410     else if (dbl_EGlpNumIsGreatZero (lp->dz[j]) && lp->vtype[col] == VBOUNDED)
01411       lp->vstat[col] = STAT_LOWER;
01412   }
01413 }
01414 
01415 /* set status values */
01416 void dbl_ILLfct_set_status_values (
01417   dbl_lpinfo * lp,
01418   int pstatus,
01419   int dstatus,
01420   int ptype,
01421   int dtype)
01422 {
01423   if (dstatus == DUAL_FEASIBLE && dtype == PHASEII)
01424   {
01425     if (!lp->ncchange)
01426     {
01427       lp->probstat.dual_feasible = 1;
01428       lp->basisstat.dual_feasible = 1;
01429       lp->basisstat.dual_infeasible = 0;
01430     }
01431   }
01432   if (dstatus == DUAL_INFEASIBLE && dtype == PHASEII)
01433   {
01434     if (!lp->ncchange)
01435     {
01436       lp->basisstat.dual_feasible = 0;
01437       lp->basisstat.dual_infeasible = 1;
01438     }
01439     if (pstatus == PRIMAL_FEASIBLE && ptype == PHASEI)
01440       if (!lp->ncchange)
01441         lp->probstat.dual_infeasible = 1;
01442   }
01443   if (pstatus == PRIMAL_FEASIBLE && ptype == PHASEII)
01444   {
01445     if (!lp->nbchange)
01446     {
01447       lp->probstat.primal_feasible = 1;
01448       lp->basisstat.primal_feasible = 1;
01449       lp->basisstat.primal_infeasible = 0;
01450     }
01451   }
01452   if (pstatus == PRIMAL_INFEASIBLE && ptype == PHASEII)
01453   {
01454     lp->basisstat.primal_feasible = 0;
01455     lp->basisstat.primal_infeasible = 1;
01456 
01457     if (dstatus == DUAL_FEASIBLE && dtype == PHASEI)
01458       lp->probstat.primal_infeasible = 1;
01459   }
01460   if (pstatus == PRIMAL_UNBOUNDED)
01461   {
01462     if (!lp->nbchange)
01463     {
01464       lp->probstat.primal_unbounded = 1;
01465       lp->basisstat.primal_unbounded = 1;
01466       lp->probstat.dual_infeasible = 1;
01467       lp->basisstat.dual_infeasible = 1;
01468       lp->basisstat.dual_feasible = 0;
01469     }
01470   }
01471   if (dstatus == DUAL_UNBOUNDED)
01472   {
01473     if (!lp->ncchange)
01474     {
01475       lp->probstat.dual_unbounded = 1;
01476       lp->basisstat.dual_unbounded = 1;
01477       lp->probstat.primal_infeasible = 1;
01478       lp->basisstat.primal_infeasible = 1;
01479       lp->basisstat.primal_feasible = 0;
01480     }
01481   }
01482   if (lp->probstat.primal_feasible && lp->probstat.dual_feasible)
01483     lp->probstat.optimal = 1;
01484 
01485   if (lp->basisstat.primal_feasible && lp->basisstat.dual_feasible)
01486     lp->basisstat.optimal = 1;
01487   else
01488     lp->basisstat.optimal = 0;
01489 }
01490 
01491 void dbl_ILLfct_init_counts (
01492   dbl_lpinfo * lp)
01493 {
01494   int i;
01495   dbl_count_struct *c = lp->cnts;
01496 
01497 #define dbl_C_VALUE(a) (1.0+(double)(a)/(PARAM_HEAP_RATIO*dbl_ILLutil_our_log2(a)))
01498   dbl_EGlpNumSet (c->y_ravg, dbl_C_VALUE (lp->nrows));
01499   dbl_EGlpNumSet (c->za_ravg, dbl_C_VALUE (lp->nnbasic));
01500   ILL_IFTRACE ("%s:%la\n", __func__, dbl_EGlpNumToLf (c->za_ravg));
01501 #undef dbl_C_VALUE
01502   c->ynz_cnt = 0;
01503   c->num_y = 0;
01504   c->znz_cnt = 0;
01505   c->num_z = 0;
01506   c->zanz_cnt = 0;
01507   c->num_za = 0;
01508   c->pnorm_cnt = 0;
01509   c->dnorm_cnt = 0;
01510   c->pinz_cnt = 0;
01511   c->num_pi = 0;
01512   c->pi1nz_cnt = 0;
01513   c->num_pi1 = 0;
01514   c->upnz_cnt = 0;
01515   c->num_up = 0;
01516   c->pupv_cnt = 0;
01517   c->dupv_cnt = 0;
01518   c->pI_iter = 0;
01519   c->pII_iter = 0;
01520   c->dI_iter = 0;
01521   c->dII_iter = 0;
01522   c->tot_iter = 0;
01523   for (i = 0; i < 10; i++)
01524   {
01525     c->pivpI[i] = 0;
01526     c->pivpII[i] = 0;
01527     c->pivdI[i] = 0;
01528     c->pivdII[i] = 0;
01529   }
01530 }
01531 
01532 static void dbl_update_piv_values (
01533   dbl_count_struct * c,
01534   int phase,
01535   const double piv2)
01536 {
01537   int i = 0;
01538   double v, piv;
01539 
01540   if (!dbl_EGlpNumIsNeqqZero(piv2))
01541     return;
01542   dbl_EGlpNumInitVar (v);
01543   dbl_EGlpNumInitVar (piv);
01544   dbl_EGlpNumCopyAbs (piv, piv2);
01545   dbl_EGlpNumOne (v);
01546   while (dbl_EGlpNumIsLess (piv, v) && i < 9)
01547   {
01548     dbl_EGlpNumDivUiTo (v, 10);
01549     i++;
01550   }
01551   switch (phase)
01552   {
01553   case PRIMAL_PHASEI:
01554     c->pivpI[i]++;
01555     break;
01556   case PRIMAL_PHASEII:
01557     c->pivpII[i]++;
01558     break;
01559   case DUAL_PHASEI:
01560     c->pivdI[i]++;
01561     break;
01562   case DUAL_PHASEII:
01563     c->pivdII[i]++;
01564     break;
01565   default:
01566     break;
01567   }
01568   dbl_EGlpNumClearVar (v);
01569   dbl_EGlpNumClearVar (piv);
01570 }
01571 
01572 void dbl_ILLfct_update_counts (
01573   dbl_lpinfo * lp,
01574   int f,
01575   int upi,
01576   const double upd)
01577 {
01578   dbl_count_struct *c = lp->cnts;
01579 
01580   switch (f)
01581   {
01582   case CNT_PPHASE1ITER:
01583     c->pI_iter++;
01584     c->tot_iter++;
01585     break;
01586   case CNT_PPHASE2ITER:
01587     c->pII_iter++;
01588     c->tot_iter++;
01589     break;
01590   case CNT_DPHASE1ITER:
01591     c->dI_iter++;
01592     c->tot_iter++;
01593     break;
01594   case CNT_DPHASE2ITER:
01595     c->dII_iter++;
01596     c->tot_iter++;
01597     break;
01598   case CNT_YNZ:
01599     c->ynz_cnt += upi;
01600     c->num_y++;
01601     break;
01602   case CNT_ZANZ:
01603     c->zanz_cnt += upi;
01604     c->num_za++;
01605     break;
01606   case CNT_PINZ:
01607     c->pinz_cnt += upi;
01608     c->num_pi++;
01609     break;
01610   case CNT_P1PINZ:
01611     c->pi1nz_cnt += upi;
01612     c->num_pi1++;
01613     break;
01614   case CNT_UPNZ:
01615     c->upnz_cnt += upi;
01616     c->num_up++;
01617     break;
01618   case CNT_PIPIV:
01619     dbl_update_piv_values (c, PRIMAL_PHASEI, upd);
01620     break;
01621   case CNT_PIIPIV:
01622     dbl_update_piv_values (c, PRIMAL_PHASEII, upd);
01623     break;
01624   case CNT_DIPIV:
01625     dbl_update_piv_values (c, DUAL_PHASEI, upd);
01626     break;
01627   case CNT_DIIPIV:
01628     dbl_update_piv_values (c, DUAL_PHASEII, upd);
01629     break;
01630   case CNT_YRAVG:
01631     dbl_EGlpNumMultUiTo (c->y_ravg, c->tot_iter);
01632     dbl_EGlpNumAddUiTo (c->y_ravg, upi);
01633     dbl_EGlpNumDivUiTo (c->y_ravg, c->tot_iter + 1);
01634     break;
01635   case CNT_ZARAVG:
01636     ILL_IFTRACE ("%s:%d:%d:%d:%la:%la", __func__, f, c->tot_iter, upi,
01637                  dbl_EGlpNumToLf (upd), dbl_EGlpNumToLf (c->za_ravg));
01638     dbl_EGlpNumMultUiTo (c->za_ravg, c->tot_iter);
01639     dbl_EGlpNumAddUiTo (c->za_ravg, upi);
01640     dbl_EGlpNumDivUiTo (c->za_ravg, c->tot_iter + 1);
01641     ILL_IFTRACE (":%la\n", dbl_EGlpNumToLf (c->za_ravg));
01642     break;
01643   }
01644 }
01645 
01646 void dbl_ILLfct_print_counts (
01647   dbl_lpinfo * lp)
01648 {
01649   int i, niter;
01650   dbl_count_struct *c = lp->cnts;
01651 
01652   c->tot_iter = c->pI_iter + c->pII_iter + c->dI_iter + c->dII_iter;
01653   niter = (c->tot_iter == 0) ? 1 : c->tot_iter;
01654   printf ("Counts for problem %s\n", lp->O->probname);
01655   if (c->num_y != 0)
01656     printf ("avg ynz = %.2f\n", (double) c->ynz_cnt / c->num_y);
01657   if (c->num_z != 0)
01658     printf ("avg znz = %.2f\n", (double) c->znz_cnt / c->num_z);
01659   if (c->num_za != 0)
01660     printf ("avg zanz = %.2f\n", (double) c->zanz_cnt / c->num_za);
01661   printf ("avg pnorm = %.2f\n", (double) c->pnorm_cnt / lp->nnbasic);
01662   printf ("avg dnorm = %.2f\n", (double) c->dnorm_cnt / lp->nrows);
01663   if (c->num_pi != 0)
01664     printf ("avg pinz = %.2f\n", (double) c->pinz_cnt / c->num_pi);
01665   if (c->num_pi1 != 0)
01666     printf ("avg piInz = %.2f\n", (double) c->pi1nz_cnt / c->num_pi1);
01667   if (c->num_up != 0)
01668     printf ("avg upnz = %.2f\n", (double) c->upnz_cnt / c->num_up);
01669 
01670   for (i = 0; i < 10; i++)
01671     printf ("piv 1.0e-%d : %d %d %d %d\n",
01672             i, c->pivpI[i], c->pivpII[i], c->pivdI[i], c->pivdII[i]);
01673 }
01674 
01675 
01676 /* c <- a + t*b */
01677 static void dbl_add_vectors (
01678   dbl_lpinfo * lp,
01679   dbl_svector * a,
01680   dbl_svector * b,
01681   dbl_svector * c,
01682   const double t)
01683 {
01684   int i, r, l;
01685   dbl_svector *w = &(lp->work);
01686 
01687   for (i = 0; i < b->nzcnt; i++)
01688   {
01689     r = b->indx[i];
01690     w->indx[i] = r;
01691     dbl_EGlpNumCopy (w->coef[r], t);
01692     dbl_EGlpNumMultTo (w->coef[r], b->coef[i]);
01693     lp->iwork[r] = 1;
01694   }
01695   l = b->nzcnt;
01696 
01697   for (i = 0; i < a->nzcnt; i++)
01698   {
01699     r = a->indx[i];
01700     if (lp->iwork[r] == 0)
01701       w->indx[l++] = r;
01702     dbl_EGlpNumAddTo (w->coef[r], a->coef[i]);
01703   }
01704   for (i = 0; i < l; i++)
01705   {
01706     r = w->indx[i];
01707     c->indx[i] = r;
01708     dbl_EGlpNumCopy (c->coef[i], w->coef[r]);
01709     dbl_EGlpNumZero (w->coef[r]);
01710     lp->iwork[r] = 0;
01711   }
01712   w->nzcnt = 0;
01713   c->nzcnt = l;
01714 }
01715 
01716 void dbl_ILLfct_update_pfeas (
01717   dbl_lpinfo * lp,
01718   int lindex,
01719   dbl_svector * srhs)
01720 {
01721   int i, k, r;
01722   int col, nz = 0;
01723   int cbnd, f;
01724   int *perm = lp->upd.perm;
01725   int *ix = lp->upd.ix;
01726   int tctr = lp->upd.tctr;
01727   double *t = lp->upd.t;
01728   double tz, *dty, ntmp;
01729   double *l, *x, *u, *pftol = &(lp->tol->ip_tol);
01730 
01731   dbl_EGlpNumInitVar (tz);
01732   dbl_EGlpNumInitVar (ntmp);
01733   dty = &(lp->upd.dty);
01734   dbl_EGlpNumZero (*dty);
01735   dbl_EGlpNumCopyAbs (tz, lp->upd.tz);
01736   dbl_EGlpNumDivUiTo (tz, 100);
01737   dbl_EGlpNumAddTo (tz, lp->upd.tz);
01738   ILL_IFTRACE ("%s:%d", __func__, tctr);
01739   for (i = 0; i < tctr && dbl_EGlpNumIsLeq (t[perm[i]], tz); i++)
01740   {
01741     cbnd = ix[perm[i]] % 10;
01742     ILL_IFTRACE (":%d", cbnd);
01743     if (cbnd == BBOUND)
01744       continue;
01745     k = ix[perm[i]] / 10;
01746     r = lp->yjz.indx[k];
01747     ILL_IFTRACE (":%d:%d:%d", k, r, lp->iwork[r]);
01748 
01749     if (lp->iwork[r] != 1)
01750     {
01751       lp->iwork[r] = 1;
01752       x = &(lp->xbz[r]);
01753       col = lp->baz[r];
01754       l = &(lp->lz[col]);
01755       u = &(lp->uz[col]);
01756 
01757       if (r != lindex)
01758       {
01759         f = 0;
01760         dbl_EGlpNumCopyDiff (ntmp, *l, *x);
01761         if (dbl_EGlpNumIsNeqq (*l, dbl_NINFTY) && dbl_EGlpNumIsLess (*pftol, ntmp))
01762           f = -1;
01763         else
01764         {
01765           dbl_EGlpNumCopyDiff (ntmp, *x, *u);
01766           if (dbl_EGlpNumIsNeqq (*u, dbl_INFTY) && dbl_EGlpNumIsLess (*pftol, ntmp))
01767             f = 1;
01768         }
01769 
01770         ILL_IFTRACE (":%d:%d", f, lp->bfeas[r]);
01771         if (f != lp->bfeas[r])
01772         {
01773           srhs->indx[nz] = r;
01774           dbl_EGlpNumSet (srhs->coef[nz], (double) (f - lp->bfeas[r]));
01775           dbl_EGlpNumAddInnProdTo (*dty, srhs->coef[nz], lp->yjz.coef[k]);
01776           nz++;
01777           lp->bfeas[r] = f;
01778         }
01779       }
01780       else
01781       {
01782         lp->bfeas[r] = 0;
01783       }
01784     }
01785   }
01786   while (--i >= 0)
01787   {
01788     cbnd = ix[perm[i]] % 10;
01789     if (cbnd == BBOUND)
01790       continue;
01791     k = ix[perm[i]] / 10;
01792     r = lp->yjz.indx[k];
01793     lp->iwork[r] = 0;
01794   }
01795   srhs->nzcnt = nz;
01796   ILL_IFTRACE (":%d\n", nz);
01797   dbl_EGlpNumClearVar (tz);
01798   dbl_EGlpNumClearVar (ntmp);
01799 }
01800 
01801 void dbl_ILLfct_compute_ppIzz (
01802   dbl_lpinfo * lp,
01803   dbl_svector * srhs,
01804   dbl_svector * ssoln)
01805 {
01806   if (srhs->nzcnt != 0)
01807   {
01808     ILL_IFTRACE ("%s:\n", __func__);
01809     dbl_ILLbasis_row_solve (lp, srhs, ssoln);
01810   }
01811 }
01812 
01813 void dbl_ILLfct_update_ppI_prices (
01814   dbl_lpinfo * lp,
01815   dbl_price_info * pinf,
01816   dbl_svector * srhs,
01817   dbl_svector * ssoln,
01818   int eindex,
01819   int lindex,
01820   const double alpha)
01821 {
01822   double ntmp;
01823 
01824   dbl_EGlpNumInitVar (ntmp);
01825   dbl_EGlpNumCopy (ntmp, alpha);
01826   ILL_IFTRACE ("%s:\n", __func__);
01827   if (lindex == -1)
01828   {
01829     if (srhs->nzcnt != 0)
01830     {
01831       dbl_ILLfct_update_pIpiz (lp, ssoln, dbl_oneLpNum);
01832       if (pinf->p_strategy == COMPLETE_PRICING)
01833       {
01834         dbl_ILLfct_compute_zA (lp, ssoln, &(lp->zA));
01835         dbl_ILLfct_update_pIdz (lp, &(lp->zA), -1, dbl_oneLpNum);
01836       }
01837     }
01838     else
01839     {
01840       if (pinf->p_strategy == COMPLETE_PRICING)
01841         dbl_ILLprice_compute_dual_inf (lp, pinf, &eindex, 1, PRIMAL_PHASEI);
01842       else
01843         dbl_ILLprice_update_mpartial_price (lp, pinf, PRIMAL_PHASEI, COL_PRICING);
01844       dbl_EGlpNumClearVar (ntmp);
01845       return;
01846     }
01847   }
01848   else
01849   {
01850     if (srhs->nzcnt == 0)
01851     {
01852       dbl_ILLfct_update_pIpiz (lp, &(lp->zz), ntmp);
01853       if (pinf->p_strategy == COMPLETE_PRICING)
01854         dbl_ILLfct_update_pIdz (lp, &(lp->zA), eindex, ntmp);
01855     }
01856     else
01857     {
01858       dbl_EGlpNumCopyFrac (ntmp, lp->upd.dty, lp->upd.piv);
01859       dbl_EGlpNumSubTo (ntmp, alpha);
01860       dbl_EGlpNumSign (ntmp);
01861       dbl_add_vectors (lp, ssoln, &(lp->zz), &(lp->zz), ntmp);
01862       dbl_ILLfct_update_pIpiz (lp, &(lp->zz), dbl_oneLpNum);
01863       if (pinf->p_strategy == COMPLETE_PRICING)
01864       {
01865         dbl_ILLfct_compute_zA (lp, &(lp->zz), &(lp->zA));
01866         dbl_ILLfct_update_pIdz (lp, &(lp->zA), eindex, dbl_oneLpNum);
01867       }
01868     }
01869     dbl_EGlpNumSet (lp->pIdz[eindex], (double) (lp->upd.fs));
01870     dbl_EGlpNumAddTo (lp->pIdz[eindex], ntmp);
01871     dbl_EGlpNumSign (lp->pIdz[eindex]);
01872   }
01873   if (pinf->p_strategy == COMPLETE_PRICING)
01874   {
01875     dbl_ILLprice_compute_dual_inf (lp, pinf, lp->zA.indx, lp->zA.nzcnt,
01876                                PRIMAL_PHASEI);
01877     if (eindex > -1)
01878       dbl_ILLprice_compute_dual_inf (lp, pinf, &eindex, 1, PRIMAL_PHASEI);
01879     dbl_ILLfct_update_counts (lp, CNT_ZARAVG, lp->zA.nzcnt, dbl_zeroLpNum);
01880   }
01881   else
01882     dbl_ILLprice_update_mpartial_price (lp, pinf, PRIMAL_PHASEI, COL_PRICING);
01883   dbl_EGlpNumClearVar (ntmp);
01884   return;
01885 }
01886 
01887 void dbl_ILLfct_update_dfeas (
01888   dbl_lpinfo * lp,
01889   int eindex,
01890   dbl_svector * srhs)
01891 {
01892   int i, j, k, c;
01893   int cbnd, col, nz = 0;
01894   int vs, vt, f;
01895   int delta;
01896   int *perm = lp->upd.perm;
01897   int *ix = lp->upd.ix;
01898   int tctr = lp->upd.tctr;
01899   int mcnt, mbeg;
01900   double *t = lp->upd.t;
01901   double *w = lp->work.coef;
01902   double tz;
01903   double *dty = &(lp->upd.dty);
01904   double *dftol = &(lp->tol->id_tol);
01905   double dj;
01906 
01907   dbl_EGlpNumInitVar (dj);
01908   dbl_EGlpNumInitVar (tz);
01909   dbl_EGlpNumZero (*dty);
01910   dbl_EGlpNumCopy (tz, lp->upd.tz);
01911   dbl_EGlpNumMultUiTo (tz, 101);
01912   dbl_EGlpNumDivUiTo (tz, 100);
01913 
01914   for (j = 0; j < tctr && dbl_EGlpNumIsLeq (t[perm[j]], tz); j++)
01915   {
01916     k = ix[perm[j]] / 10;
01917     c = lp->zA.indx[k];
01918 
01919     if (lp->iwork[c] != 1)
01920     {
01921       lp->iwork[c] = 1;
01922       cbnd = ix[perm[j]] % 10;
01923       col = lp->nbaz[c];
01924       dbl_EGlpNumCopy (dj, lp->dz[c]);
01925       vs = lp->vstat[col];
01926       vt = lp->vtype[col];
01927 
01928       if (cbnd == BSKIP)
01929       {
01930         if (!dbl_EGlpNumIsNeqZero (dj, *dftol));
01931         else if (dbl_EGlpNumIsLessZero (dj) && vs == STAT_LOWER)
01932           lp->vstat[col] = STAT_UPPER;
01933         else if (dbl_EGlpNumIsGreatZero (dj) && vs == STAT_UPPER)
01934           lp->vstat[col] = STAT_LOWER;
01935       }
01936       else if (c != eindex)
01937       {
01938         if (!dbl_EGlpNumIsNeqZero (dj, *dftol))
01939           f = 0;
01940         else if (dbl_EGlpNumIsLessZero (dj) &&
01941                  (vs == STAT_LOWER || vs == STAT_ZERO))
01942           f = -1;
01943         else if (dbl_EGlpNumIsGreatZero (dj) &&
01944                  (vs == STAT_UPPER || vs == STAT_ZERO))
01945           f = 1;
01946         else
01947           f = 0;
01948 
01949         if (f != lp->dfeas[c])
01950         {
01951           delta = f - lp->dfeas[c];
01952           mcnt = lp->matcnt[col];
01953           mbeg = lp->matbeg[col];
01954           dbl_EGlpNumSet (dj, (double) (delta));
01955           for (i = 0; i < mcnt; i++)
01956             dbl_EGlpNumAddInnProdTo (w[lp->matind[mbeg + i]], dj,
01957                                  lp->matval[mbeg + i]);
01958           dbl_EGlpNumAddInnProdTo (*dty, dj, lp->zA.coef[k]);
01959           nz = 1;
01960           lp->dfeas[c] = f;
01961         }
01962       }
01963       else
01964       {
01965         lp->dfeas[c] = 0;
01966       }
01967     }
01968   }
01969   while (--j >= 0)
01970   {
01971     k = ix[perm[j]] / 10;
01972     c = lp->zA.indx[k];
01973     lp->iwork[c] = 0;
01974   }
01975 
01976   if (nz)
01977   {
01978     for (i = 0, nz = 0; i < lp->nrows; i++)
01979       if (dbl_EGlpNumIsNeqqZero (w[i]))
01980       {
01981         dbl_EGlpNumCopy (srhs->coef[nz], w[i]);
01982         srhs->indx[nz] = i;
01983         nz++;
01984         dbl_EGlpNumZero (w[i]);
01985       }
01986   }
01987 
01988   srhs->nzcnt = nz;
01989   dbl_EGlpNumClearVar (dj);
01990   dbl_EGlpNumClearVar (tz);
01991 }
01992 
01993 void dbl_ILLfct_compute_dpIy (
01994   dbl_lpinfo * lp,
01995   dbl_svector * srhs,
01996   dbl_svector * ssoln)
01997 {
01998   if (srhs->nzcnt != 0)
01999   {
02000     dbl_ILLbasis_column_solve (lp, srhs, ssoln);
02001   }
02002 }
02003 
02004 void dbl_ILLfct_update_dpI_prices (
02005   dbl_lpinfo * lp,
02006   dbl_price_info * pinf,
02007   dbl_svector * srhs,
02008   dbl_svector * ssoln,
02009   int lindex,
02010   double alpha)
02011 {
02012   int i;
02013   double ntmp;
02014 
02015   dbl_EGlpNumInitVar (ntmp);
02016   dbl_EGlpNumZero (ntmp);
02017 
02018   if (srhs->nzcnt == 0)
02019   {
02020     dbl_ILLfct_update_xz (lp, alpha, -1, -1);
02021   }
02022   else
02023   {
02024     dbl_EGlpNumCopyFrac (ntmp, lp->upd.dty, lp->upd.piv);
02025     dbl_EGlpNumAddTo (ntmp, alpha);
02026     dbl_EGlpNumSign (ntmp);
02027     dbl_add_vectors (lp, ssoln, &(lp->yjz), &(lp->yjz), ntmp);
02028     dbl_EGlpNumSign (ntmp);
02029     for (i = 0; i < lp->yjz.nzcnt; i++)
02030       dbl_EGlpNumAddTo (lp->xbz[lp->yjz.indx[i]], lp->yjz.coef[i]);
02031   }
02032   dbl_EGlpNumSet (lp->xbz[lindex], ((double) (-lp->upd.fs)));
02033   dbl_EGlpNumAddTo (lp->xbz[lindex], ntmp);
02034 
02035   if (pinf->d_strategy == COMPLETE_PRICING)
02036   {
02037     dbl_ILLprice_compute_primal_inf (lp, pinf, lp->yjz.indx, lp->yjz.nzcnt,
02038                                  DUAL_PHASEI);
02039     dbl_ILLprice_compute_primal_inf (lp, pinf, &lindex, 1, DUAL_PHASEI);
02040     dbl_ILLfct_update_counts (lp, CNT_YRAVG, lp->yjz.nzcnt, dbl_zeroLpNum);
02041   }
02042   else
02043     dbl_ILLprice_update_mpartial_price (lp, pinf, DUAL_PHASEI, ROW_PRICING);
02044   dbl_EGlpNumClearVar (ntmp);
02045 }
02046 
02047 void dbl_ILLfct_update_dIIfeas (
02048   dbl_lpinfo * lp,
02049   int eindex,
02050   dbl_svector * srhs)
02051 {
02052   int j, k;
02053   int col, indx, vs;
02054   int *perm = lp->upd.perm;
02055   int *ix = lp->upd.ix;
02056   int tctr = lp->upd.tctr;
02057   double *zAj, *l, *u;
02058   double *dty = &(lp->upd.dty);
02059   double *t_max = &(lp->upd.tz);
02060   double *t = lp->upd.t;
02061   double delta;
02062   dbl_svector a;
02063 
02064   dbl_EGlpNumInitVar (delta);
02065   dbl_EGlpNumZero (delta);
02066   dbl_EGlpNumZero (*dty);
02067 
02068   srhs->nzcnt = 0;
02069   for (j = 0; j < tctr && dbl_EGlpNumIsLeq (t[perm[j]], *t_max); j++)
02070   {
02071     k = ix[perm[j]];
02072     indx = lp->zA.indx[k];
02073 
02074     if (indx != eindex)
02075     {
02076       zAj = &(lp->zA.coef[k]);
02077       col = lp->nbaz[indx];
02078       l = &(lp->lz[col]);
02079       u = &(lp->uz[col]);
02080       vs = lp->vstat[col];
02081       if (vs == STAT_UPPER)
02082         dbl_EGlpNumCopyDiff (delta, *l, *u);
02083       else
02084         dbl_EGlpNumCopyDiff (delta, *u, *l);
02085       dbl_EGlpNumAddInnProdTo (*dty, delta, *zAj);
02086       lp->vstat[col] = (vs == STAT_UPPER) ? STAT_LOWER : STAT_UPPER;
02087 
02088       a.nzcnt = lp->matcnt[col];
02089       a.indx = &(lp->matind[lp->matbeg[col]]);
02090       a.coef = &(lp->matval[lp->matbeg[col]]);
02091       dbl_add_vectors (lp, srhs, &a, srhs, delta);
02092     }
02093   }
02094   dbl_EGlpNumClearVar (delta);
02095 }
02096 
02097 void dbl_ILLfct_compute_dpIIy (
02098   dbl_lpinfo * lp,
02099   dbl_svector * srhs,
02100   dbl_svector * ssoln)
02101 {
02102   if (srhs->nzcnt != 0)
02103   {
02104     dbl_ILLbasis_column_solve (lp, srhs, ssoln);
02105   }
02106 }
02107 
02108 void dbl_ILLfct_update_dpII_prices (
02109   dbl_lpinfo * lp,
02110   dbl_price_info * pinf,
02111   dbl_svector * srhs,
02112   dbl_svector * ssoln,
02113   /*int eindex,*/
02114   int lindex,
02115   double eval,
02116   double alpha)
02117 {
02118   int i;
02119   dbl_svector *u;
02120 
02121   if (srhs->nzcnt == 0)
02122   {
02123     dbl_ILLfct_update_xz (lp, alpha, -1, -1);
02124     u = &(lp->yjz);
02125   }
02126   else
02127   {
02128     if (ssoln->nzcnt != 0)
02129       for (i = 0; i < ssoln->nzcnt; i++)
02130         dbl_EGlpNumSubTo (lp->xbz[ssoln->indx[i]], ssoln->coef[i]);
02131     dbl_ILLfct_update_xz (lp, alpha, -1, -1);
02132     dbl_add_vectors (lp, ssoln, &(lp->yjz), ssoln, dbl_oneLpNum);
02133     u = ssoln;
02134   }
02135   dbl_EGlpNumCopySum (lp->xbz[lindex], eval, alpha);
02136 
02137   if (pinf->d_strategy == COMPLETE_PRICING)
02138   {
02139     dbl_ILLprice_compute_primal_inf (lp, pinf, u->indx, u->nzcnt, DUAL_PHASEII);
02140     dbl_ILLprice_compute_primal_inf (lp, pinf, &lindex, 1, DUAL_PHASEII);
02141     dbl_ILLfct_update_counts (lp, CNT_YRAVG, u->nzcnt, dbl_zeroLpNum);
02142   }
02143   else
02144     dbl_ILLprice_update_mpartial_price (lp, pinf, DUAL_PHASEII, ROW_PRICING);
02145 }
02146 
02147 int dbl_ILLfct_test_pivot (
02148   dbl_lpinfo * lp,
02149   int indx,
02150   int indxtype,
02151   double piv_val)
02152 {
02153   int i;
02154   double pval, ntmp;
02155 
02156   dbl_EGlpNumInitVar (pval);
02157   dbl_EGlpNumInitVar (ntmp);
02158   dbl_EGlpNumZero (pval);
02159 
02160   if (indxtype == ROW_PIVOT)
02161   {
02162     for (i = 0; i < lp->yjz.nzcnt; i++)
02163       if (lp->yjz.indx[i] == indx)
02164       {
02165         dbl_EGlpNumCopy (pval, lp->yjz.coef[i]);
02166         break;
02167       }
02168   }
02169   else
02170   {
02171     for (i = 0; i < lp->zA.nzcnt; i++)
02172       if (lp->zA.indx[i] == indx)
02173       {
02174         dbl_EGlpNumCopy (pval, lp->zA.coef[i]);
02175         break;
02176       }
02177   }
02178   dbl_EGlpNumCopyDiff (ntmp, pval, piv_val);
02179   dbl_EGlpNumDivTo (ntmp, piv_val);
02180   if (dbl_EGlpNumIsLessZero (ntmp))
02181     dbl_EGlpNumSign (ntmp);
02182   if (dbl_EGlpNumIsLess (dbl_ALTPIV_TOLER, ntmp))
02183   {
02184 #if dbl_FCT_DEBUG > 1
02185     if (indxtype == ROW_PIVOT)
02186       printf ("y_i = %.8f, z_j = %.8f %la %la\n", dbl_EGlpNumToLf (pval),
02187               dbl_EGlpNumToLf (piv_val), dbl_EGlpNumToLf (dbl_ALTPIV_TOLER),
02188               dbl_EGlpNumToLf (ntmp));
02189     else
02190       printf ("z_j = %.8f, y_i = %.8f\n", dbl_EGlpNumToLf (pval),
02191               dbl_EGlpNumToLf (piv_val));
02192 #endif
02193     dbl_EGlpNumClearVar (ntmp);
02194     dbl_EGlpNumClearVar (pval);
02195     return 1;
02196   }
02197   dbl_EGlpNumClearVar (pval);
02198   dbl_EGlpNumClearVar (ntmp);
02199   return 0;
02200 }
02201 
02202 #if dbl_FCT_DEBUG > 0
02203 
02204 void dbl_fct_test_workvector (
02205   dbl_lpinfo * lp)
02206 {
02207   int i, err = 0;
02208 
02209   for (i = 0; i < lp->ncols; i++)
02210   {
02211     if (dbl_EGlpNumIsNeqqZero (lp->work.coef[i]))
02212     {
02213       err++;
02214       dbl_EGlpNumZero (lp->work.coef[i]);
02215     }
02216     if (lp->iwork[i] != 0)
02217     {
02218       err++;
02219       lp->iwork[i] = 0;
02220     }
02221   }
02222   if (err)
02223     printf ("bad work vector, err=%d\n", err);
02224 }
02225 
02226 void dbl_fct_test_pfeasible (
02227   dbl_lpinfo * lp)
02228 {
02229   int i, col;
02230   int err = 0;
02231   double *ftol = &(lp->tol->pfeas_tol);
02232 
02233   for (i = 0; i < lp->nrows; i++)
02234   {
02235     col = lp->baz[i];
02236 
02237     if (dbl_EGlpNumIsNeqq (lp->uz[col], dbl_INFTY)
02238         && dbl_EGlpNumIsSumLess (*ftol, lp->uz[col], lp->xbz[i]))
02239     {
02240       if (lp->bfeas[i] != 1)
02241       {
02242         err++;
02243         lp->bfeas[i] = 1;
02244       }
02245     }
02246     else if (dbl_EGlpNumIsNeqq (lp->lz[col], dbl_NINFTY)
02247              && dbl_EGlpNumIsSumLess (lp->xbz[i], *ftol, lp->lz[col]))
02248     {
02249       if (lp->bfeas[i] != -1)
02250       {
02251         err++;
02252         lp->bfeas[i] = -1;
02253       }
02254     }
02255     /* else if (lp->bfeas[i] != 0) {err++; lp->bfeas[i] = 0;} */
02256   }
02257   if (err != 0)
02258     printf ("test_pfeas err =%d\n", err);
02259 }
02260 
02261 void dbl_fct_test_dfeasible (
02262   dbl_lpinfo * lp)
02263 {
02264   int j, col;
02265   int err = 0;
02266   double *ftol = &(lp->tol->dfeas_tol);
02267   double mftol[1];
02268 
02269   dbl_EGlpNumInitVar (mftol[0]);
02270   dbl_EGlpNumCopyNeg (mftol[0], *ftol);
02271 
02272   for (j = 0; j < lp->nnbasic; j++)
02273   {
02274     col = lp->nbaz[j];
02275 
02276     if (lp->vtype[col] == VARTIFICIAL || lp->vtype[col] == VFIXED)
02277       continue;
02278     if (dbl_EGlpNumIsLess (lp->dz[j], mftol[0]) &&
02279         (lp->vstat[col] == STAT_LOWER || lp->vstat[col] == STAT_ZERO))
02280     {
02281       if (lp->dfeas[j] != -1)
02282       {
02283         err++;
02284         lp->dfeas[j] = -1;
02285       }
02286     }
02287     if (dbl_EGlpNumIsLess (*ftol, lp->dz[j]) &&
02288         (lp->vstat[col] == STAT_UPPER || lp->vstat[col] == STAT_ZERO))
02289     {
02290       if (lp->dfeas[j] != 1)
02291       {
02292         err++;
02293         lp->dfeas[j] = 1;
02294       }
02295     }
02296     /* else if (lp->dfeas[j] != 0) {err++; lp->dfeas[j] = 0;} */
02297   }
02298   if (err != 0)
02299     printf ("test_dfeas err =%d\n", err);
02300 }
02301 
02302 void dbl_fct_test_pI_x (
02303   dbl_lpinfo * lp,
02304   dbl_price_info * p)
02305 {
02306   int i;
02307   int ern = 0;
02308   double *x;
02309   double err, diff;
02310 
02311   dbl_EGlpNumInitVar (err);
02312   dbl_EGlpNumInitVar (diff);
02313   dbl_EGlpNumZero (err);
02314   x = dbl_EGlpNumAllocArray (lp->nrows);
02315 
02316   for (i = 0; i < lp->nrows; i++)
02317     dbl_EGlpNumCopy (x[i], lp->xbz[i]);
02318   dbl_ILLfct_compute_phaseI_xbz (lp);
02319   for (i = 0; i < lp->nrows; i++)
02320   {
02321     dbl_EGlpNumCopyDiff (diff, x[i], lp->xbz[i]);
02322     if (dbl_EGlpNumIsLessZero (diff))
02323       dbl_EGlpNumSign (diff);
02324     if (dbl_EGlpNumIsLess (dbl_PFEAS_TOLER, diff))
02325     {
02326       dbl_EGlpNumAddTo (err, diff);
02327       ern++;
02328       printf ("bad i = %d\n", i);
02329     }
02330   }
02331   if (dbl_EGlpNumIsNeqqZero (err))
02332     printf ("dI x err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02333   dbl_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEI);
02334   dbl_EGlpNumFreeArray (x);
02335   dbl_EGlpNumClearVar (diff);
02336   dbl_EGlpNumClearVar (err);
02337 }
02338 
02339 void dbl_fct_test_pII_x (
02340   dbl_lpinfo * lp,
02341   dbl_price_info * p)
02342 {
02343   int i;
02344   int ern = 0;
02345   double *x;
02346   double err, diff;
02347 
02348   dbl_EGlpNumInitVar (err);
02349   dbl_EGlpNumInitVar (diff);
02350   dbl_EGlpNumZero (err);
02351   x = dbl_EGlpNumAllocArray (lp->nrows);
02352 
02353   for (i = 0; i < lp->nrows; i++)
02354     dbl_EGlpNumCopy (x[i], lp->xbz[i]);
02355   dbl_ILLfct_compute_xbz (lp);
02356   for (i = 0; i < lp->nrows; i++)
02357   {
02358     dbl_EGlpNumCopyDiff (diff, x[i], lp->xbz[i]);
02359     if (dbl_EGlpNumIsLessZero (diff ))
02360       dbl_EGlpNumSign (diff);
02361     if (dbl_EGlpNumIsLess (dbl_PFEAS_TOLER, diff))
02362     {
02363       dbl_EGlpNumAddTo (err, diff);
02364       ern++;
02365       printf ("bad i = %d\n", i);
02366     }
02367   }
02368   if (dbl_EGlpNumIsNeqqZero (err))
02369     printf ("dII x err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02370   dbl_ILLprice_compute_primal_inf (lp, p, NULL, 0, DUAL_PHASEII);
02371   dbl_EGlpNumFreeArray (x);
02372   dbl_EGlpNumClearVar (diff);
02373   dbl_EGlpNumClearVar (err);
02374 }
02375 
02376 void dbl_fct_test_pI_pi_dz (
02377   dbl_lpinfo * lp,
02378   dbl_price_info * p)
02379 {
02380   int i;
02381   int ern = 0;
02382   double *pidz;
02383   double err, diff;
02384 
02385   dbl_EGlpNumInitVar (err);
02386   dbl_EGlpNumInitVar (diff);
02387   pidz = dbl_EGlpNumAllocArray (lp->ncols);
02388   dbl_EGlpNumZero (err);
02389 
02390   for (i = 0; i < lp->nrows; i++)
02391     dbl_EGlpNumCopy (pidz[i], lp->pIpiz[i]);
02392   dbl_ILLfct_compute_phaseI_piz (lp);
02393   for (i = 0; i < lp->nrows; i++)
02394   {
02395     dbl_EGlpNumCopyDiff (diff, pidz[i], lp->pIpiz[i]);
02396     if (dbl_EGlpNumIsLessZero (diff))
02397       dbl_EGlpNumSign (diff);
02398     if (dbl_EGlpNumIsLess (dbl_DFEAS_TOLER, diff))
02399     {
02400       dbl_EGlpNumAddTo (err, diff);
02401       ern++;
02402     }
02403   }
02404   if (dbl_EGlpNumIsNeqqZero (err))
02405     printf ("pI pi err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02406 
02407   dbl_EGlpNumZero (err);
02408   ern = 0;
02409   for (i = 0; i < lp->nnbasic; i++)
02410     dbl_EGlpNumCopy (pidz[i], lp->pIdz[i]);
02411   dbl_ILLfct_compute_phaseI_dz (lp);
02412   for (i = 0; i < lp->nnbasic; i++)
02413   {
02414     dbl_EGlpNumCopyDiff (diff, pidz[i], lp->pIdz[i]);
02415     if (dbl_EGlpNumIsLessZero (diff))
02416       dbl_EGlpNumSign (diff);
02417     if (dbl_EGlpNumIsLess (dbl_DFEAS_TOLER, diff))
02418     {
02419       dbl_EGlpNumAddTo (err, diff);
02420       ern++;
02421     }
02422   }
02423   if (dbl_EGlpNumIsNeqqZero (err))
02424     printf ("pI dz err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02425   dbl_ILLprice_compute_dual_inf (lp, p, NULL, 0, PRIMAL_PHASEI);
02426   dbl_EGlpNumClearVar (err);
02427   dbl_EGlpNumClearVar (diff);
02428   dbl_EGlpNumFreeArray (pidz);
02429 }
02430 
02431 void dbl_fct_test_pII_pi_dz (
02432   dbl_lpinfo * lp,
02433   dbl_price_info * p)
02434 {
02435   int i;
02436   int ern = 0;
02437   double *pidz;
02438   double err, diff;
02439 
02440   dbl_EGlpNumInitVar (err);
02441   dbl_EGlpNumInitVar (diff);
02442   dbl_EGlpNumZero (err);
02443   pidz = dbl_EGlpNumAllocArray (lp->ncols);
02444 
02445   for (i = 0; i < lp->nrows; i++)
02446     dbl_EGlpNumCopy (pidz[i], lp->piz[i]);
02447   dbl_ILLfct_compute_piz (lp);
02448   for (i = 0; i < lp->nrows; i++)
02449   {
02450     dbl_EGlpNumCopyDiff (diff, pidz[i], lp->piz[i]);
02451     if (dbl_EGlpNumIsLessZero (diff))
02452       dbl_EGlpNumSign (diff);
02453     if (dbl_EGlpNumIsLess (dbl_DFEAS_TOLER, diff))
02454     {
02455       dbl_EGlpNumAddTo (err, diff);
02456       ern++;
02457     }
02458   }
02459   if (dbl_EGlpNumIsNeqqZero (err))
02460     printf ("pII pi err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02461 
02462   dbl_EGlpNumZero (err);
02463   ern = 0;
02464   for (i = 0; i < lp->nnbasic; i++)
02465     dbl_EGlpNumCopy (pidz[i], lp->dz[i]);
02466   dbl_ILLfct_compute_dz (lp);
02467   for (i = 0; i < lp->nnbasic; i++)
02468   {
02469     dbl_EGlpNumCopyDiff (diff, pidz[i], lp->dz[i]);
02470     if (dbl_EGlpNumIsLessZero (diff))
02471       dbl_EGlpNumSign (diff);
02472     if (dbl_EGlpNumIsLess (dbl_DFEAS_TOLER, diff))
02473     {
02474       dbl_EGlpNumAddTo (err, diff);
02475       ern++;
02476     }
02477   }
02478   if (dbl_EGlpNumIsNeqqZero (err))
02479     printf ("pII dz err = %.7f, ern = %d\n", dbl_EGlpNumToLf (err), ern);
02480   /*
02481    * dbl_ILLprice_compute_dual_inf (lp, p, NULL, 0, PRIMAL_PHASEII);
02482    */
02483   dbl_EGlpNumClearVar (err);
02484   dbl_EGlpNumClearVar (diff);
02485   dbl_EGlpNumFreeArray (pidz);
02486 }
02487 
02488 #endif

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