float128_presolve.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: float128_presolve.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 static int TRACE = 0;
00025 
00026 /****************************************************************************/
00027 /*                                                                          */
00028 /*                 Presolve Routine for Simplex Method                      */
00029 /*                                                                          */
00030 /*  EXPORTED FUNCTIONS                                                      */
00031 /*                                                                          */
00032 /*    int float128_ILLlp_add_logicals (ILLlpata *lp)                                 */
00033 /*    int float128_ILLlp_presolve (float128_ILLlpdata *lp)                                    */
00034 /*    int float128_ILLlp_scale (float128_ILLlpdata *lp)                                       */
00035 /*    void float128_ILLlp_sinfo_init (float128_ILLlp_sinfo *sinfo)                            */
00036 /*    void float128_ILLlp_sinfo_free (float128_ILLlp_sinfo *sinfo)                            */
00037 /*    void float128_ILLlp_predata_init (float128_ILLlp_predata *pre)                          */
00038 /*    void float128_ILLlp_predata_free (float128_ILLlp_predata *pre)                          */
00039 /*                                                                          */
00040 /*  NOTES                                                                   */
00041 /*                                                                          */
00042 /*    presolve will assume that logicals have been added.                   */
00043 /*                                                                          */
00044 /****************************************************************************/
00045 
00046 #include "qs_config.h"
00047 #include "float128_iqsutil.h"
00048 #include "float128_lpdata.h"
00049 #include "float128_lpdefs.h"
00050 //extern float128 float128_SZERO_TOLER;
00051 
00052 #define float128_ILL_LP_STATUS_OK (0)
00053 #define float128_ILL_PRE_FEAS_TOL float128_PFEAS_TOLER  //(1e-6)
00054 #define float128_ILL_PRE_ZERO_TOL float128_PIVOT_TOLER  //(1e-10)
00055 
00056 #define float128_ILL_PRE_DELETE_EMPTY_ROW               (1)
00057 #define float128_ILL_PRE_DELETE_SINGLETON_ROW           (2)
00058 #define float128_ILL_PRE_DELETE_FIXED_VARIABLE          (3)
00059 #define float128_ILL_PRE_DELETE_FORCED_VARIABLE         (4)
00060 #define float128_ILL_PRE_DELETE_SINGLETON_VARIABLE      (5)
00061 #define float128_ILL_PRE_DELETE_FREE_SINGLETON_VARIABLE (6)
00062 #define float128_ILL_PRE_DELETE_EMPTY_COLUMN            (7)
00063 
00064 #define float128_ILL_PRE_COL_STRUC                      (0)
00065 #define float128_ILL_PRE_COL_LOGICAL                    (1)
00066 
00067 static int float128_debug = 0;
00068 
00069 typedef struct
00070 {
00071   int row;
00072   int col;
00073   char coltype;
00074   char mark;
00075   char del;
00076   float128 coef;
00077 }
00078 float128_edge;
00079 
00080 typedef struct float128_node
00081 {
00082   float128_edge **adj;
00083   float128 obj;
00084   float128 lower;
00085   float128 upper;
00086   float128 rhs;
00087   int deg;
00088   char mark;
00089   char del;
00090   char coltype;
00091   char rowsense;
00092 }
00093 float128_node;
00094 
00095 typedef struct float128_intptr
00096 {
00097   int this_val;
00098   struct float128_intptr *next;
00099 }
00100 float128_intptr;
00101 
00102 typedef struct float128_graph
00103 {
00104   float128_edge *edgelist;
00105   struct float128_node *rows;
00106   struct float128_node *cols;
00107   int ecount;
00108   int nrows;
00109   int ncols;
00110   int nzcount;
00111   float128_edge **adjspace;
00112   ILLptrworld intptrworld;
00113   int objsense;
00114 }
00115 float128_graph;
00116 
00117 void float128_ILLlp_sinfo_init (
00118   float128_ILLlp_sinfo * sinfo),
00119   float128_ILLlp_sinfo_free (
00120   float128_ILLlp_sinfo * sinfo),
00121   float128_ILLlp_predata_init (
00122   float128_ILLlp_predata * pre),
00123   float128_ILLlp_predata_free (
00124   float128_ILLlp_predata * pre),
00125   float128_ILLlp_preop_init (
00126   float128_ILLlp_preop * op),
00127   float128_ILLlp_preop_free (
00128   float128_ILLlp_preop * op),
00129   float128_ILLlp_preline_init (
00130   float128_ILLlp_preline * line),
00131   float128_ILLlp_preline_free (
00132   float128_ILLlp_preline * line);
00133 
00134 int float128_ILLlp_sinfo_print (
00135   float128_ILLlp_sinfo * s);
00136 
00137 static void float128_set_fixed_variable (
00138   float128_graph * G,
00139   int j,
00140   float128 val),
00141   float128_get_implied_rhs_bounds (
00142   float128_graph * G,
00143   int i,
00144   float128 * lb,
00145   float128 * ub),
00146   float128_get_implied_variable_bounds (
00147   float128_graph * G,
00148   int j,
00149   float128_edge * a_ij,
00150   float128 * lb,
00151   float128 * ub),
00152   float128_dump_line (
00153   float128_ILLlp_preline * line),
00154   float128_init_graph (
00155   float128_graph * G),
00156   float128_free_graph (
00157   float128_graph * G),
00158   float128_dump_graph (
00159   float128_graph * G);
00160 
00161 static int float128_simple_presolve (
00162   float128_ILLlpdata * lp,
00163   float128_ILLlp_predata * pre,
00164   float128_ILLlp_sinfo * info,
00165   int pre_types,
00166   int *status),
00167   float128_grab_lp_line (
00168   float128_graph * G,
00169   int indx,
00170   float128_ILLlp_preline * line,
00171   int row_or_col),
00172   float128_grab_lp_info (
00173   float128_graph * G,
00174   char **colnames,
00175   float128_ILLlp_sinfo * info),
00176   float128_fixed_variables (
00177   float128_graph * G,
00178   float128_ILLlp_predata * pre),
00179   float128_empty_columns (
00180   float128_graph * G,
00181   float128_ILLlp_predata * pre),
00182   float128_singleton_rows (
00183   float128_graph * G,
00184   float128_ILLlp_predata * pre,
00185   int *hit),
00186   float128_forcing_constraints (
00187   float128_graph * G,
00188   float128_ILLlp_predata * pre,
00189   int *hit),
00190   float128_singleton_columns (
00191   float128_graph * G,
00192   float128_ILLlp_predata * pre,
00193   int *hit),
00194   float128_duplicate_rows (
00195   float128_graph * G,
00196   int *hit),
00197   float128_duplicate_cols (
00198   float128_graph * G,
00199   int *hit),
00200   float128_gather_dup_lists (
00201   int *s,
00202   int count,
00203   int *duptotal,
00204   int **dupcnt,
00205   int **dupind),
00206   float128_get_next_preop (
00207   float128_ILLlp_predata * pre,
00208   float128_ILLlp_preop ** op),
00209   float128_add_to_list (
00210   ILLptrworld * world,
00211   float128_intptr ** list,
00212   int i),
00213   float128_build_graph (
00214   float128_ILLlpdata * lp,
00215   float128_graph * G);
00216 
00217 
00218 ILL_PTRWORLD_ROUTINES (float128_intptr, intptralloc, intptr_bulkalloc, intptrfree)
00219 ILL_PTRWORLD_LISTFREE_ROUTINE (float128_intptr, intptr_listfree, intptrfree)
00220 ILL_PTRWORLD_LEAKS_ROUTINE (float128_intptr, intptr_check_leaks, this_val, int)
00221      int float128_ILLlp_add_logicals (
00222   float128_ILLlpdata * lp)
00223 {
00224   int rval = 0;
00225   int ncols, nrows, nzcount, i, aindex;
00226   char *sense;
00227   float128_ILLmatrix *A;
00228 
00229   if (!lp)
00230   {
00231     fprintf (stderr, "float128_ILLlp_add_logicals called with a NULL pointer\n");
00232     rval = 1;
00233     goto CLEANUP;
00234   }
00235 
00236   printf ("float128_ILLlp_add_logicals ...\n");
00237   fflush (stdout);
00238 
00239   A = &lp->A;
00240   sense = lp->sense;
00241   ncols = lp->ncols;
00242   nrows = lp->nrows;
00243   nzcount = lp->nzcount;
00244 
00245   if (nrows == 0)
00246     goto CLEANUP;
00247   float128_EGlpNumReallocArray (&(lp->obj), lp->colsize + nrows);
00248   float128_EGlpNumReallocArray (&(lp->upper), lp->colsize + nrows);
00249   float128_EGlpNumReallocArray (&(lp->lower), lp->colsize + nrows);
00250   lp->colnames =
00251     EGrealloc (lp->colnames, sizeof (char *) * (lp->colsize + nrows));
00252   //rval = ILLutil_reallocrus_count ((void **) &(lp->colnames),
00253   //                                 lp->colsize + nrows, sizeof (char *));
00254   //ILL_CLEANUP_IF (rval);
00255   memset (lp->colnames + ncols, 0, sizeof (char *) * nrows);
00256 
00257   ILL_SAFE_MALLOC (lp->rowmap, lp->rowsize, int);
00258 
00259 
00260   A->matcnt = EGrealloc (A->matcnt, sizeof (int) * (A->matcolsize + nrows));
00261   //rval = ILLutil_reallocrus_count ((void **) &(A->matcnt),
00262   //                                 A->matcolsize + nrows, sizeof (int));
00263   //ILL_CLEANUP_IF (rval);
00264 
00265   A->matbeg = EGrealloc (A->matbeg, sizeof (int) * (A->matcolsize + nrows));
00266   //rval = ILLutil_reallocrus_count ((void **) &(A->matbeg),
00267   //                                 A->matcolsize + nrows, sizeof (int));
00268   //ILL_CLEANUP_IF (rval);
00269 
00270   A->matind = EGrealloc (A->matind, sizeof (int) * (A->matsize + nrows));
00271   //rval = ILLutil_reallocrus_count ((void **) &(A->matind),
00272   //                                 A->matsize + nrows, sizeof (int));
00273   //ILL_CLEANUP_IF (rval);
00274   float128_EGlpNumReallocArray (&(A->matval), A->matsize + nrows);
00275 
00276   for (i = 0; i < nrows; i++)
00277   {
00278     A->matind[A->matsize + i] = -1;
00279   }
00280 
00281   aindex = A->matsize - A->matfree;
00282 
00283   for (i = 0; i < nrows; i++)
00284   {
00285     lp->rowmap[i] = ncols;
00286     float128_EGlpNumZero (lp->obj[ncols]);
00287     A->matcnt[ncols] = 1;
00288     A->matbeg[ncols] = aindex;
00289     A->matind[aindex] = i;
00290     switch (sense[i])
00291     {
00292     case 'E':                 /* Arificial */
00293       float128_EGlpNumZero (lp->lower[ncols]);
00294       float128_EGlpNumZero (lp->upper[ncols]);
00295       float128_EGlpNumOne (A->matval[aindex]);
00296       break;
00297     case 'G':                 /* Surplus   */
00298       float128_EGlpNumZero (lp->lower[ncols]);
00299       float128_EGlpNumCopy (lp->upper[ncols], float128_ILL_MAXDOUBLE);
00300       float128_EGlpNumOne (A->matval[aindex]);
00301       float128_EGlpNumSign (A->matval[aindex]);
00302       break;
00303     case 'L':                 /* Slack     */
00304       float128_EGlpNumZero (lp->lower[ncols]);
00305       float128_EGlpNumCopy (lp->upper[ncols], float128_ILL_MAXDOUBLE);
00306       float128_EGlpNumOne (A->matval[aindex]);
00307       break;
00308     case 'R':                 /* Range     */
00309       float128_EGlpNumZero (lp->lower[ncols]);
00310       float128_EGlpNumCopy (lp->upper[ncols], lp->rangeval[i]);
00311       float128_EGlpNumOne (A->matval[aindex]);
00312       float128_EGlpNumSign (A->matval[aindex]);
00313       break;
00314     default:
00315       fprintf (stderr, "unknown sense %c in float128_ILLlp_add_logicals\n", sense[i]);
00316       rval = 1;
00317       goto CLEANUP;
00318     }
00319     ncols++;
00320     nzcount++;
00321     aindex++;
00322   }
00323 
00324   lp->ncols = ncols;
00325   lp->nzcount = nzcount;
00326   A->matcols = ncols;
00327 
00328   lp->colsize += nrows;
00329   A->matsize += nrows;
00330   A->matcolsize += nrows;
00331 
00332   if (lp->rA)
00333   {
00334     float128_ILLlp_rows_clear (lp->rA);
00335   }
00336   else
00337   {
00338     ILL_SAFE_MALLOC (lp->rA, 1, float128_ILLlp_rows);
00339   }
00340 
00341   rval = float128_ILLlp_rows_init (lp->rA, lp, 1);
00342   ILL_CLEANUP_IF (rval);
00343 
00344 CLEANUP:
00345 
00346   ILL_RETURN (rval, "float128_ILLlp_add_logicals");
00347 }
00348 
00349 int float128_ILLlp_scale (
00350   float128_ILLlpdata * lp)
00351 {
00352   int rval = 0;
00353   int i, j, k, col, row, nstruct, start, stop;
00354   float128_ILLmatrix *A;
00355   float128 rho;
00356   float128 *gama = 0;
00357 
00358   float128_EGlpNumInitVar (rho);
00359 
00360   /* Columns - divide by largest absolute value */
00361 
00362   if (!lp)
00363   {
00364     ILL_ERROR (rval, "float128_ILLlp_scale called with a NULL pointer");
00365   }
00366 
00367   if (lp->nrows == 0 || lp->ncols == 0)
00368     goto CLEANUP;
00369 
00370   A = &lp->A;
00371   nstruct = lp->nstruct;
00372 
00373   for (j = 0; j < nstruct; j++)
00374   {
00375     col = lp->structmap[j];
00376     float128_EGlpNumZero (rho);
00377 
00378     start = A->matbeg[col];
00379     stop = start + A->matcnt[col];
00380 
00381     for (k = start; k < stop; k++)
00382     {
00383       float128_EGlpNumSetToMaxAbs (rho, A->matval[k]);
00384     }
00385 
00386     if (float128_EGlpNumIsGreatZero (rho))
00387     {
00388       for (k = start; k < stop; k++)
00389       {
00390         float128_EGlpNumDivTo (A->matval[k], rho);
00391       }
00392       float128_EGlpNumDivTo (lp->obj[col], rho);
00393       if (float128_EGlpNumIsNeqq (lp->lower[col], float128_ILL_MINDOUBLE))
00394         float128_EGlpNumMultTo (lp->lower[col], rho);
00395       if (float128_EGlpNumIsNeqq (lp->upper[col], float128_ILL_MAXDOUBLE))
00396         float128_EGlpNumMultTo (lp->upper[col], rho);
00397     }
00398   }
00399 
00400   gama = float128_EGlpNumAllocArray (lp->nrows);
00401   for (i = 0; i < lp->nrows; i++)
00402   {
00403     float128_EGlpNumZero (gama[i]);
00404   }
00405 
00406   for (j = 0; j < nstruct; j++)
00407   {
00408     col = lp->structmap[j];
00409     start = A->matbeg[col];
00410     stop = start + A->matcnt[col];
00411 
00412     for (k = start; k < stop; k++)
00413     {
00414       row = A->matind[k];
00415       float128_EGlpNumSetToMaxAbs (gama[row], A->matval[k]);
00416     }
00417   }
00418 
00419   for (j = 0; j < nstruct; j++)
00420   {
00421     col = lp->structmap[j];
00422     start = A->matbeg[col];
00423     stop = start + A->matcnt[col];
00424 
00425     for (k = start; k < stop; k++)
00426     {
00427       row = A->matind[k];
00428       if (float128_EGlpNumIsGreatZero (gama[row]))
00429       {
00430         float128_EGlpNumDivTo (A->matval[k], gama[row]);
00431       }
00432     }
00433   }
00434 
00435   for (i = 0; i < lp->nrows; i++)
00436   {
00437     if (float128_EGlpNumIsGreatZero ( gama[i]))
00438     {
00439       float128_EGlpNumDivTo (lp->rhs[i], gama[i]);
00440       col = lp->rowmap[i];
00441       if (float128_EGlpNumIsNeqq (lp->upper[col], float128_ILL_MAXDOUBLE))
00442       {
00443         float128_EGlpNumDivTo (lp->upper[col], gama[i]);  /* Ranged row */
00444       }
00445     }
00446   }
00447 
00448   if (lp->rA)
00449   {                             /* Need to clear the row version of data */
00450     float128_ILLlp_rows_clear (lp->rA);
00451     ILL_IFFREE (lp->rA, float128_ILLlp_rows);
00452   }
00453 
00454 
00455 CLEANUP:
00456 
00457   float128_EGlpNumClearVar (rho);
00458   float128_EGlpNumFreeArray (gama);
00459   ILL_RETURN (rval, "float128_ILLlp_scale");
00460 }
00461 
00462 int float128_ILLlp_presolve (
00463   float128_ILLlpdata * lp,
00464   int pre_types)
00465 {
00466   int rval = 0;
00467   int status = float128_ILL_LP_STATUS_OK;
00468   float128_ILLlp_predata *pre = 0;
00469   float128_ILLlp_sinfo *info = 0;
00470 
00471   if (!lp)
00472   {
00473     fprintf (stderr, "float128_ILLlp_presolve called with a NULL pointer\n");
00474     rval = 1;
00475     goto CLEANUP;
00476   }
00477 
00478 
00479 /*
00480     ILLlpdata_writelp (lp, 0);
00481     printf ("\n"); fflush (stdout);
00482 */
00483 
00484   ILL_SAFE_MALLOC (pre, 1, float128_ILLlp_predata);
00485   float128_ILLlp_predata_init (pre);
00486 
00487   ILL_SAFE_MALLOC (info, 1, float128_ILLlp_sinfo);
00488   float128_ILLlp_sinfo_init (info);
00489 
00490   rval = float128_simple_presolve (lp, pre, info, pre_types, &status);
00491   ILL_CLEANUP_IF (rval);
00492   if (status != float128_ILL_LP_STATUS_OK)
00493   {
00494     printf ("float128_simple_presolve returned with bad status\n");
00495     rval = 1;
00496     goto CLEANUP;
00497   }
00498 
00499 /*
00500     rval = float128_ILLlp_sinfo_print (info);
00501     ILL_CLEANUP_IF (rval);
00502 */
00503 
00504 CLEANUP:
00505 
00506   if (rval)
00507   {
00508     if (pre)
00509     {
00510       float128_ILLlp_predata_free (pre);
00511       ILL_IFFREE (pre, float128_ILLlp_predata);
00512     }
00513 
00514     if (info)
00515     {
00516       float128_ILLlp_sinfo_free (info);
00517       ILL_IFFREE (info, float128_ILLlp_sinfo);
00518     }
00519   }
00520   else
00521   {
00522     lp->presolve = pre;
00523     lp->sinfo = info;
00524   }
00525 
00526   ILL_RETURN (rval, "float128_ILLlp_presolve");
00527 }
00528 
00529 
00530 #if 0
00531 int ILLlp_presolve_addrow (
00532   float128_lpinfo * lp,
00533   int cnt,
00534   int *ind,
00535   double *val,
00536   double rhs)
00537 {
00538   int rval = 0;
00539   float128_ILLlpdata *qslp;
00540   float128_ILLlp_sinfo *S;
00541   float128_ILLmatrix *A;
00542 
00543   /* This will need to evolve into a function that handles the task */
00544   /* of working through the presolve data to determine the new LP   */
00545   /* created when a row is added to the original LP.                */
00546 
00547   /* The copies of the obj and bound used in the simplex code are   */
00548   /* also updated in this function.                                 */
00549 
00550   if (!lp)
00551   {
00552     fprintf (stderr, "ILLlp_presolve_addrow is called without an LP\n");
00553     rval = 1;
00554     goto CLEANUP;
00555   }
00556 
00557   if (lp->presolve != 0)
00558   {
00559     fprintf (stderr, "Not yet set up to handle addrows after presolve\n");
00560     rval = 1;
00561     goto CLEANUP;
00562   }
00563 
00564   qslp = lp->O;
00565   S = qslp->sinfo;
00566   A = S->A;
00567 
00568 
00569   rval = ILLlib_matrix_addrow (A, cnt, ind, val, rhs);
00570   ILL_CLEANUP_IF (rval);
00571 
00572 
00573 CLEANUP:
00574 
00575   ILL_RETURN (rval, "ILLlp_presolve_addrow");
00576 }
00577 #endif
00578 
00579 
00580 static int float128_simple_presolve (
00581   float128_ILLlpdata * lp,
00582   float128_ILLlp_predata * pre,
00583   float128_ILLlp_sinfo * info,
00584   int pre_types,
00585   int *status)
00586 {
00587   int rval = 0;
00588   int i, hit, newhit;
00589   float128_graph G;
00590 
00591   if (status)
00592     *status = float128_ILL_LP_STATUS_OK;
00593   float128_init_graph (&G);
00594 
00595   if (!lp)
00596   {
00597     fprintf (stderr, "float128_simple_presolve called with a NULL pointer\n");
00598     rval = 1;
00599     goto CLEANUP;
00600   }
00601 
00602   printf ("Initial Rows = %d, Cols = %d, Nzcount = %d\n",
00603           lp->nrows, lp->ncols, lp->nzcount);
00604   fflush (stdout);
00605 
00606   rval = float128_build_graph (lp, &G);
00607   ILL_CLEANUP_IF (rval);
00608   if (float128_debug)
00609     float128_dump_graph (&G);
00610 
00611   if (pre_types & float128_ILL_PRE_FIXED)
00612   {
00613     rval = float128_fixed_variables (&G, pre);
00614     ILL_CLEANUP_IF (rval);
00615   }
00616 
00617   do
00618   {
00619     hit = 0;
00620     if (pre_types & float128_ILL_PRE_SINGLE_ROW)
00621     {
00622       rval = float128_singleton_rows (&G, pre, &newhit);
00623       ILL_CLEANUP_IF (rval);
00624       hit += newhit;
00625     }
00626 
00627     if (pre_types & float128_ILL_PRE_FORCING)
00628     {
00629       rval = float128_forcing_constraints (&G, pre, &newhit);
00630       ILL_CLEANUP_IF (rval);
00631       hit += newhit;
00632     }
00633 
00634     if (pre_types & float128_ILL_PRE_SINGLE_COL)
00635     {
00636       rval = float128_singleton_columns (&G, pre, &newhit);
00637       ILL_CLEANUP_IF (rval);
00638       hit += newhit;
00639     }
00640 
00641     if (pre_types & float128_ILL_PRE_DUPLICATE_ROW)
00642     {
00643       rval = float128_duplicate_rows (&G, &newhit);
00644       ILL_CLEANUP_IF (rval);
00645       hit += newhit;
00646     }
00647 
00648     if (pre_types & float128_ILL_PRE_DUPLICATE_COL)
00649     {
00650       rval = float128_duplicate_cols (&G, &newhit);
00651       ILL_CLEANUP_IF (rval);
00652       hit += newhit;
00653     }
00654 
00655 
00656 /*
00657         {
00658             int k, cnt = 0;
00659             for (i = 0; i < G.ncols; i++) {
00660                 if (G.cols[i].del == 0) {
00661                     for (k = 0; k < G.cols[i].deg; k++)  {
00662                         if (G.cols[i].adj[k]->del == 0) {
00663                             cnt++;
00664                         }
00665                     }
00666                 }
00667             }
00668             printf ("Current NZCOUNT = %d\n", cnt); fflush (stdout);
00669         }
00670 */
00671   } while (hit);
00672 
00673   if (float128_ILL_PRE_EMPTY_COL)
00674   {
00675     rval = float128_empty_columns (&G, pre);
00676     ILL_CLEANUP_IF (rval);
00677   }
00678 
00679   if (float128_debug)
00680   {
00681     printf ("Operations\n");
00682     for (i = 0; i < pre->opcount; i++)
00683     {
00684       switch (pre->oplist[i].ptype)
00685       {
00686       case float128_ILL_PRE_DELETE_EMPTY_ROW:
00687         printf ("Delete Empty Row: %d\n", pre->oplist[i].rowindex);
00688         fflush (stdout);
00689         break;
00690       case float128_ILL_PRE_DELETE_SINGLETON_ROW:
00691         printf ("Delete Singleton Row: %d (col %d)\n",
00692                 pre->oplist[i].rowindex, pre->oplist[i].colindex);
00693         fflush (stdout);
00694         float128_dump_line (&pre->oplist[i].line);
00695         break;
00696       case float128_ILL_PRE_DELETE_FIXED_VARIABLE:
00697         printf ("Delete Fixed Variable: %d\n", pre->oplist[i].colindex);
00698         fflush (stdout);
00699         float128_dump_line (&pre->oplist[i].line);
00700         break;
00701       case float128_ILL_PRE_DELETE_FORCED_VARIABLE:
00702         printf ("Delete Forced Variable: %d\n", pre->oplist[i].colindex);
00703         fflush (stdout);
00704         float128_dump_line (&pre->oplist[i].line);
00705         break;
00706       case float128_ILL_PRE_DELETE_SINGLETON_VARIABLE:
00707         printf ("Delete Singleton Variable: %d\n", pre->oplist[i].colindex);
00708         fflush (stdout);
00709         float128_dump_line (&pre->oplist[i].line);
00710         break;
00711       case float128_ILL_PRE_DELETE_FREE_SINGLETON_VARIABLE:
00712         printf ("Delete Free Singleton Variable: %d\n",
00713                 pre->oplist[i].colindex);
00714         fflush (stdout);
00715         float128_dump_line (&pre->oplist[i].line);
00716         break;
00717       case float128_ILL_PRE_DELETE_EMPTY_COLUMN:
00718         printf ("Delete Empty Column: %d\n", pre->oplist[i].colindex);
00719         fflush (stdout);
00720         float128_dump_line (&pre->oplist[i].line);
00721         break;
00722       default:
00723         fprintf (stderr, "unknon presolve operation\n");
00724         rval = 1;
00725         goto CLEANUP;
00726       }
00727     }
00728     printf ("\n");
00729   }
00730 
00731   rval = float128_grab_lp_info (&G, lp->colnames, info);
00732   ILL_CLEANUP_IF (rval);
00733 
00734 /*
00735     printf ("Final Rows = %d, Cols = %d, Nzcount = %d\n",
00736                info->nrows, info->ncols, info->nzcount);
00737     fflush (stdout);
00738 */
00739 
00740 
00741 CLEANUP:
00742 
00743   float128_free_graph (&G);
00744   ILL_RETURN (rval, "float128_simple_presolve");
00745 }
00746 
00747 static int float128_grab_lp_line (
00748   float128_graph * G,
00749   int indx,
00750   float128_ILLlp_preline * line,
00751   int row_or_col)
00752 {
00753   int rval = 0;
00754   int k, cnt;
00755   float128_node *n;
00756 
00757   if (row_or_col == 0)
00758     n = &G->rows[indx];
00759   else
00760     n = &G->cols[indx];
00761 
00762   line->count = 0;
00763 
00764   for (k = 0; k < n->deg; k++)
00765   {
00766     if (n->adj[k]->del == 0)
00767     {
00768       line->count++;
00769     }
00770   }
00771 
00772   if (line->count)
00773   {
00774     ILL_SAFE_MALLOC (line->ind, line->count, int);
00775 
00776     line->val = float128_EGlpNumAllocArray (line->count);
00777     if (!line->ind || !line->val)
00778     {
00779       fprintf (stderr, "out of memory in float128_grab_lp_line\n");
00780       rval = 1;
00781       goto CLEANUP;
00782     }
00783     for (k = 0, cnt = 0; k < n->deg; k++)
00784     {
00785       if (n->adj[k]->del == 0)
00786       {
00787         line->ind[cnt] = n->adj[k]->row;
00788         float128_EGlpNumCopy (line->val[cnt], n->adj[k]->coef);
00789         cnt++;
00790       }
00791     }
00792   }
00793 
00794   if (row_or_col == 0)
00795   {
00796     float128_EGlpNumCopy (line->rhs, n->rhs);
00797   }
00798   else
00799   {
00800     float128_EGlpNumCopy (line->obj, n->obj);
00801     float128_EGlpNumCopy (line->lower, n->lower);
00802     float128_EGlpNumCopy (line->upper, n->upper);
00803   }
00804 
00805   line->row_or_col = row_or_col;
00806 
00807 CLEANUP:
00808 
00809   ILL_RETURN (rval, "float128_grab_lp_line");
00810 }
00811 
00812 static void float128_dump_line (
00813   float128_ILLlp_preline * line)
00814 {
00815   int k;
00816 
00817   printf (" ");
00818   if (line->row_or_col == 0)
00819   {
00820     for (k = 0; k < line->count; k++)
00821     {
00822       printf (" C%d->%g", line->ind[k], float128_EGlpNumToLf (line->val[k]));
00823     }
00824     printf (" RHS->%g\n", float128_EGlpNumToLf (line->rhs));
00825   }
00826   else
00827   {
00828     for (k = 0; k < line->count; k++)
00829     {
00830       printf (" R%d->%g", line->ind[k], float128_EGlpNumToLf (line->val[k]));
00831     }
00832     printf (" Obj->%g  LB->%g  UB->%g\n", float128_EGlpNumToLf (line->obj),
00833             float128_EGlpNumToLf (line->lower), float128_EGlpNumToLf (line->upper));
00834   }
00835   fflush (stdout);
00836 }
00837 
00838 static int float128_grab_lp_info (
00839   float128_graph * G,
00840   char **colnames,
00841   float128_ILLlp_sinfo * info)
00842 {
00843   int rval = 0;
00844   int ncols = 0, nrows = 0, nzcount = 0;
00845   int i, j, k, cnt, len;
00846   float128_node *grows = G->rows;
00847   float128_node *gcols = G->cols;
00848   int *tdeg = 0;
00849   int *map = 0;
00850   char *buf = 0;
00851   float128_ILLmatrix *A = &info->A;
00852 
00853   ILL_SAFE_MALLOC (tdeg, G->ncols, int);
00854   ILL_SAFE_MALLOC (map, G->nrows, int);
00855 
00856   if (!tdeg || !map)
00857   {
00858     fprintf (stderr, "out of memory in float128_grab_lp_info\n");
00859     rval = 1;
00860     goto CLEANUP;
00861   }
00862 
00863   for (i = 0; i < G->nrows; i++)
00864   {
00865     if (grows[i].del == 0)
00866     {
00867       map[i] = nrows;
00868       nrows++;
00869     }
00870   }
00871 
00872   for (j = 0; j < G->ncols; j++)
00873   {
00874     if (gcols[j].del == 0)
00875     {
00876       tdeg[ncols] = 0;
00877       for (k = 0; k < gcols[j].deg; k++)
00878       {
00879         if (gcols[j].adj[k]->del == 0)
00880         {
00881           tdeg[ncols]++;
00882           nzcount++;
00883         }
00884       }
00885       ncols++;
00886     }
00887   }
00888 
00889   info->ncols = ncols;
00890   info->nrows = nrows;
00891   info->nzcount = nzcount;
00892 
00893   info->rowsize = nrows;
00894   info->colsize = ncols;
00895 
00896   info->rhs = float128_EGlpNumAllocArray (nrows);
00897   info->obj = float128_EGlpNumAllocArray (ncols);
00898   info->upper = float128_EGlpNumAllocArray (ncols);
00899   info->lower = float128_EGlpNumAllocArray (ncols);
00900   A->matval = float128_EGlpNumAllocArray (info->nzcount + 1);
00901   ILL_SAFE_MALLOC (A->matind, info->nzcount + 1, int);
00902   ILL_SAFE_MALLOC (A->matcnt, info->colsize, int);
00903   ILL_SAFE_MALLOC (A->matbeg, info->colsize, int);
00904 
00905   if (!info->rhs || !info->obj || !info->lower || !info->upper ||
00906       !A->matval || !A->matind || !A->matcnt || !A->matbeg)
00907   {
00908     fprintf (stderr, "out of memory in grab_lp\n");
00909     rval = 1;
00910     goto CLEANUP;
00911   }
00912 
00913   A->matind[info->nzcount] = -1;
00914   A->matsize = info->nzcount + 1;
00915   A->matcolsize = info->colsize;
00916   A->matfree = 1;
00917   A->matcols = ncols;
00918   A->matrows = nrows;
00919 
00920 
00921   nrows = 0;
00922   for (i = 0; i < G->nrows; i++)
00923   {
00924     if (grows[i].del == 0)
00925     {
00926       float128_EGlpNumCopy (info->rhs[nrows], grows[i].rhs);
00927       nrows++;
00928     }
00929   }
00930 
00931   ncols = 0;
00932   cnt = 0;
00933   for (j = 0; j < G->ncols; j++)
00934   {
00935     if (gcols[j].del == 0)
00936     {
00937       float128_EGlpNumCopy (info->obj[ncols], gcols[j].obj);
00938       float128_EGlpNumCopy (info->lower[ncols], gcols[j].lower);
00939       float128_EGlpNumCopy (info->upper[ncols], gcols[j].upper);
00940       A->matcnt[ncols] = tdeg[ncols];
00941       A->matbeg[ncols] = cnt;
00942       for (k = 0; k < gcols[j].deg; k++)
00943       {
00944         if (gcols[j].adj[k]->del == 0)
00945         {
00946           float128_EGlpNumCopy (A->matval[cnt], gcols[j].adj[k]->coef);
00947           A->matind[cnt] = map[gcols[j].adj[k]->row];
00948           cnt++;
00949         }
00950       }
00951       ncols++;
00952     }
00953   }
00954 
00955   if (colnames)
00956   {
00957     ILL_SAFE_MALLOC (info->colnames, info->colsize, char *);
00958 
00959     if (!info->colnames)
00960     {
00961       fprintf (stderr, "out of memory in grab_lp\n");
00962       rval = 1;
00963       goto CLEANUP;
00964     }
00965     for (j = 0; j < info->colsize; j++)
00966     {
00967       info->colnames[j] = 0;
00968     }
00969 
00970     ILL_SAFE_MALLOC (buf, ILL_namebufsize, char);
00971 
00972     if (!buf)
00973     {
00974       fprintf (stderr, "out of memory in grab_lp\n");
00975       rval = 1;
00976       goto CLEANUP;
00977     }
00978     ncols = 0;
00979     for (j = 0; j < G->ncols; j++)
00980     {
00981       if (gcols[j].del == 0)
00982       {
00983         if (gcols[j].coltype == float128_ILL_PRE_COL_STRUC)
00984         {
00985           len = strlen (colnames[j]) + 1;
00986           ILL_SAFE_MALLOC (info->colnames[ncols], len, char);
00987 
00988           if (!info->colnames[ncols])
00989           {
00990             fprintf (stderr, "out of memory in grab_lp\n");
00991             rval = 1;
00992             goto CLEANUP;
00993           }
00994           strcpy (info->colnames[ncols], colnames[j]);
00995         }
00996         else
00997         {
00998           for (k = 0; k < gcols[j].deg; k++)
00999           {
01000             if (gcols[j].adj[k]->del == 0)
01001             {
01002               i = gcols[j].adj[k]->row;
01003               break;
01004             }
01005           }
01006           if (k == gcols[j].deg)
01007           {
01008             fprintf (stderr, "problem with float128_graph in grab_lp\n");
01009             rval = 1;
01010             goto CLEANUP;
01011           }
01012           sprintf (buf, "s%d", i);
01013           len = strlen (buf) + 1;
01014           ILL_SAFE_MALLOC (info->colnames[ncols], len, char);
01015 
01016           if (!info->colnames[ncols])
01017           {
01018             fprintf (stderr, "out of memory in grab_lp\n");
01019             rval = 1;
01020             goto CLEANUP;
01021           }
01022           strcpy (info->colnames[ncols], buf);
01023         }
01024         ncols++;
01025       }
01026     }
01027   }
01028 
01029 /* float128_ADD STRUCT VARIABLE STUFF */
01030 
01031 
01032 CLEANUP:
01033 
01034   if (rval)
01035   {
01036     float128_ILLlp_sinfo_free (info);
01037   }
01038   ILL_IFFREE (tdeg, int);
01039   ILL_IFFREE (map, int);
01040   ILL_IFFREE (buf, char);
01041 
01042   ILL_RETURN (rval, "float128_grab_lp_info");
01043 }
01044 
01045 static int float128_fixed_variables (
01046   float128_graph * G,
01047   float128_ILLlp_predata * pre)
01048 {
01049   int rval = 0;
01050   int j;
01051   int ncols = G->ncols;
01052   float128_node *cols = G->cols;
01053   float128_ILLlp_preop *op = 0;
01054 
01055   for (j = 0; j < ncols; j++)
01056   {
01057     if (cols[j].del == 0)
01058     {
01059       if (float128_EGlpNumIsEqqual (cols[j].lower, cols[j].upper))
01060       {
01061         rval = float128_get_next_preop (pre, &op);
01062         ILL_CLEANUP_IF (rval);
01063 
01064         op->colindex = j;
01065         op->rowindex = -1;
01066         op->ptype = float128_ILL_PRE_DELETE_FIXED_VARIABLE;
01067 
01068         rval = float128_grab_lp_line (G, op->colindex, &op->line, 1);
01069         ILL_CLEANUP_IF (rval);
01070         pre->opcount++;
01071 
01072         float128_set_fixed_variable (G, j, cols[j].lower);
01073       }
01074     }
01075   }
01076 
01077 CLEANUP:
01078 
01079   ILL_RETURN (rval, "float128_fixed_variables");
01080 }
01081 
01082 static int float128_empty_columns (
01083   float128_graph * G,
01084   float128_ILLlp_predata * pre)
01085 {
01086   int rval = 0;
01087   int j, k;
01088   int ncols = G->ncols;
01089   float128_node *cols = G->cols;
01090   float128_ILLlp_preop *op = 0;
01091   float128 objtmp;
01092 
01093   float128_EGlpNumInitVar (objtmp);
01094 
01095   for (j = 0; j < ncols; j++)
01096   {
01097     if (cols[j].del == 0)
01098     {
01099       for (k = 0; k < cols[j].deg; k++)
01100       {
01101         if (cols[j].adj[k]->del == 0)
01102           break;
01103       }
01104       if (k == cols[j].deg)
01105       {
01106         rval = float128_get_next_preop (pre, &op);
01107         ILL_CLEANUP_IF (rval);
01108 
01109         op->colindex = j;
01110         op->rowindex = -1;
01111         op->ptype = float128_ILL_PRE_DELETE_EMPTY_COLUMN;
01112 
01113         rval = float128_grab_lp_line (G, op->colindex, &op->line, 1);
01114         ILL_CLEANUP_IF (rval);
01115         pre->opcount++;
01116         float128_EGlpNumCopy (objtmp, cols[j].obj);
01117         if (G->objsense < 0)
01118           float128_EGlpNumSign (objtmp);
01119         if (!float128_EGlpNumIsNeqZero (objtmp, float128_ILL_PRE_FEAS_TOL))
01120         {
01121           float128_set_fixed_variable (G, j, cols[j].lower);
01122         }
01123         else if (float128_EGlpNumIsGreatZero (objtmp))
01124         {
01125           if (float128_EGlpNumIsEqqual (cols[j].lower, float128_ILL_MINDOUBLE))
01126           {
01127             printf ("unbounded prob detected in float128_empty_columns\n");
01128             printf ("col %d, obj %g\n", j, float128_EGlpNumToLf (cols[j].obj));
01129             fflush (stdout);
01130             rval = 1;
01131             goto CLEANUP;
01132           }
01133           else
01134           {
01135             float128_set_fixed_variable (G, j, cols[j].lower);
01136           }
01137         }
01138         else if (float128_EGlpNumIsLessZero (objtmp))
01139         {
01140           if (float128_EGlpNumIsEqqual (cols[j].upper, float128_ILL_MAXDOUBLE))
01141           {
01142             printf ("unbounded prob detected in float128_empty_columns\n");
01143             printf ("col %d, obj %g\n", j, float128_EGlpNumToLf (cols[j].obj));
01144             fflush (stdout);
01145             rval = 1;
01146             goto CLEANUP;
01147           }
01148           else
01149           {
01150             float128_set_fixed_variable (G, j, cols[j].upper);
01151           }
01152         }
01153         else
01154         {
01155           float128_set_fixed_variable (G, j, cols[j].lower);
01156         }
01157       }
01158     }
01159   }
01160 
01161 CLEANUP:
01162 
01163   float128_EGlpNumClearVar (objtmp);
01164   ILL_RETURN (rval, "float128_empty_columns");
01165 }
01166 
01167 static int float128_singleton_rows (
01168   float128_graph * G,
01169   float128_ILLlp_predata * pre,
01170   int *hit)
01171 {
01172   int rval = 0;
01173   int rowindex, i, k, h;
01174   int nrows = G->nrows;
01175   float128_node *rows = G->rows;
01176   float128_node *cols = G->cols;
01177   float128_node *r, *c;
01178   float128_edge *pivot, *f;
01179   float128_intptr *next, *list = 0;
01180   int *tdeg = 0;
01181   float128 val;
01182   float128_ILLlp_preop *op = 0;
01183 
01184   float128_EGlpNumInitVar (val);
01185 
01186   *hit = 0;
01187   if (G->nrows == 0)
01188     goto CLEANUP;
01189 
01190   ILL_SAFE_MALLOC (tdeg, G->nrows, int);
01191 
01192   if (!tdeg)
01193   {
01194     fprintf (stderr, "out of memory in float128_singleton_rows\n");
01195     rval = 1;
01196     goto CLEANUP;
01197   }
01198 
01199   for (i = 0; i < nrows; i++)
01200   {
01201     if (rows[i].del == 0)
01202     {
01203       tdeg[i] = 0;
01204       for (k = 0; k < rows[i].deg; k++)
01205       {
01206         if (rows[i].adj[k]->del == 0)
01207         {
01208           tdeg[i]++;
01209         }
01210       }
01211       if (tdeg[i] <= 1)
01212       {
01213         rval = float128_add_to_list (&G->intptrworld, &list, i);
01214         ILL_CLEANUP_IF (rval);
01215       }
01216     }
01217   }
01218 
01219   while (list)
01220   {
01221     (*hit)++;
01222     rowindex = list->this_val;
01223     next = list->next;
01224     intptrfree (&G->intptrworld, list);
01225     list = next;
01226 
01227     rval = float128_get_next_preop (pre, &op);
01228     ILL_CLEANUP_IF (rval);
01229 
01230     r = &rows[rowindex];
01231 
01232     if (tdeg[rowindex] == 0)
01233     {
01234       if (float128_EGlpNumIsNeqZero (r->rhs, float128_ILL_PRE_FEAS_TOL))
01235       {
01236         printf ("infeasible row detected in singleton_row\n");
01237         printf ("empty row with rhs = %g\n", float128_EGlpNumToLf (r->rhs));
01238         fflush (stdout);
01239         rval = 1;
01240         goto CLEANUP;
01241       }
01242       op->ptype = float128_ILL_PRE_DELETE_EMPTY_ROW;
01243       op->rowindex = rowindex;
01244     }
01245     else
01246     {
01247       /*  Find the "pivot" entry and colum */
01248 
01249       for (k = 0; k < r->deg; k++)
01250       {
01251         if (r->adj[k]->del == 0)
01252           break;
01253       }
01254       if (k == r->deg)
01255       {
01256         fprintf (stderr, "lost an float128_edge in float128_singleton_rows\n");
01257         rval = 1;
01258         goto CLEANUP;
01259       }
01260 
01261       pivot = r->adj[k];
01262       c = &cols[pivot->col];
01263 
01264       /*  Store data from operation (incluing the col coefs) */
01265 
01266       op->ptype = float128_ILL_PRE_DELETE_SINGLETON_ROW;
01267       op->rowindex = rowindex;
01268       op->colindex = c - cols;
01269       float128_EGlpNumCopy (op->line.rhs, r->rhs);
01270       rval = float128_grab_lp_line (G, op->colindex, &op->line, 1);
01271       ILL_CLEANUP_IF (rval);
01272 
01273       /*  Fix the x[c] to its rhs value */
01274       /*val = r->rhs / pivot->coef; */
01275       float128_EGlpNumCopyFrac (val, r->rhs, pivot->coef);
01276       /* if (val < c->lower - float128_ILL_PRE_FEAS_TOL ||
01277        * val > c->upper + float128_ILL_PRE_FEAS_TOL) */
01278       if (float128_EGlpNumIsSumLess (val, float128_ILL_PRE_FEAS_TOL, c->lower) ||
01279           float128_EGlpNumIsSumLess (c->upper, float128_ILL_PRE_FEAS_TOL, val))
01280       {
01281         printf ("infeasible bounds detected in singleton_row %d\n", rowindex);
01282         printf ("lower->%g  upper->%g  val = %g\n",
01283                 float128_EGlpNumToLf (c->lower), float128_EGlpNumToLf (c->upper),
01284                 float128_EGlpNumToLf (val));
01285         fflush (stdout);
01286         rval = 1;
01287         goto CLEANUP;
01288       }
01289       else
01290       {
01291         float128_EGlpNumCopy (c->lower, val);
01292         float128_EGlpNumCopy (c->upper, val);
01293       }
01294 
01295       /*  Delete x[c] from other rows (and adjust their rhs) */
01296 
01297       c->del = 1;
01298 
01299       for (h = 0; h < c->deg; h++)
01300       {
01301         f = c->adj[h];
01302         if (f->del == 0)
01303         {
01304           /*rows[f->row].rhs -= (f->coef * c->lower); */
01305           float128_EGlpNumSubInnProdTo (rows[f->row].rhs, f->coef, c->lower);
01306           tdeg[f->row]--;
01307           if (tdeg[f->row] == 1)
01308           {
01309             if (f == pivot)
01310             {
01311               fprintf (stderr, "bad pivot element\n");
01312               rval = 1;
01313               goto CLEANUP;
01314             }
01315             rval = float128_add_to_list (&G->intptrworld, &list, f->row);
01316             ILL_CLEANUP_IF (rval);
01317           }
01318           f->del = 1;
01319         }
01320       }
01321     }
01322 
01323     r->del = 1;
01324     pre->opcount++;
01325   }
01326 
01327 CLEANUP:
01328 
01329   ILL_IFFREE (tdeg, int);
01330 
01331   intptr_listfree (&G->intptrworld, list);
01332   float128_EGlpNumClearVar (val);
01333   ILL_RETURN (rval, "float128_singleton_rows");
01334 }
01335 
01336 static int float128_forcing_constraints (
01337   float128_graph * G,
01338   float128_ILLlp_predata * pre,
01339   int *hit)
01340 {
01341   int rval = 0;
01342   int i, j, k, ts;
01343   float128_node *rows = G->rows;
01344   float128_node *cols = G->cols;
01345   float128_edge *e;
01346   int nrows = G->nrows;
01347   float128 ub, lb;
01348   float128_ILLlp_preop *op = 0;
01349 
01350   float128_EGlpNumInitVar (ub);
01351   float128_EGlpNumInitVar (lb);
01352 
01353   *hit = 0;
01354 
01355   for (i = 0; i < nrows; i++)
01356   {
01357     if (rows[i].del == 0)
01358     {
01359       float128_get_implied_rhs_bounds (G, i, &lb, &ub);
01360       if (float128_EGlpNumIsSumLess (rows[i].rhs, float128_ILL_PRE_FEAS_TOL, lb) ||
01361           float128_EGlpNumIsSumLess (ub, float128_ILL_PRE_FEAS_TOL, rows[i].rhs))
01362       {
01363         printf ("infeasible row detected in float128_forcing_constraints\n");
01364         printf ("Row %d:  RHS->%g  LBnd->%g  UBnd->%g\n",
01365                 i, float128_EGlpNumToLf (rows[i].rhs),
01366                 float128_EGlpNumToLf (lb), float128_EGlpNumToLf (ub));
01367         fflush (stdout);
01368         rval = 1;
01369         goto CLEANUP;
01370       }
01371       else if (float128_EGlpNumIsDiffLess (rows[i].rhs, float128_ILL_PRE_FEAS_TOL, lb) ||
01372                float128_EGlpNumIsDiffLess (ub, float128_ILL_PRE_FEAS_TOL, rows[i].rhs))
01373       {
01374         (*hit)++;
01375         ts = (float128_EGlpNumIsDiffLess (rows[i].rhs, float128_ILL_PRE_FEAS_TOL, lb) ? 0 : 1);
01376         for (k = 0; k < rows[i].deg; k++)
01377         {
01378           e = rows[i].adj[k];
01379           if (e->del == 0)
01380           {
01381             j = e->col;
01382 
01383             rval = float128_get_next_preop (pre, &op);
01384             ILL_CLEANUP_IF (rval);
01385 
01386             op->colindex = j;
01387             op->rowindex = i;
01388             op->ptype = float128_ILL_PRE_DELETE_FORCED_VARIABLE;
01389 
01390             rval = float128_grab_lp_line (G, j, &op->line, 1);
01391             ILL_CLEANUP_IF (rval);
01392             pre->opcount++;
01393 
01394             if ((ts == 0 && float128_EGlpNumIsLessZero (e->coef)) ||
01395                 (ts == 1 && float128_EGlpNumIsGreatZero (e->coef)))
01396             {
01397               float128_set_fixed_variable (G, j, cols[j].upper);
01398             }
01399             else
01400             {
01401               float128_set_fixed_variable (G, j, cols[j].lower);
01402             }
01403           }
01404         }
01405       }
01406     }
01407   }
01408 
01409 CLEANUP:
01410 
01411   float128_EGlpNumClearVar (ub);
01412   float128_EGlpNumClearVar (lb);
01413   ILL_RETURN (rval, "float128_forcing_constraints");
01414 }
01415 
01416 static int float128_singleton_columns (
01417   float128_graph * G,
01418   float128_ILLlp_predata * pre,
01419   int *hit)
01420 {
01421   int rval = 0;
01422   int ncols = G->ncols;
01423   int j, k, deg, rdeg, single = 0, irow;
01424   float128 lb, ub, b, eb;
01425   float128_node *cols = G->cols;
01426   float128_node *rows = G->rows;
01427   float128_edge *b_edge;
01428   float128_ILLlp_preop *op = 0;
01429   float128 newub, newlb;
01430   float128 a, c, l, u;
01431 
01432   float128_EGlpNumInitVar (lb);
01433   float128_EGlpNumInitVar (ub);
01434   float128_EGlpNumInitVar (eb);
01435   float128_EGlpNumInitVar (b);
01436   float128_EGlpNumInitVar (newlb);
01437   float128_EGlpNumInitVar (newub);
01438   float128_EGlpNumInitVar (a);
01439   float128_EGlpNumInitVar (c);
01440   float128_EGlpNumInitVar (l);
01441   float128_EGlpNumInitVar (u);
01442 
01443   *hit = 0;
01444   if (G->ncols == 0)
01445     goto CLEANUP;
01446 
01447   for (j = 0; j < ncols; j++)
01448   {
01449     if (cols[j].del == 0)
01450     {
01451       deg = 0;
01452       for (k = 0; k < cols[j].deg && deg <= 1; k++)
01453       {
01454         if (cols[j].adj[k]->del == 0)
01455         {
01456           single = k;
01457           deg++;
01458         }
01459       }
01460       if (deg == 1)
01461       {
01462         irow = cols[j].adj[single]->row;
01463         float128_EGlpNumCopy (b, cols[j].adj[single]->coef);
01464         b_edge = cols[j].adj[single];
01465 
01466         float128_get_implied_variable_bounds (G, j, b_edge, &lb, &ub);
01467 
01468         /*if (lb >= cols[j].lower && ub <= cols[j].upper) */
01469         if (float128_EGlpNumIsLeq (cols[j].lower, lb) &&
01470             float128_EGlpNumIsLeq (ub, cols[j].upper))
01471         {
01472           float128_edge *a_edge;
01473 
01474           /*  The jth variable can be substituted out of problem */
01475           /*        x = (c/b) - (a/b)y                           */
01476 
01477 
01478           rval = float128_get_next_preop (pre, &op);
01479           ILL_CLEANUP_IF (rval);
01480 
01481           op->colindex = j;
01482           op->rowindex = irow;
01483           op->ptype = float128_ILL_PRE_DELETE_FREE_SINGLETON_VARIABLE;
01484 
01485           rval = float128_grab_lp_line (G, irow, &op->line, 0);
01486           ILL_CLEANUP_IF (rval);
01487           pre->opcount++;
01488 
01489           /*  Adjust the objective function                      */
01490           /*     dy ==> (d - (e/b))ay   (e is obj coef of y)     */
01491           /*eb = cols[j].obj / b; */
01492           float128_EGlpNumCopyFrac (eb, cols[j].obj, b);
01493 
01494           for (k = 0; k < rows[irow].deg; k++)
01495           {
01496             a_edge = rows[irow].adj[k];
01497             if (a_edge->del == 0 && a_edge != b_edge)
01498             {
01499               /*cols[a_edge->col].obj -= (eb * a_edge->coef); */
01500               float128_EGlpNumSubInnProdTo (cols[a_edge->col].obj, eb, a_edge->coef);
01501             }
01502           }
01503 
01504 
01505           /*  Delete y from float128_graph */
01506 
01507           cols[j].del = 1;
01508 
01509           /*  Delete equation ay + bx = c */
01510 
01511           rows[irow].del = 1;
01512           for (k = 0; k < rows[irow].deg; k++)
01513           {
01514             rows[irow].adj[k]->del = 1;
01515           }
01516 
01517         }
01518         else
01519         {
01520           rdeg = 0;
01521           for (k = 0; k < rows[irow].deg && rdeg <= 2; k++)
01522           {
01523             if (rows[irow].adj[k]->del == 0)
01524             {
01525               rdeg++;
01526             }
01527           }
01528           if (rdeg == 2)
01529           {
01530             float128_edge *a_edge = 0;
01531             int col2 = 0;
01532 
01533             float128_EGlpNumCopy (newub, float128_ILL_MAXDOUBLE);
01534             float128_EGlpNumCopy (newlb, float128_ILL_MINDOUBLE);
01535             float128_EGlpNumZero (a);
01536 
01537             /*    ay + bx = c                                */
01538             /*    l <= x <= u                                */
01539             /*      x - is column singleton                  */
01540             /*      derive bounds on y and substitute out x  */
01541 
01542             float128_EGlpNumCopy (c, rows[irow].rhs);
01543             float128_EGlpNumCopy (l, cols[j].lower);
01544             float128_EGlpNumCopy (u, cols[j].upper);
01545 
01546             /* Find the ay term */
01547 
01548             for (k = 0; k < rows[irow].deg; k++)
01549             {
01550               if (rows[irow].adj[k]->del == 0 && rows[irow].adj[k]->col != j)
01551               {
01552                 a_edge = rows[irow].adj[k];
01553                 float128_EGlpNumCopy (a, rows[irow].adj[k]->coef);
01554                 col2 = rows[irow].adj[k]->col;
01555                 break;
01556               }
01557             }
01558             if (k == rows[irow].deg)
01559             {
01560               fprintf (stderr, "float128_graph error in singleton_col\n");
01561               rval = 1;
01562               goto CLEANUP;
01563             }
01564 
01565             /*  Record the operation             */
01566             /*  x is column j,  y is column col2 */
01567 
01568             rval = float128_get_next_preop (pre, &op);
01569             ILL_CLEANUP_IF (rval);
01570 
01571             op->colindex = j;
01572             op->rowindex = irow;
01573             op->ptype = float128_ILL_PRE_DELETE_SINGLETON_VARIABLE;
01574 
01575             rval = float128_grab_lp_line (G, irow, &op->line, 0);
01576             ILL_CLEANUP_IF (rval);
01577             pre->opcount++;
01578 
01579             /*  Adjust the bounds on y           */
01580             /*  Using x = c/b - (a/b)y            */
01581             /* we use eb as temporal variable here */
01582             /*if (a / b > 0) */
01583             float128_EGlpNumCopyFrac (eb, a, b);
01584             if (float128_EGlpNumIsGreatZero (eb))
01585             {
01586               /*if (l > -float128_ILL_MAXDOUBLE) */
01587               if (float128_EGlpNumIsLess (float128_ILL_MINDOUBLE, l))
01588               {
01589                 /*newub = (c / a) - (l * b) / a; */
01590                 float128_EGlpNumCopy (newub, c);
01591                 float128_EGlpNumSubInnProdTo (newub, l, b);
01592                 float128_EGlpNumDivTo (newub, a);
01593               }
01594               /*if (u < float128_ILL_MAXDOUBLE) */
01595               if (float128_EGlpNumIsLess (u, float128_ILL_MAXDOUBLE))
01596               {
01597                 /*newlb = (c / a) - (u * b) / a; */
01598                 float128_EGlpNumCopy (newlb, c);
01599                 float128_EGlpNumSubInnProdTo (newlb, u, b);
01600                 float128_EGlpNumDivTo (newlb, a);
01601               }
01602             }
01603             else
01604             {
01605               /*if (l > -float128_ILL_MAXDOUBLE) */
01606               if (float128_EGlpNumIsLess (float128_ILL_MINDOUBLE, l))
01607               {
01608                 /*newlb = (c / a) - (l * b) / a; */
01609                 float128_EGlpNumCopy (newlb, c);
01610                 float128_EGlpNumSubInnProdTo (newlb, l, b);
01611                 float128_EGlpNumDivTo (newlb, a);
01612               }
01613               /*if (u < float128_ILL_MAXDOUBLE) */
01614               if (float128_EGlpNumIsLess (u, float128_ILL_MAXDOUBLE))
01615               {
01616                 /*newub = (c / a) - (u * b) / a; */
01617                 float128_EGlpNumCopy (newub, c);
01618                 float128_EGlpNumSubInnProdTo (newub, u, b);
01619                 float128_EGlpNumDivTo (newub, a);
01620               }
01621             }
01622 
01623             if (float128_EGlpNumIsLess (cols[col2].lower, newlb))
01624               float128_EGlpNumCopy (cols[col2].lower, newlb);
01625             if (float128_EGlpNumIsLess (newub, cols[col2].upper))
01626               float128_EGlpNumCopy (cols[col2].upper, newub);
01627             float128_EGlpNumSubTo (cols[col2].obj, eb);
01628 
01629             /*  Delete x (and the bx term) from float128_graph */
01630 
01631             cols[j].del = 1;
01632             b_edge->del = 1;
01633 
01634             /*  Delete equation ay + bx = c (and the ax term) */
01635 
01636             rows[irow].del = 1;
01637             a_edge->del = 1;
01638           }
01639         }
01640       }
01641     }
01642   }
01643 
01644 
01645 CLEANUP:
01646 
01647   float128_EGlpNumClearVar (lb);
01648   float128_EGlpNumClearVar (ub);
01649   float128_EGlpNumClearVar (eb);
01650   float128_EGlpNumClearVar (b);
01651   float128_EGlpNumClearVar (newlb);
01652   float128_EGlpNumClearVar (newub);
01653   float128_EGlpNumClearVar (a);
01654   float128_EGlpNumClearVar (c);
01655   float128_EGlpNumClearVar (l);
01656   float128_EGlpNumClearVar (u);
01657   ILL_RETURN (rval, "float128_singleton_columns");
01658 }
01659 
01660 static int float128_duplicate_rows (
01661   float128_graph * G,
01662   int *hit)
01663 {
01664   int rval = 0;
01665   float128_node *cols = G->cols;
01666   float128_node *rows = G->rows;
01667   int ncols = G->ncols;
01668   int nrows = G->nrows;
01669   int *s = 0;
01670   float128 *f = 0;
01671   double szeit = ILLutil_zeit ();
01672   float128 q;
01673   int i, j, k, k2, ri, r0 = 0, n, nu = 0, got, t0, t = 1;
01674   float128_node *c;
01675 
01676   float128_EGlpNumInitVar (q);
01677 
01678 
01679   /*  Code follows J. Tomlin and J. S. Welch, OR Letters 5 (1986) 7--11 */
01680 
01681   *hit = 0;
01682   if (nrows == 0)
01683     goto CLEANUP;
01684 
01685   ILL_SAFE_MALLOC (s, nrows, int);
01686 
01687   f = float128_EGlpNumAllocArray (nrows);
01688 
01689   for (i = 0; i < nrows; i++)
01690   {
01691     if (rows[i].del || rows[i].rowsense != 'E')
01692     {
01693       s[i] = float128_ILL_MAXINT;       /* float128_ILL_MAXINT means no longer eligible    */
01694     }
01695     else
01696     {
01697       s[i] = 0;                 /* 0 means eligible, >0 means in a group */
01698       nu++;                     /* Tracks the number of eligible rows    */
01699     }
01700   }
01701 
01702   for (j = 0; j < ncols; j++)
01703   {
01704     c = &cols[j];
01705     if (c->del)
01706       continue;
01707     if (c->coltype != float128_ILL_PRE_COL_STRUC)
01708       continue;
01709 
01710     n = 0;
01711     t0 = t++;
01712 
01713     for (k = 0; k < c->deg; k++)
01714     {
01715       if (c->adj[k]->del)
01716         continue;
01717 
01718       ri = c->adj[k]->row;
01719       if (s[ri] == 0)
01720       {
01721         s[ri] = t0;
01722         float128_EGlpNumCopy (f[ri], c->adj[k]->coef);
01723         r0 = ri;
01724         n++;
01725       }
01726       else if (s[ri] < t0)
01727       {
01728         got = 0;
01729         for (k2 = k + 1; k2 < c->deg; k2++)
01730         {
01731           if (c->adj[k2]->del)
01732             continue;
01733 
01734           i = c->adj[k2]->row;
01735           if (s[i] == s[ri])
01736           {
01737             /*q = (c->adj[k]->coef * (f[i])) / (f[ri] * (c->adj[k2]->coef)); */
01738             float128_EGlpNumCopy (q, c->adj[k]->coef);
01739             float128_EGlpNumMultTo (q, f[i]);
01740             float128_EGlpNumDivTo (q, f[ri]);
01741             float128_EGlpNumDivTo (q, c->adj[k2]->coef);
01742             if (float128_EGlpNumIsEqual (q, float128_oneLpNum, float128_ILL_PRE_ZERO_TOL))
01743             {
01744               s[ri] = t;
01745               s[i] = t;
01746               got++;
01747             }
01748           }
01749         }
01750         if (got)
01751         {
01752           t++;
01753         }
01754         else
01755         {
01756           s[ri] = float128_ILL_MAXINT;
01757           if (--nu == 0)
01758             goto DONE;
01759         }
01760       }
01761     }
01762 
01763     if (n == 1)
01764     {
01765       s[r0] = float128_ILL_MAXINT;
01766       if (--nu == 0)
01767         goto DONE;
01768     }
01769   }
01770 
01771 DONE:
01772 
01773   {
01774     int idup = 0;
01775 
01776     for (i = 0; i < nrows; i++)
01777     {
01778       if (s[i] > 0 && s[i] < float128_ILL_MAXINT)
01779       {
01780         printf ("Row %d: %d\n", i, s[i]);
01781         idup++;
01782       }
01783     }
01784     printf ("Number of duplicate rows: %d\n", idup);
01785   }
01786 
01787   printf ("Time in float128_duplicate_rows: %.2f (seconds)\n", ILLutil_zeit () - szeit);
01788   fflush (stdout);
01789 
01790 CLEANUP:
01791 
01792   ILL_IFFREE (s, int);
01793 
01794   float128_EGlpNumFreeArray (f);
01795   float128_EGlpNumClearVar (q);
01796   ILL_RETURN (rval, "float128_duplicate_rows");
01797 }
01798 
01799 static int float128_duplicate_cols (
01800   float128_graph * G,
01801   int *hit)
01802 {
01803   int rval = 0;
01804   float128_node *cols = G->cols;
01805   float128_node *rows = G->rows;
01806   int ncols = G->ncols;
01807   int nrows = G->nrows;
01808   int *s = 0;
01809   float128 *f = 0;
01810   double szeit = ILLutil_zeit ();
01811   float128 q;
01812   int i, j, k, k2, ci, c0 = 0, n, nu = 0, got, t0, t = 1;
01813   float128_node *r;
01814 
01815   float128_EGlpNumInitVar (q);
01816 
01817 
01818   /*  Code follows J. Tomlin and J. S. Welch, OR Letters 5 (1986) 7--11 */
01819 
01820   *hit = 0;
01821   if (ncols == 0)
01822     goto CLEANUP;
01823 
01824   ILL_SAFE_MALLOC (s, ncols, int);
01825 
01826   f = float128_EGlpNumAllocArray (ncols);
01827 
01828   for (j = 0; j < ncols; j++)
01829   {
01830     if (cols[j].del || cols[j].coltype != float128_ILL_PRE_COL_STRUC)
01831     {
01832       s[j] = float128_ILL_MAXINT;       /* float128_ILL_MAXINT means no longer eligible    */
01833     }
01834     else
01835     {
01836       s[j] = 0;                 /* 0 means eligible, >0 means in a group */
01837       nu++;                     /* Tracks the number of eligible rows    */
01838     }
01839   }
01840 
01841   for (i = 0; i < nrows; i++)
01842   {
01843     r = &rows[i];
01844     if (r->del)
01845       continue;
01846 
01847     n = 0;
01848     t0 = t++;
01849 
01850     for (k = 0; k < r->deg; k++)
01851     {
01852       if (r->adj[k]->del)
01853         continue;
01854 
01855       ci = r->adj[k]->col;
01856       if (s[ci] == 0)
01857       {
01858         s[ci] = t0;
01859         float128_EGlpNumCopy (f[ci], r->adj[k]->coef);
01860         c0 = ci;
01861         n++;
01862       }
01863       else if (s[ci] < t0)
01864       {
01865         got = 0;
01866         for (k2 = k + 1; k2 < r->deg; k2++)
01867         {
01868           if (r->adj[k2]->del)
01869             continue;
01870 
01871           j = r->adj[k2]->col;
01872           if (s[j] == s[ci])
01873           {
01874             /*q = (r->adj[k]->coef * (f[j])) / (f[ci] * (r->adj[k2]->coef)); */
01875             float128_EGlpNumCopy (q, r->adj[k]->coef);
01876             float128_EGlpNumMultTo (q, f[j]);
01877             float128_EGlpNumDivTo (q, f[ci]);
01878             float128_EGlpNumDivTo (q, r->adj[k2]->coef);
01879             if (float128_EGlpNumIsEqual (q, float128_oneLpNum, float128_ILL_PRE_ZERO_TOL))
01880             {
01881               s[ci] = t;
01882               s[j] = t;
01883               got++;
01884             }
01885           }
01886         }
01887         if (got)
01888         {
01889           t++;
01890         }
01891         else
01892         {
01893           s[ci] = float128_ILL_MAXINT;
01894           if (--nu == 0)
01895             goto DONE;
01896         }
01897       }
01898     }
01899 
01900     if (n == 1)
01901     {
01902       s[c0] = float128_ILL_MAXINT;
01903       if (--nu == 0)
01904         goto DONE;
01905     }
01906   }
01907 
01908 DONE:
01909 
01910   {
01911     int dcount;
01912     int *dcnt;
01913     int *dlist;
01914 
01915     rval = float128_gather_dup_lists (s, ncols, &dcount, &dcnt, &dlist);
01916     ILL_CLEANUP_IF (rval);
01917   }
01918 
01919   printf ("Time in float128_duplicate_cols: %.2f (seconds)\n", ILLutil_zeit () - szeit);
01920   fflush (stdout);
01921 
01922 CLEANUP:
01923 
01924   ILL_IFFREE (s, int);
01925 
01926   float128_EGlpNumFreeArray (f);
01927   float128_EGlpNumClearVar (q);
01928   ILL_RETURN (rval, "float128_duplicate_cols");
01929 }
01930 
01931 static int float128_gather_dup_lists (
01932   /* float128_graph *G, */ int *s,
01933   /* double *f, */
01934   int count,
01935   int *duptotal,
01936   int **dupcnt,
01937   int **dupind)
01938 {
01939   int rval = 0;
01940   int *cnt = 0;
01941   int *ind = 0;
01942   int *beg = 0;
01943   int i, smax = 0, ndup = 0, total = 0;
01944 
01945   *duptotal = 0;
01946   *dupcnt = 0;
01947   *dupind = 0;
01948 
01949 
01950   for (i = 0; i < count; i++)
01951   {
01952     if (s[i] < float128_ILL_MAXINT && s[i] > smax)
01953       smax = s[i];
01954   }
01955   if (smax == 0)
01956     goto CLEANUP;
01957 
01958   ILL_SAFE_MALLOC (cnt, smax + 1, int);
01959 
01960   ILL_SAFE_MALLOC (ind, smax + 1, int);
01961 
01962   for (i = 0; i < smax + 1; i++)
01963   {
01964     cnt[i] = 0;
01965   }
01966 
01967   for (i = 0; i < count; i++)
01968   {
01969     if (s[i] < float128_ILL_MAXINT)
01970     {
01971       cnt[s[i]]++;
01972     }
01973   }
01974 
01975   if (cnt[0] > 0)
01976     printf ("%d Empty Lines\n", cnt[0]);
01977 
01978   printf ("Duplicate Classes:");
01979   fflush (stdout);
01980   for (i = 1; i < smax + 1; i++)
01981   {
01982     if (cnt[i] > 1)
01983     {
01984       ndup++;
01985       printf (" %d", cnt[i]);
01986     }
01987   }
01988   printf ("  Number %d\n", ndup);
01989   fflush (stdout);
01990 
01991   if (ndup == 0)
01992     goto CLEANUP;
01993 
01994   ILL_SAFE_MALLOC (beg, ndup, int);
01995 
01996   for (i = 1, ndup = 0; i < smax + 1; i++)
01997   {
01998     if (cnt[i] > 1)
01999     {
02000       beg[ndup] = total;
02001       total += cnt[i];
02002       ind[i] = ndup;
02003       ndup++;
02004     }
02005   }
02006 
02007   if (total == 0)
02008     goto CLEANUP;
02009 
02010   ILL_SAFE_MALLOC (*dupcnt, ndup, int);
02011 
02012   ILL_SAFE_MALLOC (*dupind, total, int);
02013 
02014   for (i = 0; i < ndup; i++)
02015   {
02016     (*dupcnt)[i] = 0;
02017   }
02018 
02019   for (i = 0; i < count; i++)
02020   {
02021     if (s[i] < float128_ILL_MAXINT && s[i] > 0)
02022     {
02023       if (cnt[s[i]] > 1)
02024       {
02025         (*dupind)[beg[ind[s[i]]] + (*dupcnt)[ind[s[i]]]] = i;
02026         (*dupcnt)[ind[s[i]]]++;
02027       }
02028     }
02029   }
02030 
02031   for (i = 0; i < ndup; i++)
02032   {
02033     int j;
02034 
02035     for (j = beg[i]; j < beg[i] + (*dupcnt)[i]; j++)
02036     {
02037       printf (" %d", (*dupind)[j]);
02038     }
02039     printf (" | ");
02040     fflush (stdout);
02041   }
02042 
02043   *duptotal = ndup;
02044 
02045 CLEANUP:
02046 
02047   ILL_IFFREE (cnt, int);
02048   ILL_IFFREE (ind, int);
02049   ILL_IFFREE (beg, int);
02050 
02051   ILL_RETURN (rval, "float128_gather_dup_lists");
02052 }
02053 
02054 static void float128_set_fixed_variable (
02055   float128_graph * G,
02056   int j,
02057   float128 val)
02058 {
02059   int k;
02060   float128_edge *e;
02061 
02062   G->cols[j].del = 1;
02063   for (k = 0; k < G->cols[j].deg; k++)
02064   {
02065     e = G->cols[j].adj[k];
02066     if (e->del == 0)
02067     {
02068       /*G->rows[e->row].rhs -= (e->coef * val); */
02069       float128_EGlpNumSubInnProdTo (G->rows[e->row].rhs, e->coef, val);
02070       e->del = 1;
02071     }
02072   }
02073 }
02074 
02075 static void float128_get_implied_rhs_bounds (
02076   float128_graph * G,
02077   int i,
02078   float128 * lb,
02079   float128 * ub)
02080 {
02081   int k;
02082   float128 l, u;
02083   float128_node *cols = G->cols;
02084   float128_node *rows = G->rows;
02085   float128_edge *e;
02086 
02087   float128_EGlpNumInitVar (u);
02088   float128_EGlpNumInitVar (l);
02089 
02090   float128_EGlpNumZero (l);
02091   for (k = 0; k < rows[i].deg; k++)
02092   {
02093     e = rows[i].adj[k];
02094     if (e->del == 0)
02095     {
02096       if (float128_EGlpNumIsLessZero (e->coef))
02097       {
02098         if (float128_EGlpNumIsEqqual (cols[e->col].upper, float128_ILL_MAXDOUBLE))
02099         {
02100           float128_EGlpNumCopy (l, float128_ILL_MINDOUBLE);
02101           break;
02102         }
02103         else
02104         {
02105           /*l += (e->coef * cols[e->col].upper); */
02106           float128_EGlpNumAddInnProdTo (l, e->coef, cols[e->col].upper);
02107         }
02108       }
02109       else if (float128_EGlpNumIsGreatZero (e->coef))
02110       {
02111         if (float128_EGlpNumIsEqqual (cols[e->col].lower, float128_ILL_MINDOUBLE))
02112         {
02113           float128_EGlpNumCopy (l, float128_ILL_MINDOUBLE);
02114           break;
02115         }
02116         else
02117         {
02118           /*l += (e->coef * cols[e->col].lower); */
02119           float128_EGlpNumAddInnProdTo (l, e->coef, cols[e->col].lower);
02120         }
02121       }
02122     }
02123   }
02124 
02125   float128_EGlpNumZero (u);
02126   for (k = 0; k < rows[i].deg; k++)
02127   {
02128     e = rows[i].adj[k];
02129     if (e->del == 0)
02130     {
02131       if (float128_EGlpNumIsLessZero (e->coef ))
02132       {
02133         if (float128_EGlpNumIsEqqual (cols[e->col].lower, float128_ILL_MINDOUBLE))
02134         {
02135           float128_EGlpNumCopy (u, float128_ILL_MAXDOUBLE);
02136         }
02137         else
02138         {
02139           /*u += (e->coef * cols[e->col].lower); */
02140           float128_EGlpNumAddInnProdTo (u, e->coef, cols[e->col].lower);
02141         }
02142       }
02143       else if (float128_EGlpNumIsGreatZero (e->coef))
02144       {
02145         if (float128_EGlpNumIsEqqual (cols[e->col].upper, float128_ILL_MAXDOUBLE))
02146         {
02147           float128_EGlpNumCopy (u, float128_ILL_MAXDOUBLE);
02148         }
02149         else
02150         {
02151           /*u += (e->coef * cols[e->col].upper); */
02152           float128_EGlpNumAddInnProdTo (u, e->coef, cols[e->col].upper);
02153         }
02154       }
02155     }
02156   }
02157 
02158   float128_EGlpNumCopy (*lb, l);
02159   float128_EGlpNumCopy (*ub, u);
02160   float128_EGlpNumClearVar (u);
02161   float128_EGlpNumClearVar (l);
02162 }
02163 
02164 static void float128_get_implied_variable_bounds (
02165   float128_graph * G,
02166   int j,
02167   float128_edge * a_ij,
02168   float128 * lb,
02169   float128 * ub)
02170 {
02171   int i = a_ij->row;
02172   float128 l, u;
02173 
02174   float128_EGlpNumInitVar (u);
02175   float128_EGlpNumInitVar (l);
02176 
02177   float128_get_implied_rhs_bounds (G, i, &l, &u);
02178   float128_EGlpNumCopy (*lb, float128_ILL_MINDOUBLE);
02179   float128_EGlpNumCopy (*ub, float128_ILL_MAXDOUBLE);
02180 
02181   if (float128_EGlpNumIsLess (float128_ILL_PRE_FEAS_TOL, a_ij->coef))
02182   {
02183     if (float128_EGlpNumIsLess (u, float128_ILL_MAXDOUBLE))
02184     {
02185       /**lb = (G->rows[i].rhs - u) / a_ij->coef + G->cols[j].upper;*/
02186       float128_EGlpNumCopyDiffRatio (*lb, G->rows[i].rhs, u, a_ij->coef);
02187       float128_EGlpNumAddTo (*lb, G->cols[j].upper);
02188     }
02189     if (float128_EGlpNumIsLess (float128_ILL_MINDOUBLE, l))
02190     {
02191       /**ub = (G->rows[i].rhs - l) / a_ij->coef + G->cols[j].lower;*/
02192       float128_EGlpNumCopyDiffRatio (*ub, G->rows[i].rhs, l, a_ij->coef);
02193       float128_EGlpNumAddTo (*ub, G->cols[j].lower);
02194     }
02195   }
02196   else if (float128_EGlpNumIsLess (a_ij->coef, float128_ILL_PRE_FEAS_TOL))
02197   {
02198     if (float128_EGlpNumIsLess (float128_ILL_MINDOUBLE, l))
02199     {
02200       /**lb = (G->rows[i].rhs - l) / a_ij->coef + G->cols[j].upper;*/
02201       float128_EGlpNumCopyDiffRatio (*lb, G->rows[i].rhs, l, a_ij->coef);
02202       float128_EGlpNumAddTo (*lb, G->cols[j].upper);
02203     }
02204     if (float128_EGlpNumIsLess (u, float128_ILL_MAXDOUBLE))
02205     {
02206       /**ub = (G->rows[i].rhs - u) / a_ij->coef + G->cols[j].lower;*/
02207       float128_EGlpNumCopyDiffRatio (*ub, G->rows[i].rhs, u, a_ij->coef);
02208       float128_EGlpNumAddTo (*ub, G->cols[j].lower);
02209     }
02210   }
02211   float128_EGlpNumClearVar (u);
02212   float128_EGlpNumClearVar (l);
02213 }
02214 
02215 static int float128_get_next_preop (
02216   float128_ILLlp_predata * pre,
02217   float128_ILLlp_preop ** op)
02218 {
02219   int rval = 0;
02220 
02221   if (pre->opcount >= pre->opsize)
02222   {
02223     pre->opsize *= 1.3;
02224     pre->opsize += 1000;
02225     if (pre->opsize < pre->opcount + 1)
02226       pre->opsize = pre->opcount + 1;
02227     pre->oplist = EGrealloc (pre->oplist, sizeof (float128_ILLlp_preop) * pre->opsize);
02228     //rval = ILLutil_reallocrus_scale ((void **) &pre->oplist,
02229     //                                 &pre->opsize, pre->opcount + 1, 1.3,
02230     //                                 sizeof (float128_ILLlp_preop));
02231     //ILL_CLEANUP_IF (rval);
02232   }
02233   *op = &pre->oplist[pre->opcount];
02234   float128_ILLlp_preop_init (*op);
02235 
02236 //CLEANUP:
02237 
02238   ILL_RETURN (rval, "float128_get_next_preop");
02239 }
02240 
02241 static int float128_add_to_list (
02242   ILLptrworld * world,
02243   float128_intptr ** list,
02244   int i)
02245 {
02246   int rval = 0;
02247   float128_intptr *ip;
02248 
02249   ip = intptralloc (world);
02250   if (!ip)
02251   {
02252     rval = 1;
02253     goto CLEANUP;
02254   }
02255   ip->this_val = i;
02256   ip->next = *list;
02257   *list = ip;
02258 
02259 CLEANUP:
02260 
02261   ILL_RETURN (rval, "float128_add_to_list");
02262 }
02263 
02264 static int float128_build_graph (
02265   float128_ILLlpdata * lp,
02266   float128_graph * G)
02267 {
02268   int rval = 0;
02269   int ncols = lp->ncols;
02270   int nrows = lp->nrows;
02271   int nzcount = lp->nzcount;
02272   int i, j, k, stop, count;
02273   float128_edge *edgelist;
02274   float128_node *rows, *cols;
02275   float128_ILLmatrix *A = &lp->A;
02276 
02277   G->objsense = lp->objsense;
02278 
02279   ILL_SAFE_MALLOC (G->rows, nrows, float128_node);
02280   if (!G->rows)
02281   {
02282     fprintf (stderr, "out of memory in float128_build_graph\n");
02283     rval = 1;
02284     goto CLEANUP;
02285   }
02286   rows = G->rows;
02287 
02288   for (i = 0; i < nrows; i++)
02289   {
02290     rows[i].rowsense = lp->sense[i];
02291     rows[i].deg = 0;
02292   }
02293 
02294   ILL_SAFE_MALLOC (G->cols, ncols, float128_node);
02295   ILL_SAFE_MALLOC (G->edgelist, nzcount, float128_edge);
02296   for (i = nzcount; i--;)
02297     float128_EGlpNumInitVar ((G->edgelist[i].coef));
02298   G->nzcount = nzcount;
02299   ILL_SAFE_MALLOC (G->adjspace, 2 * nzcount, float128_edge *);
02300 
02301   if (!G->cols || !G->edgelist || !G->adjspace)
02302   {
02303     fprintf (stderr, "out of memory in float128_build_graph\n");
02304     rval = 1;
02305     goto CLEANUP;
02306   }
02307 
02308   cols = G->cols;
02309   edgelist = G->edgelist;
02310 
02311   for (j = 0; j < ncols; j++)
02312   {
02313     stop = A->matbeg[j] + A->matcnt[j];
02314     for (k = A->matbeg[j]; k < stop; k++)
02315     {
02316       rows[A->matind[k]].deg++;
02317     }
02318   }
02319 
02320   for (i = 0, count = 0; i < nrows; i++)
02321   {
02322     rows[i].adj = G->adjspace + count;
02323     count += rows[i].deg;
02324     rows[i].deg = 0;
02325   }
02326 
02327   for (j = 0; j < ncols; j++)
02328   {
02329     cols[j].adj = G->adjspace + count;
02330     count += A->matcnt[j];
02331     cols[j].deg = 0;
02332     cols[j].coltype = float128_ILL_PRE_COL_STRUC;
02333   }
02334   for (i = 0; i < nrows; i++)
02335   {
02336     cols[lp->rowmap[i]].coltype = float128_ILL_PRE_COL_LOGICAL;
02337   }
02338 
02339   for (j = 0, count = 0; j < ncols; j++)
02340   {
02341     float128_EGlpNumCopy (cols[j].obj, lp->obj[j]);
02342     float128_EGlpNumCopy (cols[j].lower, lp->lower[j]);
02343     float128_EGlpNumCopy (cols[j].upper, lp->upper[j]);
02344     stop = A->matbeg[j] + A->matcnt[j];
02345     for (k = A->matbeg[j]; k < stop; k++)
02346     {
02347       i = A->matind[k];
02348       rows[i].adj[rows[i].deg++] = &(edgelist[count]);
02349       cols[j].adj[cols[j].deg++] = &(edgelist[count]);
02350       edgelist[count].row = i;
02351       edgelist[count].col = j;
02352       float128_EGlpNumCopy (edgelist[count].coef, A->matval[k]);
02353       edgelist[count].mark = 0;
02354       edgelist[count].del = 0;
02355       edgelist[count].coltype = cols[j].coltype;
02356       count++;
02357     }
02358   }
02359   if (count != nzcount)
02360   {
02361     fprintf (stderr, "counts are off in float128_build_graph\n");
02362     rval = 1;
02363     goto CLEANUP;
02364   }
02365 
02366   G->ecount = count;
02367   G->nrows = nrows;
02368   G->ncols = ncols;
02369 
02370   for (i = 0; i < G->nrows; i++)
02371   {
02372     G->rows[i].del = 0;
02373     float128_EGlpNumCopy (G->rows[i].rhs, lp->rhs[i]);
02374   }
02375   for (j = 0; j < G->ncols; j++)
02376   {
02377     G->cols[j].del = 0;
02378   }
02379 
02380 CLEANUP:
02381 
02382   ILL_RETURN (rval, "float128_build_graph");
02383 }
02384 
02385 static void float128_dump_graph (
02386   float128_graph * G)
02387 {
02388   int i, j, k;
02389 
02390   printf ("ecount = %d, nrows = %d, ncols = %d\n",
02391           G->ecount, G->nrows, G->ncols);
02392   fflush (stdout);
02393 
02394   for (i = 0; i < G->nrows; i++)
02395   {
02396     printf ("Row %d:", i);
02397     for (k = 0; k < G->rows[i].deg; k++)
02398     {
02399       printf (" %d", G->rows[i].adj[k]->col);
02400       if (G->rows[i].adj[k]->coltype == float128_ILL_PRE_COL_LOGICAL)
02401         printf ("S");
02402       printf ("(%g)", float128_EGlpNumToLf (G->rows[i].adj[k]->coef));
02403     }
02404     printf ("  rhs: %g", float128_EGlpNumToLf (G->rows[i].rhs));
02405     if (G->rows[i].del)
02406     {
02407       printf (" (deleted)\n");
02408     }
02409     else
02410     {
02411       printf ("\n");
02412     }
02413   }
02414 
02415   for (j = 0; j < G->ncols; j++)
02416   {
02417     if (G->cols[j].coltype == float128_ILL_PRE_COL_LOGICAL)
02418     {
02419       printf ("Slk %d:", j);
02420     }
02421     else
02422     {
02423       printf ("Col %d:", j);
02424     }
02425     for (k = 0; k < G->cols[j].deg; k++)
02426     {
02427       printf (" %d", G->cols[j].adj[k]->row);
02428     }
02429     printf ("  obj: %g  bnd: (%g, %g)", float128_EGlpNumToLf (G->cols[j].obj),
02430             float128_EGlpNumToLf (G->cols[j].lower), float128_EGlpNumToLf (G->cols[j].upper));
02431     if (G->cols[j].del)
02432     {
02433       printf (" (deleted)\n");
02434     }
02435     else
02436     {
02437       printf ("\n");
02438     }
02439   }
02440 }
02441 
02442 static void float128_init_graph (
02443   float128_graph * G)
02444 {
02445   if (G)
02446   {
02447     G->edgelist = 0;
02448     G->rows = 0;
02449     G->cols = 0;
02450     G->ecount = 0;
02451     G->nrows = 0;
02452     G->ncols = 0;
02453     G->adjspace = 0;
02454     ILLptrworld_init (&G->intptrworld);
02455   }
02456 }
02457 
02458 static void float128_free_graph (
02459   float128_graph * G)
02460 {
02461   register int i;
02462 
02463   if (G)
02464   {
02465     int total, onlist;
02466 
02467     for (i = G->nzcount; i--;)
02468       float128_EGlpNumClearVar ((G->edgelist[i].coef));
02469     ILL_IFFREE (G->edgelist, float128_edge);
02470     ILL_IFFREE (G->rows, float128_node);
02471     ILL_IFFREE (G->cols, float128_node);
02472     ILL_IFFREE (G->adjspace, float128_edge *);
02473     if (intptr_check_leaks (&G->intptrworld, &total, &onlist))
02474     {
02475       fprintf (stderr, "WARNING: %d outstanding intptrs\n", total - onlist);
02476     }
02477     ILLptrworld_delete (&G->intptrworld);
02478     float128_init_graph (G);
02479   }
02480 }
02481 
02482 int float128_ILLlp_sinfo_print (
02483   float128_ILLlp_sinfo * s)
02484 {
02485   int rval = 0;
02486   int i;
02487   float128_ILLlpdata lp;
02488   char *sense = 0;
02489 
02490   float128_ILLlpdata_init (&lp);
02491 
02492   lp.nrows = s->nrows;
02493   lp.ncols = s->ncols;
02494   lp.nzcount = s->nzcount;
02495   lp.objsense = s->objsense;
02496   lp.obj = s->obj;
02497   lp.rhs = s->rhs;
02498   lp.lower = s->lower;
02499   lp.upper = s->upper;
02500   lp.A.matval = s->A.matval;
02501   lp.A.matcnt = s->A.matcnt;
02502   lp.A.matbeg = s->A.matbeg;
02503   lp.A.matind = s->A.matind;
02504   lp.rownames = 0;
02505   lp.colnames = s->colnames;
02506   lp.objname = 0;
02507   lp.probname = 0;
02508   lp.intmarker = 0;
02509 
02510   ILL_SAFE_MALLOC (sense, s->nrows, char);
02511 
02512   if (!sense)
02513   {
02514     fprintf (stderr, "out of memory in float128_ILLlp_sinfo_print\n");
02515     rval = 1;
02516     goto CLEANUP;
02517   }
02518   for (i = 0; i < s->nrows; i++)
02519   {
02520     sense[i] = 'E';
02521   }
02522   lp.sense = sense;
02523 
02524 /*
02525     rval = ILLlpdata_writelp (&lp, 0);
02526     ILL_CLEANUP_IF (rval);
02527 */
02528 
02529 CLEANUP:
02530 
02531   ILL_IFFREE (sense, char);
02532 
02533   ILL_RETURN (rval, "float128_ILLlp_sinfo_print");
02534 }
02535 
02536 void float128_ILLlp_sinfo_init (
02537   float128_ILLlp_sinfo * sinfo)
02538 {
02539   if (sinfo)
02540   {
02541     sinfo->ncols = 0;
02542     sinfo->nrows = 0;
02543     sinfo->nzcount = 0;
02544     sinfo->rowsize = 0;
02545     sinfo->colsize = 0;
02546     sinfo->obj = 0;
02547     sinfo->rhs = 0;
02548     sinfo->lower = 0;
02549     sinfo->upper = 0;
02550     sinfo->colnames = 0;
02551     sinfo->objsense = float128_ILL_MIN;
02552     float128_ILLmatrix_init (&sinfo->A);
02553   }
02554 }
02555 
02556 void float128_ILLlp_sinfo_free (
02557   float128_ILLlp_sinfo * sinfo)
02558 {
02559   if (sinfo)
02560   {
02561     float128_EGlpNumFreeArray (sinfo->obj);
02562     float128_EGlpNumFreeArray (sinfo->lower);
02563     float128_EGlpNumFreeArray (sinfo->upper);
02564     float128_EGlpNumFreeArray (sinfo->rhs);
02565     float128_ILLmatrix_free (&sinfo->A);
02566     if (sinfo->colnames)
02567     {
02568       int i;
02569 
02570       for (i = 0; i < sinfo->ncols; i++)
02571       {
02572         ILL_IFFREE (sinfo->colnames[i], char);
02573       }
02574       ILL_IFFREE (sinfo->colnames, char *);
02575     }
02576     float128_ILLlp_sinfo_init (sinfo);
02577   }
02578 }
02579 
02580 void float128_ILLlp_predata_init (
02581   float128_ILLlp_predata * pre)
02582 {
02583   if (pre)
02584   {
02585     pre->opcount = 0;
02586     pre->opsize = 0;
02587     pre->oplist = 0;
02588     pre->r_nrows = 0;
02589     pre->r_ncols = 0;
02590     pre->colmap = 0;
02591     pre->rowmap = 0;
02592     pre->colscale = 0;
02593     pre->rowscale = 0;
02594     pre->colfixval = 0;
02595     pre->rowfixval = 0;
02596   }
02597 }
02598 
02599 void float128_ILLlp_predata_free (
02600   float128_ILLlp_predata * pre)
02601 {
02602   if (pre)
02603   {
02604     int i;
02605 
02606     for (i = 0; i < pre->opcount; i++)
02607     {
02608       float128_ILLlp_preop_free (&pre->oplist[i]);
02609     }
02610     ILL_IFFREE (pre->oplist, float128_ILLlp_preop);
02611     ILL_IFFREE (pre->colmap, int);
02612     ILL_IFFREE (pre->rowmap, int);
02613 
02614     ILL_IFFREE (pre->colscale, float128);
02615     ILL_IFFREE (pre->rowscale, float128);
02616     ILL_IFFREE (pre->colfixval, float128);
02617     ILL_IFFREE (pre->rowfixval, float128);
02618     float128_ILLlp_predata_init (pre);
02619   }
02620 }
02621 
02622 void float128_ILLlp_preop_init (
02623   float128_ILLlp_preop * op)
02624 {
02625   if (op)
02626   {
02627     op->ptype = 0;
02628     op->rowindex = -1;
02629     op->colindex = -1;
02630     float128_ILLlp_preline_init (&op->line);
02631   }
02632 }
02633 
02634 void float128_ILLlp_preop_free (
02635   float128_ILLlp_preop * op)
02636 {
02637   if (op)
02638   {
02639     float128_ILLlp_preline_free (&op->line);
02640     float128_ILLlp_preop_init (op);
02641   }
02642 }
02643 
02644 void float128_ILLlp_preline_init (
02645   float128_ILLlp_preline * line)
02646 {
02647   if (line)
02648   {
02649     float128_EGlpNumInitVar (line->rhs);
02650     float128_EGlpNumInitVar (line->obj);
02651     float128_EGlpNumInitVar (line->upper);
02652     float128_EGlpNumInitVar (line->lower);
02653     float128_EGlpNumZero (line->rhs);
02654     float128_EGlpNumZero (line->obj);
02655     float128_EGlpNumZero (line->upper);
02656     float128_EGlpNumZero (line->lower);
02657     line->count = 0;
02658     line->ind = 0;
02659     line->val = 0;
02660   }
02661 }
02662 
02663 void float128_ILLlp_preline_free (
02664   float128_ILLlp_preline * line)
02665 {
02666   if (line)
02667   {
02668     float128_EGlpNumClearVar (line->rhs);
02669     float128_EGlpNumClearVar (line->obj);
02670     float128_EGlpNumClearVar (line->upper);
02671     float128_EGlpNumClearVar (line->lower);
02672     ILL_IFFREE (line->ind, int);
02673 
02674     float128_EGlpNumFreeArray (line->val);
02675     //float128_ILLlp_preline_init (line);
02676   }
02677 }

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