factor.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: factor.c,v $ $Revision: 1.2 $ $Date: 2003/11/05 16:49:52 $"; */
00024 //static int TRACE = 0;
00025 
00026 /* implement a = max(a,abs(b)) and execute the extra code if the update is
00027  * needed */
00028 #define EGlpNumSetToMaxAbsAndDo(a,b,c) \
00029   if(EGlpNumIsGreatZero(b))\
00030   {\
00031     if(EGlpNumIsLess(a,b)){\
00032       EGlpNumCopy(a,b);\
00033       c;\
00034       }\
00035   }\
00036   else\
00037   {\
00038     EGlpNumSign(a);\
00039     if(EGlpNumIsLess(b,a)){\
00040       EGlpNumCopy(a,b);\
00041       c;\
00042       }\
00043     EGlpNumSign(a);\
00044   }
00045 
00046 
00047 #include <stdio.h>
00048 #include <stdlib.h>
00049 #include "config.h"
00050 #include "iqsutil.h"
00051 #include "lpdefs.h"
00052 #include "factor.h"
00053 #ifdef USEDMALLOC
00054 #include "dmalloc.h"
00055 #endif
00056 
00057 
00058 #undef  RECORD
00059 #undef  DEBUG_FACTOR
00060 #undef  SOLVE_DEBUG
00061 
00062 #undef  FACTOR_DEBUG
00063 #undef  UPDATE_DEBUG
00064 
00065 #undef TRACK_FACTOR
00066 #undef NOTICE_BLOWUP
00067 
00068 #undef  FACTOR_STATS
00069 #undef  UPDATE_STATS
00070 #undef  GROWTH_STATS
00071 
00072 #undef  UPDATE_STUDY
00073 
00074 #undef  SORT_RESULTS
00075 
00076 #ifdef UPDATE_STUDY
00077 int nupdate = 0;
00078 long int colspiketot = 0.0;
00079 long int rowspiketot = 0.0;
00080 long int permshifttot = 0.0;
00081 long int leftetatot = 0.0;
00082 #endif
00083 
00084 void ILLfactor_init_factor_work (
00085   factor_work * f)
00086 {
00087   f->max_k = 1000;              /* must be less than 46340 (2^15.5) */
00088   EGlpNumCopy (f->fzero_tol, SZERO_TOLER);  /* 2^-50 */
00089   EGlpNumCopy (f->szero_tol, SZERO_TOLER);  /* 2^-50 */
00090   EGlpNumCopy (f->partial_tol, OBJBND_TOLER); /* 2^-7 */
00091   f->ur_space_mul = 2.0;
00092   f->uc_space_mul = 1.1;
00093   f->lc_space_mul = 1.1;
00094   f->er_space_mul = 1000.0;
00095   f->grow_mul = 1.5;
00096   f->p = 4;
00097   f->etamax = 100;
00098   f->minmult = 1e3;
00099   f->maxmult = 1e5;
00100   f->updmaxmult = 1e7;
00101   f->dense_fract = 0.25;
00102   f->dense_min = 25;
00103   EGlpNumCopy (f->partial_cur, f->partial_tol);
00104   f->work_coef = 0;
00105   f->work_indx = 0;
00106   f->uc_inf = 0;
00107   f->ur_inf = 0;
00108   f->lc_inf = 0;
00109   f->lr_inf = 0;
00110   f->er_inf = 0;
00111   f->ucindx = 0;
00112   f->ucrind = 0;
00113   f->uccoef = 0;
00114   f->urindx = 0;
00115   f->urcind = 0;
00116   f->urcoef = 0;
00117   f->lcindx = 0;
00118   f->lccoef = 0;
00119   f->lrindx = 0;
00120   f->lrcoef = 0;
00121   f->erindx = 0;
00122   f->ercoef = 0;
00123   f->rperm = 0;
00124   f->rrank = 0;
00125   f->cperm = 0;
00126   f->crank = 0;
00127   f->dmat = 0;
00128   ILLsvector_init (&f->xtmp);
00129 }
00130 
00131 void ILLfactor_free_factor_work (
00132   factor_work * f)
00133 {
00134 #ifdef UPDATE_STUDY
00135   if (nupdate)
00136   {
00137     MESSAGE(0, "UPDATE STUDY: avg %d upd: %.2f col, %.2f row, %.2f lefteta, "
00138             "%.2f perm", nupdate, ((double) colspiketot) / nupdate, 
00139             ((double) rowspiketot) / nupdate, ((double) leftetatot) / nupdate,
00140             ((double) permshifttot) / nupdate);
00141   }
00142 #endif
00143   EGlpNumFreeArray (f->work_coef);
00144   ILL_IFFREE (f->work_indx, int);
00145 
00146   ILL_IFFREE (f->uc_inf, uc_info);
00147   if (f->dim + f->max_k > 0 && f->ur_inf)
00148   {
00149     unsigned int i = f->dim + f->max_k + 1;
00150 
00151     while (i--)
00152       EGlpNumClearVar (f->ur_inf[i].max);
00153   }
00154   ILL_IFFREE (f->ur_inf, ur_info);
00155   ILL_IFFREE (f->lc_inf, lc_info);
00156   ILL_IFFREE (f->lr_inf, lr_info);
00157   ILL_IFFREE (f->er_inf, er_info);
00158   ILL_IFFREE (f->ucindx, int);
00159   ILL_IFFREE (f->ucrind, int);
00160 
00161   EGlpNumFreeArray (f->uccoef);
00162   ILL_IFFREE (f->urindx, int);
00163   ILL_IFFREE (f->urcind, int);
00164 
00165   EGlpNumFreeArray (f->urcoef);
00166   ILL_IFFREE (f->lcindx, int);
00167 
00168   EGlpNumFreeArray (f->lccoef);
00169   ILL_IFFREE (f->lrindx, int);
00170 
00171   EGlpNumFreeArray (f->lrcoef);
00172   ILL_IFFREE (f->erindx, int);
00173 
00174   EGlpNumFreeArray (f->ercoef);
00175   ILL_IFFREE (f->rperm, int);
00176   ILL_IFFREE (f->rrank, int);
00177   ILL_IFFREE (f->cperm, int);
00178   ILL_IFFREE (f->crank, int);
00179 
00180   EGlpNumFreeArray (f->dmat);
00181   ILLsvector_free (&f->xtmp);
00182 }
00183 
00184 int ILLfactor_set_factor_iparam (
00185   factor_work * f,
00186   int param,
00187   int val)
00188 {
00189   switch (param)
00190   {
00191   case QS_FACTOR_MAX_K:
00192     f->max_k = val;
00193     break;
00194   case QS_FACTOR_P:
00195     f->p = val;
00196     break;
00197   case QS_FACTOR_ETAMAX:
00198     f->etamax = val;
00199     break;
00200   case QS_FACTOR_DENSE_MIN:
00201     f->dense_min = val;
00202     break;
00203   default:
00204     fprintf (stderr, "Invalid param %d in ILLfactor_set_factor_iparam\n",
00205              param);
00206     return 1;
00207   }
00208   return 0;
00209 }
00210 
00211 int ILLfactor_set_factor_dparam (
00212   factor_work * f,
00213   int param,
00214   EGlpNum_t val)
00215 {
00216   switch (param)
00217   {
00218   case QS_FACTOR_FZERO_TOL:
00219     EGlpNumCopy (f->fzero_tol, val);
00220     break;
00221   case QS_FACTOR_SZERO_TOL:
00222     EGlpNumCopy (f->szero_tol, val);
00223     break;
00224   case QS_FACTOR_UR_SPACE_MUL:
00225     f->ur_space_mul = EGlpNumToLf (val);
00226     break;
00227   case QS_FACTOR_UC_SPACE_MUL:
00228     f->uc_space_mul = EGlpNumToLf (val);
00229     break;
00230   case QS_FACTOR_LC_SPACE_MUL:
00231     f->lc_space_mul = EGlpNumToLf (val);
00232     break;
00233   case QS_FACTOR_LR_SPACE_MUL:
00234     f->lr_space_mul = EGlpNumToLf (val);
00235     break;
00236   case QS_FACTOR_ER_SPACE_MUL:
00237     f->er_space_mul = EGlpNumToLf (val);
00238     break;
00239   case QS_FACTOR_GROW_MUL:
00240     f->grow_mul = EGlpNumToLf (val);
00241     break;
00242   case QS_FACTOR_MAXMULT:
00243     f->maxmult = EGlpNumToLf (val);
00244     break;
00245   case QS_FACTOR_UPDMAXMULT:
00246     f->updmaxmult = EGlpNumToLf (val);
00247     break;
00248   case QS_FACTOR_DENSE_FRACT:
00249     f->dense_fract = EGlpNumToLf (val);
00250     break;
00251   case QS_FACTOR_PARTIAL_TOL:
00252     EGlpNumCopy (f->partial_tol, val);
00253     EGlpNumCopy (f->partial_cur, val);
00254     break;
00255   default:
00256     fprintf (stderr, "Invalid param %d in ILLfactor_set_factor_dparam\n",
00257              param);
00258     return 1;
00259   }
00260   return 0;
00261 }
00262 
00263 int ILLfactor_create_factor_work (
00264   factor_work * f,
00265   int dim)
00266 {
00267   int i;
00268   int rval;
00269 
00270   f->dim = dim;
00271   f->etacnt = 0;
00272   f->work_coef = EGlpNumAllocArray (dim);
00273   ILL_SAFE_MALLOC (f->work_indx, dim, int);
00274 
00275   ILL_SAFE_MALLOC (f->uc_inf, dim + (f->max_k + 1), uc_info);
00276   ILL_SAFE_MALLOC (f->ur_inf, dim + (f->max_k + 1), ur_info);
00277   ILL_SAFE_MALLOC (f->lc_inf, dim, lc_info);
00278   ILL_SAFE_MALLOC (f->lr_inf, dim, lr_info);
00279   ILL_SAFE_MALLOC (f->rperm, dim, int);
00280   ILL_SAFE_MALLOC (f->rrank, dim, int);
00281   ILL_SAFE_MALLOC (f->cperm, dim, int);
00282   ILL_SAFE_MALLOC (f->crank, dim, int);
00283 
00284   for (i = dim + f->max_k + 1; i--;)
00285     EGlpNumInitVar (f->ur_inf[i].max);
00286 
00287   for (i = 0; i < dim; i++)
00288   {
00289     EGlpNumZero (f->work_coef[i]);
00290     f->work_indx[i] = 0;
00291     f->uc_inf[i].nzcnt = 0;
00292     f->ur_inf[i].nzcnt = 0;
00293     f->lc_inf[i].nzcnt = 0;
00294     f->lr_inf[i].nzcnt = 0;
00295     f->rperm[i] = i;
00296     f->rrank[i] = i;
00297     f->cperm[i] = i;
00298     f->crank[i] = i;
00299   }
00300   for (i = 0; i <= f->max_k; i++)
00301   {
00302     f->uc_inf[dim + i].nzcnt = i;
00303     f->uc_inf[dim + i].next = dim + i;
00304     f->uc_inf[dim + i].prev = dim + i;
00305     f->ur_inf[dim + i].nzcnt = i;
00306     f->ur_inf[dim + i].next = dim + i;
00307     f->ur_inf[dim + i].prev = dim + i;
00308   }
00309 
00310   rval = ILLsvector_alloc (&f->xtmp, dim);
00311   CHECKRVALG (rval, CLEANUP);
00312 
00313   rval = 0;
00314 
00315 CLEANUP:
00316   if (rval)
00317   {
00318     ILLfactor_free_factor_work (f);
00319   }
00320   EG_RETURN (rval);
00321 }
00322 #ifdef FACTOR_DEBUG
00323 static void dump_matrix (
00324   factor_work * f,
00325   int remaining)
00326 {
00327   int dim = f->dim;
00328   ur_info *ur_inf = f->ur_inf;
00329   uc_info *uc_inf = f->uc_inf;
00330   lc_info *lc_inf = f->lc_inf;
00331   lr_info *lr_inf = f->lr_inf;
00332   er_info *er_inf = f->er_inf;
00333   int nzcnt;
00334   int beg;
00335 
00336   int i;
00337   int j;
00338 
00339   for (i = 0; i < dim; i++)
00340   {
00341     if (!remaining || ur_inf[i].next >= 0)
00342     {
00343       printf ("Row %d %d (max %.3f):", i, f->rrank[i],
00344               EGlpNumToLf (ur_inf[i].max));
00345       nzcnt = ur_inf[i].nzcnt;
00346       beg = ur_inf[i].rbeg;
00347       for (j = 0; j < nzcnt; j++)
00348       {
00349         if (j == ur_inf[i].pivcnt)
00350         {
00351           printf (" |");
00352         }
00353         printf (" %.3f*%d", EGlpNumToLf (f->urcoef[beg + j]),
00354                 f->urindx[beg + j]);
00355         if (f->urcind)
00356           printf ("@%d", f->urcind[beg + j]);
00357       }
00358       printf ("\n");
00359     }
00360   }
00361   if (f->dmat)
00362   {
00363     int start = 0;
00364 
00365     if (remaining)
00366       start = f->stage - f->dense_base;
00367     printf ("Dcols at %d %d - %d    :", f->stage - f->dense_base,
00368             f->dense_base + start, f->nstages);
00369     for (j = start; j < f->dcols; j++)
00370     {
00371       printf (" %5d", f->cperm[j + f->dense_base]);
00372     }
00373     printf ("\n");
00374     for (i = start; i < f->drows; i++)
00375     {
00376       printf ("DRow %d %d (max %.3f):", i,
00377               f->rperm[i + f->dense_base],
00378               EGlpNumToLf (ur_inf[f->rperm[i + f->dense_base]].max));
00379       for (j = start; j < f->dcols; j++)
00380       {
00381         if (j == f->drows)
00382         {
00383           printf (" |");
00384         }
00385         printf (" %.3f", EGlpNumToLf (f->dmat[i * f->dcols + j]));
00386       }
00387       printf ("\n");
00388     }
00389   }
00390 
00391   if (!remaining)
00392   {
00393     for (i = 0; i < f->stage; i++)
00394     {
00395       printf ("L col %d:", lc_inf[i].c);
00396       nzcnt = lc_inf[i].nzcnt;
00397       beg = lc_inf[i].cbeg;
00398       for (j = 0; j < nzcnt; j++)
00399       {
00400         printf (" %.3f*%d", EGlpNumToLf (f->lccoef[beg + j]),
00401                 f->lcindx[beg + j]);
00402       }
00403       printf ("\n");
00404     }
00405     for (i = f->nstages; i < f->dim; i++)
00406     {
00407       printf ("L col %d:", lc_inf[i].c);
00408       nzcnt = lc_inf[i].nzcnt;
00409       beg = lc_inf[i].cbeg;
00410       for (j = 0; j < nzcnt; j++)
00411       {
00412         printf (" %.3f*%d", EGlpNumToLf (f->lccoef[beg + j]),
00413                 f->lcindx[beg + j]);
00414       }
00415       printf ("\n");
00416     }
00417     for (i = 0; i < f->dim; i++)
00418     {
00419       if (!lr_inf[i].nzcnt)
00420         continue;
00421       printf ("L row %d:", lr_inf[i].r);
00422       nzcnt = lr_inf[i].nzcnt;
00423       beg = lr_inf[i].rbeg;
00424       for (j = 0; j < nzcnt; j++)
00425       {
00426         printf (" %.3f*%d", EGlpNumToLf (f->lrcoef[beg + j]),
00427                 f->lrindx[beg + j]);
00428       }
00429       printf ("\n");
00430     }
00431   }
00432 
00433   if (!remaining)
00434   {
00435     for (i = 0; i < f->etacnt; i++)
00436     {
00437       printf ("Eta row %d:", f->er_inf[i].r);
00438       nzcnt = er_inf[i].nzcnt;
00439       beg = er_inf[i].rbeg;
00440       for (j = 0; j < nzcnt; j++)
00441       {
00442         printf (" %.3f*%d", EGlpNumToLf (f->ercoef[beg + j]),
00443                 f->erindx[beg + j]);
00444       }
00445       printf ("\n");
00446     }
00447   }
00448 
00449   for (i = 0; i < dim; i++)
00450   {
00451     if (!remaining || uc_inf[i].next >= 0)
00452     {
00453       printf ("Col %d %d:", i, f->crank[i]);
00454       nzcnt = uc_inf[i].nzcnt;
00455       beg = uc_inf[i].cbeg;
00456       for (j = 0; j < nzcnt; j++)
00457       {
00458         if (f->uccoef != 0)
00459         {
00460           printf (" %.3f*%d", EGlpNumToLf (f->uccoef[beg + j]),
00461                   f->ucindx[beg + j]);
00462           if (f->ucrind)
00463             printf ("@%d", f->ucrind[beg + j]);
00464         }
00465         else
00466         {
00467           printf (" %d", f->ucindx[beg + j]);
00468         }
00469       }
00470       printf ("\n");
00471     }
00472   }
00473 
00474   if (!remaining)
00475   {
00476     printf ("rperm:");
00477     for (i = 0; i < dim; i++)
00478     {
00479       if (i == f->nstages)
00480         printf ("|");
00481       if (i == f->stage)
00482         printf ("|");
00483       printf (" %d", f->rperm[i]);
00484     }
00485     printf ("\n");
00486 
00487     printf ("cperm:");
00488     for (i = 0; i < dim; i++)
00489     {
00490       if (i == f->nstages)
00491         printf ("|");
00492       if (i == f->stage)
00493         printf ("|");
00494       printf (" %d", f->cperm[i]);
00495     }
00496     printf ("\n");
00497   }
00498 
00499   printf ("Rows by nzcnt:\n");
00500   for (i = 0; i <= f->max_k; i++)
00501   {
00502     if (ur_inf[dim + i].next != dim + i)
00503     {
00504       printf ("%d:", i);
00505       for (j = ur_inf[dim + i].next; j != dim + i; j = ur_inf[j].next)
00506       {
00507         printf (" %d", j);
00508       }
00509       printf ("\n");
00510     }
00511   }
00512 
00513   printf ("Cols by nzcnt:\n");
00514   for (i = 0; i <= f->max_k; i++)
00515   {
00516     if (uc_inf[dim + i].next != dim + i)
00517     {
00518       printf ("%d:", i);
00519       for (j = uc_inf[dim + i].next; j != dim + i; j = uc_inf[j].next)
00520       {
00521         printf (" %d", j);
00522       }
00523       printf ("\n");
00524     }
00525   }
00526 
00527   printf ("\n");
00528   fflush (stdout);
00529 }
00530 #endif
00531 
00532 #ifdef SORT_RESULTS
00533 static void sort_vector2 (
00534   int nzcnt,
00535   int *indx,
00536   EGlpNum_t * coef)
00537 {
00538   int i;
00539   int j;
00540   int itmp;
00541   EGlpNum_t ctmp;
00542 
00543   EGlpNumInitVar (ctmp);
00544 
00545   for (i = 1; i < nzcnt; i++)
00546   {
00547     itmp = indx[i];
00548     EGlpNumCopy (ctmp, coef[i]);
00549     for (j = i; j >= 1 && indx[j - 1] > itmp; j--)
00550     {
00551       indx[j] = indx[j - 1];
00552       EGlpNumCopy (coef[j], coef[j - 1]);
00553     }
00554     indx[j] = itmp;
00555     EGlpNumCopy (coef[j], ctmp);
00556   }
00557   EGlpNumClearVar (ctmp);
00558 }
00559 
00560 static void sort_vector (
00561   svector * x)
00562 {
00563   sort_vector2 (x->nzcnt, x->indx, x->coef);
00564 }
00565 #endif
00566 
00567 #ifdef DEBUG_FACTOR
00568 static int check_matrix (
00569   factor_work * f)
00570 {
00571   ur_info *ur_inf = f->ur_inf;
00572   uc_info *uc_inf = f->uc_inf;
00573   int rbeg;
00574   int nzcnt;
00575   int cbeg;
00576   int c;
00577   int r;
00578   int j;
00579   int nerr = 0;
00580 
00581   for (r = 0; r < f->dim; r++)
00582   {
00583     nzcnt = ur_inf[r].nzcnt;
00584     rbeg = ur_inf[r].rbeg;
00585     for (j = 0; j < nzcnt; j++)
00586     {
00587       c = f->urindx[rbeg + j];
00588       cbeg = uc_inf[c].cbeg;
00589       if (f->ucindx[cbeg + f->urcind[rbeg + j]] != r)
00590       {
00591         MESSAGE(0,"index mismatch, row %d column %d", r, c);
00592         nerr++;
00593       }
00594       if (fabs(EGlpNumToLf(f->uccoef[cbeg + f->urcind[rbeg + j]]) - EGlpNumToLf(f->urcoef[rbeg + j]))>1000*EGlpNumToLf(epsLpNum))
00595       {
00596         MESSAGE(0,"coef mismatch, row %d column %d", r, c);
00597         nerr++;
00598       }
00599     }
00600   }
00601   if (f->urindx[f->ur_space] != 0)
00602   {
00603     MESSAGE(0,"last urindx entry %d != 0", f->urindx[f->ur_space]);
00604     nerr++;
00605   }
00606 
00607   for (c = 0; c < f->dim; c++)
00608   {
00609     nzcnt = uc_inf[c].nzcnt;
00610     cbeg = uc_inf[c].cbeg;
00611     for (j = 0; j < nzcnt; j++)
00612     {
00613       r = f->ucindx[cbeg + j];
00614       rbeg = ur_inf[r].rbeg;
00615       if (f->urindx[rbeg + f->ucrind[cbeg + j]] != c)
00616       {
00617         MESSAGE(0,"index mismatch, column %d row %d", c, r);
00618         nerr++;
00619       }
00620       if (f->urcoef[rbeg + f->ucrind[cbeg + j]] != f->uccoef[cbeg + j])
00621       {
00622         MESSAGE(0,"coef mismatch, column %d row %d", c, r);
00623         nerr++;
00624       }
00625     }
00626   }
00627   if (f->ucindx[f->uc_space] != 0)
00628   {
00629     MESSAGE(0,"last ucindx entry %d != 0", f->ucindx[f->uc_space]);
00630     nerr++;
00631   }
00632   if (nerr)
00633   {
00634     dump_matrix (f, 0);
00635     return E_CHECK_FAILED;
00636   }
00637   return 0;
00638 }
00639 #endif
00640 
00641 #ifdef FACTOR_STATS
00642 static void dump_factor_stats (
00643   factor_work * f)
00644 {
00645   int dim = f->dim;
00646   int ecnt = f->etacnt;
00647   ur_info *ur_inf = f->ur_inf;
00648   lc_info *lc_inf = f->lc_inf;
00649   er_info *er_inf = f->er_inf;
00650   EGlpNum_t *urcoef = f->urcoef;
00651   EGlpNum_t *lccoef = f->lccoef;
00652   EGlpNum_t *ercoef = f->ercoef;
00653   int lnzcnt = 0;
00654   int unzcnt = 0;
00655   int enzcnt = 0;
00656   int nzcnt;
00657   int beg;
00658   EGlpNum_t umax;
00659   EGlpNum_t lmax;
00660   EGlpNum_t emax;
00661   int i;
00662   int j;
00663 
00664   EGlpNumInitVar (umax);
00665   EGlpNumInitVar (lmax);
00666   EGlpNumInitVar (emax);
00667   EGlpNumZero (umax);
00668   for (i = 0; i < dim; i++)
00669   {
00670     nzcnt = ur_inf[i].nzcnt;
00671     beg = ur_inf[i].rbeg;
00672     unzcnt += nzcnt;
00673     for (j = 0; j < nzcnt; j++)
00674     {
00675       EGlpNumSetToMaxAbs (umax, urcoef[beg + j]);
00676     }
00677   }
00678   EGlpNumZero (lmax);
00679   for (i = 0; i < dim; i++)
00680   {
00681     nzcnt = lc_inf[i].nzcnt;
00682     beg = lc_inf[i].cbeg;
00683     lnzcnt += nzcnt;
00684     for (j = 0; j < nzcnt; j++)
00685     {
00686       EGlpNumSetToMaxAbs (lmax, lccoef[beg + j]);
00687     }
00688   }
00689   EGlpNumZero (emax);
00690   for (i = 0; i < ecnt; i++)
00691   {
00692     nzcnt = er_inf[i].nzcnt;
00693     beg = er_inf[i].rbeg;
00694     enzcnt += nzcnt;
00695     for (j = 0; j < nzcnt; j++)
00696     {
00697       EGlpNumSetToMaxAbs (emax, ercoef[beg + j]);
00698     }
00699   }
00700   MESSAGE(0, "factor U %d nzs %.3e max L %d nzs %.3e max E %d nzs %.3e max",
00701           unzcnt, EGlpNumToLf (umax), lnzcnt, EGlpNumToLf (lmax), enzcnt,
00702           EGlpNumToLf (emax));
00703   fflush (stdout);
00704   EGlpNumClearVar (umax);
00705   EGlpNumClearVar (lmax);
00706   EGlpNumClearVar (emax);
00707 }
00708 #endif
00709 
00710 static void clear_work (
00711   factor_work * f)
00712 {
00713   int i;
00714   int dim = f->dim;
00715   EGlpNum_t *work_coef = f->work_coef;
00716 
00717   for (i = 0; i < dim; i++)
00718   {
00719     EGlpNumZero (work_coef[i]);
00720   }
00721 }
00722 
00723 static void load_row (
00724   factor_work * f,
00725   int r)
00726 {
00727   EGlpNum_t *prow_urcoef = f->urcoef + f->ur_inf[r].rbeg;
00728   int *prow_urindx = f->urindx + f->ur_inf[r].rbeg;
00729   int prow_nzcnt = f->ur_inf[r].nzcnt;
00730   EGlpNum_t *work_coef = f->work_coef;
00731   int *work_indx = f->work_indx;
00732   int i;
00733   int j;
00734 
00735   for (i = 0; i < prow_nzcnt; i++)
00736   {
00737     j = prow_urindx[i];
00738     EGlpNumCopy (work_coef[j], prow_urcoef[i]);
00739     work_indx[j] = 1;
00740   }
00741 }
00742 
00743 static void clear_row (
00744   factor_work * f,
00745   int r)
00746 {
00747   int *prow_urindx = f->urindx + f->ur_inf[r].rbeg;
00748   int prow_nzcnt = f->ur_inf[r].nzcnt;
00749   EGlpNum_t *work_coef = f->work_coef;
00750   int *work_indx = f->work_indx;
00751   int i;
00752   int j;
00753 
00754   for (i = 0; i < prow_nzcnt; i++)
00755   {
00756     j = prow_urindx[i];
00757     EGlpNumZero (work_coef[j]);
00758     work_indx[j] = 0;
00759   }
00760 }
00761 
00762 static int make_ur_space (
00763   factor_work * f,
00764   int space)
00765 {
00766   EGlpNum_t *new_urcoef = 0;
00767   int *new_urindx = 0;
00768   int *new_urcind = 0;
00769   EGlpNum_t *urcoef = f->urcoef;
00770   int *urindx = f->urindx;
00771   int *urcind = f->urcind;
00772   int minspace;
00773   ur_info *ur_inf = f->ur_inf;
00774   int dim = f->dim;
00775   int new_nzcnt = 0, old_nzcnt;
00776   int rbeg;
00777   int nzcnt;
00778   int i;
00779   int j;
00780   int rval;
00781 
00782   minspace = f->ur_space;
00783   nzcnt = space;
00784   for (i = 0; i < dim; i++)
00785     nzcnt += ur_inf[i].nzcnt;
00786   old_nzcnt = nzcnt;
00787   while (nzcnt * 2 >= minspace)
00788   {
00789     minspace = 1 + minspace * f->grow_mul;
00790   }
00791 
00792 #ifdef GROWTH_STATS
00793   printf ("make_ur_space growing from %d to %d...", f->ur_space, minspace);
00794   fflush (stdout);
00795 #endif
00796   new_urcoef = EGlpNumAllocArray (minspace);
00797   ILL_SAFE_MALLOC (new_urindx, minspace + 1, int);
00798 
00799   if (urcind)
00800   {
00801     ILL_SAFE_MALLOC (new_urcind, minspace, int);
00802   }
00803 
00804   if (urcind)
00805   {
00806     for (j = 0; j < dim; j++)
00807     {
00808       rbeg = ur_inf[j].rbeg;
00809       nzcnt = ur_inf[j].nzcnt;
00810       ur_inf[j].rbeg = new_nzcnt;
00811       for (i = 0; i < nzcnt; i++)
00812       {
00813         new_urindx[new_nzcnt] = urindx[rbeg + i];
00814         EGlpNumCopy (new_urcoef[new_nzcnt], urcoef[rbeg + i]);
00815         new_urcind[new_nzcnt] = urcind[rbeg + i];
00816         new_nzcnt++;
00817       }
00818     }
00819   }
00820   else
00821   {
00822     for (j = 0; j < dim; j++)
00823     {
00824       rbeg = ur_inf[j].rbeg;
00825       nzcnt = ur_inf[j].nzcnt;
00826       ur_inf[j].rbeg = new_nzcnt;
00827       for (i = 0; i < nzcnt; i++)
00828       {
00829         new_urindx[new_nzcnt] = urindx[rbeg + i];
00830         EGlpNumCopy (new_urcoef[new_nzcnt], urcoef[rbeg + i]);
00831         new_nzcnt++;
00832       }
00833     }
00834   }
00835 
00836   for (i = new_nzcnt; i < minspace; i++)
00837   {
00838     new_urindx[i] = -1;
00839   }
00840   new_urindx[minspace] = 0;
00841   EGlpNumFreeArray (f->urcoef);
00842   f->urcoef = new_urcoef;
00843   new_urcoef = 0;
00844 
00845   ILL_IFFREE (f->urindx, int);
00846 
00847   f->urindx = new_urindx;
00848   new_urindx = 0;
00849 
00850   ILL_IFFREE (f->urcind, int);
00851 
00852   f->urcind = new_urcind;
00853   new_urcind = 0;
00854 
00855   f->ur_freebeg = new_nzcnt;
00856   f->ur_space = minspace;
00857 
00858 #ifdef GROWTH_STATS
00859   MESSAGE (0,"%d/%d nonzeros", new_nzcnt, old_nzcnt);
00860   dump_factor_stats (f);
00861 #endif
00862 
00863   rval = 0;
00864 
00865 CLEANUP:
00866   ILL_IFFREE (new_urcoef, EGlpNum_t);
00867   ILL_IFFREE (new_urindx, int);
00868   ILL_IFFREE (new_urcind, int);
00869 
00870   EG_RETURN (rval);
00871 }
00872 
00873 static int make_uc_space (
00874   factor_work * f,
00875   int space)
00876 {
00877   EGlpNum_t *new_uccoef = 0;
00878   int *new_ucindx = 0;
00879   int *new_ucrind = 0;
00880   int uc_freebeg = f->uc_freebeg;
00881   EGlpNum_t *uccoef = f->uccoef;
00882   int *ucindx = f->ucindx;
00883   int *ucrind = f->ucrind;
00884   int minspace = uc_freebeg + space;
00885   uc_info *uc_inf = f->uc_inf;
00886   int dim = f->dim;
00887   int new_nzcnt = 0;
00888   int cbeg;
00889   int nzcnt;
00890   int i;
00891   int j;
00892   int rval;
00893 
00894   minspace = f->uc_space;
00895   nzcnt = space;
00896   for( i = 0 ; i < dim ; i++) nzcnt += uc_inf[i].nzcnt;
00897   while(nzcnt*2 >= minspace)
00898   {
00899     minspace = 10 + (f->grow_mul * minspace);
00900   }
00901 
00902 #ifdef GROWTH_STATS
00903   MESSAGE (0,"make_uc_space growing from %d to %d...", f->uc_space, minspace);
00904 #endif
00905 
00906   ILL_SAFE_MALLOC (new_ucindx, minspace + 1, int);
00907 
00908   if (ucrind)
00909   {
00910     new_uccoef = EGlpNumAllocArray (minspace);
00911     ILL_SAFE_MALLOC (new_ucrind, minspace, int);
00912   }
00913 
00914   if (ucrind)
00915   {
00916     for (j = 0; j < dim; j++)
00917     {
00918       cbeg = uc_inf[j].cbeg;
00919       nzcnt = uc_inf[j].nzcnt;
00920       uc_inf[j].cbeg = new_nzcnt;
00921       for (i = 0; i < nzcnt; i++)
00922       {
00923         new_ucindx[new_nzcnt] = ucindx[cbeg + i];
00924         EGlpNumCopy (new_uccoef[new_nzcnt], uccoef[cbeg + i]);
00925         new_ucrind[new_nzcnt] = ucrind[cbeg + i];
00926         new_nzcnt++;
00927       }
00928     }
00929   }
00930   else
00931   {
00932     for (j = 0; j < dim; j++)
00933     {
00934       cbeg = uc_inf[j].cbeg;
00935       nzcnt = uc_inf[j].nzcnt;
00936       uc_inf[j].cbeg = new_nzcnt;
00937       for (i = 0; i < nzcnt; i++)
00938       {
00939         new_ucindx[new_nzcnt] = ucindx[cbeg + i];
00940         new_nzcnt++;
00941       }
00942     }
00943   }
00944 
00945   for (i = new_nzcnt; i < minspace; i++)
00946   {
00947     new_ucindx[i] = -1;
00948   }
00949   new_ucindx[minspace] = 0;
00950 
00951   EGlpNumFreeArray (f->uccoef);
00952   f->uccoef = new_uccoef;
00953   new_uccoef = 0;
00954 
00955   ILL_IFFREE (f->ucindx, int);
00956 
00957   f->ucindx = new_ucindx;
00958   new_ucindx = 0;
00959 
00960   ILL_IFFREE (f->ucrind, int);
00961 
00962   f->ucrind = new_ucrind;
00963   new_ucrind = 0;
00964 
00965   f->uc_freebeg = new_nzcnt;
00966   f->uc_space = minspace;
00967 
00968 #ifdef GROWTH_STATS
00969   MESSAGE (0,"%d nonzeros", new_nzcnt);
00970   dump_factor_stats (f);
00971 #endif
00972 
00973   rval = 0;
00974 
00975 CLEANUP:
00976   ILL_IFFREE (new_uccoef, EGlpNum_t);
00977   ILL_IFFREE (new_ucindx, int);
00978   ILL_IFFREE (new_ucrind, int);
00979 
00980   EG_RETURN (rval);
00981 }
00982 
00983 static int make_lc_space (
00984   factor_work * f,
00985   int space)
00986 {
00987   EGlpNum_t *new_lccoef = 0;
00988   int *new_lcindx = 0;
00989   int lc_freebeg = f->lc_freebeg;
00990   EGlpNum_t *lccoef = f->lccoef;
00991   int *lcindx = f->lcindx;
00992   int minspace = lc_freebeg + space;
00993   int i;
00994   int rval;
00995 
00996   if (f->lc_space * f->grow_mul > minspace)
00997   {
00998     minspace = f->lc_space * f->grow_mul;
00999   }
01000 
01001 #ifdef GROWTH_STATS
01002   MESSAGE (0,"make_lc_space growing from %d to %d...", f->lc_space, minspace);
01003 #endif
01004 
01005   new_lccoef = EGlpNumAllocArray (minspace);
01006   ILL_SAFE_MALLOC (new_lcindx, minspace, int);
01007 
01008   for (i = 0; i < lc_freebeg; i++)
01009   {
01010     EGlpNumCopy (new_lccoef[i], lccoef[i]);
01011     new_lcindx[i] = lcindx[i];
01012   }
01013 
01014   EGlpNumFreeArray (lccoef);
01015   f->lccoef = new_lccoef;
01016   new_lccoef = 0;
01017 
01018   ILL_IFFREE (lcindx, int);
01019 
01020   f->lcindx = new_lcindx;
01021   new_lcindx = 0;
01022 
01023   f->lc_space = minspace;
01024 
01025 #ifdef GROWTH_STATS
01026   dump_factor_stats (f);
01027 #endif
01028 
01029   rval = 0;
01030 
01031 CLEANUP:
01032   ILL_IFFREE (new_lccoef, EGlpNum_t);
01033   ILL_IFFREE (new_lcindx, int);
01034 
01035   EG_RETURN (rval);
01036 }
01037 
01038 static void set_col_nz (
01039   factor_work * f,
01040   int c)
01041 {
01042   uc_info *uc_inf = f->uc_inf;
01043   int nzcnt = uc_inf[c].nzcnt;
01044   int max_k = f->max_k;
01045   int dim = f->dim;
01046 
01047   if (uc_inf[c].next >= 0)
01048   {
01049     uc_inf[uc_inf[c].next].prev = uc_inf[c].prev;
01050     uc_inf[uc_inf[c].prev].next = uc_inf[c].next;
01051 
01052     if (nzcnt >= max_k)
01053       nzcnt = max_k;
01054     uc_inf[c].next = uc_inf[dim + nzcnt].next;
01055     uc_inf[c].prev = dim + nzcnt;
01056     uc_inf[dim + nzcnt].next = c;
01057     uc_inf[uc_inf[c].next].prev = c;
01058   }
01059 }
01060 
01061 static void set_row_nz (
01062   factor_work * f,
01063   int r)
01064 {
01065   ur_info *ur_inf = f->ur_inf;
01066   int nzcnt = ur_inf[r].pivcnt;
01067   int max_k = f->max_k;
01068   int dim = f->dim;
01069 
01070   if (ur_inf[r].next >= 0)
01071   {
01072     ur_inf[ur_inf[r].next].prev = ur_inf[r].prev;
01073     ur_inf[ur_inf[r].prev].next = ur_inf[r].next;
01074 
01075     if (nzcnt >= max_k)
01076       nzcnt = max_k;
01077     ur_inf[r].next = ur_inf[dim + nzcnt].next;
01078     ur_inf[r].prev = dim + nzcnt;
01079     ur_inf[dim + nzcnt].next = r;
01080     ur_inf[ur_inf[r].next].prev = r;
01081   }
01082 }
01083 
01084 static void remove_col_nz (
01085   factor_work * f,
01086   int r,
01087   int c)
01088 {
01089   uc_info *uc_inf = f->uc_inf;
01090   int *ucindx = f->ucindx + uc_inf[c].cbeg;
01091   int nzcnt = uc_inf[c].nzcnt;
01092   int i;
01093 
01094   for (i = 0; i < nzcnt; i++)
01095   {
01096     if (ucindx[i] == r)
01097     {
01098       --nzcnt;
01099       ucindx[i] = ucindx[nzcnt];
01100       ucindx[nzcnt] = -1;
01101       break;
01102     }
01103   }
01104   uc_inf[c].nzcnt = nzcnt;
01105 
01106   set_col_nz (f, c);
01107 }
01108 
01109 static void remove_row_nz (
01110   factor_work * f,
01111   int r,
01112   int c)
01113 {
01114   ur_info *ur_inf = f->ur_inf;
01115   int *urindx = f->urindx + ur_inf[r].rbeg;
01116   EGlpNum_t *urcoef = f->urcoef + ur_inf[r].rbeg;
01117   int pivcnt = ur_inf[r].pivcnt;
01118   EGlpNum_t max;
01119   int tind;
01120   EGlpNum_t tcoef;
01121   int i;
01122 
01123   EGlpNumInitVar (tcoef);
01124   EGlpNumInitVar (max);
01125   EGlpNumZero (max);
01126 
01127   for (i = 0; i < pivcnt; i++)
01128   {
01129     if (urindx[i] == c)
01130     {
01131       --pivcnt;
01132       ILL_SWAP (urindx[i], urindx[pivcnt], tind);
01133       EGLPNUM_SWAP (urcoef[i], urcoef[pivcnt], tcoef);
01134       --i;
01135     }
01136     else
01137     {
01138       EGlpNumSetToMaxAbs (max, urcoef[i]);
01139     }
01140   }
01141   ur_inf[r].pivcnt = pivcnt;
01142   EGlpNumCopy (ur_inf[r].max, max);
01143   set_row_nz (f, r);
01144   EGlpNumClearVar (max);
01145   EGlpNumClearVar (tcoef);
01146 }
01147 
01148 static int add_col_nz (
01149   factor_work * f,
01150   int r,
01151   int c)
01152 {
01153   uc_info *uc_inf = f->uc_inf;
01154   int cbeg = uc_inf[c].cbeg;
01155   int nzcnt = uc_inf[c].nzcnt;
01156   int uc_freebeg = f->uc_freebeg;
01157   int *ucindx = f->ucindx;
01158   int i;
01159   int rval = 0;
01160 
01161   if (uc_inf[c].next == -1)
01162   {
01163     return 0;
01164   }
01165 
01166   if (ucindx[cbeg + nzcnt] == -1)
01167   {
01168     ucindx[cbeg + nzcnt] = r;
01169     uc_inf[c].nzcnt++;
01170     if (nzcnt + cbeg == uc_freebeg)
01171     {
01172       f->uc_freebeg = uc_freebeg + 1;
01173     }
01174   }
01175   else
01176   {
01177     if (uc_freebeg + nzcnt + 1 >= f->uc_space)
01178     {
01179       rval = make_uc_space (f, nzcnt + 1);
01180       CHECKRVALG (rval, CLEANUP);
01181       uc_freebeg = f->uc_freebeg;
01182       cbeg = uc_inf[c].cbeg;
01183       ucindx = f->ucindx;
01184     }
01185     for (i = 0; i < nzcnt; i++)
01186     {
01187       ucindx[uc_freebeg + i] = ucindx[cbeg + i];
01188       ucindx[cbeg + i] = -1;
01189     }
01190     ucindx[uc_freebeg + nzcnt] = r;
01191     uc_inf[c].cbeg = uc_freebeg;
01192     uc_inf[c].nzcnt++;
01193     f->uc_freebeg = uc_freebeg + nzcnt + 1;
01194   }
01195 
01196   set_col_nz (f, c);
01197 CLEANUP:
01198   EG_RETURN (rval);
01199 }
01200 
01201 static void disable_col (
01202   factor_work * f,
01203   int c)
01204 {
01205   uc_info *uc_inf = f->uc_inf;
01206 
01207   if (uc_inf[c].next >= 0)
01208   {
01209     uc_inf[uc_inf[c].next].prev = uc_inf[c].prev;
01210     uc_inf[uc_inf[c].prev].next = uc_inf[c].next;
01211 
01212     uc_inf[c].next = -2;
01213     uc_inf[c].prev = -2;
01214   }
01215 }
01216 
01217 static void remove_col (
01218   factor_work * f,
01219   int c)
01220 {
01221   uc_info *uc_inf = f->uc_inf;
01222   int cbeg = uc_inf[c].cbeg;
01223   int nzcnt = uc_inf[c].nzcnt;
01224   int *ucindx = f->ucindx;
01225   int i;
01226 
01227   for (i = 0; i < nzcnt; i++)
01228   {
01229     ucindx[cbeg + i] = -1;
01230   }
01231   uc_inf[c].cbeg = 0;
01232   uc_inf[c].nzcnt = 0;
01233 
01234   if (uc_inf[c].next >= 0)
01235   {
01236     uc_inf[uc_inf[c].next].prev = uc_inf[c].prev;
01237     uc_inf[uc_inf[c].prev].next = uc_inf[c].next;
01238 
01239     uc_inf[c].next = -1;
01240     uc_inf[c].prev = -1;
01241   }
01242 }
01243 
01244 static void remove_row (
01245   factor_work * f,
01246   int r)
01247 {
01248   ur_info *ur_inf = f->ur_inf;
01249 
01250   if (ur_inf[r].next >= 0)
01251   {
01252     ur_inf[ur_inf[r].next].prev = ur_inf[r].prev;
01253     ur_inf[ur_inf[r].prev].next = ur_inf[r].next;
01254 
01255     ur_inf[r].next = -1;
01256     ur_inf[r].prev = -1;
01257   }
01258 }
01259 
01260 static void find_coef (
01261   factor_work * f,
01262   int r,
01263   int c,
01264   EGlpNum_t * coef)
01265 {
01266   EGlpNum_t *prow_urcoef = f->urcoef + f->ur_inf[r].rbeg;
01267   int *prow_urindx = f->urindx + f->ur_inf[r].rbeg;
01268   int i;
01269   int prow_nzcnt = f->ur_inf[r].nzcnt;
01270 
01271   EGlpNumZero (*coef);
01272   for (i = 0; i < prow_nzcnt; i++)
01273   {
01274     if (prow_urindx[i] == c)
01275     {
01276       EGlpNumCopy (*coef, prow_urcoef[i]);
01277       return;
01278     }
01279   }
01280   fprintf (stderr, "Coefficient not found\n");
01281   return;
01282 }
01283 
01284 static int elim_row (
01285   factor_work * f,
01286   int elim_r,
01287   int r,
01288   int c,
01289   EGlpNum_t * p_pivot_coef)
01290 {
01291   ur_info *ur_inf = f->ur_inf;
01292   EGlpNum_t *work_coef = f->work_coef;
01293   int *work_indx = f->work_indx;
01294   EGlpNum_t *urcoef = f->urcoef;
01295   int *urindx = f->urindx;
01296   int prow_beg = ur_inf[r].rbeg;
01297   int prow_nzcnt = ur_inf[r].nzcnt;
01298   int prow_pivcnt = ur_inf[r].pivcnt;
01299   int fill = ur_inf[elim_r].nzcnt;
01300   int cancel = 0;
01301   EGlpNum_t max;
01302   int erow_beg;
01303   int erow_nzcnt;
01304   int erow_pivcnt;
01305   EGlpNum_t x;
01306   int i;
01307   int j;
01308   int rval = 0;
01309   EGlpNum_t elim_coef;
01310 
01311   EGlpNumInitVar (max);
01312   EGlpNumInitVar (x);
01313   EGlpNumInitVar (elim_coef);
01314   EGlpNumZero (max);
01315   find_coef (f, r, c, &elim_coef);
01316   EGlpNumDivTo (elim_coef, work_coef[c]);
01317   EGlpNumCopy (*p_pivot_coef, elim_coef);
01318 
01319   for (i = 0; i < prow_nzcnt; i++)
01320   {
01321     j = urindx[prow_beg + i];
01322     if (work_indx[j] == 1)
01323     {
01324       EGlpNumCopy (x, urcoef[prow_beg + i]);
01325       EGlpNumSubInnProdTo (x, elim_coef, work_coef[j]);
01326       if ((!(EGlpNumIsNeqZero (x, f->fzero_tol))) || j == c)
01327       {
01328         cancel++;
01329         if (j != c)
01330         {
01331           remove_col_nz (f, r, j);
01332         }
01333         if (i < prow_pivcnt)
01334         {
01335           prow_pivcnt--;
01336           prow_nzcnt--;
01337           urindx[prow_beg + i] = urindx[prow_beg + prow_pivcnt];
01338           EGlpNumCopy (urcoef[prow_beg + i], urcoef[prow_beg + prow_pivcnt]);
01339           if (prow_pivcnt != prow_nzcnt)
01340           {
01341             urindx[prow_beg + prow_pivcnt] = urindx[prow_beg + prow_nzcnt];
01342             EGlpNumCopy (urcoef[prow_beg + prow_pivcnt],
01343                          urcoef[prow_beg + prow_nzcnt]);
01344           }
01345         }
01346         else
01347         {
01348           prow_nzcnt--;
01349           urindx[prow_beg + i] = urindx[prow_beg + prow_nzcnt];
01350           EGlpNumCopy (urcoef[prow_beg + i], urcoef[prow_beg + prow_nzcnt]);
01351         }
01352         urindx[prow_beg + prow_nzcnt] = -1;
01353         i--;
01354       }
01355       else
01356       {
01357         EGlpNumCopy (urcoef[prow_beg + i], x);
01358         if (i < prow_pivcnt)
01359         {
01360           EGlpNumSetToMaxAbs (max, x);
01361         }
01362       }
01363       work_indx[j] = 0;
01364       fill--;
01365     }
01366     else
01367     {
01368       if (i < prow_pivcnt)
01369       {
01370         EGlpNumSetToMaxAbs (max, urcoef[prow_beg + i]);
01371       }
01372     }
01373   }
01374 
01375   if (fill > 0)
01376   {
01377     ur_inf[r].nzcnt = prow_nzcnt;
01378     ur_inf[r].pivcnt = prow_pivcnt;
01379     if (fill > cancel)
01380     {
01381       int ur_freebeg = f->ur_freebeg;
01382 
01383       if (ur_freebeg + prow_nzcnt + fill >= f->ur_space)
01384       {
01385         rval = make_ur_space (f, prow_nzcnt + fill);
01386         CHECKRVALG (rval, CLEANUP);
01387         urcoef = f->urcoef;
01388         urindx = f->urindx;
01389         ur_freebeg = f->ur_freebeg;
01390         prow_beg = f->ur_inf[r].rbeg;
01391       }
01392       for (i = 0; i < prow_nzcnt; i++)
01393       {
01394         urindx[ur_freebeg + i] = urindx[prow_beg + i];
01395         EGlpNumCopy (urcoef[ur_freebeg + i], urcoef[prow_beg + i]);
01396         urindx[prow_beg + i] = -1;
01397       }
01398       ur_inf[r].rbeg = ur_freebeg;
01399       f->ur_freebeg = ur_freebeg + prow_nzcnt + fill;
01400       prow_beg = ur_freebeg;
01401     }
01402 
01403     erow_beg = ur_inf[elim_r].rbeg;
01404     erow_nzcnt = ur_inf[elim_r].nzcnt;
01405     erow_pivcnt = ur_inf[elim_r].pivcnt;
01406 
01407     for (i = 0; i < erow_pivcnt; i++)
01408     {
01409       j = urindx[erow_beg + i];
01410       if (work_indx[j] == 1)
01411       {
01412         EGlpNumCopyNeg (x, elim_coef);
01413         EGlpNumMultTo (x, urcoef[erow_beg + i]);
01414         if (EGlpNumIsNeqZero (x, f->fzero_tol))
01415         {
01416           rval = add_col_nz (f, r, j);
01417           CHECKRVALG (rval, CLEANUP);
01418           if (prow_pivcnt != prow_nzcnt)
01419           {
01420             urindx[prow_beg + prow_nzcnt] = urindx[prow_beg + prow_pivcnt];
01421             EGlpNumCopy (urcoef[prow_beg + prow_nzcnt],
01422                          urcoef[prow_beg + prow_pivcnt]);
01423           }
01424           urindx[prow_beg + prow_pivcnt] = j;
01425           EGlpNumCopy (urcoef[prow_beg + prow_pivcnt], x);
01426           EGlpNumSetToMaxAbs (max, x);
01427           prow_pivcnt++;
01428           prow_nzcnt++;
01429         }
01430       }
01431       else
01432       {
01433         work_indx[j] = 1;
01434       }
01435     }
01436     for (i = erow_pivcnt; i < erow_nzcnt; i++)
01437     {
01438       j = urindx[erow_beg + i];
01439       if (work_indx[j] == 1)
01440       {
01441         EGlpNumCopyNeg (x, elim_coef);
01442         EGlpNumMultTo (x, urcoef[erow_beg + i]);
01443         if (EGlpNumIsNeqZero (x, f->fzero_tol))
01444         {
01445           rval = add_col_nz (f, r, j);
01446           CHECKRVALG (rval, CLEANUP);
01447           urindx[prow_beg + prow_nzcnt] = j;
01448           EGlpNumCopy (urcoef[prow_beg + prow_nzcnt], x);
01449           prow_nzcnt++;
01450         }
01451       }
01452       else
01453       {
01454         work_indx[j] = 1;
01455       }
01456     }
01457   }
01458   else
01459   {
01460     erow_nzcnt = ur_inf[elim_r].nzcnt;
01461     erow_beg = ur_inf[elim_r].rbeg;
01462     for (i = 0; i < erow_nzcnt; i++)
01463     {
01464       j = urindx[erow_beg + i];
01465       work_indx[j] = 1;
01466     }
01467   }
01468 
01469   ur_inf[r].nzcnt = prow_nzcnt;
01470   ur_inf[r].pivcnt = prow_pivcnt;
01471   EGlpNumCopy (ur_inf[r].max, max);
01472 
01473   set_row_nz (f, r);
01474 CLEANUP:
01475   EGlpNumClearVar (elim_coef);
01476   EGlpNumClearVar (x);
01477   EGlpNumClearVar (max);
01478   EG_RETURN (rval);
01479 }
01480 
01481 #define SETPERM(f,s,r,c) {                    \
01482         f->rperm[f->rrank[r]] = f->rperm[s];  \
01483         f->rrank[f->rperm[s]] = f->rrank[r];  \
01484         f->rperm[s] = r;                      \
01485         f->rrank[r] = s;                      \
01486                                               \
01487         f->cperm[f->crank[c]] = f->cperm[s];  \
01488         f->crank[f->cperm[s]] = f->crank[c];  \
01489         f->cperm[s] = c;                      \
01490         f->crank[c] = s;                      \
01491 }
01492 
01493 static int elim (
01494   factor_work * f,
01495   int r,
01496   int c)
01497 {
01498   uc_info *uc_inf = f->uc_inf;
01499   ur_info *ur_inf = f->ur_inf;
01500   lc_info *lc_inf = f->lc_inf;
01501   int *urindx;
01502   int *ucindx;
01503   int *lcindx;
01504   EGlpNum_t *urcoef;
01505   EGlpNum_t *lccoef;
01506   EGlpNum_t pivot_coef;
01507   int nzcnt;
01508   int lc_freebeg;
01509   int s = f->stage;
01510   int i;
01511   int j;
01512   int rval = 0;
01513 
01514   EGlpNumInitVar (pivot_coef);
01515 
01516   if (uc_inf[c].nzcnt == 1)
01517   {
01518     /* col singleton */
01519     SETPERM (f, s, r, c);
01520 
01521     lc_inf[s].cbeg = -1;
01522     lc_inf[s].c = r;
01523     lc_inf[s].nzcnt = 0;
01524     f->stage++;
01525 
01526     urindx = f->urindx + ur_inf[r].rbeg;
01527     urcoef = f->urcoef + ur_inf[r].rbeg;
01528     nzcnt = ur_inf[r].nzcnt;
01529     for (i = 0; i < nzcnt; i++)
01530     {
01531       j = urindx[i];
01532       remove_col_nz (f, r, j);
01533       if (j == c)
01534       {
01535         urindx[i] = urindx[0];
01536         urindx[0] = c;
01537         EGLPNUM_SWAP (urcoef[0], urcoef[i], pivot_coef);
01538       }
01539     }
01540     remove_row (f, r);
01541     remove_col (f, c);
01542   }
01543   else if (ur_inf[r].nzcnt == 1)
01544   {
01545     /* row singleton */
01546     --(f->nstages);
01547     SETPERM (f, f->nstages, r, c);
01548 
01549     lc_inf[f->nstages].cbeg = -1;
01550     lc_inf[f->nstages].c = r;
01551     lc_inf[f->nstages].nzcnt = 0;
01552 
01553     ucindx = f->ucindx + uc_inf[c].cbeg;
01554     nzcnt = uc_inf[c].nzcnt;
01555     for (i = 0; i < nzcnt; i++)
01556     {
01557       j = ucindx[i];
01558       remove_row_nz (f, j, c);
01559     }
01560     remove_row (f, r);
01561     remove_col (f, c);
01562   }
01563   else
01564   {
01565     SETPERM (f, s, r, c);
01566     f->stage++;
01567 
01568     nzcnt = uc_inf[c].nzcnt;
01569     if (f->lc_freebeg + nzcnt >= f->lc_space)
01570     {
01571       rval = make_lc_space (f, nzcnt);
01572       CHECKRVALG (rval, CLEANUP);
01573     }
01574     lc_freebeg = f->lc_freebeg;
01575     lc_inf[s].cbeg = lc_freebeg;
01576     lc_inf[s].c = r;
01577     lcindx = f->lcindx;
01578     lccoef = f->lccoef;
01579     load_row (f, r);
01580     ucindx = f->ucindx + uc_inf[c].cbeg;
01581     for (i = 0; i < nzcnt; i++)
01582     {
01583       j = f->ucindx[uc_inf[c].cbeg + i];
01584       if (j != r)
01585       {
01586         rval = elim_row (f, r, j, c, &pivot_coef);
01587         CHECKRVALG (rval, CLEANUP);
01588         lcindx[lc_freebeg] = j;
01589         EGlpNumCopy (lccoef[lc_freebeg], pivot_coef);
01590         lc_freebeg++;
01591 #ifdef TRACK_FACTOR
01592         EGlpNumSetToMaxAbs (f->maxelem_factor, pivot_coef);
01593         if (EGlpNumIsLess (f->maxelem_factor, ur_inf[r].max))
01594           EGlpNumCopy (f->maxelem_factor, ur_inf[r].max);
01595 #endif /* TRACK_FACTOR */
01596       }
01597     }
01598     lc_inf[s].nzcnt = lc_freebeg - lc_inf[s].cbeg;
01599     f->lc_freebeg = lc_freebeg;
01600 
01601     clear_row (f, r);
01602 
01603     urindx = f->urindx + ur_inf[r].rbeg;
01604     urcoef = f->urcoef + ur_inf[r].rbeg;
01605     nzcnt = ur_inf[r].nzcnt;
01606     for (i = 0; i < nzcnt; i++)
01607     {
01608       j = urindx[i];
01609       remove_col_nz (f, r, j);
01610       if (j == c)
01611       {
01612         urindx[i] = urindx[0];
01613         urindx[0] = c;
01614         EGLPNUM_SWAP (urcoef[0], urcoef[i], pivot_coef);
01615       }
01616     }
01617     remove_row (f, r);
01618     remove_col (f, c);
01619   }
01620 CLEANUP:
01621   EGlpNumClearVar (pivot_coef);
01622   EG_RETURN (rval);
01623 }
01624 
01625 static void find_pivot_column (
01626   factor_work * f,
01627   int c,
01628   int *p_r)
01629 {
01630   uc_info *uc_inf = f->uc_inf;
01631   ur_info *ur_inf = f->ur_inf;
01632   int *ucindx = f->ucindx;
01633   int nzcnt = uc_inf[c].nzcnt;
01634   int cbeg = uc_inf[c].cbeg;
01635   EGlpNum_t num_tmp[2];
01636   int bestnz = -1;
01637   int i;
01638   int r;
01639 
01640   EGlpNumInitVar (num_tmp[0]);
01641   EGlpNumInitVar (num_tmp[1]);
01642 
01643   *p_r = -1;
01644   for (i = 0; i < nzcnt; i++)
01645   {
01646     r = ucindx[cbeg + i];
01647     if((bestnz == -1 || ur_inf[r].pivcnt < bestnz))
01648     {
01649       find_coef (f, r, c, num_tmp);
01650       if(EGlpNumIsLessZero(num_tmp[0]))
01651         EGlpNumSign (num_tmp[0]);
01652       EGlpNumCopy (num_tmp[1], f->partial_cur);
01653       EGlpNumMultTo (num_tmp[1], ur_inf[r].max);
01654       if(EGlpNumIsLeq (num_tmp[1], num_tmp[0]))
01655       {
01656         bestnz = ur_inf[r].pivcnt;
01657         *p_r = r;
01658       }
01659     }
01660   }
01661   EGlpNumClearVar (num_tmp[0]);
01662   EGlpNumClearVar (num_tmp[1]);
01663 }
01664 
01665 static void find_pivot_row (
01666   factor_work * f,
01667   int r,
01668   int *p_c)
01669 {
01670   uc_info *uc_inf = f->uc_inf;
01671   ur_info *ur_inf = f->ur_inf;
01672   int *urindx = f->urindx;
01673   EGlpNum_t *urcoef = f->urcoef;
01674   int pivcnt = ur_inf[r].pivcnt;
01675   int rbeg = ur_inf[r].rbeg;
01676   EGlpNum_t thresh[2];
01677   int bestnz = -1;
01678   int i;
01679   int c;
01680 
01681   EGlpNumInitVar (thresh[0]);
01682   EGlpNumInitVar (thresh[1]);
01683   EGlpNumCopy (thresh[0], f->partial_cur);
01684   EGlpNumMultTo (thresh[0], ur_inf[r].max);
01685   *p_c = -1;
01686   for (i = 0; i < pivcnt; i++)
01687   {
01688     c = urindx[rbeg + i];
01689     if ((bestnz == -1 || uc_inf[c].nzcnt < bestnz))
01690     {
01691       EGlpNumCopyAbs (thresh[1], urcoef[rbeg + i]);
01692       if(EGlpNumIsLeq (thresh[0], thresh[1]))
01693       {
01694         bestnz = uc_inf[c].nzcnt;
01695         *p_c = c;
01696       }
01697     }
01698   }
01699   EGlpNumClearVar (thresh[0]);
01700   EGlpNumClearVar (thresh[1]);
01701 }
01702 
01703 static int find_pivot (
01704   factor_work * f,
01705   int *p_r,
01706   int *p_c)
01707 {
01708   uc_info *uc_inf = f->uc_inf;
01709   ur_info *ur_inf = f->ur_inf;
01710   int dim = f->dim;
01711   int max_k = f->max_k;
01712   int p = f->p;
01713   int c;
01714   int r;
01715   int mm = 0;
01716   int n = 0;
01717   int m;
01718   int k = 2;
01719 
01720   if (uc_inf[dim + 1].next != dim + 1)
01721   {
01722     c = uc_inf[dim + 1].next;
01723     r = f->ucindx[uc_inf[c].cbeg];
01724     *p_c = c;
01725     *p_r = r;
01726     return 0;
01727   }
01728   else if (ur_inf[dim + 1].next != dim + 1)
01729   {
01730     r = ur_inf[dim + 1].next;
01731     c = f->urindx[ur_inf[r].rbeg];
01732     *p_c = c;
01733     *p_r = r;
01734     return 0;
01735   }
01736   *p_r = -1;
01737   *p_c = -1;
01738   for (; k <= max_k && (mm == 0 || mm > (k - 1) * (k - 1)); k++)
01739   {
01740     if (uc_inf[dim + k].next != dim + k)
01741     {
01742       for (c = uc_inf[dim + k].next; c != dim + k; c = uc_inf[c].next)
01743       {
01744         find_pivot_column (f, c, &r);
01745         if (r >= 0)
01746         {
01747           m = (uc_inf[c].nzcnt - 1) * (ur_inf[r].pivcnt - 1);
01748           if (mm == 0 || m < mm)
01749           {
01750             mm = m;
01751             *p_c = c;
01752             *p_r = r;
01753             if (mm <= (k - 1) * (k - 1))
01754             {
01755               return 0;
01756             }
01757           }
01758         }
01759         else
01760         {
01761           c = uc_inf[c].prev;
01762           disable_col (f, uc_inf[c].next);
01763         }
01764         n++;
01765         if (n >= p && mm != 0)
01766         {
01767           return 0;
01768         }
01769       }
01770     }
01771 
01772     if (ur_inf[dim + k].next != dim + k)
01773     {
01774       for (r = ur_inf[dim + k].next; r != dim + k; r = ur_inf[r].next)
01775       {
01776         find_pivot_row (f, r, &c);
01777         if (c >= 0)
01778         {
01779           m = (uc_inf[c].nzcnt - 1) * (ur_inf[r].pivcnt - 1);
01780           if (mm == 0 || m < mm)
01781           {
01782             mm = m;
01783             *p_c = c;
01784             *p_r = r;
01785             if (mm <= k * (k - 1))
01786             {
01787               return 0;
01788             }
01789           }
01790         }
01791         n++;
01792         if (n >= p && mm != 0)
01793         {
01794           return 0;
01795         }
01796       }
01797     }
01798   }
01799   if (mm != 0)
01800   {
01801     return 0;
01802   }
01803   else
01804   {
01805     //fprintf (stderr, "No acceptable pivot found\n");
01806     return E_NO_PIVOT;
01807   }
01808 }
01809 
01810 static int create_factor_space (
01811   factor_work * f)
01812 {
01813   uc_info *uc_inf = f->uc_inf;
01814   ur_info *ur_inf = f->ur_inf;
01815   int dim = f->dim;
01816   int nzcnt;
01817   int i;
01818   int rval;
01819 
01820   nzcnt = 0;
01821   for (i = 0; i < dim; i++)
01822   {
01823     nzcnt += ur_inf[i].nzcnt;
01824   }
01825 
01826   if (f->ucindx == 0)
01827   {
01828     f->uc_space = nzcnt * f->uc_space_mul;
01829     ILL_SAFE_MALLOC (f->ucindx, f->uc_space + 1, int);
01830   }
01831 
01832   if (f->urindx == 0 || f->urcoef == 0)
01833   {
01834     ILL_IFFREE (f->urindx, int);
01835 
01836     EGlpNumFreeArray (f->urcoef);
01837     f->ur_space = nzcnt * f->ur_space_mul;
01838     ILL_SAFE_MALLOC (f->urindx, f->ur_space + 1, int);
01839 
01840     f->urcoef = EGlpNumAllocArray (f->ur_space);
01841   }
01842 
01843   if (f->lcindx == 0 || f->lccoef == 0)
01844   {
01845     ILL_IFFREE (f->lcindx, int);
01846 
01847     EGlpNumFreeArray (f->lccoef);
01848     f->lc_space = nzcnt * f->lc_space_mul;
01849     ILL_SAFE_MALLOC (f->lcindx, f->lc_space, int);
01850 
01851     f->lccoef = EGlpNumAllocArray (f->lc_space);
01852   }
01853 
01854   nzcnt = 0;
01855   for (i = 0; i < dim; i++)
01856   {
01857     ur_inf[i].rbeg = nzcnt;
01858     nzcnt += ur_inf[i].nzcnt;
01859     ur_inf[i].nzcnt = ur_inf[i].rbeg;
01860   }
01861   f->ur_freebeg = nzcnt;
01862 
01863   nzcnt = 0;
01864   for (i = 0; i < dim; i++)
01865   {
01866     uc_inf[i].cbeg = nzcnt;
01867     nzcnt += uc_inf[i].nzcnt;
01868     uc_inf[i].nzcnt = uc_inf[i].cbeg;
01869   }
01870   f->uc_freebeg = nzcnt;
01871 
01872   f->lc_freebeg = 0;
01873 
01874   rval = 0;
01875 CLEANUP:
01876   EG_RETURN (rval);
01877 }
01878 
01879 static int init_matrix (
01880   factor_work * f,
01881   int *basis,
01882   int *cbeg,
01883   int *clen,
01884   int *in_ucindx,
01885   EGlpNum_t * in_uccoef)
01886 {
01887   uc_info *uc_inf = f->uc_inf;
01888   ur_info *ur_inf = f->ur_inf;
01889   int dim = f->dim;
01890   int max_k = f->max_k;
01891   int *ucindx;
01892   int *urindx;
01893   EGlpNum_t *urcoef;
01894   int nzcnt;
01895   int beg;
01896   int i;
01897   int j;
01898   int r;
01899   int rval = 0;
01900   EGlpNum_t v;
01901   EGlpNum_t max;
01902 
01903   EGlpNumInitVar (v);
01904   EGlpNumInitVar (max);
01905 
01906   for (i = 0; i < dim; i++)
01907   {
01908     ur_inf[i].nzcnt = 0;
01909   }
01910   for (i = 0; i < dim; i++)
01911   {
01912     nzcnt = clen[basis[i]];
01913     beg = cbeg[basis[i]];
01914     uc_inf[i].nzcnt = nzcnt;
01915     for (j = 0; j < nzcnt; j++)
01916     {
01917       r = in_ucindx[beg + j];
01918       ur_inf[r].nzcnt++;
01919     }
01920   }
01921 
01922   rval = create_factor_space (f);
01923   CHECKRVALG (rval, CLEANUP);
01924 
01925   urindx = f->urindx;
01926   ucindx = f->ucindx;
01927   urcoef = f->urcoef;
01928 
01929   for (i = 0; i < dim; i++)
01930   {
01931     nzcnt = clen[basis[i]];
01932     beg = cbeg[basis[i]];
01933     for (j = 0; j < nzcnt; j++)
01934     {
01935       EGlpNumCopy (v, in_uccoef[beg + j]);
01936       if (!(EGlpNumIsNeqZero (v, f->fzero_tol)))
01937         continue;
01938       r = in_ucindx[beg + j];
01939       ucindx[uc_inf[i].nzcnt++] = r;
01940       urindx[ur_inf[r].nzcnt] = i;
01941       EGlpNumCopy (urcoef[ur_inf[r].nzcnt], v);
01942       ur_inf[r].nzcnt++;
01943     }
01944   }
01945 
01946   for (i = 0; i < dim; i++)
01947   {
01948     uc_inf[i].nzcnt -= uc_inf[i].cbeg;
01949     ur_inf[i].nzcnt -= ur_inf[i].rbeg;
01950   }
01951 
01952   j = f->uc_space;
01953   for (i = f->uc_freebeg; i < j; i++)
01954   {
01955     ucindx[i] = -1;
01956   }
01957   ucindx[j] = 0;
01958 
01959   j = f->ur_space;
01960   for (i = f->ur_freebeg; i < j; i++)
01961   {
01962     urindx[i] = -1;
01963   }
01964   urindx[j] = 0;
01965 
01966   for (i = 0; i < dim; i++)
01967   {
01968     nzcnt = ur_inf[i].nzcnt;
01969     ur_inf[i].pivcnt = nzcnt;
01970     beg = ur_inf[i].rbeg;
01971     EGlpNumZero (max);
01972     for (j = 0; j < nzcnt; j++)
01973     {
01974       EGlpNumSetToMaxAbs (max, urcoef[beg + j]);
01975     }
01976     EGlpNumCopy (ur_inf[i].max, max);
01977   }
01978 
01979   for (i = 0; i <= max_k; i++)
01980   {
01981     ur_inf[dim + i].next = dim + i;
01982     ur_inf[dim + i].prev = dim + i;
01983     uc_inf[dim + i].next = dim + i;
01984     uc_inf[dim + i].prev = dim + i;
01985   }
01986 
01987   for (i = 0; i < dim; i++)
01988   {
01989     nzcnt = uc_inf[i].nzcnt;
01990     if (nzcnt >= max_k)
01991       nzcnt = max_k;
01992     uc_inf[i].next = uc_inf[dim + nzcnt].next;
01993     uc_inf[i].prev = dim + nzcnt;
01994     uc_inf[dim + nzcnt].next = i;
01995     uc_inf[uc_inf[i].next].prev = i;
01996 
01997     nzcnt = ur_inf[i].pivcnt;
01998     if (nzcnt >= max_k)
01999       nzcnt = max_k;
02000     ur_inf[i].next = ur_inf[dim + nzcnt].next;
02001     ur_inf[i].prev = dim + nzcnt;
02002     ur_inf[dim + nzcnt].next = i;
02003     ur_inf[ur_inf[i].next].prev = i;
02004   }
02005 
02006 #ifdef TRACK_FACTOR
02007   EGlpNumZero (max);
02008   nzcnt = 0;
02009   for (i = 0; i < dim; i++)
02010   {
02011     if (EGlpNumIsLess (max, ur_inf[i].max))
02012       EGlpNumCopy (max, ur_inf[i].max);
02013     nzcnt += ur_inf[i].nzcnt;
02014   }
02015 
02016   EGlpNumCopy (f->maxelem_orig, max);
02017   f->nzcnt_orig = nzcnt;
02018   EGlpNumCopy (f->maxelem_factor, f->maxelem_orig);
02019   f->nzcnt_factor = f->nzcnt_orig;
02020 #endif /* TRACK_FACTOR */
02021 
02022   /* sentinal for column space */
02023   ucindx[f->uc_space] = 0;
02024 
02025   clear_work (f);
02026 
02027 CLEANUP:
02028   EGlpNumClearVar (max);
02029   EGlpNumClearVar (v);
02030   EG_RETURN (rval);
02031 }
02032 
02033 static int build_iteration_u_data (
02034   factor_work * f)
02035 {
02036   int dim = f->dim;
02037   ur_info *ur_inf = f->ur_inf;
02038   uc_info *uc_inf = f->uc_inf;
02039   EGlpNum_t *uccoef = 0;
02040   int *ucindx = 0;
02041   int *urindx = f->urindx;
02042   EGlpNum_t *urcoef = f->urcoef;
02043   int *ucrind = 0;
02044   int *urcind = 0;
02045   int nzcnt;
02046   int beg;
02047   int cbeg;
02048   int cnzcnt;
02049   int uc_space = f->uc_space;
02050   int er_space;
02051   int i;
02052   int j;
02053   int k;
02054   int rval;
02055 
02056   nzcnt = 0;
02057   for (i = 0; i < dim; i++)
02058   {
02059     nzcnt += ur_inf[i].nzcnt;
02060   }
02061 
02062 #ifdef TRACK_FACTOR
02063   f->nzcnt_factor = nzcnt;
02064 #endif /* TRACK_FACTOR */
02065 
02066   EGlpNumFreeArray (f->uccoef);
02067   uccoef = EGlpNumAllocArray (nzcnt);
02068   f->uccoef = uccoef;
02069 
02070   ILL_IFFREE (f->ucrind, int);
02071   ILL_SAFE_MALLOC (ucrind, nzcnt, int);
02072 
02073   f->ucrind = ucrind;
02074 
02075   ILL_IFFREE (f->urcind, int);
02076   ILL_SAFE_MALLOC (urcind, f->ur_space, int);
02077 
02078   f->urcind = urcind;
02079 
02080   if (uc_space < nzcnt)
02081   {
02082     ILL_IFFREE (f->ucindx, int);
02083     ILL_SAFE_MALLOC (f->ucindx, nzcnt + 1, int);
02084   }
02085   f->uc_space = nzcnt;
02086   uc_space = nzcnt;
02087   ucindx = f->ucindx;
02088 
02089   for (i = 0; i < dim; i++)
02090   {
02091     uc_inf[i].nzcnt = 0;
02092   }
02093 
02094   for (i = 0; i < dim; i++)
02095   {
02096     nzcnt = ur_inf[i].nzcnt;
02097     beg = ur_inf[i].rbeg;
02098     for (j = 0; j < nzcnt; j++)
02099     {
02100       uc_inf[urindx[beg + j]].nzcnt++;
02101     }
02102     ur_inf[i].delay = 0;
02103   }
02104 
02105   nzcnt = 0;
02106   for (i = 0; i < dim; i++)
02107   {
02108     uc_inf[i].cbeg = nzcnt;
02109     nzcnt += uc_inf[i].nzcnt;
02110     uc_inf[i].nzcnt = 0;
02111     uc_inf[i].delay = 0;
02112   }
02113 
02114   f->uc_freebeg = nzcnt;
02115   for (i = nzcnt; i < uc_space; i++)
02116   {
02117     ucindx[i] = -1;
02118   }
02119   ucindx[uc_space] = 0;
02120 
02121   for (i = 0; i < dim; i++)
02122   {
02123     nzcnt = ur_inf[i].nzcnt;
02124     beg = ur_inf[i].rbeg;
02125     k = urindx[beg];
02126     cbeg = uc_inf[k].cbeg;
02127     cnzcnt = uc_inf[k].nzcnt;
02128     if (cnzcnt != 0)
02129     {
02130       ucindx[cbeg + cnzcnt] = ucindx[cbeg];
02131       EGlpNumCopy (uccoef[cbeg + cnzcnt], uccoef[cbeg]);
02132       ucrind[cbeg + cnzcnt] = ucrind[cbeg];
02133       urcind[ur_inf[ucindx[cbeg]].rbeg + ucrind[cbeg]] = cnzcnt;
02134     }
02135     ucindx[cbeg] = i;
02136     EGlpNumCopy (uccoef[cbeg], urcoef[beg]);
02137     ucrind[cbeg] = 0;
02138     urcind[beg] = 0;
02139     uc_inf[k].nzcnt = cnzcnt + 1;
02140     for (j = 1; j < nzcnt; j++)
02141     {
02142       k = urindx[beg + j];
02143       cbeg = uc_inf[k].cbeg;
02144       cnzcnt = uc_inf[k].nzcnt;
02145       ucindx[cbeg + cnzcnt] = i;
02146       EGlpNumCopy (uccoef[cbeg + cnzcnt], urcoef[beg + j]);
02147       ucrind[cbeg + cnzcnt] = j;
02148       urcind[beg + j] = cnzcnt;
02149       uc_inf[k].nzcnt++;
02150     }
02151   }
02152 
02153   for (i = 0; i < dim; i++)
02154   {
02155     f->rrank[f->rperm[i]] = i;
02156   }
02157 
02158   nzcnt = f->ur_space;
02159 
02160   for (i = f->ur_freebeg; i < nzcnt; i++)
02161   {
02162     urindx[i] = -1;
02163   }
02164   urindx[nzcnt] = 0;
02165 
02166   clear_work (f);
02167 
02168   er_space = f->er_space_mul * f->etamax;
02169   ILL_SAFE_MALLOC (f->er_inf, f->etamax, er_info);
02170   ILL_SAFE_MALLOC (f->erindx, er_space, int);
02171 
02172   f->ercoef = EGlpNumAllocArray (er_space);
02173   f->etacnt = 0;
02174   f->er_freebeg = 0;
02175   f->er_space = er_space;
02176 
02177   rval = 0;
02178 
02179 CLEANUP:
02180   EG_RETURN (rval);
02181 }
02182 
02183 static int build_iteration_l_data (
02184   factor_work * f)
02185 {
02186   int dim = f->dim;
02187   lc_info *lc_inf = f->lc_inf;
02188   lr_info *lr_inf = f->lr_inf;
02189   EGlpNum_t *lrcoef = 0;
02190   int *lrindx = 0;
02191   EGlpNum_t *lccoef = f->lccoef;
02192   int *lcindx = f->lcindx;
02193   int nzcnt;
02194   int beg;
02195   int rnzcnt;
02196   int rbeg;
02197   int i;
02198   int j;
02199   int k;
02200   int c;
02201   int rval;
02202 
02203   nzcnt = 0;
02204   for (i = 0; i < dim; i++)
02205   {
02206     nzcnt += lc_inf[i].nzcnt;
02207     lr_inf[i].nzcnt = 0;
02208     lr_inf[i].delay = 0;
02209     lc_inf[lc_inf[i].c].crank = i;
02210   }
02211 
02212   EGlpNumFreeArray (f->lrcoef);
02213   if (nzcnt)
02214   {
02215     lrcoef = EGlpNumAllocArray (nzcnt);
02216     f->lrcoef = lrcoef;
02217   }
02218 
02219   ILL_IFFREE (f->lrindx, int);
02220   ILL_SAFE_MALLOC (lrindx, nzcnt + 1, int);
02221 
02222   f->lrindx = lrindx;
02223 
02224   for (i = 0; i < dim; i++)
02225   {
02226     nzcnt = lc_inf[i].nzcnt;
02227     beg = lc_inf[i].cbeg;
02228     lc_inf[i].delay = 0;
02229     for (j = 0; j < nzcnt; j++)
02230     {
02231       lr_inf[lc_inf[lcindx[beg + j]].crank].nzcnt++;
02232     }
02233   }
02234 
02235   nzcnt = 0;
02236   for (i = 0; i < dim; i++)
02237   {
02238     lr_inf[i].rbeg = nzcnt;
02239     nzcnt += lr_inf[i].nzcnt;
02240     lr_inf[i].nzcnt = 0;
02241     lr_inf[i].r = lc_inf[i].c;
02242     lr_inf[lr_inf[i].r].rrank = i;
02243   }
02244 
02245   for (i = 0; i < dim; i++)
02246   {
02247     nzcnt = lc_inf[i].nzcnt;
02248     beg = lc_inf[i].cbeg;
02249     c = lc_inf[i].c;
02250     for (j = 0; j < nzcnt; j++)
02251     {
02252       k = lc_inf[lcindx[beg + j]].crank;
02253       rbeg = lr_inf[k].rbeg;
02254       rnzcnt = lr_inf[k].nzcnt;
02255       lrindx[rbeg + rnzcnt] = c;
02256       EGlpNumCopy (lrcoef[rbeg + rnzcnt], lccoef[beg + j]);
02257       lr_inf[k].nzcnt++;
02258     }
02259   }
02260 
02261 #ifdef TRACK_FACTOR
02262   nzcnt = f->nzcnt_factor;
02263   for (i = 0; i < dim; i++)
02264   {
02265     nzcnt += lc_inf[i].nzcnt;
02266   }
02267   f->nzcnt_factor = nzcnt;
02268 
02269   EGlpNumCopy (f->maxelem_cur, f->maxelem_factor);
02270   f->nzcnt_cur = f->nzcnt_factor;
02271 
02272 /*
02273     dump_factor_stats (f);
02274     printf ("orig max  %e nzcnt %d\n", f->maxelem_orig, f->nzcnt_orig);
02275     printf ("f maxelem %e nzcnt %d\n", f->maxelem_cur, f->nzcnt_cur);
02276 */
02277 #endif /* TRACK_FACTOR */
02278 
02279   rval = 0;
02280 
02281 CLEANUP:
02282   EG_RETURN (rval);
02283 }
02284 
02285 static int handle_singularity (
02286   factor_work * f)
02287 {
02288   int rval = 0;
02289   int nsing;
02290   int *singr = 0;
02291   int *singc = 0;
02292   int i;
02293 
02294   if (f->p_nsing == 0 || f->p_singr == 0 || f->p_singc == 0)
02295   {
02296     fprintf (stderr, "singular basis, but no place for singularity data\n");
02297     return E_SING_NO_DATA;
02298   }
02299 
02300   nsing = f->nstages - f->stage;
02301   ILL_SAFE_MALLOC (singr, nsing, int);
02302   ILL_SAFE_MALLOC (singc, nsing, int);
02303 
02304   for (i = f->stage; i < f->nstages; i++)
02305   {
02306     singr[i - f->stage] = f->rperm[i];
02307     singc[i - f->stage] = f->cperm[i];
02308   }
02309   *f->p_nsing = nsing;
02310   *f->p_singr = singr;
02311   *f->p_singc = singc;
02312   singr = 0;
02313   singc = 0;
02314 
02315 CLEANUP:
02316   ILL_IFFREE (singr, int);
02317   ILL_IFFREE (singc, int);
02318 
02319   EG_RETURN (rval);
02320 }
02321 
02322 static int dense_build_matrix (
02323   factor_work * f)
02324 {
02325   EGlpNum_t *dmat = 0;
02326   int stage = f->stage;
02327   int drows = f->nstages - stage;
02328   int dcols = f->dim - stage;
02329   int dsize = drows * dcols;
02330   int *crank = f->crank;
02331   EGlpNum_t *urcoef = f->urcoef;
02332   int *urindx = f->urindx;
02333   int nzcnt;
02334   int beg;
02335   int i;
02336   int r;
02337   int j;
02338   int rval = 0;
02339 
02340   dmat = EGlpNumAllocArray (dsize);
02341 
02342   for (i = 0; i < dsize; i++)
02343     EGlpNumZero (dmat[i]);
02344 
02345   for (i = 0; i < drows; i++)
02346   {
02347     r = f->rperm[i + stage];
02348     nzcnt = f->ur_inf[r].nzcnt;
02349     beg = f->ur_inf[r].rbeg;
02350     for (j = 0; j < nzcnt; j++)
02351     {
02352       EGlpNumCopy (dmat[i * dcols - stage + crank[urindx[beg + j]]],
02353                    urcoef[beg + j]);
02354     }
02355   }
02356 
02357   f->drows = drows;
02358   f->dcols = dcols;
02359   f->dense_base = f->stage;
02360   f->dmat = dmat;
02361   dmat = 0;
02362 
02363 //CLEANUP:
02364   EGlpNumFreeArray (dmat);
02365   EG_RETURN (rval);
02366 }
02367 
02368 static int dense_find_pivot (
02369   factor_work * f,
02370   int *p_r,
02371   int *p_c)
02372 {
02373   int dcols = f->dcols;
02374   int drows = f->drows;
02375   EGlpNum_t *dmat = f->dmat;
02376   int dense_base = f->dense_base;
02377   int s = f->stage - dense_base;
02378   ur_info *ur_inf = f->ur_inf;
02379   int *rperm = f->rperm;
02380   EGlpNum_t maxval;
02381   int max_r;
02382   int max_c;
02383   int i;
02384 
02385   EGlpNumInitVar (maxval);
02386   EGlpNumZero (maxval);
02387   max_r = -1;
02388   for (i = s; i < drows; i++)
02389   {
02390     if (EGlpNumIsLess (maxval, ur_inf[rperm[dense_base + i]].max))
02391     {
02392       EGlpNumCopy (maxval, ur_inf[rperm[dense_base + i]].max);
02393       max_r = i;
02394     }
02395   }
02396   if (max_r == -1)
02397   {
02398     return E_NO_PIVOT;
02399   }
02400 
02401   EGlpNumZero (maxval);
02402   max_c = -1;
02403   for (i = s; i < drows; i++)
02404   {
02405     EGlpNumSetToMaxAbsAndDo (maxval, dmat[max_r * dcols + i], max_c = i);
02406   }
02407   if (max_c == -1)
02408   {
02409     return E_NO_PIVOT;
02410   }
02411   *p_r = max_r;
02412   *p_c = max_c;
02413 
02414   EGlpNumClearVar (maxval);
02415   return 0;
02416 }
02417 
02418 static void dense_swap (
02419   factor_work * f,
02420   int r,
02421   int c)
02422 {
02423   int dcols = f->dcols;
02424   int drows = f->drows;
02425   EGlpNum_t *dmat = f->dmat;
02426   int dense_base = f->dense_base;
02427   int s = f->stage - dense_base;
02428   int i;
02429   EGlpNum_t v;
02430 
02431   EGlpNumInitVar (v);
02432 
02433   if (r != s)
02434   {
02435     ILL_SWAP (f->rperm[dense_base + s], f->rperm[dense_base + r], i);
02436     f->rrank[f->rperm[dense_base + s]] = dense_base + s;
02437     f->rrank[f->rperm[dense_base + r]] = dense_base + r;
02438     for (i = 0; i < dcols; i++)
02439     {
02440       EGLPNUM_SWAP (dmat[s * dcols + i], dmat[r * dcols + i], v);
02441     }
02442   }
02443   if (c != s)
02444   {
02445     ILL_SWAP (f->cperm[dense_base + s], f->cperm[dense_base + c], i);
02446     f->crank[f->cperm[dense_base + s]] = dense_base + s;
02447     f->crank[f->cperm[dense_base + c]] = dense_base + c;
02448     for (i = 0; i < drows; i++)
02449     {
02450       EGLPNUM_SWAP (dmat[i * dcols + s], dmat[i * dcols + c], v);
02451     }
02452   }
02453   EGlpNumClearVar (v);
02454 }
02455 
02456 static void dense_elim (
02457   factor_work * f,
02458   int r,
02459   int c)
02460 {
02461   int dcols = f->dcols;
02462   int drows = f->drows;
02463   EGlpNum_t *dmat = f->dmat;
02464   int dense_base = f->dense_base;
02465   int s = f->stage - dense_base;
02466   ur_info *ur_inf = f->ur_inf;
02467   int *rperm = f->rperm;
02468   int i;
02469   int j;
02470   EGlpNum_t pivval;
02471   EGlpNum_t max;
02472   EGlpNum_t v;
02473   EGlpNum_t w;
02474 
02475 #ifdef TRACK_FACTOR
02476   EGlpNum_t maxelem_factor;
02477 
02478   EGlpNumInitVar (maxelem_factor);
02479   EGlpNumCopy (maxelem_factor, f->maxelem_factor);
02480 #endif
02481   EGlpNumInitVar (pivval);
02482   EGlpNumInitVar (max);
02483   EGlpNumInitVar (v);
02484   EGlpNumInitVar (w);
02485 
02486   dense_swap (f, r, c);
02487   f->stage++;
02488   EGlpNumCopyFrac (pivval, oneLpNum, dmat[s * dcols + s]);
02489   for (i = s + 1; i < drows; i++)
02490   {
02491     EGlpNumCopy (v, dmat[i * dcols + s]);
02492     if (EGlpNumIsNeqqZero (v))
02493     {
02494       EGlpNumMultTo (v, pivval);
02495       if (EGlpNumIsNeqZero (v, f->fzero_tol))
02496       {
02497         EGlpNumCopy (dmat[i * dcols + s], v);
02498 #ifdef TRACK_FACTOR
02499         EGlpNumSetToMaxAbs (maxelem_factor, v);
02500 #endif
02501         EGlpNumZero (max);
02502         for (j = s + 1; j < drows; j++)
02503         {
02504           EGlpNumCopy (w, dmat[i * dcols + j]);
02505           EGlpNumSubInnProdTo (w, v, dmat[s * dcols + j]);
02506           EGlpNumCopy (dmat[i * dcols + j], w);
02507           EGlpNumSetToMaxAbs (max, w);
02508         }
02509         for (j = drows; j < dcols; j++)
02510         {
02511           EGlpNumCopy (w, dmat[i * dcols + j]);
02512           EGlpNumSubInnProdTo (w, v, dmat[s * dcols + j]);
02513           EGlpNumCopy (dmat[i * dcols + j], w);
02514         }
02515         EGlpNumCopy (ur_inf[rperm[dense_base + i]].max, max);
02516 #ifdef TRACK_FACTOR
02517         if (EGlpNumIsLess (maxelem_factor, max))
02518           EGlpNumCopy (maxelem_factor, max);
02519 #endif
02520       }
02521       else
02522       {
02523         EGlpNumZero (dmat[i * dcols + s]);
02524       }
02525     }
02526   }
02527 #ifdef TRACK_FACTOR
02528   EGlpNumCopy (f->maxelem_factor, maxelem_factor);
02529   EGlpNumClearVar (maxelem_factor);
02530 #endif
02531   EGlpNumClearVar (pivval);
02532   EGlpNumClearVar (max);
02533   EGlpNumClearVar (v);
02534   EGlpNumClearVar (w);
02535 }
02536 
02537 static int dense_replace_row (
02538   factor_work * f,
02539   int i)
02540 {
02541   int dcols = f->dcols;
02542   int dense_base = f->dense_base;
02543   EGlpNum_t *dmat = f->dmat + i * dcols;
02544   EGlpNum_t *urcoef;
02545   ur_info *ur_inf = f->ur_inf;
02546   int *cperm = f->cperm;
02547   int r = f->rperm[dense_base + i];
02548   int *urindx;
02549   int nzcnt;
02550   int beg;
02551   int j;
02552   int rval = 0;
02553 
02554   nzcnt = 0;
02555   for (j = i; j < dcols; j++)
02556   {
02557     if (EGlpNumIsNeqZero (dmat[j], f->fzero_tol))
02558     {
02559       nzcnt++;
02560     }
02561   }
02562   if (nzcnt > ur_inf[r].nzcnt)
02563   {
02564     if (ur_inf[r].rbeg + ur_inf[r].nzcnt == f->ur_freebeg)
02565     {
02566       f->ur_freebeg = ur_inf[r].rbeg;
02567     }
02568     ur_inf[r].nzcnt = 0;
02569     if (f->ur_freebeg + nzcnt > f->ur_space)
02570     {
02571       rval = make_ur_space (f, nzcnt);
02572       CHECKRVALG (rval, CLEANUP);
02573     }
02574     ur_inf[r].rbeg = f->ur_freebeg;
02575     f->ur_freebeg += nzcnt;
02576   }
02577   beg = ur_inf[r].rbeg;
02578   urcoef = f->urcoef;
02579   urindx = f->urindx;
02580   for (j = i; j < dcols; j++)
02581   {
02582     if (EGlpNumIsNeqZero (dmat[j], f->fzero_tol))
02583     {
02584       EGlpNumCopy (urcoef[beg], dmat[j]);
02585       urindx[beg] = cperm[dense_base + j];
02586       beg++;
02587     }
02588   }
02589   ur_inf[r].nzcnt = beg - ur_inf[r].rbeg;
02590 CLEANUP:
02591   EG_RETURN (rval);
02592 }
02593 
02594 static int dense_create_col (
02595   factor_work * f,
02596   int i)
02597 {
02598   int dcols = f->dcols;
02599   int drows = f->drows;
02600   int dense_base = f->dense_base;
02601   EGlpNum_t *dmat = f->dmat;
02602   EGlpNum_t *lccoef;
02603   lc_info *lc_inf = f->lc_inf;
02604   int *rperm = f->rperm;
02605   int *lcindx;
02606   int nzcnt;
02607   int beg;
02608   int j;
02609   int rval = 0;
02610 
02611   nzcnt = 0;
02612   for (j = i + 1; j < drows; j++)
02613   {
02614     if (EGlpNumIsNeqZero (dmat[j * dcols + i], f->fzero_tol))
02615     {
02616       nzcnt++;
02617     }
02618   }
02619 
02620   if (f->lc_freebeg + nzcnt >= f->lc_space)
02621   {
02622     rval = make_lc_space (f, nzcnt);
02623     CHECKRVALG (rval, CLEANUP);
02624   }
02625   beg = f->lc_freebeg;
02626   lc_inf[dense_base + i].cbeg = beg;
02627   lc_inf[dense_base + i].c = rperm[dense_base + i];
02628   lcindx = f->lcindx;
02629   lccoef = f->lccoef;
02630 
02631   for (j = i + 1; j < drows; j++)
02632   {
02633     if (EGlpNumIsNeqZero (dmat[j * dcols + i], f->fzero_tol))
02634     {
02635       EGlpNumCopy (lccoef[beg], dmat[j * dcols + i]);
02636       lcindx[beg] = rperm[dense_base + j];
02637       beg++;
02638     }
02639   }
02640   lc_inf[dense_base + i].nzcnt = beg - lc_inf[dense_base + i].cbeg;
02641   f->lc_freebeg = beg;
02642 CLEANUP:
02643   EG_RETURN (rval);
02644 }
02645 
02646 static int dense_replace (
02647   factor_work * f)
02648 {
02649   int drows = f->drows;
02650   int rval = 0;
02651   int i;
02652 
02653   for (i = 0; i < drows; i++)
02654   {
02655     rval = dense_replace_row (f, i);
02656     CHECKRVALG (rval, CLEANUP);
02657     rval = dense_create_col (f, i);
02658     CHECKRVALG (rval, CLEANUP);
02659   }
02660   EGlpNumFreeArray (f->dmat);
02661   f->drows = 0;
02662   f->dcols = 0;
02663 CLEANUP:
02664   EG_RETURN (rval);
02665 }
02666 
02667 static int dense_factor (
02668   factor_work * f)
02669 {
02670   int r;
02671   int c;
02672   int rval = 0;
02673 
02674 #ifdef TRACK_FACTOR
02675 #ifdef NOTICE_BLOWUP
02676   double tmpsize;
02677 #endif
02678 #endif
02679 
02680 /*
02681     printf ("dense kernel, %d rows, %d  cols...\n", f->nstages - f->stage,
02682             f->dim - f->stage);
02683     fflush (stdout);
02684 */
02685 
02686   rval = dense_build_matrix (f);
02687   CHECKRVALG (rval, CLEANUP);
02688 
02689 #ifdef FACTOR_DEBUG
02690 #if (FACTOR_DEBUG+0>1)
02691   MESSAGE (0,"before Dense ILLfactor");
02692   dump_matrix (f, 1);
02693 #endif
02694 #endif
02695 
02696   while (f->stage < f->nstages)
02697   {
02698     r = f->stage - f->dense_base;
02699     rval = dense_find_pivot (f, &r, &c);
02700     if (rval == E_NO_PIVOT)
02701     {
02702       rval = handle_singularity (f);
02703       CHECKRVALG (rval, CLEANUP);
02704       return E_SINGULAR_INTERNAL;
02705     }
02706     else
02707     {
02708       CHECKRVALG (rval, CLEANUP);
02709     }
02710 #ifdef FACTOR_DEBUG
02711 #if (FACTOR_DEBUG+0>2)
02712     MESSAGE (0,"dense pivot elem: %d %d", r, c);
02713 #endif
02714 #endif /* FACTOR_DEBUG */
02715     dense_elim (f, r, c);
02716 
02717 #ifdef TRACK_FACTOR
02718 #ifdef NOTICE_BLOWUP
02719     tmpsize = f->maxmult * EGlpNumToLf (f->maxelem_orig);
02720     if (tmpsize < EGlpNumToLf (f->maxelem_factor) &&
02721         EGlpNumIsLess (f->partial_cur, oneLpNum))
02722     {
02723       return E_FACTOR_BLOWUP;
02724     }
02725 #endif /* NOTICE_BLOWUP */
02726 #endif /* TRACK_FACTOR */
02727 
02728 #ifdef FACTOR_DEBUG
02729 #if (FACTOR_DEBUG+0>1)
02730     MESSAGE (0,"After dense pivot stage %d (%d) of %d (%d)",
02731             f->stage - f->dense_base, f->stage,
02732             f->nstages - f->dense_base, f->nstages);
02733 #endif
02734 #if (FACTOR_DEBUG+0>2)
02735     dump_matrix (f, 1);
02736 #endif
02737 #endif /* FACTOR_DEBUG */
02738   }
02739 
02740 #ifdef FACTOR_DEBUG
02741   MESSAGE (0,"After dense ILLfactor:\n");
02742   dump_matrix (f, 0);
02743 #endif /* FACTOR_DEBUG */
02744 
02745   rval = dense_replace (f);
02746   CHECKRVALG (rval, CLEANUP);
02747 
02748 #ifdef FACTOR_DEBUG
02749   MESSAGE (0,"After replacement:\n");
02750   dump_matrix (f, 0);
02751 #endif /* FACTOR_DEBUG */
02752 
02753 CLEANUP:
02754   EG_RETURN (rval);
02755 }
02756 
02757 #ifdef RECORD
02758 EGioFile_t *fsave = 0;
02759 int fsavecnt = 0;
02760 #endif /* RECORD */
02761 
02762 static int ILLfactor_try (
02763   factor_work * f,
02764   int *basis,
02765   int *cbeg,
02766   int *clen,
02767   int *cindx,
02768   EGlpNum_t * ccoef)
02769 {
02770   int rval = 0;
02771   int r;
02772   int c;
02773 
02774 #ifdef TRACK_FACTOR
02775 #ifdef NOTICE_BLOWUP
02776   EGlpNum_t tmpsize;
02777 
02778   EGlpNumInitVar (tmpsize);
02779 #endif
02780 #endif
02781 
02782 #ifdef RECORD
02783   {
02784     int ncol = 0;
02785     int nzcnt = 0;
02786     int dim = f->dim;
02787     int i;
02788     int j;
02789     char fnambuf[40];
02790 
02791     for (i = 0; i < dim; i++)
02792     {
02793       if (basis[i] > ncol)
02794         ncol = basis[i];
02795     }
02796     ncol++;
02797     for (i = 0; i < ncol; i++)
02798     {
02799       nzcnt += clen[i];
02800     }
02801     if (fsave)
02802       EGioClose (fsave);
02803     #if HAVE_ZLIB
02804     sprintf (fnambuf, "prob.mat.%d.gz", fsavecnt);
02805     #elif HAVE_BZLIB
02806     sprintf (fnambuf, "prob.mat.%d.bz2", fsavecnt);
02807     #else
02808     sprintf (fnambuf, "prob.mat.%d", fsavecnt);
02809     #endif
02810     fsavecnt++;
02811     fsave = EGioOpen (fnambuf, "w");
02812     EGioPrintf (fsave, "%d %d %d\n", f->dim, ncol, nzcnt);
02813     for (i = 0; i < dim; i++)
02814     {
02815       EGioPrintf (fsave, "%d ", basis[i]);
02816     }
02817     EGioPrintf (fsave, "\n");
02818     for (i = 0; i < ncol; i++)
02819     {
02820       EGioPrintf (fsave, "%d", clen[i]);
02821       for (j = 0; j < clen[i]; j++)
02822       {
02823         EGioPrintf (fsave, " %d %.16lg", cindx[cbeg[i] + j],
02824                  EGlpNumToLf (ccoef[cbeg[i] + j]));
02825       }
02826       EGioPrintf (fsave, "\n");
02827     }
02828     EGioPrintf (fsave, "\n");
02829     EGioFlush (fsave);
02830   }
02831 #endif /* RECORD */
02832 
02833   rval = init_matrix (f, basis, cbeg, clen, cindx, ccoef);
02834   CHECKRVALG (rval, CLEANUP);
02835 
02836   f->stage = 0;
02837   f->nstages = f->dim;
02838 
02839 #ifdef FACTOR_DEBUG
02840   MESSAGE (0,"Initial matrix:");
02841 #if (FACTOR_DEBUG+0>1)
02842   dump_matrix (f, 0);
02843 #endif
02844 #endif /* FACTOR_DEBUG */
02845 #ifdef FACTOR_STATS
02846   printf ("Initial matrix: ");
02847   dump_factor_stats (f);
02848 #endif /* FACTOR_STATS */
02849 
02850   while (f->stage < f->nstages)
02851   {
02852     rval = find_pivot (f, &r, &c);
02853     if (rval == E_NO_PIVOT)
02854     {
02855       rval = handle_singularity (f);
02856       CHECKRVALG (rval, CLEANUP);
02857       return 0;
02858     }
02859     else
02860     {
02861       CHECKRVALG (rval, CLEANUP);
02862     }
02863     if (f->ur_inf[r].pivcnt > f->dense_fract * (f->nstages - f->stage) &&
02864         f->uc_inf[c].nzcnt > f->dense_fract * (f->nstages - f->stage) &&
02865         f->nstages - f->stage > f->dense_min)
02866     {
02867       rval = dense_factor (f);
02868       if (rval == E_SINGULAR_INTERNAL)
02869         return 0;
02870       if (rval)
02871         return rval;
02872       break;
02873     }
02874 #ifdef FACTOR_DEBUG
02875     MESSAGE (0,"pivot elem: %d %d", r, c);
02876 #endif /* FACTOR_DEBUG */
02877     rval = elim (f, r, c);
02878     CHECKRVALG (rval, CLEANUP);
02879 
02880 #ifdef TRACK_FACTOR
02881 #ifdef NOTICE_BLOWUP
02882     EGlpNumSet (tmpsize, f->maxmult);
02883     EGlpNumMultTo (tmpsize, f->maxelem_orig);
02884     if (EGlpNumIsLess (tmpsize, f->maxelem_factor) &&
02885         EGlpNumIsLess (f->partial_cur, oneLpNum))
02886     {
02887       return E_FACTOR_BLOWUP;
02888     }
02889 #endif /* NOTICE_BLOWUP */
02890 #endif /* TRACK_FACTOR */
02891 
02892 #ifdef FACTOR_DEBUG
02893 #if (FACTOR_DEBUG+0>3)
02894     MESSAGE (0,"After pivot stage %d of %d", f->stage, f->nstages);
02895     dump_matrix (f, 0);
02896 #endif
02897 #endif /* FACTOR_DEBUG */
02898   }
02899 
02900   rval = build_iteration_u_data (f);
02901   CHECKRVALG (rval, CLEANUP);
02902 
02903   rval = build_iteration_l_data (f);
02904   CHECKRVALG (rval, CLEANUP);
02905 
02906 #ifdef TRACK_FACTOR
02907 #ifdef NOTICE_BLOWUP
02908   EGlpNumSet (tmpsize, f->minmult);
02909   EGlpNumMultTo (tmpsize, f->maxelem_orig);
02910   if (EGlpNumIsLess (f->maxelem_factor, tmpsize) &&
02911       EGlpNumIsLess (f->partial_tol, f->partial_cur))
02912   {
02913     if (EGlpNumIsGreaDbl (f->partial_cur, 0.5))
02914     {
02915       EGlpNumSet (f->partial_cur, 0.5);
02916     }
02917     else if (EGlpNumIsGreaDbl (f->partial_cur, 0.25))
02918     {
02919       EGlpNumSet (f->partial_cur, 0.25);
02920     }
02921     else if (EGlpNumIsGreaDbl (f->partial_cur, 0.1))
02922     {
02923       EGlpNumSet (f->partial_cur, 0.1);
02924     }
02925     else
02926     {
02927       EGlpNumDivUiTo (f->partial_cur, 10);
02928     }
02929     if (EGlpNumIsLess (f->partial_cur, f->partial_tol))
02930     {
02931       EGlpNumCopy (f->partial_cur, f->partial_tol);
02932     }
02933 /*  Bico - comment out for dist 
02934         fprintf (stderr, "factor good, lowering partial tolerance to %.2f\n",
02935                  f->partial_cur);
02936 */
02937   }
02938 #endif /* NOTICE_BLOWUP */
02939 #endif /* TRACK_FACTOR */
02940 
02941 #ifdef FACTOR_DEBUG
02942   MESSAGE(0,"Factored matrix:");
02943 #if (FACTOR_DEBUG+0>1)
02944   dump_matrix (f, 0);
02945 #endif
02946 #endif /* FACTOR_DEBUG */
02947 
02948 #ifdef FACTOR_STATS
02949   printf ("Factored matrix: ");
02950   dump_factor_stats (f);
02951 #endif /* FACTOR_STATS */
02952 CLEANUP:
02953 #ifdef TRACK_FACTOR
02954 #ifdef NOTICE_BLOWUP
02955   EGlpNumClearVar (tmpsize);
02956 #endif
02957 #endif
02958   EG_RETURN (rval);
02959 }
02960 
02961 int ILLfactor (
02962   factor_work * f,
02963   int *basis,
02964   int *cbeg,
02965   int *clen,
02966   int *cindx,
02967   EGlpNum_t * ccoef,
02968   int *p_nsing,
02969   int **p_singr,
02970   int **p_singc)
02971 {
02972   int rval;
02973 
02974   f->p_nsing = p_nsing;
02975   f->p_singr = p_singr;
02976   f->p_singc = p_singc;
02977   *p_nsing = 0;
02978 
02979 AGAIN:
02980   rval = ILLfactor_try (f, basis, cbeg, clen, cindx, ccoef);
02981   if (rval == E_FACTOR_BLOWUP)
02982   {
02983     if (EGlpNumIsLessDbl (f->partial_cur, 0.1))
02984     {
02985       EGlpNumMultUiTo (f->partial_cur, 10);
02986     }
02987     else if (EGlpNumIsLessDbl (f->partial_cur, 0.25))
02988     {
02989       EGlpNumSet (f->partial_cur, 0.25);
02990     }
02991     else if (EGlpNumIsLessDbl (f->partial_cur, 0.5))
02992     {
02993       EGlpNumSet (f->partial_cur, 0.5);
02994     }
02995     else if (EGlpNumIsLess (f->partial_cur, oneLpNum))
02996     {
02997       EGlpNumOne (f->partial_cur);
02998     }
02999     else
03000     {
03001       EG_RETURN (rval);
03002     }
03003 /* Bico - comment out for dist
03004         fprintf (stderr, "factor blowup, changing partial tolerance to %.2f\n",
03005                  f->partial_cur);
03006 */
03007     goto AGAIN;
03008   }
03009   EG_RETURN (rval);
03010 }
03011 
03012 static void ILLfactor_ftranl (
03013   factor_work * f,
03014   EGlpNum_t * a)
03015 {
03016   int *lcindx = f->lcindx;
03017   lc_info *lc_inf = f->lc_inf;
03018   EGlpNum_t *lccoef = f->lccoef;
03019   int dim = f->dim;
03020   int beg;
03021   int nzcnt;
03022   int i;
03023   int j;
03024   EGlpNum_t v;
03025 
03026   EGlpNumInitVar (v);
03027 
03028   for (i = 0; i < dim; i++)
03029   {
03030     EGlpNumCopy (v, a[lc_inf[i].c]);
03031     if (EGlpNumIsNeqqZero (v))
03032     {
03033       nzcnt = lc_inf[i].nzcnt;
03034       beg = lc_inf[i].cbeg;
03035       for (j = 0; j < nzcnt; j++)
03036       {
03037         EGlpNumSubInnProdTo (a[lcindx[beg + j]], v, lccoef[beg + j]);
03038       }
03039     }
03040 #ifdef SOLVE_DEBUG
03041 #if (SOLVE_DEBUG+0 > 1)
03042     fpritf (stderr,"ILLfactor_ftran a after l %d:", i);
03043     for (j = 0; j < f->dim; j++)
03044     {
03045       fprintf (stderr," %.3f", EGlpNumToLf (a[j]));
03046     }
03047     MESSAGE(0," ");
03048 #endif
03049 #endif /* SOLVE_DEBUG */
03050   }
03051 #ifdef SOLVE_DEBUG
03052 #if (SOLVE_DEBUG+0 <= 1)
03053   printf ("ILLfactor_ftran a after l:");
03054   for (j = 0; j < f->dim; j++)
03055   {
03056     printf (" %.3f", EGlpNumToLf (a[j]));
03057   }
03058   printf ("\n");
03059   fflush (stdout);
03060 #endif
03061 #endif /* SOLVE_DEBUG */
03062   EGlpNumClearVar (v);
03063 }
03064 
03065 #if 0
03066 static void ftranl3_delay (
03067   factor_work * f,
03068   int c)
03069 {
03070   lc_info *lc_inf = f->lc_inf;
03071   int nzcnt;
03072   int *indx;
03073   int i;
03074 
03075   c = lc_inf[c].crank;
03076   nzcnt = lc_inf[c].nzcnt;
03077   indx = f->lcindx + lc_inf[c].cbeg;
03078   for (i = 0; i < nzcnt; i++)
03079   {
03080     c = indx[i];
03081     if (lc_inf[c].delay++ == 0)
03082     {
03083       ftranl3_delay (f, c);
03084     }
03085   }
03086 }
03087 #endif
03088 
03089 static void ftranl3_delay2 (
03090   factor_work * f,
03091   int c)
03092 {
03093   lc_info *lc_inf = f->lc_inf;
03094   int nzcnt;
03095   int *indx;
03096   int i;
03097   int last;
03098 
03099   do
03100   {
03101     c = lc_inf[c].crank;
03102     nzcnt = lc_inf[c].nzcnt;
03103     indx = f->lcindx + lc_inf[c].cbeg;
03104     last = -1;
03105     for (i = 0; i < nzcnt; i++)
03106     {
03107       c = indx[i];
03108       if (lc_inf[c].delay++ == 0)
03109       {
03110         if (last >= 0)
03111         {
03112           ftranl3_delay2 (f, last);
03113         }
03114         last = c;
03115       }
03116     }
03117     c = last;
03118   } while (c >= 0);
03119 }
03120 
03121 #if 0
03122 static void ftranl3_process (
03123   factor_work * f,
03124   int c,
03125   svector * x)
03126 {
03127   lc_info *lc_inf = f->lc_inf;
03128   EGlpNum_t *work = f->work_coef;
03129   int nzcnt;
03130   int *indx;
03131   int i;
03132   EGlpNum_t *coef;
03133   EGlpNum_t v;
03134 
03135   EGlpNumInitVar (v);
03136 
03137   EGlpNumCopy (v, work[c]);
03138   EGlpNumZero (work[c]);
03139   if (EGlpNumIsNeqqZero (v))
03140   {
03141     x->indx[x->nzcnt] = c;
03142     EGlpNumCopy (x->coef[x->nzcnt], v);
03143     x->nzcnt++;
03144   }
03145   c = lc_inf[c].crank;
03146   nzcnt = lc_inf[c].nzcnt;
03147   indx = f->lcindx + lc_inf[c].cbeg;
03148   coef = f->lccoef + lc_inf[c].cbeg;
03149   for (i = 0; i < nzcnt; i++)
03150   {
03151     c = indx[i];
03152     EGlpNumSubInnProdTo (work[c], v, coef[i]);
03153     if (--lc_inf[c].delay == 0)
03154     {
03155       ftranl3_process (f, c, x);
03156     }
03157   }
03158   EGlpNumClearVar (v);
03159 }
03160 #endif
03161 
03162 static void ftranl3_process2 (
03163   factor_work * f,
03164   int c,
03165   svector * x)
03166 {
03167   lc_info *lc_inf = f->lc_inf;
03168   EGlpNum_t *work = f->work_coef;
03169   int nzcnt;
03170   int *indx;
03171   EGlpNum_t *coef;
03172   int i;
03173   int last;
03174   EGlpNum_t v;
03175 
03176   EGlpNumInitVar (v);
03177 
03178   do
03179   {
03180     EGlpNumCopy (v, work[c]);
03181     EGlpNumZero (work[c]);
03182     if (EGlpNumIsNeqqZero (v))
03183     {
03184       x->indx[x->nzcnt] = c;
03185       EGlpNumCopy (x->coef[x->nzcnt], v);
03186       x->nzcnt++;
03187     }
03188     c = lc_inf[c].crank;
03189     nzcnt = lc_inf[c].nzcnt;
03190     indx = f->lcindx + lc_inf[c].cbeg;
03191     coef = f->lccoef + lc_inf[c].cbeg;
03192     last = -1;
03193     for (i = 0; i < nzcnt; i++)
03194     {
03195       c = indx[i];
03196       EGlpNumSubInnProdTo (work[c], v, coef[i]);
03197       if (--lc_inf[c].delay == 0)
03198       {
03199         if (last >= 0)
03200         {
03201           ftranl3_process2 (f, last, x);
03202         }
03203         last = c;
03204       }
03205     }
03206     c = last;
03207   } while (c >= 0);
03208   EGlpNumClearVar (v);
03209 }
03210 
03211 static void ILLfactor_ftranl3 (
03212   factor_work * f,
03213   svector * a,
03214   svector * x)
03215 {
03216   EGlpNum_t *work = f->work_coef;
03217   int anzcnt = a->nzcnt;
03218   int *aindx = a->indx;
03219   EGlpNum_t *acoef = a->coef;
03220   lc_info *lc_inf = f->lc_inf;
03221   int i;
03222 
03223   for (i = 0; i < anzcnt; i++)
03224   {
03225     if (lc_inf[aindx[i]].delay++ == 0)
03226     {
03227       ftranl3_delay2 (f, aindx[i]);
03228     }
03229     EGlpNumCopy (work[aindx[i]], acoef[i]);
03230   }
03231   x->nzcnt = 0;
03232   for (i = 0; i < anzcnt; i++)
03233   {
03234     if (--lc_inf[aindx[i]].delay == 0)
03235     {
03236       ftranl3_process2 (f, aindx[i], x);
03237     }
03238   }
03239 #ifdef SOLVE_DEBUG
03240   printf ("ILLfactor_ftran x after l3:");
03241   for (i = 0; i < x->nzcnt; i++)
03242   {
03243     printf (" %.3f*%d", EGlpNumToLf (x->coef[i]), x->indx[i]);
03244   }
03245   printf ("\n");
03246   fflush (stdout);
03247 #endif /* SOLVE_DEBUG */
03248 }
03249 
03250 static void ILLfactor_ftrane (
03251   factor_work * f,
03252   EGlpNum_t * a)
03253 {
03254   int *erindx = f->erindx;
03255   EGlpNum_t *ercoef = f->ercoef;
03256   er_info *er_inf = f->er_inf;
03257   int etacnt = f->etacnt;
03258   int beg;
03259   int nzcnt;
03260   int i;
03261   int j;
03262   EGlpNum_t v;
03263 
03264   EGlpNumInitVar (v);
03265 
03266   for (i = 0; i < etacnt; i++)
03267   {
03268     EGlpNumCopy (v, a[er_inf[i].r]);
03269     nzcnt = er_inf[i].nzcnt;
03270     beg = er_inf[i].rbeg;
03271     for (j = 0; j < nzcnt; j++)
03272     {
03273       EGlpNumSubInnProdTo (v, ercoef[beg + j], a[erindx[beg + j]]);
03274     }
03275     EGlpNumCopy (a[er_inf[i].r], v);
03276 #ifdef SOLVE_DEBUG
03277 #if (SOLVE_DEBUG+0 > 1)
03278     printf ("ILLfactor_ftran a after eta %d:", i);
03279     for (j = 0; j < f->dim; j++)
03280     {
03281       printf (" %.3f", EGlpNumToLf (a[j]));
03282     }
03283     printf ("\n");
03284     fflush (stdout);
03285 #endif
03286 #endif /* SOLVE_DEBUG */
03287   }
03288 #ifdef SOLVE_DEBUG
03289 #if (SOLVE_DEBUG+0 <= 1)
03290   printf ("ILLfactor_ftran a after eta:");
03291   for (j = 0; j < f->dim; j++)
03292   {
03293     printf (" %.3f", EGlpNumToLf (a[j]));
03294   }
03295   printf ("\n");
03296   fflush (stdout);
03297 #endif
03298 #endif /* SOLVE_DEBUG */
03299   EGlpNumClearVar (v);
03300 }
03301 
03302 static void ILLfactor_ftrane2 (
03303   factor_work * f,
03304   svector * a)
03305 {
03306   int *erindx = f->erindx;
03307   EGlpNum_t *ercoef = f->ercoef;
03308   er_info *er_inf = f->er_inf;
03309   int etacnt = f->etacnt;
03310   int beg;
03311   int nzcnt;
03312   int anzcnt = a->nzcnt;
03313   int *aindx = a->indx;
03314   EGlpNum_t *acoef = a->coef;
03315   EGlpNum_t *work_coef = f->work_coef;
03316   int *work_indx = f->work_indx;
03317   int i;
03318   int j;
03319   int r;
03320   EGlpNum_t v;
03321 
03322   EGlpNumInitVar (v);
03323 
03324   for (i = 0; i < anzcnt; i++)
03325   {
03326     EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03327     work_indx[aindx[i]] = i + 1;
03328   }
03329   for (i = 0; i < etacnt; i++)
03330   {
03331     r = er_inf[i].r;
03332     EGlpNumCopy (v, work_coef[r]);
03333     nzcnt = er_inf[i].nzcnt;
03334     beg = er_inf[i].rbeg;
03335     for (j = 0; j < nzcnt; j++)
03336     {
03337       EGlpNumSubInnProdTo (v, ercoef[beg + j], work_coef[erindx[beg + j]]);
03338     }
03339     if (EGlpNumIsNeqqZero (v))
03340     {
03341       EGlpNumCopy (work_coef[r], v);
03342       if (work_indx[r] == 0)
03343       {
03344         EGlpNumCopy (acoef[anzcnt], v);
03345         aindx[anzcnt] = r;
03346         work_indx[r] = anzcnt + 1;
03347         anzcnt++;
03348       }
03349       else
03350       {
03351         EGlpNumCopy (acoef[work_indx[r] - 1], v);
03352       }
03353     }
03354     else
03355     {
03356       EGlpNumZero (work_coef[r]);
03357       if (work_indx[r])
03358       {
03359         EGlpNumZero (acoef[work_indx[r] - 1]);
03360       }
03361     }
03362 #ifdef SOLVE_DEBUG
03363 #if (SOLVE_DEBUG+0 > 1)
03364     printf ("ILLfactor_ftran a after eta2 %d:", i);
03365     for (j = 0; j < anzcnt; j++)
03366     {
03367       printf (" %.3f*%d", EGlpNumToLf (acoef[j]), aindx[j]);
03368     }
03369     printf ("\n");
03370     fflush (stdout);
03371 #endif
03372 #endif /* SOLVE_DEBUG */
03373   }
03374   i = 0;
03375   while (i < anzcnt)
03376   {
03377     EGlpNumZero (work_coef[aindx[i]]);
03378     work_indx[aindx[i]] = 0;
03379     if (EGlpNumIsNeqZero (acoef[i], f->fzero_tol))
03380     {
03381       /*if (acoef[i] > fzero_tol || acoef[i] < -fzero_tol) */
03382       i++;
03383     }
03384     else
03385     {
03386       --anzcnt;
03387       EGlpNumCopy (acoef[i], acoef[anzcnt]);
03388       aindx[i] = aindx[anzcnt];
03389     }
03390   }
03391   a->nzcnt = anzcnt;
03392 
03393 #ifdef SOLVE_DEBUG
03394 #if (SOLVE_DEBUG+0 <= 1)
03395   printf ("ILLfactor_ftran a after eta2:");
03396   for (j = 0; j < anzcnt; j++)
03397   {
03398     printf (" %.3f*%d", EGlpNumToLf (acoef[j]), aindx[j]);
03399   }
03400   printf ("\n");
03401   fflush (stdout);
03402 #endif
03403 #endif /* SOLVE_DEBUG */
03404   EGlpNumClearVar (v);
03405 }
03406 
03407 static void ILLfactor_ftranu (
03408   factor_work * f,
03409   EGlpNum_t * a,
03410   svector * x)
03411 {
03412   int *ucindx = f->ucindx;
03413   EGlpNum_t *uccoef = f->uccoef;
03414   uc_info *uc_inf = f->uc_inf;
03415   int *cperm = f->cperm;
03416   int *rperm = f->rperm;
03417   int dim = f->dim;
03418   int xnzcnt = 0;
03419   int *xindx = x->indx;
03420   EGlpNum_t *xcoef = x->coef;
03421   int nzcnt;
03422   int beg;
03423   int i;
03424   int j;
03425   EGlpNum_t v;
03426 
03427   EGlpNumInitVar (v);
03428 
03429   for (i = dim - 1; i >= 0; i--)
03430   {
03431     EGlpNumCopy (v, a[rperm[i]]);
03432     if (EGlpNumIsNeqqZero (v))  /*((v = a[rperm[i]]) != 0.0) */
03433     {
03434       j = cperm[i];
03435       beg = uc_inf[j].cbeg;
03436       EGlpNumDivTo (v, uccoef[beg]);
03437       if (EGlpNumIsNeqZero (v, f->szero_tol))
03438       {
03439         /*if (v > szero_tol || v < -szero_tol) */
03440         xindx[xnzcnt] = j;
03441         EGlpNumCopy (xcoef[xnzcnt], v);
03442         xnzcnt++;
03443       }
03444       nzcnt = uc_inf[j].nzcnt;
03445       for (j = 1; j < nzcnt; j++)
03446       {
03447         EGlpNumSubInnProdTo (a[ucindx[beg + j]], v, uccoef[beg + j]);
03448       }
03449       EGlpNumZero (a[rperm[i]]);
03450     }
03451   }
03452   x->nzcnt = xnzcnt;
03453 #ifdef SOLVE_DEBUG
03454   printf ("ILLfactor_ftran x after u:");
03455   for (j = 0; j < x->nzcnt; j++)
03456   {
03457     printf (" %.3f*%d", EGlpNumToLf (x->coef[j]), x->indx[j]);
03458   }
03459   printf ("\n");
03460   fflush (stdout);
03461 #endif /* SOLVE_DEBUG */
03462   EGlpNumClearVar (v);
03463 }
03464 
03465 
03466 #if 0
03467 static void ftranu3_delay (
03468   factor_work * f,
03469   int c)
03470 {
03471   uc_info *uc_inf = f->uc_inf;
03472   int nzcnt;
03473   int *indx;
03474   int i;
03475 
03476   c = f->cperm[f->rrank[c]];
03477   nzcnt = uc_inf[c].nzcnt;
03478   indx = f->ucindx + uc_inf[c].cbeg;
03479   for (i = 1; i < nzcnt; i++)
03480   {
03481     c = indx[i];
03482     if (uc_inf[c].delay++ == 0)
03483     {
03484       ftranu3_delay (f, c);
03485     }
03486   }
03487 }
03488 #endif
03489 
03490 static void ftranu3_delay2 (
03491   factor_work * f,
03492   int c)
03493 {
03494   uc_info *uc_inf = f->uc_inf;
03495   int nzcnt;
03496   int *indx;
03497   int i;
03498   int last;
03499 
03500   do
03501   {
03502     c = f->cperm[f->rrank[c]];
03503     nzcnt = uc_inf[c].nzcnt;
03504     indx = f->ucindx + uc_inf[c].cbeg;
03505     last = -1;
03506     for (i = 1; i < nzcnt; i++)
03507     {
03508       c = indx[i];
03509       if (uc_inf[c].delay++ == 0)
03510       {
03511         if (last >= 0)
03512         {
03513           ftranu3_delay2 (f, last);
03514         }
03515         last = c;
03516       }
03517     }
03518     c = last;
03519   } while (c >= 0);
03520 }
03521 
03522 #if 0
03523 static void ftranu3_process (
03524   factor_work * f,
03525   int c,
03526   svector * x)
03527 {
03528   uc_info *uc_inf = f->uc_inf;
03529   EGlpNum_t *work = f->work_coef;
03530   int nzcnt;
03531   int *indx;
03532   EGlpNum_t *coef;
03533   int i;
03534   EGlpNum_t v;
03535 
03536   EGlpNumInitVar (v);
03537 
03538   EGlpNumCopy (v, work[c]);
03539   EGlpNumZero (work[c]);
03540   c = f->cperm[f->rrank[c]];
03541   nzcnt = uc_inf[c].nzcnt;
03542   indx = f->ucindx + uc_inf[c].cbeg;
03543   coef = f->uccoef + uc_inf[c].cbeg;
03544   EGlpNumDivTo (v, coef[0]);
03545   if (EGlpNumIsNeqZero (v, f->szero_tol))
03546     /*if (v > szero_tol || v < -szero_tol) */
03547   {
03548     x->indx[x->nzcnt] = c;
03549     EGlpNumCopy (x->coef[x->nzcnt], v);
03550     x->nzcnt++;
03551   }
03552   for (i = 1; i < nzcnt; i++)
03553   {
03554     c = indx[i];
03555     EGlpNumSubInnProdTo (work[c], v, coef[i]);
03556     if (--uc_inf[c].delay == 0)
03557     {
03558       ftranu3_process (f, c, x);
03559     }
03560   }
03561   EGlpNumClearVar (v);
03562 }
03563 #endif
03564 
03565 static void ftranu3_process2 (
03566   factor_work * f,
03567   int c,
03568   svector * x)
03569 {
03570   uc_info *uc_inf = f->uc_inf;
03571   EGlpNum_t *work = f->work_coef;
03572   int nzcnt;
03573   int *indx;
03574   EGlpNum_t *coef;
03575   int i;
03576   int last;
03577   EGlpNum_t v;
03578 
03579   EGlpNumInitVar (v);
03580 
03581   do
03582   {
03583     EGlpNumCopy (v, work[c]);
03584     EGlpNumZero (work[c]);
03585     c = f->cperm[f->rrank[c]];
03586     nzcnt = uc_inf[c].nzcnt;
03587     indx = f->ucindx + uc_inf[c].cbeg;
03588     coef = f->uccoef + uc_inf[c].cbeg;
03589     EGlpNumDivTo (v, coef[0]);
03590     if (EGlpNumIsNeqZero (v, f->szero_tol))
03591       /*if (v > szero_tol || v < -szero_tol) */
03592     {
03593       x->indx[x->nzcnt] = c;
03594       EGlpNumCopy (x->coef[x->nzcnt], v);
03595       x->nzcnt++;
03596     }
03597     last = -1;
03598     for (i = 1; i < nzcnt; i++)
03599     {
03600       c = indx[i];
03601       EGlpNumSubInnProdTo (work[c], v, coef[i]);
03602       if (--uc_inf[c].delay == 0)
03603       {
03604         if (last >= 0)
03605         {
03606           ftranu3_process2 (f, last, x);
03607         }
03608         last = c;
03609       }
03610     }
03611     c = last;
03612   } while (c >= 0);
03613   EGlpNumClearVar (v);
03614 }
03615 
03616 static void ILLfactor_ftranu3 (
03617   factor_work * f,
03618   svector * a,
03619   svector * x)
03620 {
03621   EGlpNum_t *work = f->work_coef;
03622   int anzcnt = a->nzcnt;
03623   int *aindx = a->indx;
03624   EGlpNum_t *acoef = a->coef;
03625   uc_info *uc_inf = f->uc_inf;
03626   int i;
03627 
03628   for (i = 0; i < anzcnt; i++)
03629   {
03630     if (uc_inf[aindx[i]].delay++ == 0)
03631     {
03632       ftranu3_delay2 (f, aindx[i]);
03633     }
03634     EGlpNumCopy (work[aindx[i]], acoef[i]);
03635   }
03636   x->nzcnt = 0;
03637   for (i = 0; i < anzcnt; i++)
03638   {
03639     if (--uc_inf[aindx[i]].delay == 0)
03640     {
03641       ftranu3_process2 (f, aindx[i], x);
03642     }
03643   }
03644 #ifdef SOLVE_DEBUG
03645   printf ("ILLfactor_ftran x after u3:");
03646   for (i = 0; i < x->nzcnt; i++)
03647   {
03648     printf (" %.3f*%d", EGlpNumToLf (x->coef[i]), x->indx[i]);
03649   }
03650   printf ("\n");
03651   fflush (stdout);
03652 #endif /* SOLVE_DEBUG */
03653 }
03654 
03655 /* ILLfactor_ftran solves Bx=a for x */
03656 void ILLfactor_ftran (
03657   factor_work * f,
03658   svector * a,
03659   svector * x)
03660 {
03661   int i;
03662   int nzcnt;
03663   int sparse;
03664   int *aindx;
03665   EGlpNum_t *acoef;
03666   EGlpNum_t *work_coef = f->work_coef;
03667 
03668 #ifdef RECORD
03669   {
03670     EGioPrintf (fsave, "f %d", a->nzcnt);
03671     for (i = 0; i < a->nzcnt; i++)
03672     {
03673       EGioPrintf (fsave, " %d %.16e", a->indx[i], EGlpNumToLf (a->coef[i]));
03674     }
03675     EGioPrintf (fsave, "\n");
03676     EGioFlush (fsave);
03677   }
03678 #endif /* RECORD */
03679 #ifdef DEBUG_FACTOR
03680   {
03681     printf ("ILLfactor_ftran a:");
03682     for (i = 0; i < a->nzcnt; i++)
03683     {
03684       printf (" %d %la", a->indx[i], EGlpNumToLf (a->coef[i]));
03685     }
03686     printf ("\n");
03687     fflush (stdout);
03688   }
03689 #endif /* DEBUG_FACTOR */
03690 
03691   if (a->nzcnt >= SPARSE_FACTOR * f->dim)
03692   {
03693     nzcnt = a->nzcnt;
03694     aindx = a->indx;
03695     acoef = a->coef;
03696     for (i = 0; i < nzcnt; i++)
03697     {
03698       EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03699     }
03700     sparse = 0;
03701   }
03702   else
03703   {
03704     sparse = 1;
03705   }
03706 
03707   if (sparse)
03708   {
03709     ILLfactor_ftranl3 (f, a, &f->xtmp);
03710     if (f->xtmp.nzcnt >= SPARSE_FACTOR * f->dim)
03711     {
03712       nzcnt = f->xtmp.nzcnt;
03713       aindx = f->xtmp.indx;
03714       acoef = f->xtmp.coef;
03715 
03716       for (i = 0; i < nzcnt; i++)
03717       {
03718         EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03719       }
03720       sparse = 0;
03721     }
03722   }
03723   else
03724   {
03725     ILLfactor_ftranl (f, work_coef);
03726   }
03727 
03728   if (sparse)
03729   {
03730     ILLfactor_ftrane2 (f, &f->xtmp);
03731     if (f->xtmp.nzcnt >= SPARSE_FACTOR * f->dim)
03732     {
03733       nzcnt = f->xtmp.nzcnt;
03734       aindx = f->xtmp.indx;
03735       acoef = f->xtmp.coef;
03736 
03737       for (i = 0; i < nzcnt; i++)
03738       {
03739         EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03740       }
03741       sparse = 0;
03742     }
03743   }
03744   else
03745   {
03746     ILLfactor_ftrane (f, work_coef);
03747   }
03748 
03749   if (sparse)
03750   {
03751     ILLfactor_ftranu3 (f, &f->xtmp, x);
03752   }
03753   else
03754   {
03755     ILLfactor_ftranu (f, work_coef, x);
03756   }
03757 
03758 #ifdef SORT_RESULTS
03759   sort_vector (x);
03760 #endif
03761 
03762 #ifdef DEBUG_FACTOR
03763   {
03764     printf ("ILLfactor_ftran x:");
03765     for (i = 0; i < x->nzcnt; i++)
03766     {
03767       printf (" %d %la", x->indx[i], EGlpNumToLf (x->coef[i]));
03768     }
03769     printf ("\n");
03770     fflush (stdout);
03771   }
03772 #endif /* DEBUG_FACTOR */
03773   return;
03774 }
03775 
03776 /* ILLfactor_ftran_update solves Bx=a for x, and also returns upd, where Ux=upd */
03777 void ILLfactor_ftran_update (
03778   factor_work * f,
03779   svector * a,
03780   svector * upd,
03781   svector * x)
03782 {
03783   int i;
03784   int nzcnt;
03785   int dim;
03786   int sparse;
03787   int *aindx;
03788   EGlpNum_t *acoef;
03789   EGlpNum_t *work_coef = f->work_coef;
03790 
03791 #ifdef RECORD
03792   {
03793     EGioPrintf (fsave, "F %d", a->nzcnt);
03794     for (i = 0; i < a->nzcnt; i++)
03795     {
03796       EGioPrintf (fsave, " %d %.16e", a->indx[i], EGlpNumToLf (a->coef[i]));
03797     }
03798     EGioPrintf (fsave, "\n");
03799     EGioFlush (fsave);
03800   }
03801 #endif /* RECORD */
03802 #ifdef DEBUG_FACTOR
03803   {
03804     printf ("ILLfactor_ftran_update a:");
03805     for (i = 0; i < a->nzcnt; i++)
03806     {
03807       printf (" %d %.3f", a->indx[i], EGlpNumToLf (a->coef[i]));
03808     }
03809     printf ("\n");
03810     fflush (stdout);
03811   }
03812 #endif /* DEBUG_FACTOR */
03813 
03814   if (a->nzcnt >= SPARSE_FACTOR * f->dim)
03815   {
03816     aindx = a->indx;
03817     acoef = a->coef;
03818     nzcnt = a->nzcnt;
03819 
03820     for (i = 0; i < nzcnt; i++)
03821     {
03822       EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03823     }
03824     sparse = 0;
03825   }
03826   else
03827   {
03828     sparse = 1;
03829   }
03830 
03831   if (sparse)
03832   {
03833     ILLfactor_ftranl3 (f, a, upd);
03834     if (upd->nzcnt >= SPARSE_FACTOR * f->dim)
03835     {
03836       nzcnt = upd->nzcnt;
03837       aindx = upd->indx;
03838       acoef = upd->coef;
03839 
03840       for (i = 0; i < nzcnt; i++)
03841       {
03842         EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03843       }
03844       sparse = 0;
03845     }
03846   }
03847   else
03848   {
03849     ILLfactor_ftranl (f, work_coef);
03850   }
03851 
03852   if (sparse)
03853   {
03854     ILLfactor_ftrane2 (f, upd);
03855     if (upd->nzcnt >= SPARSE_FACTOR * f->dim)
03856     {
03857       nzcnt = upd->nzcnt;
03858       aindx = upd->indx;
03859       acoef = upd->coef;
03860 
03861       for (i = 0; i < nzcnt; i++)
03862       {
03863         EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
03864       }
03865       sparse = 0;
03866     }
03867   }
03868   else
03869   {
03870     ILLfactor_ftrane (f, work_coef);
03871     nzcnt = 0;
03872     dim = f->dim;
03873     aindx = upd->indx;
03874     acoef = upd->coef;
03875     for (i = 0; i < dim; i++)
03876     {
03877       if (EGlpNumIsNeqqZero (work_coef[i]))
03878       {
03879         if (EGlpNumIsNeqZero (work_coef[i], f->szero_tol))
03880           /*if(work_coef[i] > szero_tol || work_coef[i] < -szero_tol) */
03881         {
03882           aindx[nzcnt] = i;
03883           EGlpNumCopy (acoef[nzcnt], work_coef[i]);
03884           nzcnt++;
03885         }
03886       }
03887     }
03888     upd->nzcnt = nzcnt;
03889   }
03890 
03891   if (sparse)
03892   {
03893     ILLfactor_ftranu3 (f, upd, x);
03894   }
03895   else
03896   {
03897     ILLfactor_ftranu (f, work_coef, x);
03898   }
03899 
03900 #ifdef SORT_RESULTS
03901   sort_vector (upd);
03902   sort_vector (x);
03903 #endif
03904 
03905 #ifdef DEBUG_FACTOR
03906   {
03907     printf ("ILLfactor_ftran update x:");
03908     for (i = 0; i < x->nzcnt; i++)
03909     {
03910       printf (" %d %.3f", x->indx[i], EGlpNumToLf (x->coef[i]));
03911     }
03912     printf ("\n");
03913     fflush (stdout);
03914   }
03915 #endif /* DEBUG_FACTOR */
03916 }
03917 
03918 
03919 static void ILLfactor_btranl2 (
03920   factor_work * f,
03921   EGlpNum_t * x)
03922 {
03923   int *lrindx = f->lrindx;
03924   EGlpNum_t *lrcoef = f->lrcoef;
03925   lr_info *lr_inf = f->lr_inf;
03926   int dim = f->dim;
03927   int nzcnt;
03928   int beg;
03929   int i;
03930   int j;
03931   EGlpNum_t v;
03932 
03933   EGlpNumInitVar (v);
03934 
03935   for (i = dim - 1; i >= 0; i--)
03936   {
03937 #ifdef SOLVE_DEBUG
03938 #if (SOLVE_DEBUG+0 > 1)
03939     printf ("ILLfactor_btran x before l2 %d:", i);
03940     for (j = 0; j < f->dim; j++)
03941     {
03942       printf (" %.3f", EGlpNumToLf (x[j]));
03943     }
03944     printf ("\n");
03945     fflush (stdout);
03946 #endif
03947 #endif /* SOLVE_DEBUG */
03948     EGlpNumCopy (v, x[lr_inf[i].r]);
03949     if (EGlpNumIsNeqqZero (v))
03950     {
03951       nzcnt = lr_inf[i].nzcnt;
03952       beg = lr_inf[i].rbeg;
03953       for (j = 0; j < nzcnt; j++)
03954       {
03955         EGlpNumSubInnProdTo (x[lrindx[beg + j]], v, lrcoef[beg + j]);
03956       }
03957     }
03958   }
03959 #ifdef SOLVE_DEBUG
03960 #if (SOLVE_DEBUG+0 <= 1)
03961   printf ("ILLfactor_btran x after l2:");
03962   for (j = 0; j < f->dim; j++)
03963   {
03964     printf (" %.3f", EGlpNumToLf (x[j]));
03965   }
03966   printf ("\n");
03967   fflush (stdout);
03968 #endif
03969 #endif /* SOLVE_DEBUG */
03970   EGlpNumClearVar (v);
03971 }
03972 
03973 #if 0
03974 static void btranl3_delay (
03975   factor_work * f,
03976   int r)
03977 {
03978   lr_info *lr_inf = f->lr_inf;
03979   int nzcnt;
03980   int *indx;
03981   int i;
03982 
03983   r = lr_inf[r].rrank;
03984   nzcnt = lr_inf[r].nzcnt;
03985   indx = f->lrindx + lr_inf[r].rbeg;
03986   for (i = 0; i < nzcnt; i++)
03987   {
03988     r = indx[i];
03989     if (lr_inf[r].delay++ == 0)
03990     {
03991       btranl3_delay (f, r);
03992     }
03993   }
03994 }
03995 #endif
03996 
03997 static void btranl3_delay2 (
03998   factor_work * f,
03999   int r)
04000 {
04001   lr_info *lr_inf = f->lr_inf;
04002   int nzcnt;
04003   int *indx;
04004   int i;
04005   int last;
04006 
04007   do
04008   {
04009     r = lr_inf[r].rrank;
04010     nzcnt = lr_inf[r].nzcnt;
04011     indx = f->lrindx + lr_inf[r].rbeg;
04012     last = -1;
04013     for (i = 0; i < nzcnt; i++)
04014     {
04015       r = indx[i];
04016       if (lr_inf[r].delay++ == 0)
04017       {
04018         if (last >= 0)
04019         {
04020           btranl3_delay2 (f, last);
04021         }
04022         last = r;
04023       }
04024     }
04025     r = last;
04026   } while (r >= 0);
04027 }
04028 
04029 #if 0
04030 static void btranl3_process (
04031   factor_work * f,
04032   int r,
04033   svector * x)
04034 {
04035   lr_info *lr_inf = f->lr_inf;
04036   EGlpNum_t *work = f->work_coef;
04037   int nzcnt;
04038   int *indx;
04039   EGlpNum_t *coef;
04040   int i;
04041   EGlpNum_t v;
04042 
04043   EGlpNumInitVar (v);
04044 
04045   EGlpNumCopy (v, work[r]);
04046   EGlpNumZero (work[r]);
04047   if (EGlpNumIsNeqZero (v, f->szero_tol))
04048     /*if (v > szero_tol || v < -szero_tol) */
04049   {
04050     x->indx[x->nzcnt] = r;
04051     EGlpNumCopy (x->coef[x->nzcnt], v);
04052     x->nzcnt++;
04053   }
04054   r = lr_inf[r].rrank;
04055   nzcnt = lr_inf[r].nzcnt;
04056   indx = f->lrindx + lr_inf[r].rbeg;
04057   coef = f->lrcoef + lr_inf[r].rbeg;
04058   for (i = 0; i < nzcnt; i++)
04059   {
04060     r = indx[i];
04061     EGlpNumSubInnProdTo (work[r], v, coef[i]);
04062     if (--lr_inf[r].delay == 0)
04063     {
04064       btranl3_process (f, r, x);
04065     }
04066   }
04067   EGlpNumClearVar (v);
04068 }
04069 #endif
04070 
04071 static void btranl3_process2 (
04072   factor_work * f,
04073   int r,
04074   svector * x)
04075 {
04076   lr_info *lr_inf = f->lr_inf;
04077   EGlpNum_t *work = f->work_coef;
04078   int nzcnt;
04079   int *indx;
04080   EGlpNum_t *coef;
04081   int i;
04082   int last;
04083   EGlpNum_t v;
04084 
04085   EGlpNumInitVar (v);
04086 
04087   do
04088   {
04089     EGlpNumCopy (v, work[r]);
04090     EGlpNumZero (work[r]);
04091     if (EGlpNumIsNeqZero (v, f->szero_tol))
04092       /*if (v > szero_tol || v < -szero_tol) */
04093     {
04094       x->indx[x->nzcnt] = r;
04095       EGlpNumCopy (x->coef[x->nzcnt], v);
04096       x->nzcnt++;
04097     }
04098     r = lr_inf[r].rrank;
04099     nzcnt = lr_inf[r].nzcnt;
04100     indx = f->lrindx + lr_inf[r].rbeg;
04101     coef = f->lrcoef + lr_inf[r].rbeg;
04102     last = -1;
04103     for (i = 0; i < nzcnt; i++)
04104     {
04105       r = indx[i];
04106       EGlpNumSubInnProdTo (work[r], v, coef[i]);
04107       if (--lr_inf[r].delay == 0)
04108       {
04109         if (last >= 0)
04110         {
04111           btranl3_process2 (f, last, x);
04112         }
04113         last = r;
04114       }
04115     }
04116     r = last;
04117   } while (r >= 0);
04118   EGlpNumClearVar (v);
04119 }
04120 
04121 static void ILLfactor_btranl3 (
04122   factor_work * f,
04123   svector * a,
04124   svector * x)
04125 {
04126   EGlpNum_t *work = f->work_coef;
04127   int anzcnt = a->nzcnt;
04128   int *aindx = a->indx;
04129   EGlpNum_t *acoef = a->coef;
04130   lr_info *lr_inf = f->lr_inf;
04131   int i;
04132 
04133   for (i = 0; i < anzcnt; i++)
04134   {
04135     if (lr_inf[aindx[i]].delay++ == 0)
04136     {
04137       btranl3_delay2 (f, aindx[i]);
04138     }
04139     EGlpNumCopy (work[aindx[i]], acoef[i]);
04140   }
04141   x->nzcnt = 0;
04142   for (i = 0; i < anzcnt; i++)
04143   {
04144     if (--lr_inf[aindx[i]].delay == 0)
04145     {
04146       btranl3_process2 (f, aindx[i], x);
04147     }
04148   }
04149 #ifdef SOLVE_DEBUG
04150   printf ("ILLfactor_btran x after l3:");
04151   for (i = 0; i < x->nzcnt; i++)
04152   {
04153     printf (" %.3f*%d", EGlpNumToLf (x->coef[i]), x->indx[i]);
04154   }
04155   printf ("\n");
04156   fflush (stdout);
04157 #endif /* SOLVE_DEBUG */
04158 }
04159 
04160 static void ILLfactor_btrane (
04161   factor_work * f,
04162   EGlpNum_t * x)
04163 {
04164   int *erindx = f->erindx;
04165   EGlpNum_t *ercoef = f->ercoef;
04166   er_info *er_inf = f->er_inf;
04167   int etacnt = f->etacnt;
04168   int beg;
04169   int nzcnt;
04170   int i;
04171   int j;
04172   EGlpNum_t v;
04173 
04174   EGlpNumInitVar (v);
04175 
04176   for (i = etacnt - 1; i >= 0; i--)
04177   {
04178 #ifdef SOLVE_DEBUG
04179 #if (SOLVE_DEBUG+0 > 1)
04180     printf ("ILLfactor_btran x before eta %d:", i);
04181     for (j = 0; j < f->dim; j++)
04182     {
04183       printf (" %.3f", EGlpNumToLf (x[j]));
04184     }
04185     printf ("\n");
04186     fflush (stdout);
04187 #endif
04188 #endif /* SOLVE_DEBUG */
04189     EGlpNumCopy (v, x[er_inf[i].r]);
04190     if (EGlpNumIsNeqqZero (v))
04191     {
04192       nzcnt = er_inf[i].nzcnt;
04193       beg = er_inf[i].rbeg;
04194       for (j = 0; j < nzcnt; j++)
04195       {
04196         EGlpNumSubInnProdTo (x[erindx[beg + j]], v, ercoef[beg + j]);
04197       }
04198     }
04199   }
04200 #ifdef SOLVE_DEBUG
04201 #if (SOLVE_DEBUG+0 <= 1)
04202   printf ("ILLfactor_btran x after eta:");
04203   for (j = 0; j < f->dim; j++)
04204   {
04205     printf (" %.3f", EGlpNumToLf (x[j]));
04206   }
04207   printf ("\n");
04208   fflush (stdout);
04209 #endif
04210 #endif /* SOLVE_DEBUG */
04211   EGlpNumClearVar (v);
04212 }
04213 
04214 static void ILLfactor_btrane2 (
04215   factor_work * f,
04216   svector * x)
04217 {
04218   int *erindx = f->erindx;
04219   EGlpNum_t *ercoef = f->ercoef;
04220   er_info *er_inf = f->er_inf;
04221   int etacnt = f->etacnt;
04222   int beg;
04223   int nzcnt;
04224   int xnzcnt = x->nzcnt;
04225   int *xindx = x->indx;
04226   EGlpNum_t *xcoef = x->coef;
04227   EGlpNum_t *work_coef = f->work_coef;
04228   int *work_indx = f->work_indx;
04229   int i;
04230   int j;
04231   EGlpNum_t v;
04232 
04233   EGlpNumInitVar (v);
04234 
04235   for (i = 0; i < xnzcnt; i++)
04236   {
04237     EGlpNumCopy (work_coef[xindx[i]], xcoef[i]);
04238     work_indx[xindx[i]] = i + 1;
04239   }
04240   for (i = etacnt - 1; i >= 0; i--)
04241   {
04242 #ifdef SOLVE_DEBUG
04243 #if (SOLVE_DEBUG+0 > 1)
04244     printf ("ILLfactor_btran x before eta2 %d:", i);
04245     for (j = 0; j < xnzcnt; j++)
04246     {
04247       printf (" %.3f*%d", EGlpNumToLf (work_coef[xindx[j]]), xindx[j]);
04248     }
04249     printf ("\n");
04250     fflush (stdout);
04251 #endif
04252 #endif /* SOLVE_DEBUG */
04253     EGlpNumCopy (v, work_coef[er_inf[i].r]);
04254     if (EGlpNumIsNeqqZero (v))
04255     {
04256       nzcnt = er_inf[i].nzcnt;
04257       beg = er_inf[i].rbeg;
04258       for (j = 0; j < nzcnt; j++)
04259       {
04260         if (work_indx[erindx[beg + j]] == 0)
04261         {
04262           work_indx[erindx[beg + j]] = xnzcnt;
04263           xindx[xnzcnt++] = erindx[beg + j];
04264         }
04265         EGlpNumSubInnProdTo (work_coef[erindx[beg + j]], v, ercoef[beg + j]);
04266       }
04267     }
04268   }
04269 
04270   j = 0;
04271   while (j < xnzcnt)
04272   {
04273     EGlpNumCopy (xcoef[j], work_coef[xindx[j]]);
04274     EGlpNumZero (work_coef[xindx[j]]);
04275     work_indx[xindx[j]] = 0;
04276     if (!EGlpNumIsNeqqZero (xcoef[j]))
04277     {
04278       --xnzcnt;
04279       xindx[j] = xindx[xnzcnt];
04280     }
04281     else
04282     {
04283       j++;
04284     }
04285   }
04286   x->nzcnt = xnzcnt;
04287 
04288 #ifdef SOLVE_DEBUG
04289 #if (SOLVE_DEBUG+0 <= 1)
04290   printf ("ILLfactor_btran x after eta2:");
04291   for (j = 0; j < xnzcnt; j++)
04292   {
04293     printf (" %.3f*%d", EGlpNumToLf (xcoef[j]), xindx[j]);
04294   }
04295   printf ("\n");
04296   fflush (stdout);
04297 #endif
04298 #endif /* SOLVE_DEBUG */
04299   EGlpNumClearVar (v);
04300 }
04301 
04302 static void ILLfactor_btranu (
04303   factor_work * f,
04304   EGlpNum_t * a,
04305   svector * x)
04306 {
04307   int *urindx = f->urindx;
04308   EGlpNum_t *urcoef = f->urcoef;
04309   ur_info *ur_inf = f->ur_inf;
04310   int *rperm = f->rperm;
04311   int *cperm = f->cperm;
04312   int dim = f->dim;
04313   int xnzcnt = 0;
04314   int *xindx = x->indx;
04315   EGlpNum_t *xcoef = x->coef;
04316   int nzcnt;
04317   int beg;
04318   int i;
04319   int j;
04320   EGlpNum_t v;
04321 
04322   EGlpNumInitVar (v);
04323 
04324   for (i = 0; i < dim; i++)
04325   {
04326     EGlpNumCopy (v, a[cperm[i]]);
04327     if (EGlpNumIsNeqqZero (v))
04328     {
04329       j = rperm[i];
04330       beg = ur_inf[j].rbeg;
04331       EGlpNumDivTo (v, urcoef[beg]);
04332       if (EGlpNumIsNeqZero (v, f->szero_tol)) /*
04333                                                * if (v > szero_tol || v < -szero_tol) */
04334       {
04335         xindx[xnzcnt] = j;
04336         EGlpNumCopy (xcoef[xnzcnt], v);
04337         xnzcnt++;
04338       }
04339       nzcnt = ur_inf[j].nzcnt;
04340       for (j = 1; j < nzcnt; j++)
04341       {
04342         EGlpNumSubInnProdTo (a[urindx[beg + j]], v, urcoef[beg + j]);
04343       }
04344       EGlpNumZero (a[cperm[i]]);
04345     }
04346   }
04347   x->nzcnt = xnzcnt;
04348 #ifdef SOLVE_DEBUG
04349   printf ("ILLfactor_btran x after u:");
04350   for (i = 0; i < x->nzcnt; i++)
04351   {
04352     printf (" %.3f*%d", EGlpNumToLf (x->coef[i]), x->indx[i]);
04353   }
04354   printf ("\n");
04355   fflush (stdout);
04356 #endif /* SOLVE_DEBUG */
04357   EGlpNumClearVar (v);
04358 }
04359 
04360 
04361 #if 0
04362 static void btranu3_delay (
04363   factor_work * f,
04364   int r)
04365 {
04366   ur_info *ur_inf = f->ur_inf;
04367   int nzcnt;
04368   int *indx;
04369   int i;
04370 
04371   r = f->rperm[f->crank[r]];
04372   nzcnt = ur_inf[r].nzcnt;
04373   indx = f->urindx + ur_inf[r].rbeg;
04374   for (i = 1; i < nzcnt; i++)
04375   {
04376     r = indx[i];
04377     if (ur_inf[r].delay++ == 0)
04378     {
04379       btranu3_delay (f, r);
04380     }
04381   }
04382 }
04383 #endif
04384 
04385 static void btranu3_delay2 (
04386   factor_work * f,
04387   int r)
04388 {
04389   ur_info *ur_inf = f->ur_inf;
04390   int nzcnt;
04391   int *indx;
04392   int i;
04393   int last;
04394 
04395   do
04396   {
04397     r = f->rperm[f->crank[r]];
04398     nzcnt = ur_inf[r].nzcnt;
04399     indx = f->urindx + ur_inf[r].rbeg;
04400     last = -1;
04401     for (i = 1; i < nzcnt; i++)
04402     {
04403       r = indx[i];
04404       if (ur_inf[r].delay++ == 0)
04405       {
04406         if (last >= 0)
04407         {
04408           btranu3_delay2 (f, last);
04409         }
04410         last = r;
04411       }
04412     }
04413     r = last;
04414   } while (r >= 0);
04415 }
04416 
04417 #if 0
04418 static void btranu3_process (
04419   factor_work * f,
04420   int r,
04421   svector * x)
04422 {
04423   ur_info *ur_inf = f->ur_inf;
04424   EGlpNum_t *work = f->work_coef;
04425   int nzcnt;
04426   int *indx;
04427   EGlpNum_t *coef;
04428   int i;
04429   EGlpNum_t v;
04430 
04431   EGlpNumInitVar (v);
04432 
04433   EGlpNumCopy (v, work[r]);
04434   EGlpNumZero (work[r]);
04435   r = f->rperm[f->crank[r]];
04436   nzcnt = ur_inf[r].nzcnt;
04437   indx = f->urindx + ur_inf[r].rbeg;
04438   coef = f->urcoef + ur_inf[r].rbeg;
04439   EGlpNumDivTo (v, coef[0]);
04440   if (EGlpNumIsNeqqZero (v))
04441   {
04442     x->indx[x->nzcnt] = r;
04443     EGlpNumCopy (x->coef[x->nzcnt], v);
04444     x->nzcnt++;
04445   }
04446   for (i = 1; i < nzcnt; i++)
04447   {
04448     r = indx[i];
04449     EGlpNumSubInnProdTo (work[r], v, coef[i]);
04450     if (--ur_inf[r].delay == 0)
04451     {
04452       btranu3_process (f, r, x);
04453     }
04454   }
04455   EGlpNumClearVar (v);
04456 }
04457 #endif
04458 
04459 static void btranu3_process2 (
04460   factor_work * f,
04461   int r,
04462   svector * x)
04463 {
04464   ur_info *ur_inf = f->ur_inf;
04465   EGlpNum_t *work = f->work_coef;
04466   int nzcnt;
04467   int *indx;
04468   EGlpNum_t *coef;
04469   int i;
04470   int last;
04471   EGlpNum_t v;
04472 
04473   EGlpNumInitVar (v);
04474 
04475   do
04476   {
04477     EGlpNumCopy (v, work[r]);
04478     EGlpNumZero (work[r]);
04479     r = f->rperm[f->crank[r]];
04480     nzcnt = ur_inf[r].nzcnt;
04481     indx = f->urindx + ur_inf[r].rbeg;
04482     coef = f->urcoef + ur_inf[r].rbeg;
04483     EGlpNumDivTo (v, coef[0]);
04484     if (EGlpNumIsNeqqZero (v))
04485     {
04486       x->indx[x->nzcnt] = r;
04487       EGlpNumCopy (x->coef[x->nzcnt], v);
04488       x->nzcnt++;
04489     }
04490     last = -1;
04491     for (i = 1; i < nzcnt; i++)
04492     {
04493       r = indx[i];
04494       EGlpNumSubInnProdTo (work[r], v, coef[i]);
04495       if (--ur_inf[r].delay == 0)
04496       {
04497         if (last >= 0)
04498         {
04499           btranu3_process2 (f, last, x);
04500         }
04501         last = r;
04502       }
04503     }
04504     r = last;
04505   } while (r >= 0);
04506   EGlpNumClearVar (v);
04507 }
04508 
04509 static void ILLfactor_btranu3 (
04510   factor_work * f,
04511   svector * a,
04512   svector * x)
04513 {
04514   EGlpNum_t *work = f->work_coef;
04515   int anzcnt = a->nzcnt;
04516   int *aindx = a->indx;
04517   EGlpNum_t *acoef = a->coef;
04518   ur_info *ur_inf = f->ur_inf;
04519   int i;
04520 
04521   for (i = 0; i < anzcnt; i++)
04522   {
04523     if (ur_inf[aindx[i]].delay++ == 0)
04524     {
04525       btranu3_delay2 (f, aindx[i]);
04526     }
04527     EGlpNumCopy (work[aindx[i]], acoef[i]);
04528   }
04529   x->nzcnt = 0;
04530   for (i = 0; i < anzcnt; i++)
04531   {
04532     if (--ur_inf[aindx[i]].delay == 0)
04533     {
04534       btranu3_process2 (f, aindx[i], x);
04535     }
04536   }
04537 #ifdef SOLVE_DEBUG
04538   printf ("ILLfactor_btran x after u3:");
04539   for (i = 0; i < x->nzcnt; i++)
04540   {
04541     printf (" %.3f*%d", EGlpNumToLf (x->coef[i]), x->indx[i]);
04542   }
04543   printf ("\n");
04544   fflush (stdout);
04545 #endif /* SOLVE_DEBUG */
04546 }
04547 
04548 /* ILLfactor_btran solves x^tB=a^t (or, B^t x = a) for x */
04549 void ILLfactor_btran (
04550   factor_work * f,
04551   svector * a,
04552   svector * x)
04553 {
04554   int i;
04555   int nzcnt;
04556   int sparse;
04557   int *aindx = a->indx;
04558   EGlpNum_t *acoef = a->coef;
04559   EGlpNum_t *work_coef = f->work_coef;
04560   int dim = f->dim;
04561 
04562 #ifdef RECORD
04563   {
04564     EGioPrintf (fsave, "b %d", a->nzcnt);
04565     for (i = 0; i < a->nzcnt; i++)
04566     {
04567       EGioPrintf (fsave, " %d %.16e", a->indx[i], EGlpNumToLf (a->coef[i]));
04568     }
04569     EGioPrintf (fsave, "\n");
04570     EGioFlush (fsave);
04571   }
04572 #endif /* RECORD */
04573 #ifdef DEBUG_FACTOR
04574   {
04575     printf ("ILLfactor_btran a:");
04576     for (i = 0; i < a->nzcnt; i++)
04577     {
04578       printf (" %d %.3f", a->indx[i], EGlpNumToLf (a->coef[i]));
04579     }
04580     printf ("\n");
04581     fflush (stdout);
04582   }
04583 #endif /* DEBUG_FACTOR */
04584 
04585   if (a->nzcnt >= SPARSE_FACTOR * f->dim)
04586   {
04587     aindx = a->indx;
04588     acoef = a->coef;
04589     work_coef = f->work_coef;
04590     nzcnt = a->nzcnt;
04591     for (i = 0; i < nzcnt; i++)
04592     {
04593       EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
04594     }
04595     sparse = 0;
04596   }
04597   else
04598   {
04599     sparse = 1;
04600   }
04601 
04602   if (sparse)
04603   {
04604     ILLfactor_btranu3 (f, a, &f->xtmp);
04605   }
04606   else
04607   {
04608     ILLfactor_btranu (f, work_coef, &f->xtmp);
04609   }
04610 
04611   if (f->xtmp.nzcnt >= SPARSE_FACTOR * f->dim)
04612   {
04613     aindx = f->xtmp.indx;
04614     acoef = f->xtmp.coef;
04615     work_coef = f->work_coef;
04616     nzcnt = f->xtmp.nzcnt;
04617     for (i = 0; i < nzcnt; i++)
04618     {
04619       EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
04620     }
04621     sparse = 0;
04622   }
04623   else
04624   {
04625     sparse = 1;
04626   }
04627 
04628   if (sparse)
04629   {
04630     ILLfactor_btrane2 (f, &f->xtmp);
04631     if (f->xtmp.nzcnt >= SPARSE_FACTOR * f->dim)
04632     {
04633       aindx = f->xtmp.indx;
04634       acoef = f->xtmp.coef;
04635       work_coef = f->work_coef;
04636       nzcnt = f->xtmp.nzcnt;
04637       for (i = 0; i < nzcnt; i++)
04638       {
04639         EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
04640       }
04641       sparse = 0;
04642     }
04643   }
04644   else
04645   {
04646     ILLfactor_btrane (f, work_coef);
04647   }
04648 
04649   if (sparse)
04650   {
04651     ILLfactor_btranl3 (f, &f->xtmp, x);
04652   }
04653   else
04654   {
04655     ILLfactor_btranl2 (f, work_coef);
04656     dim = f->dim;
04657     nzcnt = 0;
04658     aindx = x->indx;
04659     acoef = x->coef;
04660     for (i = 0; i < dim; i++)
04661     {
04662       if (EGlpNumIsNeqqZero (work_coef[i]))
04663       {
04664         if (EGlpNumIsNeqZero (work_coef[i], f->szero_tol))
04665           /*if (work_coef[i] > szero_tol || work_coef[i] < -szero_tol) */
04666         {
04667           aindx[nzcnt] = i;
04668           EGlpNumCopy (acoef[nzcnt], work_coef[i]);
04669           nzcnt++;
04670         }
04671         EGlpNumZero (work_coef[i]);
04672       }
04673     }
04674     x->nzcnt = nzcnt;
04675   }
04676 
04677 #ifdef SORT_RESULTS
04678   sort_vector (x);
04679 #endif
04680 
04681 #ifdef DEBUG_FACTOR
04682   {
04683     printf ("ILLfactor_btran x:");
04684     for (i = 0; i < x->nzcnt; i++)
04685     {
04686       printf (" %d %.3f", x->indx[i], EGlpNumToLf (x->coef[i]));
04687     }
04688     printf ("\n");
04689     fflush (stdout);
04690   }
04691 #endif /* DEBUG_FACTOR */
04692   return;
04693 }
04694 
04695 static int expand_col (
04696   factor_work * f,
04697   int col)
04698 {
04699   uc_info *uc_inf = f->uc_inf + col;
04700   int uc_freebeg = f->uc_freebeg;
04701   int nzcnt = uc_inf->nzcnt;
04702   int cbeg;
04703   EGlpNum_t *uccoef;
04704   int *ucindx;
04705   int *ucrind;
04706   int i;
04707   int rval = 0;
04708 
04709   if (uc_freebeg + nzcnt + 1 >= f->uc_space)
04710   {
04711     rval = make_uc_space (f, nzcnt + 1);
04712     CHECKRVALG (rval, CLEANUP);
04713     uc_freebeg = f->uc_freebeg;
04714   }
04715   cbeg = uc_inf->cbeg;
04716   uccoef = f->uccoef;
04717   ucindx = f->ucindx;
04718   ucrind = f->ucrind;
04719 
04720   for (i = 0; i < nzcnt; i++)
04721   {
04722     EGlpNumCopy (uccoef[uc_freebeg + i], uccoef[cbeg + i]);
04723     ucindx[uc_freebeg + i] = ucindx[cbeg + i];
04724     ucrind[uc_freebeg + i] = ucrind[cbeg + i];
04725     ucindx[cbeg + i] = -1;
04726   }
04727 
04728   uc_inf->cbeg = uc_freebeg;
04729   f->uc_freebeg = uc_freebeg + nzcnt;
04730 CLEANUP:
04731   EG_RETURN (rval);
04732 }
04733 
04734 static int expand_row (
04735   factor_work * f,
04736   int row)
04737 {
04738   ur_info *ur_inf = f->ur_inf + row;
04739   int ur_freebeg = f->ur_freebeg;
04740   int nzcnt = ur_inf->nzcnt;
04741   int rbeg;
04742   EGlpNum_t *urcoef;
04743   int *urindx;
04744   int *urcind;
04745   int i;
04746   int rval = 0;
04747 
04748   if (ur_freebeg + nzcnt + 1 >= f->ur_space)
04749   {
04750     rval = make_ur_space (f, nzcnt + 1);
04751     CHECKRVALG (rval, CLEANUP);
04752     ur_freebeg = f->ur_freebeg;
04753   }
04754   rbeg = ur_inf->rbeg;
04755   urcoef = f->urcoef;
04756   urindx = f->urindx;
04757   urcind = f->urcind;
04758 
04759   for (i = 0; i < nzcnt; i++)
04760   {
04761     EGlpNumCopy (urcoef[ur_freebeg + i], urcoef[rbeg + i]);
04762     urindx[ur_freebeg + i] = urindx[rbeg + i];
04763     urcind[ur_freebeg + i] = urcind[rbeg + i];
04764     urindx[rbeg + i] = -1;
04765   }
04766 
04767   ur_inf->rbeg = ur_freebeg;
04768   f->ur_freebeg = ur_freebeg + nzcnt;
04769 CLEANUP:
04770   EG_RETURN (rval);
04771 }
04772 
04773 static int add_nonzero (
04774   factor_work * f,
04775   int row,
04776   int col,
04777   EGlpNum_t val)
04778 {
04779   ur_info *ur_inf = f->ur_inf + row;
04780   uc_info *uc_inf = f->uc_inf + col;
04781   int cnzcnt = uc_inf->nzcnt;
04782   int rnzcnt = ur_inf->nzcnt;
04783   int cloc = uc_inf->cbeg + cnzcnt;
04784   int rloc = ur_inf->rbeg + rnzcnt;
04785   int rval = 0;
04786 
04787   if (f->ucindx[cloc] != -1)
04788   {
04789     rval = expand_col (f, col);
04790     CHECKRVALG (rval, CLEANUP);
04791     cloc = uc_inf->cbeg + cnzcnt;
04792   }
04793   TESTG ((rval = (rloc < 0 || rloc > f->ur_space)), CLEANUP,
04794          "rloc %d outside boundaries [0:%d]", rloc, f->ur_space);
04795   if (f->urindx[rloc] != -1)
04796   {
04797     rval = expand_row (f, row);
04798     CHECKRVALG (rval, CLEANUP);
04799     rloc = ur_inf->rbeg + rnzcnt;
04800   }
04801   f->ucindx[cloc] = row;
04802   EGlpNumCopy (f->uccoef[cloc], val);
04803   f->ucrind[cloc] = rnzcnt;
04804   f->urindx[rloc] = col;
04805   EGlpNumCopy (f->urcoef[rloc], val);
04806   f->urcind[rloc] = cnzcnt;
04807 
04808   if (cloc == f->uc_freebeg)
04809     f->uc_freebeg++;
04810   if (rloc == f->ur_freebeg)
04811     f->ur_freebeg++;
04812 
04813   uc_inf->nzcnt = cnzcnt + 1;
04814   ur_inf->nzcnt = rnzcnt + 1;
04815 CLEANUP:
04816   EG_RETURN (rval);
04817 }
04818 
04819 static int delete_nonzero_row (
04820   factor_work * f,
04821   int row,
04822   int ind)
04823 {
04824   ur_info *ur_inf = f->ur_inf;
04825   EGlpNum_t *urcoef = f->urcoef;
04826   int *urindx = f->urindx;
04827   int *urcind = f->urcind;
04828   int *ucrind = f->ucrind;
04829   int rbeg = ur_inf[row].rbeg;
04830   int nzcnt = ur_inf[row].nzcnt - 1;
04831   int cbeg, rval = 0;
04832   #ifdef DEBUG_FACTOR
04833   TESTG((rval=(nzcnt<0)),CLEANUP,"Deleting empty row %d ind %d!",row, ind);
04834   #endif
04835 
04836   if (ind != nzcnt)
04837   {
04838     EGlpNumCopy (urcoef[rbeg + ind], urcoef[rbeg + nzcnt]);
04839     urindx[rbeg + ind] = urindx[rbeg + nzcnt];
04840     urcind[rbeg + ind] = urcind[rbeg + nzcnt];
04841     cbeg = f->uc_inf[urindx[rbeg + nzcnt]].cbeg;
04842     ucrind[cbeg + urcind[rbeg + nzcnt]] = ind;
04843     urindx[rbeg + nzcnt] = -1;
04844   }
04845   ur_inf[row].nzcnt = nzcnt;
04846   #ifdef DEBUG_FACTOR
04847   CLEANUP:
04848   #endif
04849   return rval;
04850 }
04851 
04852 static void delete_nonzero_col (
04853   factor_work * f,
04854   int col,
04855   int ind)
04856 {
04857   uc_info *uc_inf = f->uc_inf;
04858   EGlpNum_t *uccoef = f->uccoef;
04859   int *ucindx = f->ucindx;
04860   int *ucrind = f->ucrind;
04861   int *urcind = f->urcind;
04862   int cbeg = uc_inf[col].cbeg;
04863   int nzcnt = uc_inf[col].nzcnt - 1;
04864   int rbeg;
04865 
04866   if (ind != nzcnt)
04867   {
04868     EGlpNumCopy (uccoef[cbeg + ind], uccoef[cbeg + nzcnt]);
04869     ucindx[cbeg + ind] = ucindx[cbeg + nzcnt];
04870     ucrind[cbeg + ind] = ucrind[cbeg + nzcnt];
04871     rbeg = f->ur_inf[ucindx[cbeg + nzcnt]].rbeg;
04872     urcind[rbeg + ucrind[cbeg + nzcnt]] = ind;
04873     ucindx[cbeg + nzcnt] = -1;
04874   }
04875   uc_inf[col].nzcnt = nzcnt;
04876 }
04877 
04878 static int delete_column (
04879   factor_work * f,
04880   int col)
04881 {
04882   uc_info *uc_inf = f->uc_inf;
04883   int beg = uc_inf[col].cbeg;
04884   int nzcnt = uc_inf[col].nzcnt;
04885   int *ucindx = f->ucindx + beg;
04886   int *ucrind = f->ucrind + beg;
04887   int i, rval = 0;
04888   #ifdef DEBUG_FACTOR
04889   rval = check_matrix(f);
04890   TESTG(rval,CLEANUP,"Corrupted Matrix");
04891   #endif
04892 
04893   for (i = 0; i < nzcnt; i++)
04894   {
04895     rval = delete_nonzero_row (f, ucindx[i], ucrind[i]);
04896     CHECKRVALG(rval,CLEANUP);
04897     ucindx[i] = -1;
04898   }
04899   uc_inf[col].nzcnt = 0;
04900 
04901 #ifdef TRACK_FACTOR
04902   f->nzcnt_cur -= nzcnt;
04903 #endif
04904 
04905   #ifdef DEBUG_FACTOR
04906   rval = check_matrix(f);
04907   TESTG(rval,CLEANUP,"Corrupted Matrix");
04908   #endif /* DEBUG_FACTOR */
04909   CLEANUP:
04910   EG_RETURN(rval);
04911 }
04912 
04913 static int delete_row (
04914   factor_work * f,
04915   int row,
04916   svector * x)
04917 {
04918   ur_info *ur_inf = f->ur_inf;
04919   int beg = ur_inf[row].rbeg;
04920   int nzcnt = ur_inf[row].nzcnt;
04921   int *urindx = f->urindx + beg;
04922   EGlpNum_t *urcoef = f->urcoef + beg;
04923   int *urcind = f->urcind + beg;
04924   int i,rval=0;
04925 
04926   #ifdef DEBUG_FACTOR
04927   rval = check_matrix(f);
04928   TESTG(rval,CLEANUP,"Corrupted Matrix");
04929   #endif /* DEBUG_FACTOR */
04930 
04931   for (i = 0; i < nzcnt; i++)
04932   {
04933     x->indx[i] = urindx[i];
04934     EGlpNumCopy (x->coef[i], urcoef[i]);
04935     delete_nonzero_col (f, urindx[i], urcind[i]);
04936     urindx[i] = -1;
04937   }
04938   x->nzcnt = nzcnt;
04939   ur_inf[row].nzcnt = 0;
04940 
04941 #ifdef TRACK_FACTOR
04942   f->nzcnt_cur -= nzcnt;
04943 #endif
04944 
04945   #ifdef DEBUG_FACTOR
04946   rval = check_matrix(f);
04947   TESTG(rval,CLEANUP,"Corrupted Matrix");
04948   CLEANUP:
04949   #endif /* DEBUG_FACTOR */
04950   return rval;
04951 }
04952 
04953 static int create_column (
04954   factor_work * f,
04955   svector * a,
04956   int col,
04957   int *p_last_rank)
04958 {
04959   int *rrank = f->rrank;
04960   int nzcnt = a->nzcnt;
04961   int *aindx = a->indx;
04962   EGlpNum_t *acoef = a->coef;
04963   int i;
04964   int j;
04965   int rval = 0;
04966   int last_rank = -1;
04967 
04968 #ifdef TRACK_FACTOR
04969   EGlpNum_t max;
04970 
04971   EGlpNumInitVar (max);
04972   EGlpNumCopy (max, f->maxelem_cur);
04973 #endif /* TRACK_FACTOR */
04974 
04975   last_rank = 0;
04976 
04977   for (i = 0; i < nzcnt; i++)
04978   {
04979     rval = add_nonzero (f, aindx[i], col, acoef[i]);
04980     CHECKRVALG (rval, CLEANUP);
04981 #ifdef TRACK_FACTOR
04982     EGlpNumSetToMaxAbs (max, acoef[i]);
04983 #endif /* TRACK_FACTOR */
04984     j = rrank[aindx[i]];
04985     if (j > last_rank)
04986       last_rank = j;
04987   }
04988   *p_last_rank = last_rank;
04989 
04990 #ifdef TRACK_FACTOR
04991   f->nzcnt_cur += nzcnt;
04992   EGlpNumCopy (f->maxelem_cur, max);
04993   EGlpNumClearVar (max);
04994 #endif /* TRACK_FACTOR */
04995 
04996   #ifdef DEBUG_FACTOR
04997   rval = check_matrix(f);
04998   TESTG(rval,CLEANUP,"Corrupted Matrix");
04999   #endif /* DEBUG_FACTOR */
05000 
05001   CLEANUP:
05002   EG_RETURN (rval);
05003 }
05004 
05005 #ifdef UPDATE_STUDY
05006 static int column_rank (
05007   factor_work * f,
05008   int col)
05009 {
05010   int *cperm = f->cperm;
05011   int dim = f->dim;
05012   int i;
05013 
05014   for (i = 0; i < dim; i++)
05015   {
05016     if (cperm[i] == col)
05017     {
05018       return i;
05019     }
05020   }
05021   return 0;
05022 }
05023 #endif
05024 
05025 static void shift_permutations (
05026   factor_work * f,
05027   int rank_p,
05028   int rank_r)
05029 {
05030   int *cperm = f->cperm;
05031   int *crank = f->crank;
05032   int *rperm = f->rperm;
05033   int *rrank = f->rrank;
05034   int col_p = cperm[rank_p];
05035   int row_p = rperm[rank_p];
05036   int i;
05037 
05038   for (i = rank_p; i < rank_r; i++)
05039   {
05040     cperm[i] = cperm[i + 1];
05041     crank[cperm[i]] = i;
05042     rperm[i] = rperm[i + 1];
05043     rrank[rperm[i]] = i;
05044   }
05045   cperm[rank_r] = col_p;
05046   crank[col_p] = rank_r;
05047   rperm[rank_r] = row_p;
05048   rrank[row_p] = rank_r;
05049 }
05050 
05051 static int eliminate_row (
05052   factor_work * f,
05053   int rank_p,
05054   int rank_r)
05055 {
05056   ur_info *ur_inf = f->ur_inf;
05057   int *rperm = f->rperm;
05058   int *cperm = f->cperm;
05059   int *urindx = f->urindx;
05060   EGlpNum_t *urcoef = f->urcoef;
05061   int *erindx = f->erindx;
05062   EGlpNum_t *ercoef = f->ercoef;
05063   EGlpNum_t *work_coef = f->work_coef;
05064   int er_freebeg = f->er_freebeg;
05065   int er_space = f->er_space;
05066   int beg;
05067   int nzcnt;
05068   int i;
05069   int j;
05070   int c;
05071   int r;
05072   EGlpNum_t pivot_mul;
05073 
05074 #ifdef TRACK_FACTOR
05075   EGlpNum_t max;
05076 
05077   EGlpNumInitVar (max);
05078   EGlpNumCopy (max, f->maxelem_cur);
05079 #endif /* TRACK_FACTOR */
05080   EGlpNumInitVar (pivot_mul);
05081 
05082   for (i = rank_p; i < rank_r; i++)
05083   {
05084     c = cperm[i];
05085     if (EGlpNumIsNeqZero (work_coef[c], f->fzero_tol))  /*
05086                                                          * if (work_coef[c] > fzero_tol || work_coef[c] < -fzero_tol) */
05087     {
05088       r = rperm[i];
05089       beg = ur_inf[r].rbeg;
05090       nzcnt = ur_inf[r].nzcnt;
05091       EGlpNumCopyFrac (pivot_mul, work_coef[c], urcoef[beg]);
05092       EGlpNumZero (work_coef[c]);
05093       for (j = 1; j < nzcnt; j++)
05094       {
05095         EGlpNumSubInnProdTo (work_coef[urindx[beg + j]], pivot_mul, urcoef[beg + j]); /* 0.85 */
05096       }
05097       if (er_freebeg >= er_space)
05098       {
05099         /* fprintf (stderr, "no space in eliminate_row\n"); */
05100 #ifdef TRACK_FACTOR
05101         EGlpNumClearVar (max);
05102 #endif
05103         EGlpNumClearVar (pivot_mul);
05104         return E_UPDATE_NOSPACE;
05105       }
05106       erindx[er_freebeg] = r;
05107       EGlpNumCopy (ercoef[er_freebeg], pivot_mul);
05108 #ifdef TRACK_FACTOR
05109       EGlpNumSetToMaxAbs (max, pivot_mul);
05110 #endif /* TRACK_FACTOR */
05111       er_freebeg++;
05112     }
05113     else
05114     {
05115       EGlpNumZero (work_coef[c]);
05116     }
05117   }
05118   f->er_freebeg = er_freebeg;
05119 #ifdef TRACK_FACTOR
05120   EGlpNumCopy (f->maxelem_cur, max);
05121   EGlpNumClearVar (max);
05122 #endif /* TRACK_FACTOR */
05123   EGlpNumClearVar (pivot_mul);
05124   return 0;
05125 }
05126 
05127 static int create_row (
05128   factor_work * f,
05129   EGlpNum_t * a,
05130   int row,
05131   int minrank)
05132 {
05133   int *cperm = f->cperm;
05134   int dim = f->dim;
05135   int i;
05136   int j;
05137   int rval = 0;
05138 
05139 #ifdef TRACK_FACTOR
05140   EGlpNum_t max;
05141 
05142   EGlpNumInitVar (max);
05143   EGlpNumCopy (max, f->maxelem_cur);
05144 #endif /* TRACK_FACTOR */
05145 
05146   for (i = minrank; i < dim; i++)
05147   {
05148     if (EGlpNumIsNeqqZero (a[cperm[i]]))
05149     {
05150       j = cperm[i];
05151       if (EGlpNumIsNeqZero (a[j], f->fzero_tol))  /*
05152                                                    * if (a[j] > fzero_tol || a[j] < -fzero_tol) */
05153       {
05154         rval = add_nonzero (f, row, j, a[j]);
05155         CHECKRVALG (rval, CLEANUP);
05156 #ifdef TRACK_FACTOR
05157         EGlpNumSetToMaxAbs (max, a[j]);
05158 #endif /* TRACK_FACTOR */
05159       }
05160       EGlpNumZero (a[j]);
05161     }
05162   }
05163 
05164 #ifdef TRACK_FACTOR
05165   f->nzcnt_cur += f->ur_inf[row].nzcnt;
05166   EGlpNumCopy (f->maxelem_cur, max);
05167   EGlpNumClearVar (max);
05168 #endif /* TRACK_FACTOR */
05169 
05170   #ifdef DEBUG_FACTOR
05171   rval = check_matrix(f);
05172   TESTG(rval,CLEANUP,"Corrupted Matrix");
05173   #endif /* DEBUG_FACTOR */
05174   CLEANUP:
05175   EG_RETURN (rval);
05176 }
05177 
05178 static void serow_delay (
05179   factor_work * f,
05180   int r,
05181   int rank_r)
05182 {
05183   ur_info *ur_inf = f->ur_inf;
05184   int *crank = f->crank;
05185   int nzcnt;
05186   int *indx;
05187   int i;
05188   int last;
05189 
05190   do
05191   {
05192     r = f->rperm[crank[r]];
05193     nzcnt = ur_inf[r].nzcnt;
05194     indx = f->urindx + ur_inf[r].rbeg;
05195     last = -1;
05196     for (i = 1; i < nzcnt; i++)
05197     {
05198       r = indx[i];
05199       if (ur_inf[r].delay++ == 0 && crank[r] < rank_r)
05200       {
05201         if (last >= 0)
05202         {
05203           serow_delay (f, last, rank_r);
05204         }
05205         last = r;
05206       }
05207     }
05208     r = last;
05209   } while (r >= 0);
05210 }
05211 
05212 static int serow_process (
05213   factor_work * f,
05214   int r,
05215   svector * newr,
05216   int rank_r)
05217 {
05218   ur_info *ur_inf = f->ur_inf;
05219   EGlpNum_t *work = f->work_coef;
05220   int nzcnt;
05221   int *indx;
05222   EGlpNum_t *coef;
05223   int i;
05224   EGlpNum_t v;
05225   int last;
05226   int rval;
05227 
05228   EGlpNumInitVar (v);
05229 
05230   do
05231   {
05232     EGlpNumCopy (v, work[r]);
05233     EGlpNumZero (work[r]);
05234     if (f->crank[r] >= rank_r)
05235     {
05236       if (EGlpNumIsNeqZero (v, f->fzero_tol)) /*
05237                                                * if (v > fzero_tol || v < -fzero_tol) */
05238       {
05239         /* stash this nonzero in the resulting row */
05240 #ifdef TRACK_FACTOR
05241         EGlpNumSetToMaxAbs (f->maxelem_cur, v);
05242 #endif /* TRACK_FACTOR */
05243         newr->indx[newr->nzcnt] = r;
05244         EGlpNumCopy (newr->coef[newr->nzcnt], v);
05245         newr->nzcnt++;
05246         EGlpNumClearVar (v);
05247         return 0;
05248       }
05249       else
05250       {
05251         EGlpNumClearVar (v);
05252         return 0;
05253       }
05254     }
05255     r = f->rperm[f->crank[r]];
05256     nzcnt = ur_inf[r].nzcnt;
05257     indx = f->urindx + ur_inf[r].rbeg;
05258     coef = f->urcoef + ur_inf[r].rbeg;
05259     EGlpNumDivTo (v, coef[0]);
05260     if (EGlpNumIsNeqZero (v, f->fzero_tol)) /*
05261                                              * if (v > fzero_tol || v < -fzero_tol) */
05262     {
05263       /* stash v in eta */
05264       if (f->er_freebeg >= f->er_space)
05265       {
05266         /* fprintf (stderr, "no space in eliminate_row\n"); */
05267         EGlpNumClearVar (v);
05268         return E_UPDATE_NOSPACE;
05269       }
05270       f->erindx[f->er_freebeg] = r;
05271       EGlpNumCopy (f->ercoef[f->er_freebeg], v);
05272 #ifdef TRACK_FACTOR
05273       EGlpNumSetToMaxAbs (f->maxelem_cur, v);
05274 #endif /* TRACK_FACTOR */
05275       f->er_freebeg++;
05276     }
05277     last = -1;
05278     for (i = 1; i < nzcnt; i++)
05279     {
05280       r = indx[i];
05281       EGlpNumSubInnProdTo (work[r], v, coef[i]);
05282       if (--ur_inf[r].delay == 0)
05283       {
05284         if (last >= 0)
05285         {
05286           rval = serow_process (f, last, newr, rank_r);
05287           if (rval)
05288           {
05289             EGlpNumClearVar (v);
05290             return rval;
05291           }
05292         }
05293         last = r;
05294       }
05295     }
05296     r = last;
05297   } while (r >= 0);
05298   EGlpNumClearVar (v);
05299   return 0;
05300 }
05301 
05302 static int sparse_eliminate_row (
05303   factor_work * f,
05304   svector * x,
05305   int row_p,
05306   int rank_r)
05307 {
05308   EGlpNum_t *work = f->work_coef;
05309   int xnzcnt = x->nzcnt;
05310   int *xindx = x->indx;
05311   EGlpNum_t *xcoef = x->coef;
05312   ur_info *ur_inf = f->ur_inf;
05313   int *crank = f->crank;
05314   int i;
05315   int j;
05316   int rval = 0;
05317   svector newr;
05318 
05319   newr.indx = 0;
05320   newr.coef = 0;
05321 
05322   for (i = 0; i < xnzcnt; i++)
05323   {
05324     j = xindx[i];
05325     if (ur_inf[j].delay++ == 0 && crank[j] < rank_r)
05326     {
05327       serow_delay (f, j, rank_r);
05328     }
05329     EGlpNumCopy (work[j], xcoef[i]);
05330   }
05331 
05332   newr.nzcnt = 0;
05333   ILL_SAFE_MALLOC (newr.indx, f->dim, int);
05334 
05335   newr.coef = EGlpNumAllocArray (f->dim);
05336 
05337   for (i = 0; i < xnzcnt; i++)
05338   {
05339     j = xindx[i];
05340     if (--ur_inf[j].delay == 0)
05341     {
05342       rval = serow_process (f, j, &newr, rank_r);
05343       CHECKRVALG (rval, CLEANUP);
05344     }
05345   }
05346 
05347   for (i = 0; i < newr.nzcnt; i++)
05348   {
05349     rval = add_nonzero (f, row_p, newr.indx[i], newr.coef[i]);
05350     CHECKRVALG (rval, CLEANUP);
05351   }
05352 
05353 #ifdef TRACK_FACTOR
05354   f->nzcnt_cur += newr.nzcnt;
05355 #endif /* TRACK_FACTOR */
05356 
05357 CLEANUP:
05358   EGlpNumFreeArray (newr.coef);
05359   ILL_IFFREE (newr.indx, int);
05360 
05361   /* Bico 031210 - chg from ILL_RETURN */
05362   EG_RETURN (rval);
05363 }
05364 
05365 static int move_pivot_row (
05366   factor_work * f,
05367   int r,
05368   int c)
05369 {
05370   ur_info *ur_inf = f->ur_inf + r;
05371   uc_info *uc_inf = f->uc_inf;
05372   int beg = ur_inf->rbeg;
05373   int nzcnt = ur_inf->nzcnt;
05374   int *urindx = f->urindx;
05375   int *urcind = f->urcind;
05376   int *ucrind = f->ucrind;
05377   EGlpNum_t *urcoef = f->urcoef;
05378   EGlpNum_t dt;
05379   int it;
05380   int i;
05381 
05382   if (urindx[beg] == c)
05383     return 0;
05384   EGlpNumInitVar (dt);
05385 
05386   for (i = 1; i < nzcnt; i++)
05387   {
05388     if (urindx[beg + i] == c)
05389     {
05390       EGLPNUM_SWAP (urcoef[beg], urcoef[beg + i], dt);
05391       ILL_SWAP (urcind[beg], urcind[beg + i], it);
05392       urindx[beg + i] = urindx[beg];
05393       urindx[beg] = c;
05394       ucrind[uc_inf[c].cbeg + urcind[beg]] = 0;
05395       ucrind[uc_inf[urindx[beg + i]].cbeg + urcind[beg + i]] = i;
05396       EGlpNumClearVar (dt);
05397       return 0;
05398     }
05399   }
05400   MESSAGE (__QS_SB_VERB, "pivot row nonzero not found");
05401   EGlpNumClearVar (dt);
05402   return E_UPDATE_SINGULAR_ROW;
05403 }
05404 
05405 static int move_pivot_col (
05406   factor_work * f,
05407   int c,
05408   int r)
05409 {
05410   uc_info *uc_inf = f->uc_inf + c;
05411   ur_info *ur_inf = f->ur_inf;
05412   int beg = uc_inf->cbeg;
05413   int nzcnt = uc_inf->nzcnt;
05414   int *ucindx = f->ucindx;
05415   int *ucrind = f->ucrind;
05416   int *urcind = f->urcind;
05417   EGlpNum_t *uccoef = f->uccoef;
05418   EGlpNum_t dt;
05419   int i, it;
05420 
05421   if (ucindx[beg] == r)
05422     return 0;
05423   EGlpNumInitVar (dt);
05424 
05425   for (i = 1; i < nzcnt; i++)
05426   {
05427     if (ucindx[beg + i] == r)
05428     {
05429       EGLPNUM_SWAP (uccoef[beg], uccoef[beg + i], dt);
05430       ILL_SWAP (ucrind[beg], ucrind[beg + i], it);
05431       ucindx[beg + i] = ucindx[beg];
05432       ucindx[beg] = r;
05433       urcind[ur_inf[r].rbeg + ucrind[beg]] = 0;
05434       urcind[ur_inf[ucindx[beg + i]].rbeg + ucrind[beg + i]] = i;
05435       EGlpNumClearVar (dt);
05436       return 0;
05437     }
05438   }
05439   MESSAGE(__QS_SB_VERB, "pivot col nonzero not found");
05440   EGlpNumClearVar (dt);
05441   return E_UPDATE_SINGULAR_COL;
05442 }
05443 
05444 static int move_pivot (
05445   factor_work * f,
05446   int rank_r)
05447 {
05448   int r = f->rperm[rank_r];
05449   int c = f->cperm[rank_r];
05450   int rval = 0;
05451 
05452   rval = move_pivot_row (f, r, c);
05453   CHECKRVALG (rval, CLEANUP);
05454 
05455   #ifdef DEBUG_FACTOR
05456   rval = check_matrix(f);
05457   TESTG(rval,CLEANUP,"Corrupted Matrix");
05458   #endif /* DEBUG_FACTOR */
05459 
05460   rval = move_pivot_col (f, c, r);
05461   if(rval != E_UPDATE_SINGULAR_COL) CHECKRVALG (rval, CLEANUP);
05462   else goto CLEANUP;
05463 
05464   #ifdef DEBUG_FACTOR
05465   rval = check_matrix(f);
05466   TESTG(rval,CLEANUP,"Corrupted Matrix");
05467   #endif /* DEBUG_FACTOR */
05468 
05469   CLEANUP:
05470   if(rval != E_UPDATE_SINGULAR_COL) EG_RETURN (rval);             /* Bico 031209 - chg from RETURN */
05471   return rval;
05472 }
05473 
05474 int ILLfactor_update (
05475   factor_work * f,
05476   svector * a,
05477   int col_p,
05478   int *p_refact)
05479 {
05480   int row_p;
05481   int rank_r = 0;
05482   int rank_p = 0;
05483   int rval = 0;
05484   int nzcnt;
05485   int *aindx;
05486   EGlpNum_t *acoef;
05487   EGlpNum_t *work_coef = f->work_coef;
05488 
05489 #ifdef TRACK_FACTOR
05490 #ifdef NOTICE_BLOWUP
05491   EGlpNum_t tmpsize;
05492 #endif
05493 #endif
05494   int i;
05495 
05496 #ifdef RECORD
05497   {
05498     EGioPrintf (fsave, "u %d %d", col_p, a->nzcnt);
05499     for (i = 0; i < a->nzcnt; i++)
05500     {
05501       EGioPrintf (fsave, " %d %.16e", a->indx[i], EGlpNumToLf (a->coef[i]));
05502     }
05503     EGioPrintf (fsave, "\n");
05504     EGioFlush (fsave);
05505   }
05506 #endif /* RECORD */
05507 
05508 #ifdef DEBUG_FACTOR
05509   {
05510     printf ("ILLfactor_update col %d:", col_p);
05511     for (i = 0; i < a->nzcnt; i++)
05512     {
05513       printf (" %.3f*%d", EGlpNumToLf (a->coef[i]), a->indx[i]);
05514     }
05515     printf ("\n");
05516     fflush (stdout);
05517   }
05518 #endif /* DEBUG_FACTOR */
05519 
05520   #ifdef DEBUG_FACTOR
05521   rval = check_matrix(f);
05522   TESTG(rval,CLEANUP,"Corrupted Matrix");
05523   #endif /* DEBUG_FACTOR */
05524 
05525   if (f->etacnt >= f->etamax)
05526   {
05527     *p_refact = 1;
05528     return 0;
05529   }
05530 
05531 #ifdef UPDATE_STUDY
05532   nupdate++;
05533 #endif
05534 
05535   row_p = f->ucindx[f->uc_inf[col_p].cbeg];
05536 
05537   rval = delete_column (f, col_p);
05538   CHECKRVALG (rval, CLEANUP);
05539 
05540   rval = create_column (f, a, col_p, &rank_r);
05541   /* if (rval) fprintf (stderr, "create_column failed\n"); */
05542   CHECKRVALG (rval, CLEANUP);
05543 
05544   rank_p = f->crank[col_p];
05545 #ifdef UPDATE_STUDY
05546   if (rank_p != f->rrank[row_p] || rank_p != column_rank (f, col_p))
05547   {
05548     printf ("rank_p %d rrank[row_p] %d column_rank(f,col_p) %d\n",
05549             rank_p, f->rrank[row_p], column_rank (f, col_p));
05550   }
05551   if (rank_r > rank_p)
05552   {
05553     permshifttot += rank_r - rank_p;
05554   }
05555   for (i = 0; i < a->nzcnt; i++)
05556   {
05557     if (f->rrank[a->indx[i]] > rank_p)
05558       colspiketot++;
05559   }
05560   for (i = 0; i < f->ur_inf[row_p].nzcnt; i++)
05561   {
05562     if (f->crank[f->urindx[f->ur_inf[row_p].rbeg + i]] <= rank_r &&
05563         f->crank[f->urindx[f->ur_inf[row_p].rbeg + i]] != rank_p)
05564     {
05565       rowspiketot++;
05566     }
05567   }
05568 #endif
05569 
05570   shift_permutations (f, rank_p, rank_r);
05571 
05572   rval = delete_row (f, row_p, &f->xtmp);
05573   CHECKRVALG(rval,CLEANUP);
05574 
05575   f->er_inf[f->etacnt].rbeg = f->er_freebeg;
05576   f->er_inf[f->etacnt].r = row_p;
05577 
05578   if (f->xtmp.nzcnt >= SPARSE_FACTOR * f->dim)
05579   {
05580     nzcnt = f->xtmp.nzcnt;
05581     aindx = f->xtmp.indx;
05582     acoef = f->xtmp.coef;
05583 
05584     for (i = 0; i < nzcnt; i++)
05585     {
05586       EGlpNumCopy (work_coef[aindx[i]], acoef[i]);
05587     }
05588 
05589     rval = eliminate_row (f, rank_p, rank_r);
05590     /* if (rval) fprintf (stderr, "eliminate_row failed\n"); */
05591     CHECKRVALG (rval, CLEANUP);
05592 
05593     rval = create_row (f, f->work_coef, row_p, rank_r);
05594     /* if (rval) fprintf (stderr, "create_row failed\n"); */
05595     CHECKRVALG (rval, CLEANUP);
05596   }
05597   else
05598   {
05599     rval = sparse_eliminate_row (f, &f->xtmp, row_p, rank_r);
05600     /* if (rval) fprintf (stderr, "sparse_eliminate_row failed\n"); */
05601     CHECKRVALG (rval, CLEANUP);
05602   }
05603 
05604   if (f->er_freebeg - f->er_inf[f->etacnt].rbeg > 0)
05605   {
05606     f->er_inf[f->etacnt].nzcnt = f->er_freebeg - f->er_inf[f->etacnt].rbeg;
05607 #ifdef TRACK_FACTOR
05608     f->nzcnt_cur += f->er_inf[f->etacnt].nzcnt;
05609 #endif /* TRACK_FACTOR */
05610 #ifdef UPDATE_STUDY
05611     leftetatot += f->er_inf[f->etacnt].nzcnt;
05612 #endif
05613 
05614 #ifdef SORT_RESULTS
05615     sort_vector2 (f->er_inf[f->etacnt].nzcnt,
05616                   f->erindx + f->er_inf[f->etacnt].rbeg,
05617                   f->ercoef + f->er_inf[f->etacnt].rbeg);
05618 #endif
05619 
05620     f->etacnt++;
05621   }
05622 
05623   rval = move_pivot (f, rank_r);
05624   /* if (rval) fprintf (stderr, "move_pivot failed\n"); */
05625   if(rval != E_UPDATE_SINGULAR_COL) CHECKRVALG (rval, CLEANUP);
05626   else goto CLEANUP;
05627 
05628 #ifdef UPDATE_DEBUG
05629   printf ("Updated factorization:\n");
05630 #if (UPDATE_DEBUG+0>1)
05631   dump_matrix (f, 0);
05632 #endif
05633   fflush (stdout);
05634 #endif /* UPDATE_DEBUG */
05635 
05636 #ifdef TRACK_FACTOR
05637 #ifdef NOTICE_BLOWUP
05638   EGlpNumInitVar (tmpsize);
05639   EGlpNumSet (tmpsize, f->updmaxmult);
05640   EGlpNumMultTo (tmpsize, f->maxelem_orig);
05641   if (EGlpNumIsLess (tmpsize, f->maxelem_cur))
05642   {
05643 /* Bico - comment out for dist 
05644         fprintf (stderr, "factor_update blowup max cur %e max orig %e\n",
05645                  f->maxelem_cur, f->maxelem_orig);
05646 */
05647     EGlpNumClearVar (tmpsize);
05648     return E_FACTOR_BLOWUP;
05649   }
05650   EGlpNumClearVar (tmpsize);
05651 #endif /* NOTICE_BLOWUP */
05652 #endif
05653 #ifdef UPDATE_STATS
05654   dump_factor_stats (f);
05655 #endif
05656 CLEANUP:
05657   if(rval != E_UPDATE_SINGULAR_COL) EG_RETURN (rval);             /* Bico 031209 - chg from RETURN */
05658   return rval;
05659 }

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