dbl_basis.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_basis.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 static int TRACE = 0;
00025 
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <math.h>
00029 #include "qs_config.h"
00030 #include "config.h"
00031 #include "dbl_sortrus.h"
00032 #include "dbl_iqsutil.h"
00033 #include "dbl_lpdefs.h"
00034 #include "dbl_qstruct.h"
00035 #include "dbl_qsopt.h"
00036 #include "dbl_basis.h"
00037 #include "dbl_fct.h"
00038 #include "dbl_lp.h"
00039 #include "dbl_lib.h"
00040 #ifdef USEDMALLOC
00041 #include "dmalloc.h"
00042 #endif
00043 
00044 //#define DJZERO_TOLER dbl_PFEAS_TOLER
00045 #define dbl_BASIS_STATS 0
00046 //#define dbl_BASIS_DEBUG 10
00047 #define dbl_BASIS_DEBUG 0
00048 
00049 void dbl_ILLbasis_init_vardata (
00050   dbl_var_data * vd)
00051 {
00052   memset (vd, 0, sizeof (dbl_var_data));
00053   dbl_EGlpNumInitVar (vd->cmax);
00054 }
00055 
00056 void dbl_ILLbasis_clear_vardata (
00057   dbl_var_data * vd)
00058 {
00059   dbl_EGlpNumClearVar (vd->cmax);
00060   memset (vd, 0, sizeof (dbl_var_data));
00061 }
00062 
00063 static void dbl_get_var_info (
00064   dbl_lpinfo * lp,
00065   dbl_var_data * v);
00066 
00067 static int dbl_init_slack_basis (
00068   dbl_lpinfo * lp,
00069   int *vstat,
00070   int *irow,
00071   int *rrow,
00072   int *unitcol,
00073   int *icol,
00074   int *rcol),
00075   dbl_get_initial_basis1 (
00076   dbl_lpinfo * lp,
00077   int *vstat),
00078   dbl_get_initial_basis2 (
00079   dbl_lpinfo * lp,
00080   int *vstat),
00081   dbl_set_basis_indices (
00082   dbl_lpinfo * lp,
00083   int *vstat),
00084   dbl_choose_basis (
00085   int algorithm,
00086   double pinf1,
00087   double dinf1,
00088   double pinf2,
00089   double dinf2);
00090 
00091 void dbl_ILLbasis_init_basisinfo (
00092   dbl_lpinfo * lp)
00093 {
00094   lp->baz = 0;
00095   lp->nbaz = 0;
00096   lp->vstat = 0;
00097   lp->vindex = 0;
00098   lp->f = 0;
00099 }
00100 
00101 void dbl_ILLbasis_free_basisinfo (
00102   dbl_lpinfo * lp)
00103 {
00104   ILL_IFFREE (lp->baz, int);
00105   ILL_IFFREE (lp->nbaz, int);
00106   ILL_IFFREE (lp->vstat, int);
00107   ILL_IFFREE (lp->vindex, int);
00108 
00109   if (lp->f)
00110   {
00111     dbl_ILLfactor_free_factor_work (lp->f);
00112     dbl_EGlpNumClearVar (lp->f->fzero_tol);
00113     dbl_EGlpNumClearVar (lp->f->szero_tol);
00114     dbl_EGlpNumClearVar (lp->f->partial_tol);
00115     dbl_EGlpNumClearVar (lp->f->maxelem_orig);
00116     dbl_EGlpNumClearVar (lp->f->maxelem_factor);
00117     dbl_EGlpNumClearVar (lp->f->maxelem_cur);
00118     dbl_EGlpNumClearVar (lp->f->partial_cur);
00119     ILL_IFFREE (lp->f, dbl_factor_work);
00120   }
00121 }
00122 
00123 int dbl_ILLbasis_build_basisinfo (
00124   dbl_lpinfo * lp)
00125 {
00126   int rval = 0;
00127 
00128   ILL_SAFE_MALLOC (lp->baz, lp->O->nrows, int);
00129   ILL_SAFE_MALLOC (lp->nbaz, lp->O->ncols, int);
00130   ILL_SAFE_MALLOC (lp->vstat, lp->O->ncols, int);
00131   ILL_SAFE_MALLOC (lp->vindex, lp->O->ncols, int);
00132 
00133   lp->fbasisid = -1;
00134 
00135 CLEANUP:
00136   if (rval)
00137     dbl_ILLbasis_free_basisinfo (lp);
00138   EG_RETURN (rval);
00139 }
00140 
00141 int dbl_ILLbasis_load (
00142   dbl_lpinfo * lp,
00143   dbl_ILLlp_basis * B)
00144 {
00145   int rval = 0;
00146   char *cstat = B->cstat;
00147   char *rstat = B->rstat;
00148   int *structmap = lp->O->structmap;
00149   int *rowmap = lp->O->rowmap;
00150   char *sense = lp->O->sense;
00151   int i, j, ncols = lp->O->ncols, nrows = lp->O->nrows, nstruct = lp->O->nstruct;
00152   int basic = 0, nonbasic = 0;
00153 
00154   dbl_ILLbasis_free_basisinfo (lp);
00155   dbl_ILLbasis_init_basisinfo (lp);
00156   rval = dbl_ILLbasis_build_basisinfo (lp);
00157   CHECKRVALG (rval, CLEANUP);
00158 
00159   for (i = 0; i < nstruct; i++)
00160   {
00161     j = structmap[i];
00162     if (cstat[i] == QS_COL_BSTAT_BASIC)
00163     {
00164       lp->vstat[j] = STAT_BASIC;
00165       lp->baz[basic] = j;
00166       lp->vindex[j] = basic;
00167       basic++;
00168     }
00169     else
00170     {
00171       lp->nbaz[nonbasic] = j;
00172       lp->vindex[j] = nonbasic;
00173       nonbasic++;
00174       switch (cstat[i])
00175       {
00176       case QS_COL_BSTAT_LOWER:
00177         lp->vstat[j] = STAT_LOWER;
00178         break;
00179       case QS_COL_BSTAT_UPPER:
00180         lp->vstat[j] = STAT_UPPER;
00181         break;
00182       case QS_COL_BSTAT_FREE:
00183         lp->vstat[j] = STAT_ZERO;
00184         break;
00185       default:
00186         fprintf (stderr, "unknown col basis stat 1: %c\n", cstat[i]);
00187         rval = 1;
00188         goto CLEANUP;
00189       }
00190     }
00191   }
00192 
00193   for (i = 0; i < nrows; i++)
00194   {
00195     j = rowmap[i];
00196     if (sense[i] == 'R')
00197     {
00198       if (rstat[i] == QS_ROW_BSTAT_BASIC)
00199       {
00200         lp->vstat[j] = STAT_BASIC;
00201         lp->baz[basic] = j;
00202         lp->vindex[j] = basic;
00203         basic++;
00204       }
00205       else
00206       {
00207         lp->nbaz[nonbasic] = j;
00208         lp->vindex[j] = nonbasic;
00209         nonbasic++;
00210         switch (rstat[i])
00211         {
00212         case QS_ROW_BSTAT_LOWER:
00213           lp->vstat[j] = STAT_LOWER;
00214           break;
00215         case QS_ROW_BSTAT_UPPER:
00216           lp->vstat[j] = STAT_UPPER;
00217           break;
00218         default:
00219           fprintf (stderr, "unknown range basis stat 2\n");
00220           rval = 1;
00221           goto CLEANUP;
00222         }
00223       }
00224     }
00225     else
00226     {
00227       switch (rstat[i])
00228       {
00229       case QS_ROW_BSTAT_BASIC:
00230         lp->vstat[j] = STAT_BASIC;
00231         lp->baz[basic] = j;
00232         lp->vindex[j] = basic;
00233         basic++;
00234         break;
00235       case QS_ROW_BSTAT_LOWER:
00236         lp->vstat[j] = STAT_LOWER;
00237         lp->nbaz[nonbasic] = j;
00238         lp->vindex[j] = nonbasic;
00239         nonbasic++;
00240         break;
00241       default:
00242         fprintf (stderr, "unknown row basis stat 3\n");
00243         rval = 1;
00244         goto CLEANUP;
00245       }
00246     }
00247   }
00248 
00249   if (basic + nonbasic != ncols)
00250   {
00251     fprintf (stderr, "error in counts in ILLopt_load_basis\n");
00252     rval = 1;
00253     goto CLEANUP;
00254   }
00255 
00256   if (lp->fbasisid != 0)
00257     lp->basisid = 0;
00258   else
00259     lp->basisid = 1;
00260 
00261 CLEANUP:
00262 
00263   EG_RETURN (rval);
00264 }
00265 
00266 int dbl_ILLbasis_tableau_row (
00267   dbl_lpinfo * lp,
00268   int row,
00269   double * brow,
00270   double * trow,
00271   double * rhs,
00272   int strict)
00273 {
00274   int rval = 0;
00275   int i;
00276   int singular = 0;
00277   int indx;
00278   double coef;
00279   double sum;
00280   dbl_svector z, zA;
00281 
00282   dbl_EGlpNumInitVar (coef);
00283   dbl_EGlpNumInitVar (sum);
00284   dbl_EGlpNumZero (sum);
00285 
00286   dbl_ILLsvector_init (&z);
00287   dbl_ILLsvector_init (&zA);
00288 
00289   if (lp->basisid == -1)
00290   {
00291     fprintf (stderr, "dbl_ILLbasis_tableau_row: no basis\n");
00292     rval = E_GENERAL_ERROR;
00293     ILL_CLEANUP;
00294   }
00295   if (lp->fbasisid != lp->basisid)
00296   {                             /* Needs to be changed */
00297     rval = dbl_ILLbasis_factor (lp, &singular);
00298     CHECKRVALG (rval, CLEANUP);
00299     if (singular)
00300     {
00301       MESSAGE (__QS_SB_VERB, "Singular Basis found!");
00302       rval = E_BASIS_SINGULAR;
00303       ILL_CLEANUP;
00304     }
00305   }
00306   if (brow == NULL)
00307   {
00308     fprintf (stderr, "No array for basis inverse row\n");
00309     rval = E_GENERAL_ERROR;
00310     ILL_CLEANUP;
00311   }
00312 
00313   rval = dbl_ILLsvector_alloc (&z, lp->nrows);
00314   CHECKRVALG (rval, CLEANUP);
00315   dbl_ILLfct_compute_zz (lp, &z, row);
00316 
00317   for (i = 0; i < lp->O->nrows; i++)
00318     dbl_EGlpNumZero (brow[i]);
00319   for (i = 0; i < z.nzcnt; i++)
00320   {
00321     indx = z.indx[i];
00322     dbl_EGlpNumCopy (coef, z.coef[i]);
00323     dbl_EGlpNumCopy (brow[indx], coef);
00324     dbl_EGlpNumAddInnProdTo (sum, coef, lp->bz[indx]);
00325   }
00326 
00327   if (rhs != NULL)
00328     dbl_EGlpNumCopy (*rhs, sum);
00329   if (trow != NULL)
00330   {
00331     if (!strict)
00332     {
00333       rval = dbl_ILLsvector_alloc (&zA, lp->ncols);
00334       if (rval)
00335         ILL_CLEANUP;
00336       ILL_IFTRACE ("%s:\n", __func__);
00337       rval = dbl_ILLfct_compute_zA (lp, &z, &zA);
00338       CHECKRVALG (rval, CLEANUP);
00339 
00340       for (i = 0; i < lp->ncols; i++)
00341         dbl_EGlpNumZero (trow[i]);
00342       for (i = 0; i < zA.nzcnt; i++)
00343         dbl_EGlpNumCopy (trow[lp->nbaz[zA.indx[i]]], zA.coef[i]);
00344       dbl_EGlpNumOne (trow[lp->baz[row]]);
00345     }
00346     else
00347     {
00348       dbl_ILLfct_compute_vA (lp, &z, trow);
00349     }
00350   }
00351 
00352 #if dbl_BASIS_DEBUG > 0
00353   if (rhs != NULL && trow != NULL)
00354   {
00355     double *tr = NULL;
00356 
00357     dbl_EGlpNumZero (sum);
00358     if (strict)
00359       tr = trow;
00360     else
00361     {
00362       tr = dbl_EGlpNumAllocArray (lp->ncols);
00363       dbl_ILLfct_compute_vA (lp, &z, tr);
00364     }
00365     for (i = 0; i < lp->nrows; i++)
00366       if (dbl_EGlpNumIsGreatZero (tr[lp->baz[i]]))
00367         dbl_EGlpNumAddTo (sum, tr[lp->baz[i]]);
00368       else
00369         dbl_EGlpNumSubTo (sum, tr[lp->baz[i]]);
00370     dbl_EGlpNumCopy (coef, dbl_oneLpNum);
00371     dbl_EGlpNumSubTo (coef, sum);
00372     if (dbl_EGlpNumIsLessZero (coef))
00373       dbl_EGlpNumSign (coef);
00374     if (dbl_EGlpNumIsLess (dbl_PIVZ_TOLER, coef))
00375       fprintf (stderr, "tableau: bas computed = %.12f\n", dbl_EGlpNumToLf (sum));
00376     if (!strict)
00377       dbl_EGlpNumFreeArray (tr);
00378 #if dbl_BASIS_DEBUG > 1
00379     dbl_EGlpNumZero (sum);
00380     for (i = 0; i < lp->ncols; i++)
00381     {
00382       if (lp->vstat[i] == STAT_BASIC)
00383         dbl_EGlpNumAddInnProdTo (sum, lp->xbz[lp->vindex[i]], trow[i]);
00384       else if (lp->vstat[i] == STAT_UPPER)
00385         dbl_EGlpNumAddInnProdTo (sum, lp->uz[i], trow[i]);
00386       else if (lp->vstat[i] == STAT_LOWER)
00387         dbl_EGlpNumAddInnProdTo (sum, lp->lz[i], trow[i]);
00388     }
00389     dbl_EGlpNumSet (coef, 1e-10);
00390     if (dbl_EGlpNumIsNeq (sum, *rhs, coef))
00391       fprintf (stderr, "tableau rhs = %.9f, computed = %.9f\n",
00392                dbl_EGlpNumToLf (*rhs), dbl_EGlpNumToLf (sum));
00393 #endif
00394   }
00395 #endif
00396 
00397 CLEANUP:
00398   dbl_ILLsvector_free (&z);
00399   dbl_ILLsvector_free (&zA);
00400   dbl_EGlpNumClearVar (coef);
00401   dbl_EGlpNumClearVar (sum);
00402   return rval;
00403 }
00404 
00405 static void dbl_get_var_info (
00406   dbl_lpinfo * lp,
00407   dbl_var_data * v)
00408 {
00409   int i = 0;
00410 
00411   v->nartif = 0;
00412   v->nslacks = 0;
00413   v->nfree = 0;
00414   v->nbndone = 0;
00415   v->nbounded = 0;
00416   v->nfixed = 0;
00417   dbl_EGlpNumCopy (v->cmax, dbl_NINFTY);
00418 
00419   for (i = 0; i < lp->ncols; i++)
00420   {
00421     switch (lp->vtype[i])
00422     {
00423     case VARTIFICIAL:
00424       v->nartif++;
00425       break;
00426     case VFREE:
00427       v->nfree++;
00428       break;
00429     case VLOWER:
00430     case VUPPER:
00431       if (lp->vclass[i] == CLASS_LOGICAL)
00432         v->nslacks++;
00433       else
00434         v->nbndone++;
00435       break;
00436 
00437     case VFIXED:
00438       v->nfixed++;
00439     case VBOUNDED:
00440       if (lp->vclass[i] == CLASS_LOGICAL)
00441         v->nslacks++;
00442       else
00443         v->nbounded++;
00444       break;
00445     }
00446     dbl_EGlpNumSetToMaxAbs (v->cmax, lp->cz[i]);
00447   }
00448 
00449 #if dbl_BASIS_STATS > 0
00450   printf ("cols = %d, acols = %d, total  = %d, nrows = %d, nlog = %d\n",
00451           lp->ncols, lp->ncols - lp->nrows,
00452           v->nartif + v->nfree + v->nslacks + v->nbndone + v->nbounded,
00453           lp->nrows, v->nartif + v->nslacks);
00454 #endif
00455 }
00456 
00457 static int dbl_init_slack_basis (
00458   dbl_lpinfo * lp,
00459   int *vstat,
00460   int *irow,
00461   int *rrow,
00462   int *unitcol,
00463   int *icol,
00464   int *rcol)
00465 {
00466   int j, r, vt;
00467   int nslacks = 0;
00468 
00469   for (j = 0; j < lp->ncols; j++)
00470   {
00471     r = lp->matind[lp->matbeg[j]];
00472     vt = lp->vtype[j];
00473 
00474     if ((vt == VUPPER || vt == VLOWER || vt == VBOUNDED || vt == VFIXED) &&
00475         lp->vclass[j] == CLASS_LOGICAL)
00476     {
00477 
00478       vstat[j] = STAT_BASIC;
00479       irow[r] = 1;
00480       rrow[r] = 1;
00481       unitcol[r] = j;
00482       if (icol != NULL)
00483       {
00484         icol[j] = 1;
00485         rcol[j] = 1;
00486       }
00487       nslacks++;
00488     }
00489     else if (vt == VARTIFICIAL)
00490     {
00491       unitcol[r] = j;
00492       vstat[j] = STAT_UPPER;
00493     }
00494     else if (vt == VFREE)
00495       vstat[j] = STAT_ZERO;
00496     else if (vt == VFIXED || vt == VUPPER)
00497       vstat[j] = STAT_UPPER;
00498     else if (vt == VLOWER)
00499       vstat[j] = STAT_LOWER;
00500     else if (vt == VBOUNDED)
00501     {
00502       if (fabs (dbl_EGlpNumToLf (lp->lz[j])) < fabs (dbl_EGlpNumToLf (lp->uz[j])))
00503         vstat[j] = STAT_LOWER;
00504       else
00505         vstat[j] = STAT_UPPER;
00506     }
00507   }
00508   return nslacks;
00509 }
00510 
00511 static int dbl_primal_col_select (
00512   dbl_lpinfo * lp,
00513   int *vstat,
00514   int *irow,
00515   int *rrow,
00516   int *unitcol,
00517   double * v,
00518   int *perm,
00519   int *porder,
00520   int nbelem,
00521   int pcols)
00522 {
00523   int i, j, k, tr, r = 0;
00524   int mcnt, mbeg;
00525   int *matbeg = lp->matbeg;
00526   int *matcnt = lp->matcnt;
00527   int *matind = lp->matind;
00528   double *matval = lp->matval;
00529   double alpha, val, maxelem;
00530 
00531   dbl_EGlpNumInitVar (alpha);
00532   dbl_EGlpNumInitVar (val);
00533   dbl_EGlpNumInitVar (maxelem);
00534 
00535   for (k = 0; k < pcols; k++)
00536   {
00537     j = porder[perm[k]];
00538     mcnt = matcnt[j];
00539     mbeg = matbeg[j];
00540 
00541     dbl_EGlpNumCopy (alpha, dbl_NINFTY);
00542     dbl_EGlpNumCopy (maxelem, dbl_NINFTY);
00543 
00544     for (i = 0; i < mcnt; i++)
00545     {
00546       dbl_EGlpNumCopyAbs (val, matval[mbeg + i]);
00547       if (dbl_EGlpNumIsLess (maxelem, val))
00548         dbl_EGlpNumCopy (maxelem, val);
00549       if (rrow[matind[mbeg + i]] == 0 && dbl_EGlpNumIsLess (alpha, val))
00550       {
00551         dbl_EGlpNumCopy (alpha, val);
00552         r = matind[mbeg + i];
00553       }
00554     }
00555     dbl_EGlpNumCopy (val, maxelem);
00556     dbl_EGlpNumMultTo (val, dbl_PARAM_IBASIS_RPIVOT);
00557     if (dbl_EGlpNumIsLess (val, alpha))
00558     {
00559       vstat[j] = STAT_BASIC;
00560       nbelem++;
00561       irow[r] = 1;
00562       dbl_EGlpNumCopy (v[r], alpha);
00563       for (i = 0; i < mcnt; i++)
00564         if (dbl_EGlpNumIsNeqqZero (matval[mbeg + i]))
00565           rrow[matind[mbeg + i]]++;
00566     }
00567     else
00568     {
00569       dbl_EGlpNumCopy (alpha, dbl_NINFTY);
00570       for (i = 0; i < mcnt; i++)
00571       {
00572         tr = matind[mbeg + i];
00573         dbl_EGlpNumCopyAbs (val, matval[mbeg + i]);
00574         dbl_EGlpNumDivTo (val, dbl_PARAM_IBASIS_RTRIANG);
00575         if (dbl_EGlpNumIsNeqq (v[tr], dbl_INFTY) && dbl_EGlpNumIsLess (v[tr], val))
00576         {
00577           dbl_EGlpNumZero (alpha);
00578           break;
00579         }
00580         dbl_EGlpNumCopyAbs (val, matval[mbeg + i]);
00581         if (irow[tr] == 0 && dbl_EGlpNumIsLess (alpha, val))
00582         {
00583           dbl_EGlpNumCopy (alpha, val);
00584           r = tr;
00585         }
00586       }
00587       if (dbl_EGlpNumIsNeqqZero (alpha) && dbl_EGlpNumIsNeqq (alpha, dbl_NINFTY))
00588       {
00589         vstat[j] = STAT_BASIC;
00590         nbelem++;
00591         irow[r] = 1;
00592         dbl_EGlpNumCopy (v[r], alpha);
00593         for (i = 0; i < mcnt; i++)
00594           if (dbl_EGlpNumIsNeqqZero (matval[mbeg + i]))
00595             rrow[matind[mbeg + i]]++;
00596       }
00597     }
00598   }
00599 #if dbl_BASIS_STATS > 0
00600   printf ("nartifs = %d\n", lp->nrows - nbelem);
00601 #endif
00602 
00603   if (nbelem < lp->nrows)
00604   {
00605     for (i = 0; i < lp->nrows; i++)
00606     {
00607       if (irow[i] == 0)
00608       {
00609         if (unitcol[i] != -1)
00610         {
00611           vstat[unitcol[i]] = STAT_BASIC;
00612           nbelem++;
00613         }
00614         else
00615         {
00616           fprintf (stderr, "Error: Not enough artificials\n");
00617           return -1;
00618         }
00619       }
00620     }
00621   }
00622   dbl_EGlpNumClearVar (alpha);
00623   dbl_EGlpNumClearVar (val);
00624   dbl_EGlpNumClearVar (maxelem);
00625   return nbelem;
00626 }
00627 
00628 /* This is an implementation of the initial basis procedure
00629    in: "Implementing the simplex method: the initial basis", by
00630    Bob Bixby.
00631    Goals: choose initial variables to go into basis which satisfy:
00632    1) vars are slacks, 2) vars have freedom to move
00633    3) initial submatrix is nonsingular, 4) low objective function
00634    contribution.
00635 */
00636 static int dbl_get_initial_basis1 (
00637   dbl_lpinfo * lp,
00638   int *vstat)
00639 {
00640   int rval = 0;
00641   int i, j, tot1 = 0, tot2 = 0;
00642   int nbelem = 0, nslacks = 0;
00643   int tfree = 0, tbndone = 0;
00644   int tbounded = 0;
00645   int *irow = NULL, *rrow = NULL;
00646   int *perm = NULL, *porder = NULL;
00647   int *unitcol = NULL;
00648   double cmax;
00649   double *v = NULL;
00650   double *qpenalty = NULL;
00651   dbl_var_data vd;
00652 
00653   dbl_ILLbasis_init_vardata (&vd);
00654   dbl_EGlpNumInitVar (cmax);
00655 
00656   dbl_get_var_info (lp, &vd);
00657   if (!dbl_EGlpNumIsNeqqZero (vd.cmax))
00658     dbl_EGlpNumOne (cmax);
00659   else
00660   {
00661     dbl_EGlpNumCopy (cmax, vd.cmax);
00662     dbl_EGlpNumMultUiTo (cmax, 1000);
00663   }
00664 
00665   ILL_SAFE_MALLOC (irow, lp->nrows, int);
00666   ILL_SAFE_MALLOC (rrow, lp->nrows, int);
00667 
00668   v = dbl_EGlpNumAllocArray (lp->nrows);
00669   ILL_SAFE_MALLOC (unitcol, lp->nrows, int);
00670 
00671   for (i = 0; i < lp->nrows; i++)
00672   {
00673     unitcol[i] = -1;
00674     dbl_EGlpNumCopy (v[i], dbl_INFTY);
00675     irow[i] = 0;
00676     rrow[i] = 0;
00677   }
00678 
00679   nslacks = dbl_init_slack_basis (lp, vstat, irow, rrow, unitcol, NULL, NULL);
00680   if (nslacks != vd.nslacks)
00681   {
00682     printf ("complain: incorrect basis info(slacks)\n");
00683     rval = E_SIMPLEX_ERROR;
00684     ILL_CLEANUP;
00685   }
00686   if (nslacks == lp->nrows)
00687     ILL_CLEANUP;
00688   nbelem = nslacks;
00689   if (nbelem < lp->nrows)
00690   {
00691     for (i = 0; i < lp->nrows; i++)
00692     {
00693       if (irow[i] == 0)
00694       {
00695         if (unitcol[i] != -1)
00696         {
00697           vstat[unitcol[i]] = STAT_BASIC;
00698           nbelem++;
00699         }
00700         else
00701         {
00702           fprintf (stderr, "Error: Not enough artificials\n");
00703           return -1;
00704         }
00705       }
00706     }
00707   }
00708   ILL_CLEANUP;
00709 
00710   tot1 = vd.nfree + vd.nbndone;
00711   tot2 = vd.nfree + vd.nbndone + vd.nbounded;
00712   ILL_SAFE_MALLOC (perm, tot2, int);
00713   ILL_SAFE_MALLOC (porder, tot2, int);
00714 
00715   qpenalty = dbl_EGlpNumAllocArray (tot2);
00716 
00717   for (j = 0; j < lp->ncols; j++)
00718   {
00719     if (vstat[j] == STAT_BASIC)
00720       continue;
00721 
00722     switch (lp->vtype[j])
00723     {
00724     case VFREE:
00725       porder[tfree] = j;
00726       perm[tfree] = tfree;
00727       dbl_EGlpNumCopyFrac (qpenalty[tfree], lp->cz[j], cmax);
00728       tfree++;
00729       break;
00730 
00731     case VLOWER:
00732     case VUPPER:
00733       porder[vd.nfree + tbndone] = j;
00734       perm[vd.nfree + tbndone] = tbndone;
00735       dbl_EGlpNumCopyFrac (qpenalty[vd.nfree + tbndone], lp->cz[j], cmax);
00736       if (lp->vtype[j] == VLOWER)
00737         dbl_EGlpNumAddTo (qpenalty[vd.nfree + tbndone], lp->lz[j]);
00738       else
00739         dbl_EGlpNumSubTo (qpenalty[vd.nfree + tbndone], lp->uz[j]);
00740       tbndone++;
00741       break;
00742 
00743     case VFIXED:
00744     case VBOUNDED:
00745       porder[tot1 + tbounded] = j;
00746       perm[tot1 + tbounded] = tbounded;
00747       dbl_EGlpNumCopyFrac (qpenalty[tot1 + tbndone], lp->cz[j], cmax);
00748       dbl_EGlpNumAddTo (qpenalty[tot1 + tbndone], lp->lz[j]);
00749       dbl_EGlpNumSubTo (qpenalty[tot1 + tbndone], lp->uz[j]);
00750       tbounded++;
00751       break;
00752     }
00753   }
00754   if (tfree != vd.nfree || tbndone != vd.nbndone || tbounded != vd.nbounded)
00755   {
00756     printf ("complain: incorrect basis info \n");
00757     rval = E_SIMPLEX_ERROR;
00758     ILL_CLEANUP;
00759   }
00760 
00761   dbl_ILLutil_EGlpNum_perm_quicksort (perm, qpenalty, vd.nfree);
00762   dbl_ILLutil_EGlpNum_perm_quicksort (perm + vd.nfree, qpenalty + vd.nfree,
00763                                   vd.nbndone);
00764   dbl_ILLutil_EGlpNum_perm_quicksort (perm + tot1, qpenalty + tot1, vd.nbounded);
00765 
00766   for (i = 0; i < vd.nbndone; i++)
00767     perm[vd.nfree + i] += vd.nfree;
00768   for (i = 0; i < vd.nbounded; i++)
00769     perm[tot1 + i] += tot1;
00770 
00771   nbelem =
00772     dbl_primal_col_select (lp, vstat, irow, rrow, unitcol, v, perm, porder, nbelem,
00773                        tot2);
00774   if (nbelem != lp->nrows)
00775   {
00776     printf ("complain: incorrect final basis size\n");
00777     rval = E_SIMPLEX_ERROR;
00778     ILL_CLEANUP;
00779   }
00780 
00781 CLEANUP:
00782   dbl_EGlpNumClearVar (cmax);
00783   if (rval)
00784     dbl_ILLbasis_free_basisinfo (lp);
00785   ILL_IFFREE (irow, int);
00786   ILL_IFFREE (rrow, int);
00787 
00788   dbl_EGlpNumFreeArray (v);
00789   ILL_IFFREE (perm, int);
00790   ILL_IFFREE (porder, int);
00791   ILL_IFFREE (unitcol, int);
00792 
00793   dbl_EGlpNumFreeArray (qpenalty);
00794   dbl_ILLbasis_clear_vardata (&vd);
00795   EG_RETURN (rval);
00796 }
00797 
00798 static int dbl_get_initial_basis2 (
00799   dbl_lpinfo * lp,
00800   int *vstat)
00801 {
00802   int rval = 0;
00803   int i, j, k, tot1, tot2;
00804   int rbeg, rcnt, mcnt;
00805   int nbelem = 0, nslacks = 0;
00806   int tfree = 0, tbndone = 0;
00807   int tbounded = 0;
00808   int *irow = NULL, *rrow = NULL;
00809   int *perm = NULL, *porder = NULL;
00810   int *unitcol = NULL;
00811   double *v = NULL;
00812   double *qpenalty = NULL;
00813   int col = 0, s_i = 0, selc = 0;
00814   int *icol = NULL, *rcol = NULL;
00815   int *plen = NULL;
00816   double *dj = NULL;
00817   dbl_var_data vd;
00818   double seldj;
00819   double selv;
00820   double c_dj;
00821   double cmax;
00822 
00823   dbl_EGlpNumInitVar (seldj);
00824   dbl_EGlpNumInitVar (selv);
00825   dbl_EGlpNumInitVar (c_dj);
00826   dbl_EGlpNumInitVar (cmax);
00827   dbl_EGlpNumZero (c_dj);
00828   dbl_EGlpNumZero (selv);
00829   dbl_EGlpNumZero (seldj);
00830   dbl_ILLbasis_init_vardata (&vd);
00831 
00832   dbl_get_var_info (lp, &vd);
00833 
00834   ILL_SAFE_MALLOC (irow, lp->nrows, int);
00835   ILL_SAFE_MALLOC (rrow, lp->nrows, int);
00836 
00837   v = dbl_EGlpNumAllocArray (lp->nrows);
00838   ILL_SAFE_MALLOC (unitcol, lp->nrows, int);
00839   ILL_SAFE_MALLOC (icol, lp->ncols, int);
00840   ILL_SAFE_MALLOC (rcol, lp->ncols, int);
00841 
00842   dj = dbl_EGlpNumAllocArray (lp->ncols);
00843 
00844   for (i = 0; i < lp->nrows; i++)
00845   {
00846     unitcol[i] = -1;
00847     dbl_EGlpNumCopy (v[i], dbl_INFTY);
00848     irow[i] = 0;
00849     rrow[i] = 0;
00850   }
00851   /* assign all d_j */
00852   for (i = 0; i < lp->ncols; i++)
00853   {
00854     icol[i] = 0;
00855     rcol[i] = 0;
00856     dbl_EGlpNumCopy (dj[i], lp->cz[i]);
00857   }
00858 
00859   nslacks = dbl_init_slack_basis (lp, vstat, irow, rrow, unitcol, icol, rcol);
00860   if (nslacks != vd.nslacks)
00861   {
00862     printf ("complain: incorrect basis info\n");
00863     rval = E_SIMPLEX_ERROR;
00864     ILL_CLEANUP;
00865   }
00866   if (nslacks == lp->nrows)
00867     ILL_CLEANUP;
00868   nbelem = nslacks;
00869 
00870   /* allocate maximum required space for perm etc. */
00871   ILL_SAFE_MALLOC (perm, lp->ncols, int);
00872   ILL_SAFE_MALLOC (porder, lp->ncols, int);
00873   ILL_SAFE_MALLOC (plen, lp->nrows, int);
00874 
00875   qpenalty = dbl_EGlpNumAllocArray (lp->ncols);
00876 
00877   /* find all unit rows and record lengths */
00878   for (i = 0; i < lp->nrows; i++)
00879   {
00880     if (irow[i] != 1)
00881     {
00882       rbeg = lp->rowbeg[i];
00883       rcnt = lp->rowcnt[i];
00884       for (j = 0; j < rcnt; j++)
00885       {
00886         dbl_EGlpNumCopyAbs (cmax, lp->rowval[rbeg + j]);
00887         if (dbl_EGlpNumIsNeqq (cmax, dbl_oneLpNum))
00888           break;
00889       }
00890       if (j == rcnt)
00891       {
00892         perm[s_i] = s_i;
00893         porder[s_i] = i;
00894         plen[s_i] = rcnt;
00895         s_i++;
00896       }
00897     }
00898   }
00899 
00900   /*sort all unit rows */
00901   ILLutil_int_perm_quicksort (perm, plen, s_i);
00902 
00903   /* now go through the unit rows */
00904   for (k = 0; k < s_i; k++)
00905   {
00906     i = porder[perm[k]];
00907     rbeg = lp->rowbeg[i];
00908     rcnt = lp->rowcnt[i];
00909     selc = -1;
00910     dbl_EGlpNumCopy (seldj, dbl_INFTY);
00911     dbl_EGlpNumZero (selv);
00912 
00913     /* for every row s_i, compute min {d_j : d_j <0 , j is u or l or fr} */
00914     for (j = 0; j < rcnt; j++)
00915     {
00916       col = lp->rowind[rbeg + j];
00917       if (rcol[col] == 1)
00918         break;
00919       if (dbl_EGlpNumIsLessZero (dj[col]))
00920       {
00921         if (dbl_EGlpNumIsLess (dj[col], seldj))
00922         {
00923           selc = col;
00924           dbl_EGlpNumCopy (seldj, dj[col]);
00925           dbl_EGlpNumCopy (selv, lp->rowval[rbeg + j]);
00926         }
00927       }
00928     }
00929     /* select pivot element and update all d_j's */
00930     if (selc != -1)
00931     {
00932       nbelem++;
00933       irow[i] = 1;
00934       rrow[i] = 1;
00935       icol[selc] = 1;
00936       dbl_EGlpNumCopyFrac (c_dj, dj[selc], selv);
00937       vstat[selc] = STAT_BASIC;
00938       for (j = 0; j < rcnt; j++)
00939       {
00940         col = lp->rowind[rbeg + j];
00941         dbl_EGlpNumSubInnProdTo (dj[col], lp->rowval[rbeg + j], c_dj);
00942         rcol[col] = 1;
00943       }
00944     }
00945   }
00946 #if dbl_BASIS_STATS > 0
00947   printf ("unit rows = %d\n", s_i);
00948   printf ("nslacks %d, unit rows selected = %d\n", nslacks, nbelem - nslacks);
00949 #endif
00950   /* now go through remaining cols with dj = 0 */
00951   tot1 = vd.nfree + vd.nbndone;
00952 
00953   if (!dbl_EGlpNumIsNeqqZero (vd.cmax))
00954     dbl_EGlpNumOne (cmax);
00955   else
00956   {
00957     dbl_EGlpNumCopy (cmax, vd.cmax);
00958     dbl_EGlpNumMultUiTo (cmax, 1000);
00959   }
00960   for (j = 0; j < lp->ncols; j++)
00961   {
00962     if (vstat[j] == STAT_BASIC)
00963       continue;
00964     if (icol[j] == 1 || dbl_EGlpNumIsNeqZero (dj[j], dbl_BD_TOLER))
00965       continue;
00966     mcnt = lp->matcnt[j];
00967 
00968     dbl_EGlpNumSet (c_dj, (double) mcnt);
00969     switch (lp->vtype[j])
00970     {
00971     case VFREE:
00972       porder[tfree] = j;
00973       perm[tfree] = tfree;
00974       dbl_EGlpNumCopyFrac (qpenalty[tfree], lp->cz[j], cmax);
00975       dbl_EGlpNumAddTo (qpenalty[tfree], c_dj);
00976       tfree++;
00977       break;
00978 
00979     case VLOWER:
00980     case VUPPER:
00981       porder[vd.nfree + tbndone] = j;
00982       perm[vd.nfree + tbndone] = tbndone;
00983       dbl_EGlpNumCopyFrac (qpenalty[vd.nfree + tbndone], lp->cz[j], cmax);
00984       dbl_EGlpNumAddTo (qpenalty[vd.nfree + tbndone], c_dj);
00985       if (lp->vtype[j] == VLOWER)
00986         dbl_EGlpNumAddTo (qpenalty[vd.nfree + tbndone], lp->lz[j]);
00987       else
00988         dbl_EGlpNumSubTo (qpenalty[vd.nfree + tbndone], lp->uz[j]);
00989       tbndone++;
00990       break;
00991 
00992     case VFIXED:
00993     case VBOUNDED:
00994       porder[tot1 + tbounded] = j;
00995       perm[tot1 + tbounded] = tbounded;
00996       dbl_EGlpNumCopyFrac (qpenalty[tot1 + tbounded], lp->cz[j], cmax);
00997       dbl_EGlpNumAddTo (qpenalty[tot1 + tbounded], lp->lz[j]);
00998       dbl_EGlpNumSubTo (qpenalty[tot1 + tbounded], lp->uz[j]);
00999       dbl_EGlpNumAddTo (qpenalty[tot1 + tbounded], c_dj);
01000       tbounded++;
01001       break;
01002     }
01003   }
01004 #if dbl_BASIS_STATS > 0
01005   printf ("bfree %d, bone %d, bbnd %d\n", tfree, tbndone, tbounded);
01006 #endif
01007 
01008   dbl_ILLutil_EGlpNum_perm_quicksort (perm, qpenalty, tfree);
01009   dbl_ILLutil_EGlpNum_perm_quicksort (perm + vd.nfree, qpenalty + vd.nfree,
01010                                   tbndone);
01011   dbl_ILLutil_EGlpNum_perm_quicksort (perm + tot1, qpenalty + tot1, tbounded);
01012 
01013   tot2 = tfree + tbndone;
01014   for (i = 0; i < tbndone; i++)
01015   {
01016     perm[tfree + i] = perm[vd.nfree + i] + tfree;
01017     porder[tfree + i] = porder[vd.nfree + i];
01018   }
01019   for (i = 0; i < tbounded; i++)
01020   {
01021     perm[tot2 + i] = perm[tot1 + i] + tot2;
01022     porder[tot2 + i] = porder[tot1 + i];
01023   }
01024   tot2 += tbounded;
01025 
01026   nbelem =
01027     dbl_primal_col_select (lp, vstat, irow, rrow, unitcol, v, perm, porder, nbelem,
01028                        tot2);
01029   if (nbelem != lp->nrows)
01030   {
01031     printf ("complain: incorrect final basis size\n");
01032     rval = E_SIMPLEX_ERROR;
01033     ILL_CLEANUP;
01034   }
01035 
01036 CLEANUP:
01037   if (rval)
01038     dbl_ILLbasis_free_basisinfo (lp);
01039 
01040   ILL_IFFREE (irow, int);
01041   ILL_IFFREE (rrow, int);
01042 
01043   dbl_EGlpNumFreeArray (v);
01044   ILL_IFFREE (unitcol, int);
01045   ILL_IFFREE (icol, int);
01046   ILL_IFFREE (rcol, int);
01047 
01048   dbl_EGlpNumFreeArray (dj);
01049   ILL_IFFREE (perm, int);
01050   ILL_IFFREE (porder, int);
01051   ILL_IFFREE (plen, int);
01052 
01053   dbl_EGlpNumFreeArray (qpenalty);
01054   dbl_EGlpNumClearVar (seldj);
01055   dbl_EGlpNumClearVar (selv);
01056   dbl_EGlpNumClearVar (c_dj);
01057   dbl_EGlpNumClearVar (cmax);
01058   dbl_ILLbasis_clear_vardata (&vd);
01059   EG_RETURN (rval);
01060 }
01061 
01062 static int dbl_set_basis_indices (
01063   dbl_lpinfo * lp,
01064   int *vstat)
01065 {
01066   int i, b = 0, nb = 0;
01067   int vs;
01068 
01069   for (i = 0; i < lp->ncols; i++)
01070   {
01071     vs = vstat[i];
01072     lp->vstat[i] = vs;
01073 
01074     if (vs == STAT_BASIC)
01075     {
01076       lp->baz[b] = i;
01077       lp->vindex[i] = b;
01078       b++;
01079     }
01080     else if (vs == STAT_UPPER || vs == STAT_LOWER || vs == STAT_ZERO)
01081     {
01082       lp->nbaz[nb] = i;
01083       lp->vindex[i] = nb;
01084       nb++;
01085     }
01086     else
01087     {
01088       fprintf (stderr, "Error in basis creation\n");
01089       return E_SIMPLEX_ERROR;
01090     }
01091   }
01092   if (b != lp->nrows)
01093   {
01094     fprintf (stderr, "Error 2 in basis creation\n");
01095     return E_SIMPLEX_ERROR;
01096   }
01097   else if (nb != lp->nnbasic)
01098   {
01099     fprintf (stderr, "Error 3 in basis creation\n");
01100     return E_SIMPLEX_ERROR;
01101   }
01102   return 0;
01103 }
01104 
01105 int dbl_ILLbasis_get_initial (
01106   dbl_lpinfo * lp,
01107   int algorithm)
01108 {
01109   int rval = 0;
01110   int *vstat = NULL;
01111 
01112   dbl_ILLbasis_free_basisinfo (lp);
01113   dbl_ILLbasis_init_basisinfo (lp);
01114   rval = dbl_ILLbasis_build_basisinfo (lp);
01115   CHECKRVALG (rval, CLEANUP);
01116 
01117   ILL_SAFE_MALLOC (vstat, lp->ncols, int);
01118 
01119   if (algorithm == PRIMAL_SIMPLEX)
01120     rval = dbl_get_initial_basis1 (lp, vstat);
01121   else
01122     rval = dbl_get_initial_basis2 (lp, vstat);
01123 
01124   if (rval == E_SIMPLEX_ERROR)
01125   {
01126     #ifdef HAVE_ZLIB_H
01127     EGioFile_t *f = EGioOpen ("bad.lp.gz", "w");
01128     #else
01129     #ifdef HAVE_BZLIB_H
01130     EGioFile_t *f = EGioOpen ("bad.lp.bz2", "w");
01131     #else
01132     EGioFile_t *f = EGioOpen ("bad.lp", "w");
01133     #endif
01134     #endif
01135     int tval = dbl_ILLwrite_lp_file (lp->O, f, NULL);
01136     if (tval)
01137     {
01138       fprintf (stderr, "Error writing bad lp\n");
01139     }
01140     if (f != NULL)
01141       EGioClose (f);
01142   }
01143   CHECKRVALG (rval, CLEANUP);
01144 
01145   rval = dbl_set_basis_indices (lp, vstat);
01146   CHECKRVALG (rval, CLEANUP);
01147   lp->basisid = 0;
01148 
01149 CLEANUP:
01150   ILL_IFFREE (vstat, int);
01151 
01152   EG_RETURN (rval);
01153 }
01154 
01155 static int dbl_choose_basis (
01156   int algorithm,
01157   double pinf1,
01158   double dinf1,
01159   double pinf2,
01160   double dinf2)
01161 {
01162 /* We changed the constant definitions outside here, the actual numbers are
01163  * asigned in dbl_lpdata.c. the values are as follows:
01164  * dbl_CB_EPS = 0.001;
01165  * dbl_CB_PRI_RLIMIT = 0.25;
01166  * dbl_CB_INF_RATIO = 10.0; 
01167  * */
01168   int choice = 1;
01169   double rp, rd;
01170 
01171   if (algorithm == PRIMAL_SIMPLEX)
01172   {
01173     dbl_EGlpNumInitVar (rp);
01174     dbl_EGlpNumInitVar (rd);
01175     dbl_EGlpNumCopyDiff (rp, pinf1, pinf2);
01176     dbl_EGlpNumCopyDiff (rd, dinf1, dinf2);
01177     if (dbl_EGlpNumIsLeq (rp, dbl_CB_EPS) && dbl_EGlpNumIsLeq (rd, dbl_CB_EPS))
01178       choice = 1;
01179     else
01180     {
01181       dbl_EGlpNumSign (rp);
01182       dbl_EGlpNumSign (rd);
01183       if (dbl_EGlpNumIsLeq (rp, dbl_CB_EPS) && dbl_EGlpNumIsLeq (rd, dbl_CB_EPS))
01184         choice = 2;
01185       else if (dbl_EGlpNumIsLess (pinf1, pinf2) && dbl_EGlpNumIsLess (dinf2, dinf1))
01186       {
01187         choice = 1;
01188         dbl_EGlpNumCopyFrac (rp, pinf1, pinf2);
01189         dbl_EGlpNumCopyFrac (rd, dinf2, dinf1);
01190         dbl_EGlpNumMultTo (rd, dbl_CB_INF_RATIO);
01191         if (dbl_EGlpNumIsLess (dbl_CB_PRI_RLIMIT, rp) && (dbl_EGlpNumIsLess (rd, rp)))
01192           choice = 2;
01193       }
01194       else if (dbl_EGlpNumIsLess (pinf2, pinf1) && dbl_EGlpNumIsLess (dinf1, dinf2))
01195       {
01196         choice = 2;
01197         dbl_EGlpNumCopyFrac (rp, pinf2, pinf1);
01198         dbl_EGlpNumCopyFrac (rd, dinf1, dinf2);
01199         dbl_EGlpNumMultTo (rd, dbl_CB_INF_RATIO);
01200         if (dbl_EGlpNumIsLess (dbl_CB_PRI_RLIMIT, rp) && dbl_EGlpNumIsLess (rd, rp))
01201           choice = 1;
01202       }
01203       else
01204         choice = 1;
01205     }
01206     dbl_EGlpNumClearVar (rp);
01207     dbl_EGlpNumClearVar (rd);
01208   }
01209   ILL_IFTRACE ("%s:%d\n", __func__, choice);
01210   return choice;
01211 }
01212 
01213 int dbl_ILLbasis_get_cinitial (
01214   dbl_lpinfo * lp,
01215   int algorithm)
01216 {
01217   int rval = 0;
01218   int *vstat1 = NULL;
01219   int *vstat2 = NULL;
01220   int singular;
01221   int choice = 0;
01222 
01223 #if dbl_BASIS_STATS > 0
01224   int i, nz1 = 0, nz2 = 0;
01225 #endif
01226   double pinf1, pinf2, dinf1, dinf2;
01227   dbl_feas_info fi;
01228 
01229   dbl_EGlpNumInitVar (pinf1);
01230   dbl_EGlpNumInitVar (pinf2);
01231   dbl_EGlpNumInitVar (dinf1);
01232   dbl_EGlpNumInitVar (dinf2);
01233   dbl_EGlpNumInitVar (fi.totinfeas);
01234 
01235   dbl_ILLbasis_free_basisinfo (lp);
01236   dbl_ILLbasis_init_basisinfo (lp);
01237   rval = dbl_ILLbasis_build_basisinfo (lp);
01238   CHECKRVALG (rval, CLEANUP);
01239 
01240   ILL_SAFE_MALLOC (vstat1, lp->ncols, int);
01241   ILL_SAFE_MALLOC (vstat2, lp->ncols, int);
01242 
01243   if (algorithm != PRIMAL_SIMPLEX)
01244   {
01245     rval = dbl_get_initial_basis2 (lp, vstat2);
01246     CHECKRVALG (rval, CLEANUP);
01247     rval = dbl_set_basis_indices (lp, vstat2);
01248     lp->basisid = 0;
01249     ILL_CLEANUP;
01250   }
01251 
01252   rval = dbl_get_initial_basis1 (lp, vstat1);
01253   CHECKRVALG (rval, CLEANUP);
01254   rval = dbl_get_initial_basis2 (lp, vstat2);
01255   CHECKRVALG (rval, CLEANUP);
01256   lp->basisid = 0;
01257 
01258   /* handle first basis */
01259   rval = dbl_set_basis_indices (lp, vstat1);
01260   CHECKRVALG (rval, CLEANUP);
01261 #if dbl_BASIS_STATS > 0
01262   for (i = 0; i < lp->nrows; i++)
01263     nz1 += lp->matcnt[lp->baz[i]];
01264 #endif
01265   rval = dbl_ILLbasis_factor (lp, &singular);
01266   if (singular)
01267     MESSAGE (__QS_SB_VERB, "Singular Basis found!");
01268   CHECKRVALG (rval, CLEANUP);
01269 
01270   dbl_ILLfct_compute_piz (lp);
01271   dbl_ILLfct_compute_dz (lp);
01272   dbl_ILLfct_dual_adjust (lp, dbl_zeroLpNum);
01273   dbl_ILLfct_compute_xbz (lp);
01274 
01275   dbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
01276   dbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01277   dbl_EGlpNumCopy (pinf1, lp->pinfeas);
01278   dbl_EGlpNumCopy (dinf1, lp->dinfeas);
01279   /*
01280    * dbl_ILLfct_compute_pobj (lp);  obj1p = lp->objval;
01281    * dbl_ILLfct_compute_dobj (lp);  obj1d = lp->objval;
01282    */
01283 
01284   /* handle second basis */
01285   rval = dbl_set_basis_indices (lp, vstat2);
01286   CHECKRVALG (rval, CLEANUP);
01287 #if dbl_BASIS_STATS > 0
01288   for (i = 0; i < lp->nrows; i++)
01289     nz2 += lp->matcnt[lp->baz[i]];
01290 #endif
01291   rval = dbl_ILLbasis_factor (lp, &singular);
01292   if (singular)
01293     MESSAGE (__QS_SB_VERB, "Singular Basis found!");
01294   CHECKRVALG (rval, CLEANUP);
01295 
01296   dbl_ILLfct_compute_piz (lp);
01297   dbl_ILLfct_compute_dz (lp);
01298   dbl_ILLfct_dual_adjust (lp, dbl_zeroLpNum);
01299   dbl_ILLfct_compute_xbz (lp);
01300 
01301   dbl_ILLfct_check_pfeasible (lp, &fi, lp->tol->pfeas_tol);
01302   dbl_ILLfct_check_dfeasible (lp, &fi, lp->tol->dfeas_tol);
01303   dbl_EGlpNumCopy (pinf2, lp->pinfeas);
01304   dbl_EGlpNumCopy (dinf2, lp->dinfeas);
01305 
01306 #if dbl_BASIS_STATS > 0
01307   printf ("b1: nz %d pinf %.2f dinf %.2f\n", nz1, dbl_EGlpNumToLf (pinf1),
01308           dbl_EGlpNumToLf (dinf1));
01309   printf ("b2: nz %d pinf %.2f dinf %.2f\n", nz2, dbl_EGlpNumToLf (pinf2),
01310           dbl_EGlpNumToLf (dinf2));
01311 #endif
01312   choice = dbl_choose_basis (algorithm, pinf1, dinf1, pinf2, dinf2);
01313   if (choice == 1)
01314   {
01315     lp->fbasisid = -1;
01316     rval = dbl_set_basis_indices (lp, vstat1);
01317     CHECKRVALG (rval, CLEANUP);
01318   }
01319 
01320 CLEANUP:
01321   if (rval == E_SIMPLEX_ERROR)
01322   {
01323     #ifdef HAVE_ZLIB_H
01324     EGioFile_t *fil = EGioOpen ("bad.lp.gz", "w");
01325     #else
01326     #ifdef HAVE_BZLIB_H
01327     EGioFile_t *fil = EGioOpen ("bad.lp.bz2", "w");
01328     #else
01329     EGioFile_t *fil = EGioOpen ("bad.lp", "w");
01330     #endif
01331     #endif
01332     int tval = dbl_ILLwrite_lp_file (lp->O, fil, NULL);
01333 
01334     if (tval)
01335     {
01336       fprintf (stderr, "Error writing bad lp\n");
01337     }
01338     if (fil != NULL)
01339       EGioClose (fil);
01340   }
01341   ILL_IFFREE (vstat1, int);
01342   ILL_IFFREE (vstat2, int);
01343 
01344   dbl_EGlpNumClearVar (pinf1);
01345   dbl_EGlpNumClearVar (pinf2);
01346   dbl_EGlpNumClearVar (dinf1);
01347   dbl_EGlpNumClearVar (dinf2);
01348   dbl_EGlpNumClearVar (fi.totinfeas);
01349   EG_RETURN (rval);
01350 }
01351 
01352 int dbl_ILLbasis_factor (
01353   dbl_lpinfo * lp,
01354   int *singular)
01355 {
01356   int rval = 0;
01357   int i;
01358   int eindex;
01359   int lindex;
01360   int ltype;
01361   int lvstat;
01362   int nsing = 0;
01363   int *singr = 0;
01364   int *singc = 0;
01365 
01366   *singular = 0;
01367   do
01368   {
01369     if (lp->f)
01370     {
01371       dbl_ILLfactor_free_factor_work (lp->f);
01372     }
01373     else
01374     {
01375       ILL_SAFE_MALLOC (lp->f, 1, dbl_factor_work);
01376       dbl_EGlpNumInitVar (lp->f->fzero_tol);
01377       dbl_EGlpNumInitVar (lp->f->szero_tol);
01378       dbl_EGlpNumInitVar (lp->f->partial_tol);
01379       dbl_EGlpNumInitVar (lp->f->maxelem_orig);
01380       dbl_EGlpNumInitVar (lp->f->maxelem_factor);
01381       dbl_EGlpNumInitVar (lp->f->maxelem_cur);
01382       dbl_EGlpNumInitVar (lp->f->partial_cur);
01383       dbl_ILLfactor_init_factor_work (lp->f);
01384     }
01385     rval = dbl_ILLfactor_create_factor_work (lp->f, lp->O->nrows);
01386     CHECKRVALG (rval, CLEANUP);
01387 
01388     rval = dbl_ILLfactor (lp->f, lp->baz, lp->matbeg, lp->matcnt,
01389                       lp->matind, lp->matval, &nsing, &singr, &singc);
01390     CHECKRVALG (rval, CLEANUP);
01391 
01392     if (nsing != 0)
01393     {
01394       *singular = 1;
01395       MESSAGE (__QS_SB_VERB, "Found singular basis!");
01396       for (i = 0; i < nsing; i++)
01397       {
01398         eindex = lp->vindex[lp->O->rowmap[singr[i]]];
01399         lindex = singc[i];
01400         ltype = lp->vtype[lp->baz[lindex]];
01401 
01402         if (ltype == VBOUNDED || ltype == VLOWER || ltype == VARTIFICIAL)
01403           lvstat = STAT_LOWER;
01404         else if (ltype == VUPPER)
01405           lvstat = STAT_UPPER;
01406         else
01407           lvstat = STAT_ZERO;
01408 
01409         dbl_ILLfct_update_basis_info (lp, eindex, lindex, lvstat);
01410         lp->basisid++;
01411       }
01412       ILL_IFFREE (singr, int);
01413       ILL_IFFREE (singc, int);
01414     }
01415 
01416   } while (nsing != 0);
01417 
01418   lp->fbasisid = lp->basisid;
01419 
01420 CLEANUP:
01421   ILL_IFFREE (singr, int);
01422   ILL_IFFREE (singc, int);
01423 
01424   if (rval)
01425     fprintf (stderr, "Error: unknown in %s\n", __func__);
01426   EG_RETURN (rval);
01427 }
01428 
01429 int dbl_ILLbasis_refactor (
01430   dbl_lpinfo * lp)
01431 {
01432   int sing = 0;
01433   int rval = 0;
01434 
01435   rval = dbl_ILLbasis_factor (lp, &sing);
01436   if (sing)
01437   {
01438     MESSAGE (__QS_SB_VERB, "Singular Basis found!");
01439     rval = QS_LP_CHANGE_PREC;
01440     return rval;
01441   }
01442   EG_RETURN (rval);
01443 }
01444 
01445 void dbl_ILLbasis_column_solve (
01446   dbl_lpinfo * lp,
01447   dbl_svector * rhs,
01448   dbl_svector * soln)
01449 {
01450   dbl_ILLfactor_ftran (lp->f, rhs, soln);
01451 }
01452 
01453 void dbl_ILLbasis_column_solve_update (
01454   dbl_lpinfo * lp,
01455   dbl_svector * rhs,
01456   dbl_svector * upd,
01457   dbl_svector * soln)
01458 {
01459   dbl_ILLfactor_ftran_update (lp->f, rhs, upd, soln);
01460 }
01461 
01462 void dbl_ILLbasis_row_solve (
01463   dbl_lpinfo * lp,
01464   dbl_svector * rhs,
01465   dbl_svector * soln)
01466 {
01467   dbl_ILLfactor_btran (lp->f, rhs, soln);
01468 }
01469 
01470 int dbl_ILLbasis_update (
01471   dbl_lpinfo * lp,
01472   dbl_svector * y,
01473   int lindex,
01474   int *refactor,
01475   int *singular)
01476 {
01477 #if 0                           /* To always refactor, change 0 to 1 */
01478   *refactor = 1;
01479   return dbl_ILLbasis_factor (lp, singular);
01480 #else
01481 
01482   int rval = 0;
01483 
01484   *refactor = 0;
01485   rval = dbl_ILLfactor_update (lp->f, y, lindex, refactor);
01486   if (rval == E_FACTOR_BLOWUP || rval == E_UPDATE_SINGULAR_ROW
01487       || rval == E_UPDATE_SINGULAR_COL)
01488   {
01489 /* Bico - comment out for dist
01490        fprintf(stderr, "Warning: numerically bad basis in dbl_ILLfactor_update\n");
01491 */
01492     *refactor = 1;
01493     rval = 0;
01494   }
01495   if (rval == E_UPDATE_NOSPACE)
01496   {
01497     *refactor = 1;
01498     rval = 0;
01499   }
01500 
01501   if (*refactor)
01502   {
01503     rval = dbl_ILLbasis_factor (lp, singular);
01504     if (*singular)
01505       MESSAGE (__QS_SB_VERB, "Singular Basis found!");
01506   }
01507   if (rval)
01508   {
01509     EGioFile_t *eout = 0;
01510     int tval;
01511 
01512     printf ("write bad lp to factor.lp\n");
01513     fflush (stdout);
01514     #ifdef HAVE_ZLIB_H
01515     eout = EGioOpen ("factor.lp.gz", "w");
01516     #else
01517     #ifdef HAVE_BZLIB_H
01518     eout = EGioOpen ("factor.lp.bz2", "w");
01519     #else
01520     eout = EGioOpen ("factor.lp", "w");
01521     #endif
01522     #endif
01523     if (!eout)
01524     {
01525       fprintf (stderr, "could not open file to write bad factor lp\n");
01526     }
01527     else
01528     {
01529       tval = dbl_ILLwrite_lp_file (lp->O, eout, NULL);
01530       if (tval)
01531       {
01532         fprintf (stderr, "error while writing bad factor lp\n");
01533       }
01534       EGioClose (eout);
01535     }
01536 
01537     printf ("write bad basis to factor.bas\n");
01538     fflush (stdout);
01539     tval = dbl_ILLlib_writebasis (lp, 0, "factor.bas");
01540     if (tval)
01541     {
01542       fprintf (stderr, "error while writing factor basis\n");
01543     }
01544   }
01545 
01546   EG_RETURN (rval);
01547 #endif
01548 }

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