Logo Search packages:      
Sourcecode: rsymphony version File versions  Download package

CoinOslFactorization3.cpp

/* $Id: CoinOslFactorization3.cpp 1202 2009-09-25 15:27:15Z forrest $ */
/* Copyright (C) 1987, 2009, International Business Machines
   Corporation and others.  All Rights Reserved. */
#include "CoinOslC.h"
#define GO_DENSE 70
#define GO_DENSE_RATIO 1.8
int c_ekkclco(const EKKfactinfo *fact,int *hcoli,
          int *mrstrt, int *hinrow, int xnewro);

void c_ekkclcp(const int *hcol, const double *dels, const int * mrstrt,
           int *hrow, double *dels2, int *mcstrt,
           int *hincol, int itype, int nnrow, int nncol,
           int ninbas);

int c_ekkcmfc(EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          EKKHlink *mwork, void *maction_void,
          int nnetas,
          int *nsingp, int *xrejctp,
          int *xnewrop, int xnewco, 
          int *ncompactionsp);

int c_ekkcmfy(EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          EKKHlink *mwork, void *maction_void,
          int nnetas,
          int *nsingp, int *xrejctp,
          int *xnewrop, int xnewco, 
          int *ncompactionsp);

int c_ekkcmfd(EKKfactinfo *fact,
          int *mcol, 
          EKKHlink *rlink, EKKHlink *clink,
          int *maction,
          int nnetas,
          int *nnentlp, int *nnentup,
          int *nsingp);
int c_ekkford(const EKKfactinfo *fact,const int *hinrow, const int *hincol,
          int *hpivro, int *hpivco,
          EKKHlink *rlink, EKKHlink *clink);
void c_ekkrowq(int *hrow, int *hcol, double *dels, 
           int *mrstrt,
           const int *hinrow, int nnrow, int ninbas);
int c_ekkrwco(const EKKfactinfo *fact,double *dluval, int *hcoli, int *
          mrstrt, int *hinrow, int xnewro);

int c_ekkrwcs(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
          const int *hinrow, const EKKHlink *mwork,
          int nfirst);

void c_ekkrwct(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
           const int *hinrow, const EKKHlink *mwork,
           const EKKHlink *rlink,
           const short *msort, double *dsort,
           int nlast, int xnewro);

int c_ekkshff(EKKfactinfo *fact,
          EKKHlink *clink, EKKHlink *rlink, 
          int xnewro);

void c_ekkshfv(EKKfactinfo *fact, EKKHlink *rlink, EKKHlink *clink,
           int xnewro);
int c_ekktria(EKKfactinfo *fact,
            EKKHlink * rlink,
            EKKHlink * clink,
          int *nsingp, 
          int *xnewcop, int *xnewrop,
          int *nlrowtp,
          const int ninbas);
#if 0
static void c_ekkafpv(int *hentry, int *hcoli,
           double *dluval, int *mrstrt,
           int *hinrow, int nentry)
{
  int j;
  int nel, krs;
  int koff;
  int irow;
  int ientry;
  int * index;
  
  for (ientry = 0; ientry < nentry; ++ientry) {
#ifdef INTEL
    int * els_long,maxaij_long;
#endif
    double * els;
    irow = UNSHIFT_INDEX(hentry[ientry]);
    nel = hinrow[irow];
    krs = mrstrt[irow];
    index=&hcoli[krs];
    els=&dluval[krs];
#ifdef INTEL
    els_long=reinterpret_cast<int *> (els);
    maxaij_long=0;
#else
    double maxaij = 0.f;
#endif
    koff = 0;
    j=0;
    if ((nel&1)!=0) {
#ifdef INTEL
      maxaij_long = els_long[1] & 0x7fffffff;
#else
      maxaij=fabs(els[0]);
#endif
      j=1;
    }
    
    while (j<nel) {
#ifdef INTEL
      UNROLL_LOOP_BODY2({
        int d_long = els_long[1+(j<<1)] & 0x7fffffff;
        if (maxaij_long < d_long) {
          maxaij_long = d_long;
          koff=j;
        }
        j++;
      });
#else
      UNROLL_LOOP_BODY2({
        double d = fabs(els[j]);
        if (maxaij < d) {
          maxaij = d;
          koff=j;
        }
        j++;
      });
#endif
    }
    
    SWAP(int, index[koff], index[0]);
    SWAP(double, els[koff], els[0]);
  }
} /* c_ekkafpv */
#endif

/*     Uwe H. Suhl, March 1987 */
/*     This routine processes col singletons during the LU-factorization. */
/*     Return codes (checked version 1.11): */
/*    0: ok */
/*      6: pivot element too small */
/*     -43: system error at label 420 (ipivot not found) */
/*
 * This routine processes singleton columns during factorization of the
 * nucleus.  It is very similar to the first part of c_ekktria,
 * but is more complex, because we now have to maintain the length
 * lists.
 * The differences are:
 * (1) here we use the length list for length 1 rather than a queue.
 * This routine is only called if it is known that there is a singleton
 * column.
 *
 * (2) here we maintain hrowi by moving the last entry into the pivot
 * column entry; that means we don't have to search for the pivot row
 * entry like we do in c_ekktria.
 *
 * (3) here the hlink data structure is in use for the length lists,
 * so we maintain it as we shorten rows and removing columns altogether.
 *
 */
int c_ekkcsin(EKKfactinfo *fact,
            EKKHlink *rlink, EKKHlink *clink,
            
            int *nsingp)
{
#if 1
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  //double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  const int nrow  = fact->nrow;
  const double drtpiv   = fact->drtpiv;
  
  
  int j, k, kc, kce, kcs, nzj;
  double pivot;
  int kipis, kipie;
  int jpivot;
#ifndef NDEBUG
  int kpivot=-1;
#else
  int kpivot=-1;
#endif
  
  bool small_pivot = false;
  
  
  /* next singleton column.
   * Note that when the pivot column itself was removed from the
   * list, the column in the list after it (if any) moves to the
   * head of the list.
   * Also, if any column from the pivot row was reduced to length 1,
   * then it will have been added to the list and now be in front.
   */
  for (jpivot = hpivco[1]; jpivot > 0; jpivot = hpivco[1]) {
    const int ipivot = hrowi[mcstrt[jpivot]]; /* (2) */
    assert(ipivot);
    /* The pivot row is being eliminated (3) */
    C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, ipivot);
    
    /* Loop over nonzeros in pivot row: */
    kipis = mrstrt[ipivot];
    kipie = kipis + hinrow[ipivot] - 1;
    for (k = kipis; k <= kipie; ++k) {
      j = hcoli[k];
      
      /*
       * We're eliminating column jpivot,
       * so we're eliminating the row it occurs in,
       * so every column in this row is becoming one shorter.
       *
       * I don't know why we don't do the same for rejected columns.
       *
       * if xrejct is false, then no column has ever been rejected
       * and this test wouldn't have to be made.
       * However, that means this whole loop would have to be copied.
       */
      if (! (clink[j].pre > nrow)) {
      C_EKK_REMOVE_LINK(hpivco, hincol, clink, j); /* (3) */
      }
      --hincol[j];
      
      kcs = mcstrt[j];
      kce = kcs + hincol[j];
      for (kc = kcs; kc <= kce; ++kc) {
      if (ipivot == hrowi[kc]) {
        break;
      }
      }
      /* ASSERT !(kc>kce) */
      
      /* (2) */
      hrowi[kc] = hrowi[kce];
      hrowi[kce] = 0;
      
      if (j == jpivot) {
      /* remember the slot corresponding to the pivot column */
      kpivot = k;
      }
      else {
      /*
       * We just reduced the length of the column.
       * If we haven't eliminated all of its elements completely,
       * then we have to put it back in its new length list.
       *
       * If the column was rejected, we only put it back in a length
       * list when it has been reduced to a singleton column,
       * because it would just be rejected again.
       */
      nzj = hincol[j];
      if (! (nzj <= 0) &&
          ! (clink[j].pre > nrow && nzj != 1)) {
        C_EKK_ADD_LINK(hpivco, nzj, clink, j); /* (3) */
      }
      }
    }
    assert (kpivot>0);
    
    /* store pivot sequence number */
    ++fact->npivots;
    rlink[ipivot].pre = -fact->npivots;
    clink[jpivot].pre = -fact->npivots;
    
    /* compute how much room we'll need later */
    fact->nuspike += hinrow[ipivot];
    
    /* check the pivot */
    pivot = dluval[kpivot];
    if (fabs(pivot) < drtpiv) {
      /* pivot element too small */
      small_pivot = true;
      rlink[ipivot].pre = -nrow - 1;
      clink[jpivot].pre = -nrow - 1;
      ++(*nsingp);
    }
    
    /* swap the pivoted column entry with the first entry in the row */
    dluval[kpivot] = dluval[kipis];
    dluval[kipis] = pivot;
    hcoli[kpivot] = hcoli[kipis];
    hcoli[kipis] = jpivot;
  }
  
  return (small_pivot);
} /* c_ekkcsin */


/*     Uwe H. Suhl, March 1987 */
/*     This routine processes row singletons during the computation of */
/*     an LU-decomposition for the nucleus. */
/*     Return codes (checked version 1.16): */
/*      -5: not enough space in row file */
/*      -6: not enough space in column file */
/*      7: pivot element too small */
/*     -52: system error at label 220 (ipivot not found) */
/*     -53: system error at label 400 (jpivot not found) */
int c_ekkrsin(EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          EKKHlink *mwork, int nfirst,
          int *nsingp, 
          int *xnewcop, int *xnewrop,
          int *nnentup, 
          int *kmxetap, int *ncompactionsp,
          int *nnentlp)
  
{
#if 1
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  //double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  const int nrow  = fact->nrow;
  const double drtpiv   = fact->drtpiv;
  
  int xnewro      = *xnewrop;
  int xnewco      = *xnewcop;
  int kmxeta      = *kmxetap;
  int nnentu      = *nnentup;
  int ncompactions      = *ncompactionsp;
  int nnentl      = *nnentlp;
  
  int i, j, k, kc, kr, npr, nzi;
  double pivot;
  int kjpis, kjpie, knprs, knpre;
  double elemnt, maxaij;
  int ipivot, epivco, lstart;
#ifndef NDEBUG
  int kpivot=-1;
#else
  int kpivot=-1;
#endif
  int irtcod = 0;
  const int nnetas      = fact->nnetas;
  
  lstart = nnetas - nnentl + 1;
  
  
  for (ipivot = hpivro[1]; ipivot > 0; ipivot = hpivro[1]) {
    const int jpivot = hcoli[mrstrt[ipivot]];
    
    kjpis = mcstrt[jpivot];
    kjpie = kjpis + hincol[jpivot] ;
    for (k = kjpis; k < kjpie; ++k) {
      i = hrowi[k];
      
      /*
       * We're eliminating row ipivot,
       * so we're eliminating the column it occurs in,
       * so every row in this column is becoming one shorter.
       *
       * No exception is made for rejected rows.
       */
      C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
    }
    
    /* The pivot column is being eliminated */
    /* I don't know why there is an exception for rejected columns */
    if (! (clink[jpivot].pre > nrow)) {
      C_EKK_REMOVE_LINK(hpivco, hincol, clink, jpivot);
    }
    
    epivco = hincol[jpivot] - 1;
    kjpie = kjpis + epivco;
    for (kc = kjpis; kc <= kjpie; ++kc) {
      if (ipivot == hrowi[kc]) {
      break;
      }
    }
    /* ASSERT !(kc>kjpie) */
    
    /* move the last column entry into this deleted one to keep */
    /* the entries compact */
    hrowi[kc] = hrowi[kjpie];
    
    hrowi[kjpie] = 0;
    
    /* store pivot sequence number */
    ++fact->npivots;
    rlink[ipivot].pre = -fact->npivots;
    clink[jpivot].pre = -fact->npivots;
    
    /* Check if row or column files have to be compressed */
    if (! (xnewro + epivco < lstart)) {
      if (! (nnentu + epivco < lstart)) {
      return (-5);
      }
      {
      int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
      kmxeta += xnewro - iput ;
      xnewro = iput - 1;
      ++ncompactions;
      }
    }
    if (! (xnewco + epivco < lstart)) {
      if (! (nnentu + epivco < lstart)) {
      return (-5);
      }
      xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
      ++ncompactions;
    }
    
    /* This column has no more entries in it */
    hincol[jpivot] = 0;
    
    /* Perform numerical part of elimination. */
    pivot = dluval[mrstrt[ipivot]];
    if (fabs(pivot) < drtpiv) {
      irtcod = 7;
      rlink[ipivot].pre = -nrow - 1;
      clink[jpivot].pre = -nrow - 1;
      ++(*nsingp);
    }
    
    /* If epivco is 0, then we can treat this like a singleton column (?)*/
    if (! (epivco <= 0)) {
      ++fact->xnetal;
      mcstrt[fact->xnetal] = lstart - 1;
      hpivco[fact->xnetal] = ipivot;
      
      /* Loop over nonzeros in pivot column. */
      kjpis = mcstrt[jpivot];
      kjpie = kjpis + epivco ;
      nnentl+=epivco;
      nnentu-=epivco;
      for (kc = kjpis; kc < kjpie; ++kc) {
      npr = hrowi[kc];
      /* zero out the row entries as we go along */
      hrowi[kc] = 0;
      
      /* each row in the column is getting shorter */
      --hinrow[npr];
      
      /* find the entry in this row for the pivot column */
      knprs = mrstrt[npr];
      knpre = knprs + hinrow[npr];
      for (kr = knprs; kr <= knpre; ++kr) {
        if (jpivot == hcoli[kr])
          break;
      }
      /* ASSERT !(kr>knpre) */
      
      elemnt = dluval[kr];
      /* move the last pivot column entry into this one */
      /* to keep entries compact */
      dluval[kr] = dluval[knpre];
      hcoli[kr] = hcoli[knpre];
      
      /*
       * c_ekkmltf put the largest entries in front, and 
       * we want to maintain that property.
       * There is only a problem if we just pivoted out the first
       * entry, and there is more than one entry in the list.
       */
      if (! (kr != knprs || hinrow[npr] <= 1)) {
        maxaij = 0.f;
        for (k = knprs; k <= knpre; ++k) {
          if (! (fabs(dluval[k]) <= maxaij)) {
            maxaij = fabs(dluval[k]);
            kpivot = k;
          }
        }
        assert (kpivot>0);
        maxaij = dluval[kpivot];
        dluval[kpivot] = dluval[knprs];
        dluval[knprs] = maxaij;
        
        j = hcoli[kpivot];
        hcoli[kpivot] = hcoli[knprs];
        hcoli[knprs] = j;
      }
      
      /* store elementary row transformation */
      --lstart;
      dluval[lstart] = -elemnt / pivot;
      hrowi[lstart] = SHIFT_INDEX(npr);
      
      /* Only add the row back in a length list if it isn't empty */
      nzi = hinrow[npr];
      if (! (nzi <= 0)) {
        C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
      }
      }
      ++fact->nuspike;
    }
  }
  
  *xnewrop = xnewro;
  *xnewcop = xnewco;
  *kmxetap = kmxeta;
  *nnentup = nnentu;
  *ncompactionsp = ncompactions;
  *nnentlp = nnentl;
  
  return (irtcod);
} /* c_ekkrsin */


int c_ekkfpvt(const EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          int *nsingp, int *xrejctp,
          int *xipivtp, int *xjpivtp)
{
  double zpivlu = fact->zpivlu;
#if 1
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  //double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, j, k, ke, kk, ks, nz, nz1, kce, kcs, kre, krs;
  double minsze;
  int marcst, mincst, mincnt, trials, nentri;
  int jpivot=-1;
  bool rjectd;
  int ipivot;
  const int nrow  = fact->nrow;
  int irtcod = 0;
  
  /* this used to be initialized in c_ekklfct */
  const int xtrial = 1;
  
  trials = 0;
  ipivot = 0;
  mincst = COIN_INT_MAX;
  mincnt = COIN_INT_MAX;
  for (nz = 2; nz <= nrow; ++nz) {
    nz1 = nz - 1;
    if (mincnt <= nz) {
      goto L900;
    }
    
    /* Search rows for a pivot */
    for (i = hpivro[nz]; ! (i <= 0); i = rlink[i].suc) {
      
      ks = mrstrt[i];
      ke = ks + nz - 1;
      /* Determine magnitude of minimal acceptable element */
      minsze = fabs(dluval[ks]) * zpivlu;
      for (k = ks; k <= ke; ++k) {
      /* Consider a column only if it passes the stability test */
      if (! (fabs(dluval[k]) < minsze)) {
        j = hcoli[k];
        marcst = nz1 * hincol[j];
        if (! (marcst >= mincst)) {
          mincst = marcst;
          mincnt = hincol[j];
          ipivot = i;
          jpivot = j;
          if (mincnt <= nz + 1) {
            goto L900;
          }
        }
      }
      }
      ++trials;
      
      if (trials >= xtrial) {
      goto L900;
      }
    }
    
    /* Search columns for a pivot */
    j = hpivco[nz];
    while (! (j <= 0)) {
      /* XSEARD = XSEARD + 1 */
      rjectd = false;
      kcs = mcstrt[j];
      kce = kcs + nz - 1;
      for (k = kcs; k <= kce; ++k) {
      i = hrowi[k];
      nentri = hinrow[i];
      marcst = nz1 * nentri;
      if (! (marcst >= mincst)) {
        /* Determine magnitude of minimal acceptable element */
        minsze = fabs(dluval[mrstrt[i]]) * zpivlu;
        krs = mrstrt[i];
        kre = krs + nentri - 1;
        for (kk = krs; kk <= kre; ++kk) {
          if (hcoli[kk] == j)
            break;
        }
        /* ASSERT (kk <= kre) */
        
        /* perform stability test */
        if (! (fabs(dluval[kk]) < minsze)) {
          mincst = marcst;
          mincnt = nentri;
          ipivot = i;
          jpivot = j;
          rjectd = false;
          if (mincnt <= nz) {
            goto L900;
          }
        }
        else {
          if (ipivot == 0) {
            rjectd = true;
          }
        }
      }
      }
      ++trials;
      if (trials >= xtrial && ipivot > 0) {
      goto L900;
      }
      if (rjectd) {
      int jsuc = clink[j].suc;
      ++(*xrejctp);
      C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
      clink[j].pre = nrow + 1;
      j = jsuc;
      }
      else {
      j = clink[j].suc;
      }
    }
  }
  
  /* FLAG REJECTED ROWS (should this be columns ?) */
  for (j = 1; j <= nrow; ++j) {
    if (hinrow[j] == 0) {
      rlink[j].pre = -nrow - 1;
      ++(*nsingp);
    }
  }
  irtcod = 10;
  
 L900:
  *xipivtp = ipivot;
  *xjpivtp = jpivot;
  return (irtcod);
} /* c_ekkfpvt */
void c_ekkprpv(EKKfactinfo *fact,
           EKKHlink *rlink, EKKHlink *clink,
           
           int xrejct,
           int ipivot, int jpivot)
{
#if 1
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  //double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, k;
  int kc;
  double pivot;
  int kipis = mrstrt[ipivot];
  int kipie = kipis + hinrow[ipivot] - 1;
  
#ifndef NDEBUG
  int kpivot=-1;
#else
  int kpivot=-1;
#endif
  const int nrow  = fact->nrow;
  
  /*     Update data structures */
  {
    int kjpis = mcstrt[jpivot];
    int kjpie = kjpis + hincol[jpivot] ;
    for (k = kjpis; k < kjpie; ++k) {
      i = hrowi[k];
      C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
    }
  }
  
  for (k = kipis; k <= kipie; ++k) {
    int j = hcoli[k];
    
    if ((xrejct == 0) ||
      ! (clink[j].pre > nrow)) {
      C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
    }
    
    --hincol[j];
    int kcs = mcstrt[j];
    int kce = kcs + hincol[j];
    
    for (kc = kcs; kc < kce ; kc ++) {
      if (hrowi[kc] == ipivot) 
      break;
    }
    assert (kc<kce||hrowi[kce]==ipivot);
    hrowi[kc] = hrowi[kce];
    hrowi[kce] = 0;
    if (j == jpivot) {
      kpivot = k;
    }
  }
  assert (kpivot>0);
  
  /* Store the pivot sequence number */
  ++fact->npivots;
  rlink[ipivot].pre = -fact->npivots;
  clink[jpivot].pre = -fact->npivots;
  
  pivot = dluval[kpivot];
  dluval[kpivot] = dluval[kipis];
  dluval[kipis] = pivot;
  hcoli[kpivot] = hcoli[kipis];
  hcoli[kipis] = jpivot;
  
} /* c_ekkprpv */

/*
 * c_ekkclco is almost exactly like c_ekkrwco.
 */
int c_ekkclco(const EKKfactinfo *fact,int *hcoli, int *mrstrt, int *hinrow, int xnewro)
{
#if 0
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, k, nz, kold;
  int kstart;
  const int nrow  = fact->nrow;
  
  for (i = 1; i <= nrow; ++i) {
    nz = hinrow[i];
    if (0 < nz) {
      /* save the last column entry of row i in hinrow */
      /* and replace that entry with -i */
      k = mrstrt[i] + nz - 1;
      hinrow[i] = hcoli[k];
      hcoli[k] = -i;
    }
  }
  
  kstart = 0;
  kold = 0;
  for (k = 1; k <= xnewro; ++k) {
    if (hcoli[k] != 0) {
      ++kstart;
      
      /* if this is the last entry for the row... */
      if (hcoli[k] < 0) {
      /* restore the entry */
      i = -hcoli[k];
      hcoli[k] = hinrow[i];
      
      /* update mrstart and hinrow */
      mrstrt[i] = kold + 1;
      hinrow[i] = kstart - kold;
      kold = kstart;
      }
      
      hcoli[kstart] = hcoli[k];
    }
  }
  
  /* INSERTED INCASE CALLED FROM YTRIAN JJHF */
  mrstrt[nrow + 1] = kstart + 1;
  
  return (kstart);
} /* c_ekkclco */



#undef MACTION_T
#define COIN_OSL_CMFC
#define MACTION_T short int
int c_ekkcmfc(EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          EKKHlink *mwork, void *maction_void,
          int nnetas,
          int *nsingp, int *xrejctp,
          int *xnewrop, int xnewco, 
          int *ncompactionsp)
  
#include "CoinOslC.h"
#undef COIN_OSL_CMFC
#undef MACTION_T
  static int c_ekkidmx(int n, const double *dx)
{
  int ret_val;
  int i;
  double dmax;
  --dx;
  
  /* Function Body */
  
  if (n < 1) {
    return (0);
  }
  
  if (n == 1) {
    return (1);
  }
  
  ret_val = 1;
  dmax = fabs(dx[1]);
  for (i = 2; i <= n; ++i) {
    if (fabs(dx[i]) > dmax) {
      ret_val = i;
      dmax = fabs(dx[i]);
    }
  }
  
  return ret_val;
} /* c_ekkidmx */
/*     Return codes in IRTCOD/IRTCOD are */
/*     4: numerical problems */
/*     5: not enough space in row file */
/*     6: not enough space in column file */
int c_ekkcmfd(EKKfactinfo *fact,
          int *mcol, 
          EKKHlink *rlink, EKKHlink *clink,
          int *maction,
          int nnetas,
          int *nnentlp, int *nnentup,
          int *nsingp)
{
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
  int nnentl      = *nnentlp;
  int nnentu      = *nnentup;
  int storeZero = fact->ndenuc;
  
  int mkrs[8];
  double dpivyy[8];
  
  /* Local variables */
  int i, j;
  double d0, dx;
  int nz, ndo, krs;
  int kend, jcol;
  int irow, iput, jrow, krxs;
  int mjcol[8];
  double pivot;
  int count;
  int ilast, isort;
  double dpivx, dsave;
  double dpivxx[8];
  double multip;
  int lstart, ndense, krlast, ifirst, krfirst, kcount, idense, ipivot,
    jdense, kchunk, jpivot;
  
  const int nrow  = fact->nrow;
  
  int irtcod = 0;
  
  lstart = nnetas - nnentl + 1;
  
  /* put list of columns in last HROWI */
  /* fix row order once for all */
  ndense = nrow - fact->npivots;
  iput = ndense + 1;
  for (i = 1; i <= nrow; ++i) {
    if (hpivro[i] > 0) {
      irow = hpivro[i];
      for (j = 1; j <= nrow; ++j) {
      --iput;
      maction[iput] = irow;
      irow = rlink[irow].suc;
      if (irow == 0) {
        break;
      }
      }
    }
  }
  if (iput != 1) {
    ++(*nsingp);
  }
  else {
    /*     Use HCOLI just for last row */
    ilast = maction[1];
    krlast = mrstrt[ilast];
    ifirst = maction[ndense];
    krfirst = mrstrt[ifirst];
    /*     put list of columns in last HCOLI */
    iput = 0;
    for (i = 1; i <= nrow; ++i) {
      if (clink[i].pre >= 0) {
      hcoli[krlast + iput] = i;
      ++iput;
      }
    }
    if (iput != ndense) {
      ++(*nsingp);
    }
    else {
      ndo = ndense / 8;
      /*     do most */
      for (kcount = 1; kcount <= ndo; ++kcount) {
        idense = ndense;
        isort = 8;
        for (count = ndense; count >= ndense - 7; --count) {
        ipivot = maction[count];
        krs = mrstrt[ipivot];
        --isort;
        mkrs[isort] = krs;
        }
        isort = 8;
        for (count = ndense; count >= ndense - 7; --count) {
        /* Find a pivot element */
        --isort;
        ipivot = maction[count];
        krs = mkrs[isort];
        jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
        pivot = dluval[krs + jcol];
        --idense;
        mcol[count] = jcol;
        mjcol[isort] = mcol[count];
        dluval[krs + jcol] = dluval[krs + idense];
        if (fabs(pivot) < fact->zeroTolerance) {
          pivot = 0.;
          dpivx = 0.;
        } else {
          dpivx = 1. / pivot;
        }
        dluval[krs + idense] = pivot;
        dpivxx[isort] = dpivx;
        for (j = isort - 1; j >= 0; --j) {
          krxs = mkrs[j];
          multip = -dluval[krxs + jcol] * dpivx;
          dluval[krxs + jcol] = dluval[krxs + idense];
          /*           for moment skip if zero */
          if (fabs(multip) > fact->zeroTolerance) {
            for (i = 0; i < idense; ++i) {
            dluval[krxs + i] += multip * dluval[krs + i];
            }
          } else {
            multip = 0.;
          }
          dluval[krxs + idense] = multip;
        }
        }
      /*       sort all U in rows already done */
        for (i = 7; i >= 0; --i) {
        /* ****     this is important bit */
        krs = mkrs[i];
        for (j = i - 1; j >= 0; --j) {
          jcol = mjcol[j];
          dsave = dluval[krs + jcol];
          dluval[krs + jcol] = dluval[krs + idense + j];
          dluval[krs + idense + j] = dsave;
        }
        }
      /*       leave IDENSE as it is */
        if (ndense <= 400) {
        for (jrow = ndense - 8; jrow >= 1; --jrow) {
          irow = maction[jrow];
          krxs = mrstrt[irow];
          for (j = 7; j >= 0; --j) {
            jcol = mjcol[j];
            dsave = dluval[krxs + jcol];
            dluval[krxs + jcol] = dluval[krxs + idense + j];
            dluval[krxs + idense + j] = dsave;
          }
          for (j = 7; j >= 0; --j) {
            krs = mkrs[j];
            jdense = idense + j;
            dpivx = dpivxx[j];
            multip = -dluval[krxs + jdense] * dpivx;
            if (fabs(multip) <= fact->zeroTolerance) {
            multip = 0.;
            }
            dpivyy[j] = multip;
            dluval[krxs + jdense] = multip;
            for (i = idense; i < jdense; ++i) {
            dluval[krxs + i] += multip * dluval[krs + i];
            }
          }
          for (i = 0; i < idense; ++i) {
            dx = dluval[krxs + i];
            d0 = dpivyy[0] * dluval[mkrs[0] + i];
            dx += dpivyy[1] * dluval[mkrs[1] + i];
            d0 += dpivyy[2] * dluval[mkrs[2] + i];
            dx += dpivyy[3] * dluval[mkrs[3] + i];
            d0 += dpivyy[4] * dluval[mkrs[4] + i];
            dx += dpivyy[5] * dluval[mkrs[5] + i];
            d0 += dpivyy[6] * dluval[mkrs[6] + i];
            dx += dpivyy[7] * dluval[mkrs[7] + i];
            dluval[krxs + i] = d0 + dx;
          }
        }
        } else {
        for (jrow = ndense - 8; jrow >= 1; --jrow) {
          irow = maction[jrow];
          krxs = mrstrt[irow];
          for (j = 7; j >= 0; --j) {
            jcol = mjcol[j];
            dsave = dluval[krxs + jcol];
            dluval[krxs + jcol] = dluval[krxs + idense + j];
            dluval[krxs + idense + j] = dsave;
          }
          for (j = 7; j >= 0; --j) {
            krs = mkrs[j];
            jdense = idense + j;
            dpivx = dpivxx[j];
            multip = -dluval[krxs + jdense] * dpivx;
            if (fabs(multip) <= fact->zeroTolerance) {
            multip = 0.;
            }
            dluval[krxs + jdense] = multip;
            for (i = idense; i < jdense; ++i) {
            dluval[krxs + i] += multip * dluval[krs + i];
            }
          }
        }
        for (kchunk = 0; kchunk < idense; kchunk += 400) {
          kend = CoinMin(idense - 1, kchunk + 399);
          for (jrow = ndense - 8; jrow >= 1; --jrow) {
            irow = maction[jrow];
            krxs = mrstrt[irow];
            for (j = 7; j >= 0; --j) {
            dpivyy[j] = dluval[krxs + idense + j];
            }
            for (i = kchunk; i <= kend; ++i) {
            dx = dluval[krxs + i];
            d0 = dpivyy[0] * dluval[mkrs[0] + i];
            dx += dpivyy[1] * dluval[mkrs[1] + i];
            d0 += dpivyy[2] * dluval[mkrs[2] + i];
            dx += dpivyy[3] * dluval[mkrs[3] + i];
            d0 += dpivyy[4] * dluval[mkrs[4] + i];
            dx += dpivyy[5] * dluval[mkrs[5] + i];
            d0 += dpivyy[6] * dluval[mkrs[6] + i];
            dx += dpivyy[7] * dluval[mkrs[7] + i];
            dluval[krxs + i] = d0 + dx;
            }
          }
        }
        }
      /*       resort all U in rows already done */
        for (i = 7; i >= 0; --i) {
        krs = mkrs[i];
        for (j = 0; j < i; ++j) {
          jcol = mjcol[j];
          dsave = dluval[krs + jcol];
          dluval[krs + jcol] = dluval[krs + idense + j];
          dluval[krs + idense + j] = dsave;
        }
        }
        ndense += -8;
      }
      idense = ndense;
      /*     do remainder */
      for (count = ndense; count >= 1; --count) {
      /*        Find a pivot element */
        ipivot = maction[count];
        krs = mrstrt[ipivot];
        jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
        pivot = dluval[krs + jcol];
        --idense;
        mcol[count] = jcol;
        dluval[krs + jcol] = dluval[krs + idense];
        if (fabs(pivot) < fact->zeroTolerance) {
        dluval[krs + idense] = 0.;
        } else {
        dpivx = 1. / pivot;
        dluval[krs + idense] = pivot;
        for (jrow = idense; jrow >= 1; --jrow) {
          irow = maction[jrow];
          krxs = mrstrt[irow];
          multip = -dluval[krxs + jcol] * dpivx;
          dluval[krxs + jcol] = dluval[krxs + idense];
          /*           for moment skip if zero */
          if (fabs(multip) > fact->zeroTolerance) {
            dluval[krxs + idense] = multip;
            for (i = 0; i < idense; ++i) {
            dluval[krxs + i] += multip * dluval[krs + i];
            }
          } else {
            dluval[krxs + idense] = 0.;
          }
        }
        }
      }
      /*     now create in form for OSL */
      ndense = nrow - fact->npivots;
      idense = ndense;
      for (count = ndense; count >= 1; --count) {
      /*        Find a pivot element */
        ipivot = maction[count];
        krs = mrstrt[ipivot];
        --idense;
        jcol = mcol[count];
        jpivot = hcoli[krlast + jcol];
        ++fact->npivots;
        pivot = dluval[krs + idense];
        if (pivot == 0.) {
        hinrow[ipivot] = 0;
        rlink[ipivot].pre = -nrow - 1;
        ++(*nsingp);
        irtcod = 10;
        } else {
        rlink[ipivot].pre = -fact->npivots;
        clink[jpivot].pre = -fact->npivots;
        hincol[jpivot] = 0;
        ++fact->xnetal;
        mcstrt[fact->xnetal] = lstart - 1;
        hpivco[fact->xnetal] = ipivot;
        for (jrow = idense; jrow >= 1; --jrow) {
          irow = maction[jrow];
          krxs = mrstrt[irow];
          multip = dluval[krxs + idense];
          /*           for moment skip if zero */
          if (multip != 0.||storeZero) {
            /* Store elementary row transformation */
            ++nnentl;
            --nnentu;
            --lstart;
            dluval[lstart] = multip;
            hrowi[lstart] = SHIFT_INDEX(irow);
          }
        }
        hcoli[krlast + jcol] = hcoli[krlast + idense];
        /*         update pivot row and last row HCOLI */
        dluval[krs + idense] = dluval[krs];
        hcoli[krlast + idense] = hcoli[krlast];
        nz = 1;
        dluval[krs] = pivot;
        hcoli[krs] = jpivot;
        if (!storeZero) {
          for (i = 1; i <= idense; ++i) {
            if (fabs(dluval[krs + i]) > fact->zeroTolerance) {
            ++nz;
            hcoli[krs + nz - 1] = hcoli[krlast + i];
            dluval[krs + nz - 1] = dluval[krs + i];
            }
          }
          hinrow[ipivot] = nz;
        } else {
          for (i = 1; i <= idense; ++i) {
            ++nz;
            hcoli[krs + nz - 1] = hcoli[krlast + i];
            dluval[krs + nz - 1] = dluval[krs + i];
          }
          hinrow[ipivot] = nz;
        }
        }
      }
    }
  }
  
  *nnentlp = nnentl;
  *nnentup = nnentu;
  
  return (irtcod);
} /* c_ekkcmfd */
/* ***C_EKKCMFC */

/*
 * Generate a variant of c_ekkcmfc that uses an maction array of type
 * int rather than short.
 */
#undef MACTION_T
#define C_EKKCMFY
#define COIN_OSL_CMFC
#define MACTION_T int
int c_ekkcmfy(EKKfactinfo *fact,
          EKKHlink *rlink, EKKHlink *clink,
          EKKHlink *mwork, void *maction_void,
          int nnetas,
          int *nsingp, int *xrejctp,
          int *xnewrop, int xnewco, 
          int *ncompactionsp)
  
#include "CoinOslC.h"
#undef COIN_OSL_CMFC
#undef C_EKKCMFY
#undef MACTION_T
int c_ekkford(const EKKfactinfo *fact,const int *hinrow, const int *hincol,
          int *hpivro, int *hpivco,
          EKKHlink *rlink, EKKHlink *clink)
{
  int i, iri, nzi;
  const int nrow  = fact->nrow;
  int nsing = 0;
  
  /*     Uwe H. Suhl, August 1986 */
  /*     Builds linked lists of rows and cols of nucleus for efficient */
  /*     pivot searching. */

  memset(hpivro+1,0,nrow*sizeof(int));
  memset(hpivco+1,0,nrow*sizeof(int));
  for (i = 1; i <= nrow; ++i) {
    //hpivro[i] = 0;
    //hpivco[i] = 0;
    assert(rlink[i].suc == 0);
    assert(clink[i].suc == 0);
  }
  
  /*     Generate double linked list of rows having equal numbers of */
  /*     nonzeros in each row. Skip pivotal rows. */
  for (i = 1; i <= nrow; ++i) {
    if (! (rlink[i].pre < 0)) {
      nzi = hinrow[i];
      if (nzi <= 0) {
      ++nsing;
      rlink[i].pre = -nrow - 1;
      }
      else {
      iri = hpivro[nzi];
      hpivro[nzi] = i;
      rlink[i].suc = iri;
      rlink[i].pre = 0;
      if (iri != 0) {
        rlink[iri].pre = i;
      }
      }
    }
  }
  
  /*     Generate double linked list of cols having equal numbers of */
  /*     nonzeros in each col. Skip pivotal cols. */
  for (i = 1; i <= nrow; ++i) {
    if (! (clink[i].pre < 0)) {
      nzi = hincol[i];
      if (nzi <= 0) {
      ++nsing;
      clink[i].pre = -nrow - 1;
      }
      else {
      iri = hpivco[nzi];
      hpivco[nzi] = i;
      clink[i].suc = iri;
      clink[i].pre = 0;
      if (iri != 0) {
        clink[iri].pre = i;
      }
      }
    }
  }
  
  return (nsing);
} /* c_ekkford */

/*    c version of OSL from 36100 */

/*     Assumes that a basis exists in correct form */
/*     Calls Uwe's routines (approximately) */
/*     Then if OK shuffles U into column order */
/*     Return codes: */

/*     0: everything ok */
/*     1: everything ok but performance would be better if more space */
/*        would be make available */
/*     4: growth rate of element in U too big */
/*     5: not enough space in row file */
/*     6: not enough space in column file */
/*     7: pivot too small - col sing */
/*     8: pivot too small - row sing */
/*    10: matrix is singular */

/* I suspect c_ekklfct never returns 1 */
/*
 * layout of data
 *
 * dluval/hcoli:  (L^-1)B     - hole - L factors
 *
 * The L factors are written from high to low, starting from nnetas.
 * There are nnentl factors in L.  lstart the next entry to use for the
 * L factors.  Eventually, (L^-1)B turns into U.
 
 * The ninbas coefficients of matrix B are originally in the start of
 * dluval/hcoli.  As L transforms it, rows may have to be expanded.
 * If there is room, they are copied to the start of the hole,
 * otherwise the first part of this area is compacted, and hopefully
 * there is then room.
 * There are nnentu coefficients in (L^-1)B.
 * nnentu + nnentl >= ninbas.
 * nnentu + nnentl == ninbas if there has been no fill-in.
 * nnentu is decreased when the pivot eliminates elements
 * (in which case there is a corresponding increase in nnentl),
 * and if pivoting happens to cancel out factors (in which case
 * there is no corresponding increase in L).
 * nnentu is increased if there is fill-in (no decrease in L).
 * If nnentu + nnentl >= nnetas, then we've run out of room.
 * It is not the case that the elements of (L^-1)B are all in the
 * first nnentu positions of dluval/hcoli, but that is of course
 * the lower bound on the number of positions needed to store it.
 * nuspik is roughly the sum of the row lengths of the rows that were pivoted
 * out.  singleton rows in c_ekktria do not change nuspik, but
 * c_ekkrsin does increment it for each singleton row.
 * That is, there are nuspik elements that in the upper part of (L^-1)B,
 * and (nnentu - nuspik) elements left in B.
 */
/*
 * As part of factorization, we test candidate pivots for numerical
 * stability; if the largest element in a row/col is much larger than
 * the smallest, this generally causes problems.  To easily determine
 * what the largest element is, we ensure that it is always in front.
 * This establishes this property; later on we take steps to preserve it.
 */
static void c_ekkmltf(const EKKfactinfo *fact,double *dluval, int *hcoli,
                const int *mrstrt, const int *hinrow,
                const EKKHlink *rlink)
{
#if 0
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, j, k;
#ifndef NDEBUG
  int koff=-1;
#else
  int koff;
#endif
  const int nrow  = fact->nrow;
  
  
  for (i = 1; i <= nrow; ++i) {
    /* ignore rows that have already been pivoted */
    /* if it is a singleton row, the property trivially holds */
    if (! (rlink[i].pre < 0 || hinrow[i] <= 1)) {
      const int krs = mrstrt[i];
      const int kre = krs + hinrow[i] - 1;
      
      double maxaij = 0.f;
      
      /* this assumes that at least one of the dluvals is non-zero. */
      for (k = krs; k <= kre; ++k) {
      if (! (fabs(dluval[k]) <= maxaij)) {
        maxaij = fabs(dluval[k]);
        koff = k;
      }
      }
      assert (koff>0);
      maxaij = dluval[koff];
      j = hcoli[koff];
      
      dluval[koff] = dluval[krs];
      hcoli[koff] = hcoli[krs];
      
      dluval[krs] = maxaij;
      hcoli[krs] = j;
    }
  }
} /* c_ekkmltf */
int c_ekklfct( register EKKfactinfo *fact)
{
  const int nrow  = fact->nrow;
  int ninbas = fact->xcsadr[nrow+1]-1;
  int ifvsol = fact->ifvsol;
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr;
  int *hpivco     = fact->kcpadr;
  
  
  EKKHlink *rlink = fact->kp1adr;
  EKKHlink *clink = fact->kp2adr;
  EKKHlink *mwork = (reinterpret_cast<EKKHlink*>(fact->kw1adr))-1;
  
  int nsing, kdnspt, xnewro, xnewco;
  int i;
  int xrejct;
  int irtcod;
  const int nnetas      = fact->nnetas;
  
  int ncompactions;
  double save_drtpiv = fact->drtpiv;
  double save_zpivlu = fact->zpivlu;
  if (ifvsol > 0 && fact->invok < 0) {
    fact->zpivlu =  CoinMin(0.9, fact->zpivlu * 10.);
    fact->drtpiv=1.0e-8;
  }
  
  rlink --;
  clink --;
  
  /* Function Body */
  hcoli[nnetas] = 1;
  hrowi[nnetas] = 1;
  dluval[nnetas] = 0.0;
  /*     set amount of work */
  xrejct = 0;
  nsing = 0;
  kdnspt = nnetas + 1;
  fact->ndenuc = 0;
  /*     Triangularize */
  irtcod = c_ekktria(fact,rlink,clink,
               &nsing, 
               &xnewco, &xnewro,
               &ncompactions, ninbas);
  fact->nnentl = ninbas - fact->nnentu;
  
  if (irtcod < 0) {
    /* no space or system error */
    goto L8000;
  }
  
  if (irtcod != 0 && fact->invok >= 0) {
    goto L8500;   /* 7 or 8 - pivot too small */
  }
#if 0
  /* is this necessary ? */
  lstart = nnetas - fact->nnentl + 1;
  for (i = lstart; i <= nnetas; ++i) {
    hrowi[i] = (hcoli[i] << 3);
  }
#endif
  
  /* See if finished */
  if (! (fact->npivots >= nrow)) {
    int nsing1;
    
    /*     No - do nucleus */
    
    nsing1 = c_ekkford(fact,hinrow, hincol, hpivro, hpivco, rlink, clink);
    nsing+= nsing1;
    if (nsing1 != 0 && fact->invok >= 0) {
      irtcod=7;
      goto L8500;
    }
    c_ekkmltf(fact,dluval, hcoli, mrstrt, hinrow, rlink);
    
    {
      bool callcmfy = false;
      
      if (nrow > 32767) {
      int count = 0;
      for (i = 1; i <= nrow; ++i) {
        count = CoinMax(count,hinrow[i]);
      }
      if (count + nrow - fact->npivots > 32767) {
        /* will have to use I*4 version of CMFC */
        /* no changes to pointer params */
        callcmfy = true;
      }
      }
      
      irtcod = (callcmfy ? c_ekkcmfy : c_ekkcmfc)
      (fact, 
       rlink, clink,
       mwork, &mwork[nrow + 1],
       nnetas,
       &nsing, &xrejct,
       &xnewro, xnewco,
       &ncompactions);
      
      /* irtcod one of 0,-5,7,10 */
    }
    
    if (irtcod < 0) {
      goto L8000;
    }
    kdnspt = nnetas - fact->nnentl;
  }
  
  /*     return if error */
  if (nsing > 0 || irtcod == 10) {
    irtcod = 99;
  }
  /* irtcod one of 0,7,99 */
  
  if (irtcod != 0) {
    goto L8500;
  }
  ++fact->xnetal;
  mcstrt[fact->xnetal] = nnetas - fact->nnentl;
  
  /* give message if tight on memory */
  if (ncompactions > 2 ) {
    if (1) {
      int etasize =CoinMax(4*fact->nnentu+(nnetas-fact->nnentl)+1000,fact->eta_size);
      fact->eta_size=CoinMin(static_cast<int>(1.2*fact->eta_size),etasize);
      if (fact->maxNNetas>0&&fact->eta_size>
        fact->maxNNetas) {
      fact->eta_size=fact->maxNNetas;
      }
    } /* endif */
  }
  /*       Shuffle U and multiply L by 8 (if assembler) */
  {
    int jrtcod = c_ekkshff(fact, clink, rlink,
                     xnewro);
    
    /* nR_etas is the number of R transforms;
     * it is incremented only in c_ekketsj.
     */
    fact->nR_etas = 0;
    /*fact->R_etas_start = mcstrt+nrow+fact->nnentl+3;*/
    fact->R_etas_start[1] = /*kdnspt - 1*/0;    /* magic */
    fact->R_etas_index = &fact->xeradr[kdnspt - 1];
    fact->R_etas_element = &fact->xeeadr[kdnspt - 1];
    
    
    
    if (jrtcod != 0) {
      irtcod = jrtcod;
      /* irtcod == 2 */
    }
  }
  goto L8500;
  
  /* Fatal error */
 L8000:
  
  if (1) {
    if (fact->maxNNetas != fact->eta_size &&
      nnetas) {
      /* return and get more space */
      
      /* double eta_size, unless that exceeds max (if there is one) */
      fact->eta_size = fact->eta_size<<1;
      if (fact->maxNNetas > 0 &&
        fact->eta_size > fact->maxNNetas) {
      fact->eta_size = fact->maxNNetas;
      }
      return (5);
    }
  }
  /*c_ekkmesg_no_i1(121, -irtcod);*/
  irtcod = 3;
  
 L8500:
  /* restore pivot tolerance */
  fact->drtpiv=save_drtpiv;
  fact->zpivlu=save_zpivlu;
#ifndef NDEBUG
  if (fact->rows_ok) {
    int * hinrow=fact->xrnadr;
    if (!fact->xe2adr) {
      for (int i=1;i<=fact->nrow;i++) {
      assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
      }
    }
  }
#endif
  return (irtcod);
} /* c_ekklfct */


/*
  summary of return codes
  
  c_ekktria:
  7 small pivot
  -5 no memory
  
  c_ekkcsin:
  returns true if small pivot
  
  c_ekkrsin:
  -5 no memory
  7 small pivot
  
  c_ekkfpvt:
  10: no pivots found (singular)
  
  c_ekkcmfd:
  10: zero pivot (not just small)
  
  c_ekkcmfc:
  -5:  no memory
  any non-zero code from c_ekkcsin, c_ekkrsin, c_ekkfpvt, c_ekkprpv, c_ekkcmfd
  
  c_ekkshff:
  2:  singular
  
  c_ekklfct:
  any positive code from c_ekktria, c_ekkcmfc, c_ekkshff (2,7,10)
  *except* 10, which is changed to 99.
  all negative return codes are changed to 5 or 3
  (5 == ran out of memory but could get more,
  3 == ran out of memory, no luck)
  so:  2,3,5,7,99
  
  c_ekklfct1:
  1: c_ekksmem_invert failed
  2: c_ekkslcf/c_ekkslct ran out of room
  any return code from c_ekklfct, except 2 and 5
*/

void c_ekkrowq(int *hrow, int *hcol, double *dels, 
           int *mrstrt,
           const int *hinrow, int nnrow, int ninbas)
{
  int i, k, iak, jak;
  double daik;
  int iloc;
  double dsave;
  int isave, jsave;
  
  /* Order matrix rowwise using MRSTRT, DELS, HCOL */
  
  k = 1;
  /* POSITION AFTER END OF ROW */
  for (i = 1; i <= nnrow; ++i) {
    k += hinrow[i];
    mrstrt[i] = k;
  }
  
  for (k = ninbas; k >= 1; --k) {
    iak = hrow[k];
    if (iak != 0) {
      daik = dels[k];
      jak = hcol[k];
      hrow[k] = 0;
      while (1) {
      --mrstrt[iak];
      
      iloc = mrstrt[iak];
      
      dsave = dels[iloc];
      isave = hrow[iloc];
      jsave = hcol[iloc];
      dels[iloc] = daik;
      hrow[iloc] = 0;
      hcol[iloc] = jak;
      
      if (isave == 0)
        break;
      daik = dsave;
      iak = isave;
      jak = jsave;
      }
      
    }
  }
} /* c_ekkrowq */



int c_ekkrwco(const EKKfactinfo *fact,double *dluval, 
          int *hcoli, int *mrstrt, int *hinrow, int xnewro)
{
  int i, k, nz, kold;
  int kstart;
  const int nrow  = fact->nrow;
  
  for (i = 1; i <= nrow; ++i) {
    nz = hinrow[i];
    if (0 < nz) {
      /* save the last column entry of row i in hinrow */
      /* and replace that entry with -i */
      k = mrstrt[i] + nz - 1;
      hinrow[i] = hcoli[k];
      hcoli[k] = -i;
    }
  }
  
  kstart = 0;
  kold = 0;
  for (k = 1; k <= xnewro; ++k) {
    if (hcoli[k] != 0) {
      ++kstart;
      
      /* if this is the last entry for the row... */
      if (hcoli[k] < 0) {
      /* restore the entry */
      i = -hcoli[k];
      hcoli[k] = hinrow[i];
      
      /* update mrstart and hinrow */
      /* ACTUALLY, hinrow should already be accurate */
      mrstrt[i] = kold + 1;
      hinrow[i] = kstart - kold;
      kold = kstart;
      }
      
      /* move the entry */
      dluval[kstart] = dluval[k];
      hcoli[kstart] = hcoli[k];
    }
  }
  
  return (kstart);
} /* c_ekkrwco */



int c_ekkrwcs(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
          const int *hinrow, const EKKHlink *mwork,
          int nfirst)
{
#if 0
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, k, k1, k2, nz;
  int irow, iput;
  const int nrow  = fact->nrow;
  
  /*     Compress row file */
  
  iput = 1;
  irow = nfirst;
  for (i = 1; i <= nrow; ++i) {
    nz = hinrow[irow];
    k1 = mrstrt[irow];
    if (k1 != iput) {
      mrstrt[irow] = iput;
      k2 = k1 + nz - 1;
      for (k = k1; k <= k2; ++k) {
      dluval[iput] = dluval[k];
      hcoli[iput] = hcoli[k];
      ++iput;
      }
    } else {
      iput += nz;
    }
    irow = mwork[irow].suc;
  }
  
  return (iput);
} /* c_ekkrwcs */
void c_ekkrwct(const EKKfactinfo *fact,double *dluval, int *hcoli, int *mrstrt,
           const int *hinrow, const EKKHlink *mwork,
           const EKKHlink *rlink,
           const short *msort, double *dsort,
           int nlast, int xnewro)
{
#if 0
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
#endif
  int i, k, k1, nz, icol;
  int kmax;
  int irow, iput;
  int ilook;
  const int nrow  = fact->nrow;
  
  iput = xnewro;
  irow = nlast;
  kmax = nrow - fact->npivots;
  for (i = 1; i <= nrow; ++i) {
    nz = hinrow[irow];
    k1 = mrstrt[irow] - 1;
    if (rlink[irow].pre < 0) {
      /* pivoted on already */
      iput -= nz;
      if (k1 != iput) {
      mrstrt[irow] = iput + 1;
      for (k = nz; k >= 1; --k) {
        dluval[iput + k] = dluval[k1 + k];
        hcoli[iput + k] = hcoli[k1 + k];
      }
      }
    } else {
      /* not pivoted - going dense */
      iput -= kmax;
      mrstrt[irow] = iput + 1;
      c_ekkdzero( kmax, &dsort[1]);
      for (k = 1; k <= nz; ++k) {
      icol = hcoli[k1 + k];
      ilook = msort[icol];
      dsort[ilook] = dluval[k1 + k];
      }
      c_ekkdcpy(kmax,
            (dsort+1), (dluval+iput + 1));
    }
    irow = mwork[irow].pre;
  }
} /* c_ekkrwct */
/*     takes Uwe's modern structures and puts them back 20 years */
int c_ekkshff(EKKfactinfo *fact,
          EKKHlink *clink, EKKHlink *rlink, 
          int xnewro)
{
  int *hpivro     = fact->krpadr; 
  
  int i, j;
  int nbas, icol;
  int ipiv;
  const int nrow  = fact->nrow;
  int nsing;
  
  for (i = 1; i <= nrow; ++i) {
    j = -rlink[i].pre;
    rlink[i].pre = j;
    if (j > 0 && j <= nrow) {
      hpivro[j] = i;
    }
    j = -clink[i].pre;
    clink[i].pre = j;
  }
  /* hpivro[j] is now (hopefully) the row that was pivoted on step j */
  /* rlink[i].pre is the step in which row i was pivoted */
  
  nbas = 0;
  nsing = 0;
  /* Decide if permutation wanted */
  fact->first_dense=nrow-fact->ndenuc+1+1;
  fact->last_dense=nrow;
  
  /* rlink[].suc is dead at this point */
  
  /*
   * replace the the basis index 
   * with the pivot (or permuted) index generated by factorization.
   * This eventually goes into mpermu.
   */
  for (icol = 1; icol <= nrow; ++icol) {
    int ibasis = icol;
    ipiv = clink[ibasis].pre;
    
    if (0 < ipiv && ipiv <= nrow) {
      rlink[ibasis].suc = ipiv;
      ++nbas;
    }
  }
  
  nsing = nrow - nbas;
  if (nsing > 0) {
    abort();
  }
  
  /* if we reach here, then rlink[1..nrow].suc == clink[1..nrow].pre */
  
  /* switch off sparse update if any dense section */
  {
    const int notMuchRoom = (fact->nnentu + xnewro + 10 > fact->nnetas - fact->nnentl);
    
    /* must be same as in c_ekkshfv */
    if (fact->ndenuc || notMuchRoom||nrow<C_EKK_GO_SPARSE) {
#if PRINT_DEBUG
      if (fact->if_sparse_update) {
      printf("**** Switching off sparse update - dense - c_ekkshff\n");
      }
#endif
      fact->if_sparse_update=0;
    }
  }
  
  /* hpivro[1..nrow] is not read by c_ekkshfv */
  c_ekkshfv(fact,
          rlink, clink, 
        xnewro);
  
  return (0);
} /* c_ekkshff */
/* sorts on indices dragging elements with */
static void c_ekk_sort2(int * key , double * array2,int number)
{
  int minsize=10;
  int n = number;
  int sp;
  int *v = key;
  int *m, t;
  int * ls[32] , * rs[32];
  int *l , *r , c;
  double it;
  int j;
  /*check already sorted  */
#ifndef LONG_MAX
#define LONG_MAX 0x7fffffff;
#endif
  int last=-LONG_MAX;
  for (j=0;j<number;j++) {
    if (key[j]>=last) {
      last=key[j];
    } else {
      break;
    } /* endif */
  } /* endfor */
  if (j==number) {
    return;
  } /* endif */
  sp = 0 ; ls[sp] = v ; rs[sp] = v + (n-1) ;
  while( sp >= 0 )
    {
      if ( rs[sp] - ls[sp] > minsize )
      {
        l = ls[sp] ; r = rs[sp] ; m = l + (r-l)/2 ;
        if ( *l > *m )
          {
            t = *l ; *l = *m ; *m = t ;
            it = array2[l-v] ; array2[l-v] = array2[m-v] ; array2[m-v] = it ;
          }
        if ( *m > *r )
          {
            t = *m ; *m = *r ; *r = t ;
            it = array2[m-v] ; array2[m-v] = array2[r-v] ; array2[r-v] = it ;
            if ( *l > *m )
            {
              t = *l ; *l = *m ; *m = t ;
              it = array2[l-v] ; array2[l-v] = array2[m-v] ; array2[m-v] = it ;
            }
          }
        c = *m ;
        while ( r - l > 1 )
          {
            while ( *(++l) < c ) ;
            while ( *(--r) > c ) ;
            t = *l ; *l = *r ; *r = t ;
            it = array2[l-v] ; array2[l-v] = array2[r-v] ; array2[r-v] = it ;
          }
        l = r - 1 ;
        if ( l < m )
          {  ls[sp+1] = ls[sp] ;
            rs[sp+1] = l      ;
            ls[sp  ] = r      ;
          }
        else
          {  ls[sp+1] = r      ;
            rs[sp+1] = rs[sp] ;
            rs[sp  ] = l      ;
          }
        sp++ ;
      }
      else sp-- ;
    }
  for ( l = v , m = v + (n-1) ; l < m ; l++ )
    {  if ( *l > *(l+1) )
      {
        c = *(l+1) ;
        it = array2[(l-v)+1] ;
        for ( r = l ; r >= v && *r > c ; r-- )
          {
            *(r+1) = *r ;
            array2[(r-v)+1] = array2[(r-v)] ;
          }
        *(r+1) = c ;
        array2[(r-v)+1] = it ;
      }
    }
}
/*     For each row compute reciprocal of pivot element and  take out of */
/*     Also use HLINK(1 to permute column numbers */
/*     and HPIVRO to permute row numbers */
/*     Sort into column order as was stored by row */
/*     If Assembler then shift row numbers in L by 3 */
/*     Put column numbers in U for L-U update */
/*     and multiply U elements by - reciprocal of pivot element */
/*     and set up backward pointers for pivot rows */
void c_ekkshfv(EKKfactinfo *fact,
             EKKHlink *rlink, EKKHlink *clink,
           int xnewro)
{
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  double *dvalpv = fact->kw3adr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *hpivro     = fact->krpadr; 
  int *hpivco     = fact->kcpadr;
  double *dpermu = fact->kadrpm;
  double * de2val = fact->xe2adr ? fact->xe2adr-1: 0;
  int nnentu      = fact->nnentu;
  int xnetal      = fact->xnetal;
  
  int numberSlacks; /* numberSlacks not read */
  
  int i, j, k, kk, nel;
  int nroom;
  bool need_more_space;
  int ndenuc=fact->ndenuc;
  int if_sparse_update=fact->if_sparse_update; 
  int nnentl = fact->nnentl;
  int nnetas = fact->nnetas;
  
  int *ihlink     = (reinterpret_cast<int*> (clink))+1;     /* can't use rlink for simple loop below */
  
  const int nrow        = fact->nrow;
  const int maxinv      = fact->maxinv;
  
  /* this is not just a temporary - c_ekkbtrn etc use this */
  int *mpermu     = (reinterpret_cast<int*> (dpermu+nrow))+1;
  
  int * temp = ihlink+nrow;
  int * temp2 = temp+nrow;
  const int notMuchRoom = (nnentu + xnewro + 10 > nnetas - nnentl);
  
  /* compress hlink and make simpler */
  for (i = 1; i <= nrow; ++i) {
    mpermu[i] = rlink[i].pre;
    ihlink[i] = rlink[i].suc;
  }
  /* mpermu[i] == the step in which row i was pivoted */
  /* ihlink[i] == the step in which col i was pivoted */
  
  /* must be same as in c_ekkshff */
  if (fact->ndenuc||notMuchRoom||nrow<C_EKK_GO_SPARSE) {
    int ninbas;
    
    /* CHANGE COLUMN NUMBERS AND FILL IN RECIPROCALS */
    /* ALSO RECOMPUTE NUMBER IN COLUMN */
    /* initialize with a fake pivot in each column */
    c_ekkscpy_0_1(nrow, 1, &hincol[1]);
    
    if (notMuchRoom) {
      fact->eta_size=static_cast<int>(1.05*fact->eta_size);
      
      /* eta_size can be no larger than maxNNetas */
      if (fact->maxNNetas > 0 &&
        fact->eta_size > fact->maxNNetas) {
      fact->eta_size=fact->maxNNetas;
      }
    } /* endif */
    
    /* For each row compute reciprocal of pivot element and take out of U */
    /* Also use ihlink to permute column numbers */
    /* the rows are not stored compactly or in order,
     * so we have to find out where the last one is stored */
    ninbas=0;
    for (i = 1; i <= nrow; ++i) {
      int jpiv=mpermu[i];
      int nin=hinrow[i];
      int krs = mrstrt[i];
      int kre = krs + nin;
      
      temp[jpiv]=krs;
      temp2[jpiv]=nin;
      
      ninbas = CoinMax(kre, ninbas);
      
      /* c_ekktria etc ensure that the first row entry is the pivot */
      dvalpv[jpiv] = 1. / dluval[krs];
      hcoli[krs] = 0;   /* probably needed for c_ekkrowq */
      /* room for the pivot has already been allocated, so hincol ok */
      
      for (kk = krs + 1; kk < kre; ++kk) {
      int j = ihlink[hcoli[kk]];
      hcoli[kk] = j;          /* permute the col index */
      hrowi[kk] = jpiv; /* permute the row index */
      ++hincol[j];
      }
    }
    /* temp [mpermu[i]] == mrstrt[i] */
    /* temp2[mpermu[i]] == hinrow[i] */
    
    ninbas--;     /* ???? */
    c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
    c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
    
    /* now mrstrt, hinrow, hcoli and hrowi have been permuted */
    
    /* Sort into column order as was stored by row */
    /* There will be an empty entry in front of each each column,
     * because we initialized hincol to 1s, and c_ekkrowq fills in
     * entries from the back */
    c_ekkrowq(hcoli, hrowi, dluval, mcstrt, hincol, nrow, ninbas);
    
    
    /* The shuffle zeroed out column pointers */
    /* Put them back for L-U update */
    /* Also multiply U elements by - reciprocal of pivot element */
    /* Also decrement mcstrt/hincol to give "real" sizes */
    for (i = 1; i <= nrow; ++i) {
      int kx = --mcstrt[i];
      nel = --hincol[i];
      hrowi[kx] = nel;
      dluval[kx] = dvalpv[i];
#ifndef NO_SHIFT
      for (int j=kx+1;j<=kx+nel;j++)
      hrowi[j] = SHIFT_INDEX(hrowi[j]);
#endif
    }
    
    /* sort dense part */
    for (i=nrow-ndenuc+1; i<=nrow; i++) {
      int kx = mcstrt[i]+1;   /* "real" entries start after pivot */
      int nel = hincol[i];
      c_ekk_sort2(&hrowi[kx],&dluval[kx],nel);
    }
    
    /* Recompute number in U */
    nnentu = mcstrt[nrow] + hincol[nrow];
    mcstrt[nrow + 4] = nnentu + 1;  /* magic - AND DEAD */
    
    /* as not much room switch off fast etas */
    mrstrt[1] = 0;                  /* magic */
    fact->rows_ok = false;
    i = nrow + maxinv + 5;    /* DEAD */
  } else {
    /* *************************************** */
    /*       enough memory to do a bit faster */
    /*       For each row compute reciprocal of pivot element and */
    /*       take out of U */
    /*       Also use HLINK(1 to permute column numbers */
    int ninbas=0;
    int ilast; /* last available entry */
    int spareSpace;
    int * hcoli2;
    double * dluval2;
    /*int * hlink2 = ihlink+nrow;
      int * mrstrt2 = hlink2+nrow;*/
    int extraSpace=10000;
    /* mwork has order of row copy */
    EKKHlink *mwork     = (reinterpret_cast<EKKHlink*>(fact->kw1adr))-1;
    fact->rows_ok = true;
    
    if (if_sparse_update) {
      ilast=nnetas-nnentl;
    } else {
      /* missing out nnentl stuff */
      ilast=nnetas;
    }
    spareSpace=ilast-nnentu;
    need_more_space=false;
    hcoli2 = hcoli+spareSpace;
    /*     save clean row copy if enough room */
    nroom = (spareSpace) / nrow;
    if ((nnentu<<3)>150*maxinv) {
      extraSpace=150*maxinv;
    } else {
      extraSpace=nnentu<<3;
    }
    if (nrow<10000) {
      if (nroom < 10) {
      need_more_space=true;
      }
    } else {
      if (nroom < 5&&!if_sparse_update) {
      need_more_space=true;
      }
    }
    if (nroom > CoinMin(50,maxinv)) {
      need_more_space=false;
    }
    if (need_more_space) {
      if (if_sparse_update) {
      int i1=fact->eta_size+10*nrow;
      fact->eta_size=static_cast<int>(1.2*fact->eta_size);
      if (i1>fact->eta_size) {
        fact->eta_size=i1;
      }
      } else {
      fact->eta_size=static_cast<int>(1.05*fact->eta_size);
      }
    } else {
      if (nroom<11) {
      if (if_sparse_update) {
        int i1=fact->eta_size+(11-nroom)*nrow;
        fact->eta_size=static_cast<int>(1.2*fact->eta_size);
        if (i1>fact->eta_size) {
          fact->eta_size=i1;
        }
      }
      }
    }
    if (fact->maxNNetas>0&&fact->eta_size>
      fact->maxNNetas) {
      fact->eta_size=fact->maxNNetas;
    }
    {
      /* we can swap de2val and dluval to save copying */
      int * eta_last=mpermu+nrow*2+3;
      int * eta_next=eta_last+nrow+2;
      int last=0;
      eta_last[0]=-1;
      if (nnentl) {
      /* went into c_ekkcmfc - if not then in order */
      int next;
      /*next=mwork[((nrow+1)<<1)+1];*/
      next=mwork[nrow+1].pre;
#ifdef DEBUG
      j=mrstrt[next];
#endif
      for (i = 1; i <= nrow; ++i) {
        int iperm=mpermu[next];
        eta_next[last]=iperm;
        eta_last[iperm]=last;
        temp[iperm] = mrstrt[next];
        temp2[iperm] = hinrow[next];
#ifdef DEBUG
        if (mrstrt[next]!=j) abort();
        j=mrstrt[next]+hinrow[next];
#endif
        /*next= mwork[(next<<1)+2];*/
        next= mwork[next].suc;
        last=iperm;
      }
      } else {
#ifdef DEBUG
      j=0;
#endif
      for (i = 1; i <= nrow; ++i) {
        int iperm=mpermu[i];
        eta_next[last]=iperm;
        eta_last[iperm]=last;
        temp[iperm] = mrstrt[i];
        temp2[iperm] = hinrow[i];
        last=iperm;
#ifdef DEBUG
        if (mrstrt[i]<=j) abort();
        if (i>1&&mrstrt[i]!=j+hinrow[i-1]) abort();
        j=mrstrt[i];
#endif
      }
      }
      eta_next[last]=nrow+1;
      eta_last[nrow+1]=last;
      eta_next[nrow+1]=nrow+2;
      c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
      c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
      i=eta_last[nrow+1];
      ninbas=mrstrt[i]+hinrow[i]-1;
#ifdef DEBUG
      if (spareSpace<ninbas) {
        abort();
      }
#endif
      c_ekkizero( nrow, &hincol[1]);
#ifdef DEBUG
      for (i=nrow; i>0; i--) {
      int krs = mrstrt[i];
      int jpiv = hcoli[krs];
      if (ihlink[jpiv]!=i) abort();
      }
#endif
      for (i = 1; i <= ninbas; ++i) {
      k = hcoli[i];
      k = ihlink[k];
#ifdef DEBUG
        if (k<=0||k>nrow) abort();
#endif
      hcoli[i]=k;
      hincol[k]++;
      }
#ifdef DEBUG
      for (i=nrow; i>0; i--) {
      int krs = mrstrt[i];
      int jpiv = hcoli[krs];
      if (jpiv!=i) abort();
      if (krs>ninbas) abort();
      }
#endif
      /*       Sort into column order as was stored by row */
      k = 1;
      /*        Position */
      for (kk = 1; kk <= nrow; ++kk) {
      nel=hincol[kk];
        mcstrt[kk] = k;
      hrowi[k]=nel-1;
        k += hincol[kk];
      hincol[kk]=0;
      }
      if (de2val) {
      dluval2=de2val;
      } else {
      dluval2=dluval+ninbas;
      }
      nnentu = k-1;
      mcstrt[nrow + 4] = nnentu + 1;
      /* create column copy */
      for (i=nrow; i>0; i--) {
      int krs = mrstrt[i];
      int kre = krs + hinrow[i];
      hinrow[i]--;
      mrstrt[i]++;
        {
        int kx = mcstrt[i];
        /*nel = hincol[i];
          if (hrowi[kx]!=nel) abort();
          hrowi[kx] = nel-1;*/
        dluval2[kx] = 1.0 /dluval[krs];
        /*hincol[i]=0;*/
        for (kk = krs + 1; kk < kre; ++kk) {
          int j = hcoli[kk];
          int iput = hincol[j]+1;
          hincol[j]=iput;
          iput+= mcstrt[j];
          hrowi[iput] = SHIFT_INDEX(i);
          dluval2[iput] = dluval[kk];
        }
      }
      }
      if (de2val) {
      double * a=dluval;
      double * address;
      /* move first down */
      i=eta_next[0];
      {
        int krs=mrstrt[i];
        nel=hinrow[i];
        for (j=1;j<=nel;j++) {
          hcoli[j]=hcoli[j+krs-1];
          dluval[j]=dluval[j+krs-1];
        }
      }
      mrstrt[i]=1;
      /****** swap dluval and de2val !!!! ******/
      /* should work even for dspace */
      /* move L part across */
      address=fact->xeeadr+1;
      fact->xeeadr=fact->xe2adr-1;
      fact->xe2adr=address;
      if (nnentl) {
        int n=xnetal-nrow-maxinv-5;
        int j1,j2;
        int * mcstrt2=mcstrt+nrow+maxinv+4;
        j2 = mcstrt2[1];
        j1 = mcstrt2[n+1]+1;
#if 0
        memcpy(de2val+j1,dluval+j1,(j2-j1+1)*sizeof(double));
#else
        c_ekkdcpy(j2-j1+1,
              (dluval+j1),(de2val+j1));
#endif
      }
      dluval = de2val; 
      de2val = a;
      } else {
      /* copy down dluval */
#if 0
      memcpy(&dluval[1],&dluval2[1],ninbas*sizeof(double));
#else
      c_ekkdcpy(ninbas,
            (dluval2+1),(dluval+1));
#endif
      }
      /* sort dense part */
      for (i=nrow-ndenuc+1;i<=nrow;i++) {
      int kx = mcstrt[i]+1;
      int nel = hincol[i];
      c_ekk_sort2(&hrowi[kx],&dluval[kx],nel);
      }
    }
    mrstrt[nrow + 1] = ilast + 1;
  }
  /* Find first non slack */
  for (i = 1; i <= nrow; ++i) {
    int kcs = mcstrt[i];
    if (hincol[i] != 0 || dluval[kcs] != SLACK_VALUE) {
      break;
    }
  }
  numberSlacks = i - 1;
  {
    /* set slacks to 1 */
    int * array = fact->krpadr + ( fact->nrowmx+2);
    int nSet = (numberSlacks)>>5;
    int n2 = (fact->nrowmx+32)>>5;
    int i;
    memset(array,0xff,nSet*sizeof(int));
    memset(array+nSet,0,(n2-nSet)*sizeof(int));
    for (i=nSet<<5;i<=numberSlacks;i++)
      c_ekk_Set(array,i);
    c_ekk_Unset(array,fact->nrow+1); /* make sure off end not slack */
#ifndef NDEBUG
    for (i=1;i<=numberSlacks;i++)
      assert (c_ekk_IsSet(array,i));
    for (;i<=fact->nrow;i++)
      assert (!c_ekk_IsSet(array,i));
#endif
  }
  
  /* and set up backward pointers */
  /* clean up HPIVCO for fancy assembler stuff */
  /* xnetal was initialized to nrow + maxinv + 4 in c_ekktria, and grows */
  c_ekkscpy_0_1(maxinv + 1, 1, &hpivco[nrow+4]);      /* magic */
  
  hpivco[xnetal] = 1;
  /* shuffle down for gaps so can get rid of hpivco for L */
  {
    const int lstart    = nrow + maxinv + 5;
    int n=xnetal-lstart ;     /* number of L entries */
    int add,iel;
    int * hpivco_L = &hpivco[lstart];
    int * mcstrt_L = &mcstrt[lstart];
    if (nnentl) {
      /* elements of L were stored in descending order in dluval/hcoli */
      int kle = mcstrt_L[0];
      int kls = mcstrt_L[n]+1;
      
      if(if_sparse_update) {
      int i2,iel;
      int * mrstrt2 = &mrstrt[nrow];
      
      /* need row copy of L */
      /* hpivro is spare for counts; just used as a temp buffer */
      c_ekkizero( nrow, &hpivro[1]);
      
      /* permute L indices; count L row lengths */
      for (iel = kls; iel <= kle; ++iel) {
        int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
        hpivro[jrow]++;
        hrowi[iel] = SHIFT_INDEX(jrow);
      }
      {
        int ibase=nnetas-nnentl+1;
        int firstDoRow=0;
        for (i=1;i<=nrow;i++) {
          mrstrt2[i]=ibase;
          if (hpivro[i]&&!firstDoRow) {
            firstDoRow=i;
          }
          ibase+=hpivro[i];
          hpivro[i]=mrstrt2[i];
        }
        if (!firstDoRow) {
          firstDoRow=nrow+1;
        }
        mrstrt2[i]=ibase;
        fact->firstDoRow = firstDoRow;
      }
      i2=mcstrt_L[n];
      for (i = n-1; i >= 0; --i) {
        int i1 = mcstrt_L[i];
        int ipiv=hpivco_L[i];
        ipiv=mpermu[ipiv];
        hpivco_L[i]=ipiv;
        for (iel=i2 ; iel < i1; iel++) {
          int irow = UNSHIFT_INDEX(hrowi[iel+1]);
          int iput=hpivro[irow];
          hpivro[irow]=iput+1;
          hcoli[iput]=ipiv;
          de2val[iput]=dluval[iel+1];
        }
        i2=i1;
      }
      } else {
      /* just permute row numbers */
      
      for (j = 0; j < n; ++j) {
        hpivco_L[j] = mpermu[hpivco_L[j]];
      }
      for (iel = kls; iel <= kle; ++iel) {
        int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
        hrowi[iel] = SHIFT_INDEX(jrow);
      }
      }
      
      add=hpivco_L[n-1]-hpivco_L[0]-n+1;
      if (add) {
      int i;
      int last = hpivco_L[n-1];
      int laststart = mcstrt_L[n];
      int base=hpivco_L[0]-1;
      /* adjust so numbers match */
      mcstrt_L-=base;
      hpivco_L-=base;
      mcstrt_L[last]=laststart;
      for (i=n-1;i>=0;i--) {
        int ipiv=hpivco_L[i+base];
        while (ipiv<last) {
          mcstrt_L[last-1]=laststart;
          hpivco_L[last-1]=last;
          last--;
        }
        laststart=mcstrt_L[i+base];
        mcstrt_L[last-1]=laststart;
        hpivco_L[last-1]=last;
        last--;
      }
      xnetal+=add;
      }
    }
    //int lstart=fact->lstart;
    //const int * COIN_RESTRICT hpivco    = fact->kcpadr;
    fact->firstLRow = hpivco[lstart];
  }
  fact->nnentu = nnentu;
  fact->xnetal = xnetal;
  /* now we have xnetal * we can set up pointers */
  clp_setup_pointers(fact);
  
  /* this is the array used in c_ekkbtrn; it is passed to c_ekkbtju as hpivco.
   * this gets modified by F-T as we pivot columns in and out.
   */
  {
    /* do new hpivco */
    int * hpivco_new = fact->kcpadr+1;
    int * back = &fact->kcpadr[2*nrow+maxinv+4];
    /* set zeroth to stop illegal read */
    back[0]=1;
    
    hpivco_new[nrow+1]=nrow+1; /* deliberate loop for dense tests */
    hpivco_new[0]=1;
    
    for (i=1;i<=nrow;i++) {
      hpivco_new[i]=i+1;
      back[i+1]=i;
    }
    back[1]=0;
    
    fact->first_dense = CoinMax(fact->first_dense,4);
    fact->numberSlacks=numberSlacks;
    fact->lastSlack=numberSlacks;
    fact->firstNonSlack=hpivco_new[numberSlacks];
  }
  
  /* also zero out permute region and nonzero */
  c_ekkdzero( nrow, (dpermu+1));
  
  if (if_sparse_update) {
    char * nonzero = reinterpret_cast<char *> (&mpermu[nrow+1]);  /* used in c_ekkbtrn */
    /*c_ekkizero(nrow,(int *)nonzero);*/
    c_ekkczero(nrow,nonzero);
    /*memset(nonzero,0,nrow*sizeof(int));*/ /* for faster method */
  }
  for (i = 1; i <= nrow; ++i) {
    hpivro[mpermu[i]] = i;
  }
  
} /* c_ekkshfv */


static void c_ekkclcp1(const int *hcol, const int * mrstrt,
                 int *hrow, int *mcstrt,
                 int *hincol, int nnrow, int nncol,
                 int ninbas)
{
  int i, j, kc, kr, kre, krs, icol; 
  int iput;
  
  /* Create columnwise storage of row indices */
  
  kc = 1;
  for (j = 1; j <= nncol; ++j) {
    mcstrt[j] = kc;
    kc += hincol[j];
    hincol[j] = 0;
  }
  mcstrt[nncol + 1] = ninbas + 1;
  
  for (i = 1; i <= nnrow; ++i) {
    krs = mrstrt[i];
    kre = mrstrt[i + 1] - 1;
    for (kr = krs; kr <= kre; ++kr) {
      icol = hcol[kr];
      iput = hincol[icol];
      hincol[icol] = iput + 1;
      iput += mcstrt[icol];
      hrow[iput] = i;
    }
  }
} /* c_ekkclcp */
inline void c_ekkclcp2(const int *hcol, const double *dels, const int * mrstrt,
                 int *hrow, double *dels2, int *mcstrt,
                 int *hincol, int nnrow, int nncol,
                 int ninbas)
{
  int i, j, kc, kr, kre, krs, icol; 
  int iput;
  
  /* Create columnwise storage of row indices */
  
  kc = 1;
  for (j = 1; j <= nncol; ++j) {
    mcstrt[j] = kc;
    kc += hincol[j];
    hincol[j] = 0;
  }
  mcstrt[nncol + 1] = ninbas + 1;
  
  for (i = 1; i <= nnrow; ++i) {
    krs = mrstrt[i];
    kre = mrstrt[i + 1] - 1;
    for (kr = krs; kr <= kre; ++kr) {
      icol = hcol[kr];
      iput = hincol[icol];
      hincol[icol] = iput + 1;
      iput += mcstrt[icol];
      hrow[iput] = i;
      dels2[iput] = dels[kr];
    }
  }
} /* c_ekkclcp */
int c_ekkslcf( register const EKKfactinfo *fact)
{
  int * hrow = fact->xeradr;
  int * hcol = fact->xecadr;
  double * dels = fact->xeeadr;
  int * hinrow = fact->xrnadr;
  int * hincol = fact->xcnadr;
  int * mrstrt = fact->xrsadr;
  int * mcstrt = fact->xcsadr;
  const int nrow = fact->nrow;
  int ninbas;
  /* space for etas */
  const int nnetas      = fact->nnetas;
  ninbas=mcstrt[nrow+1]-1;
  
  /* Now sort */
  if (ninbas << 1 > nnetas) {
    /* Put it in row order */
    int i,k;
    c_ekkrowq(hrow, hcol, dels, mrstrt, hinrow, nrow, ninbas);
    k = 1;
    for (i = 1; i <= nrow; ++i) {
      mrstrt[i] = k;
      k += hinrow[i];
    }
    mrstrt[nrow + 1] = k;
    
    /* make a column copy without the extra values */
    c_ekkclcp1(hcol, mrstrt, hrow, mcstrt, hincol, nrow, nrow, ninbas); 
  } else {
    /* Move elements up memory */
    c_ekkdcpy(ninbas,
          (dels+1), (dels+ninbas + 1));
    
    /* make a row copy with the extra values */
    c_ekkclcp2(hrow, &dels[ninbas], mcstrt, hcol, dels, mrstrt, hinrow, nrow, nrow, ninbas); 
  }
  return (ninbas);
} /* c_ekkslcf */
/*     Uwe H. Suhl, September 1986 */
/*     Removes lower and upper triangular factors from the matrix. */
/*     Code for routine: 102 */
/*     Return codes: */
/*    0: ok */
/*      -5: not enough space in row file */
/*      7: pivot too small - col sing */
/*
 * This selects singleton columns and rows for the LU factorization.
 * Singleton columns require no 
 *
 * (1) Note that columns are processed using a queue, not a stack;
 * this produces better pivots.
 *
 * (2) At most nrows elements are ever entered into the queue.
 *
 * (3) When pivoting singleton columns, every column that is part of
 * the pivot row is shortened by one, including the singleton column
 * itself; the hincol entries are updated appropriately.
 * Thus, pivoting on a singleton column may create other singleton columns
 * (but not singleton rows).
 * The dual property is true for rows.
 *
 * (4) Row entries (hrowi) are not changed when pivoting singleton columns.
 * Singleton columns that are created as a result of pivoting the
 * rows of other singleton columns will therefore have row entries
 * corresponding to those pivoted rows.  Since we need to find the
 * row entry for the row being pivoted, we have to
 * search its row entries for the one whose hlink entry indicates
 * that it has not yet been pivoted.
 *
 * (9) As a result of pivoting columns, sections in hrowi corresponding to
 * pivoted columns are no longer needed, and entries in sections
 * for non-pivoted columns may have entries corresponding to pivoted rows.
 * This is why hrowi needs to be compacted.
 *
 * (5) When the row_pre and col_pre fields of the hlink struct contain
 * negative values, they indicate that the row has been pivoted, and
 * the negative of that value is the pivot order.
 * That is the only use for these fields in this routine.
 *
 * (6) This routine assumes that hlink is initialized to zeroes.
 * Under this assumption, the following is an invariant in this routine:
 *
 *    (clink[i].pre < 0) ==> (hincol[i]==0)
 *
 * The converse is not true; see (15).
 *
 * The dual is also true, but only while pivoting singletong rows,
 * since we don't update hinrow while pivoting columns;
 * THESE VALUES ARE USED LATER, BUT I DON'T UNDERSTAND HOW YET.
 *
 * (7) hpivco is used for two purposes.  The low end is used to implement the
 * queue when pivoting columns; the high end is used to hold eta-matrix
 * entries. 
 *
 * (8) As a result of pivoting columns, for all i:1<=i<=nrow, either
 *    hinrow[i] has not changed
 * or
 *    hinrow[i] = 0
 * This is another way of saying that pivoting singleton columns cannot
 * create singleton rows.
 * The dual holds for hincol after pivoting rows.
 *
 * (10) In constrast to (4), while pivoting rows we
 * do not let the hcoli get out-of-date.  That is because as part of
 * the process of numerical pivoting we have to find the row entries
 * for all the rows in the pivot column, so we may as well keep the
 * entries up to date.  This is done by moving the last column entry
 * for each row into the entry that was used for the pivot column.
 *
 * (11) When pivoting a column, we must find the pivot row entry in
 * its row table.  Sometimes we search for other things at the same time.
 * The same is true for pivoting columns.  This search should never
 * fail.
 *
 * (12) Information concerning the eta matrices is stored in the high
 * ends of arrays that are also used to store information concerning
 * the basis; these arrays are: hpivco, mcstrt, dluval and hcoli.
 * Information is only stored in these arrays as a part of pivoting
 * singleton rows, since the only thing that needs to be saved as
 * a part of pivoting singleton columns is which rows and columns were chosen,
 * and this is stored in hlink.
 * Since they have to share the same array, the eta information grows
 * downward instead of upward. Eventually, eta information may grow
 * down to the top of the basis information.  As pivoting proceeds,
 * more and more of this information is no longer needed, so when this
 * happens we can try compacting the arrays to see if we can recover
 * enough space.  lstart points at the bottom entry in the arrays,
 * xnewro/xnewco at the top of the basis information, and each time we
 * pivot a singleton row we know that we will need exactly as many new
 * entries as there are rows in the pivot column, so we can easily
 * determine if we need more room.  The variable maxinv may be used
 * to reserve extra room when inversion starts.
 *
 * (13) Eta information is stored in a fashion that is similar to how
 * matrices are stored.  There is one entry in hpivco and mcstrt for
 * each eta (other than the initial ones for singleton columns and
 * for singleton rows that turn out to be singleton columns),
 * in the order they were chosen.  hpivco records the pivot row,
 * and mcstrt points at the first entry in the other two arrays
 * for this row.  dluval contains the actual eta values for the column,
 * and hcoli the rows these values were in.
 * These entries in mcstrt and hpivco grow upward; they start above
 * the entries used to store basis information.
 * (Actually, I don't see why they need to start maxinv+4 entries past the top).
 *
 * (14) c_ekkrwco assumes that invalidated hrowi/hcoli entries contain 0.
 *
 * (15) When pivoting singleton columns, it may possibly happen
 * that a row with all singleton column entries is created.
 * In this case, all of the columns will be enqueued, and pivoting
 * on any of them eliminates the rest, without their being chosen
 * as pivots.  The dual holds for singleton rows.
 * DOES THIS INDICATE A SINGULARITY?
 *
 * (15) There are some aspects of the implementation that I find odd.
 * hrowi is not set to 0 for pivot rows while pivoting singleton columns,
 * which would make sense to me.  Things don't work if this isn't done,
 * so the information is used somehow later on.  Also, the information
 * for the pivot column is shifted to the front of the pivot row
 * when pivoting singleton columns; this is also necessary for reasons
 * I don't understand.
 */
int c_ekktria(EKKfactinfo *fact,
            EKKHlink * rlink,
            EKKHlink * clink,
          int *nsingp, 
          int *xnewcop, int *xnewrop,
          int *ncompactionsp,
          const int ninbas)
{
  const int nrow  = fact->nrow;
  const int maxinv      = fact->maxinv;
  int *hcoli      = fact->xecadr;
  double *dluval  = fact->xeeadr;
  int *mrstrt     = fact->xrsadr;
  int *hrowi      = fact->xeradr;
  int *mcstrt     = fact->xcsadr;
  int *hinrow     = fact->xrnadr;
  int *hincol     = fact->xcnadr;
  int *stack      = fact->krpadr; /* normally hpivro */
  int *hpivco     = fact->kcpadr;
  const double drtpiv   = fact->drtpiv;
  CoinZeroN(reinterpret_cast<int *>(rlink+1),static_cast<int>(nrow*(sizeof(EKKHlink)/sizeof(int))));
  CoinZeroN(reinterpret_cast<int *>(clink+1),static_cast<int>(nrow*(sizeof(EKKHlink)/sizeof(int))));
  
  fact->npivots   = 0;
  /*      Use NUSPIK to keep sum of deactivated row counts */
  fact->nuspike   = 0;
  int xnetal      = nrow + maxinv + 4;
  int xnewro      = mrstrt[nrow] + hinrow[nrow] - 1;
  int xnewco      = xnewro;
  int kmxeta      = ninbas;
  int ncompactions      = 0;
  
  int i, j, k, kc, kce, kcs, npr;
  double pivot;
  int kipis, kipie, kjpis, kjpie, knprs, knpre;
  int ipivot, jpivot, stackc, stackr;
#ifndef NDEBUG
  int kpivot=-1;
#else
  int kpivot=-1;
#endif
  int epivco, kstart, maxstk;
  int irtcod = 0;
  int lastSlack=0;
  
  int lstart = fact->nnetas + 1;
  /*int nnentu    = ninbas; */
  int lstart_minus_nnentu=lstart-ninbas;
  /* do initial column singletons - as can do faster */
  for (jpivot = 1; jpivot <= nrow; ++jpivot) {
    if (hincol[jpivot] == 1) {
      ipivot = hrowi[mcstrt[jpivot]];
      if (ipivot>lastSlack) {
      lastSlack=ipivot;
      } else {
        /* so we can't put a structural over a slack */
      break;
      }
      kipis = mrstrt[ipivot];
#if 1
      assert (hcoli[kipis]==jpivot);
#else
      if (hcoli[kipis]!=jpivot) {
      kpivot=kipis+1;
      while(hcoli[kpivot]!=jpivot) kpivot++;
#ifdef DEBUG
      kipie = kipis + hinrow[ipivot] ;
      if (kpivot>=kipie) {
        abort();
      }
#endif
        pivot=dluval[kpivot];
      dluval[kpivot] = dluval[kipis];
      dluval[kipis] = pivot;
      hcoli[kpivot] = hcoli[kipis];
      hcoli[kipis] = jpivot;
      }
#endif
      if (dluval[kipis]==SLACK_VALUE) {
      /* record the new pivot row and column */
      ++fact->npivots;
      rlink[ipivot].pre = -fact->npivots;
      clink[jpivot].pre = -fact->npivots;
      hincol[jpivot]=0;
      fact->nuspike += hinrow[ipivot];
      } else {
      break;
      }
    } else {
      break;
    }
  }
  /* Fill queue with other column singletons and clean up */
  maxstk = 0;
  for (j = 1; j <= nrow; ++j) {
    if (hincol[j]) {
      int n=0;
      kcs = mcstrt[j];
      kce = mcstrt[j + 1];
      for (k = kcs; k < kce; ++k) {
      if (! (rlink[hrowi[k]].pre < 0)) {
        n++;
      }
      }
      hincol[j] = n;
      if (n == 1) {
      /* we just created a new singleton column - enqueue it */
      ++maxstk;
      stack[maxstk] = j;
      }
    }
  }
  stackc = 0; /* (1) */
  
  while (! (stackc >= maxstk)) {    /* (1) */
    /* dequeue the next entry */
    ++stackc;
    jpivot = stack[stackc];
    
    /* (15) */
    if (hincol[jpivot] != 0) {
      
      for (k = mcstrt[jpivot]; rlink[hrowi[k]].pre < 0; k++) {
      /* (4) */
      }
      ipivot = hrowi[k];
      
      /* All the columns in this row are being shortened. */
      kipis = mrstrt[ipivot];
      kipie = kipis + hinrow[ipivot] ;
      for (k = kipis; k < kipie; ++k) {
      j = hcoli[k];
      --hincol[j];      /* (3) (6) */
      
      if (j == jpivot) {
        kpivot = k;           /* (11) */
      } else if (hincol[j] == 1) {
        /* we just created a new singleton column - enqueue it */
        ++maxstk;
        stack[maxstk] = j;
      }
      }
      
      /* record the new pivot row and column */
      ++fact->npivots;
      rlink[ipivot].pre = -fact->npivots;
      clink[jpivot].pre = -fact->npivots;
      
      fact->nuspike += hinrow[ipivot];
      
      /* check the pivot */
      assert (kpivot>0);
      pivot = dluval[kpivot];
      if (fabs(pivot) < drtpiv) {
      irtcod = 7;
      ++(*nsingp);
      rlink[ipivot].pre = -nrow - 1;
      clink[jpivot].pre = -nrow - 1;
      }
      
      /* swap the pivot column entry with the first one. */
      /* I don't know why. */
      dluval[kpivot] = dluval[kipis];
      dluval[kipis] = pivot;
      hcoli[kpivot] = hcoli[kipis];
      hcoli[kipis] = jpivot;
    }
  }
  /* (8) */
  
  /* The entire basis may already be triangular */
  if (fact->npivots < nrow) {
    
    /* (9) */
    kstart = 0;
    for (j = 1; j <= nrow; ++j) {
      if (! (clink[j].pre < 0)) {
      kcs = mcstrt[j];
      kce = mcstrt[j + 1];
      
      mcstrt[j] = kstart + 1;
      
      for (k = kcs; k < kce; ++k) {
        if (! (rlink[hrowi[k]].pre < 0)) {
          ++kstart;
          hrowi[kstart] = hrowi[k];
        }
      }
      hincol[j] = kstart - mcstrt[j] + 1;
      }
    }
    xnewco = kstart;
    
    
    /* Fill stack with initial row singletons that haven't been pivoted away */
    stackr = 0;
    for (i = 1; i <= nrow; ++i) {
      if (! (rlink[i].pre < 0) &&
        (hinrow[i] == 1)) {
      ++stackr;
      stack[stackr] = i;
      }
    }
    
    while (! (stackr <= 0)) {
      ipivot = stack[stackr];
      assert (ipivot);
      --stackr;
      
#if 1
      assert (rlink[ipivot].pre>=0);
#else
      /* This test is probably unnecessary:  rlink[i].pre < 0 ==> hinrow[i]==0 */
      if (rlink[ipivot].pre < 0) {
      continue;
      }
#endif
      
      /* (15) */
      if (hinrow[ipivot] != 0) {
      
      /* This is a singleton row, which means it has exactly one column */
      jpivot = hcoli[mrstrt[ipivot]];
      
      kjpis = mcstrt[jpivot];
      epivco = hincol[jpivot] - 1;
      hincol[jpivot] = 0;     /* this column is being pivoted away */
      
      /* (11) */
      kjpie = kjpis + epivco;
      for (k = kjpis; k <= kjpie; ++k) {
        if (ipivot == hrowi[k])
          break;
      }
      /* ASSERT (k <= kjpie) */
      
      /* move the last row entry for the pivot column into the pivot row's entry */
      /* I don't know why */
      hrowi[k] = hrowi[kjpie];
      
      /* invalidate the (old) last row entry of the pivot column */
      /* I don't know why */
      hrowi[kjpie] = 0;
      
      /* (12) */
      if (! (xnewro + epivco < lstart)) {
        int kstart;
        
        if (! (epivco < lstart_minus_nnentu)) {
          irtcod = -5;
          break;
        }
        kstart = c_ekkrwco(fact,dluval, hcoli, mrstrt, hinrow, xnewro);
        ++ncompactions;
        kmxeta += (xnewro - kstart) << 1;
        xnewro = kstart;
      }
      if (! (xnewco + epivco < lstart)) {
        if (! (epivco < lstart_minus_nnentu)) {
          irtcod = -5;
          break;
        }
        xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
        ++ncompactions;
        
        /*     HINCOL MAY HAVE CHANGED ??? (JJHF) */
        epivco = hincol[jpivot];
      }
      
      /* record the new pivot row and column */
      ++fact->npivots;
      rlink[ipivot].pre = -fact->npivots;
      clink[jpivot].pre = -fact->npivots;
      
      /* no update for nuspik */
      
      /* check the pivot */
      pivot = dluval[mrstrt[ipivot]];
      if (fabs(pivot) < drtpiv) {
        /* If the pivot is too small, reject it, but keep going */
        irtcod = 7;
        rlink[ipivot].pre = -nrow - 1;
        clink[jpivot].pre = -nrow - 1;
      }
      
      /* Perform numerical part of elimination. */
      if (! (epivco <= 0)) {
        ++xnetal;
        mcstrt[xnetal] = lstart - 1;
        hpivco[xnetal] = ipivot;
        pivot = -1.f / pivot;
        
        kcs = mcstrt[jpivot];
        kce = kcs + epivco - 1;
        hincol[jpivot] = 0;
        
        for (kc = kcs; kc <= kce; ++kc) {
          npr = hrowi[kc];
          
          /* why bother? */
          hrowi[kc] = 0;
          
          --hinrow[npr];      /* (3) */
          if (hinrow[npr] == 1) {
            /* this may create new singleton rows */
            ++stackr;
            stack[stackr] = npr;
          }
          
          /* (11) */
          knprs = mrstrt[npr];
          knpre = knprs + hinrow[npr];
          for (k = knprs; k <= knpre; ++k) {
            if (jpivot == hcoli[k]) {
            kpivot = k;
            break;
            }
          }
          /* ASSERT (kpivot <= knpre) */
          
          {
            /* (10) */
            double elemnt = dluval[kpivot];
            dluval[kpivot] = dluval[knpre];
            hcoli[kpivot] = hcoli[knpre];
            
            hcoli[knpre] = 0; /* (14) */
            
            /* store elementary row transformation */
            --lstart;
            dluval[lstart] = elemnt * pivot;
            hcoli[lstart] = npr;
          }
        }
      }
      }
    }
  }
  /* (8) */
  
  *xnewcop = xnewco;
  *xnewrop = xnewro;
  fact->xnetal = xnetal;
  fact->nnentu = lstart - lstart_minus_nnentu;
  fact->kmxeta = kmxeta;
  *ncompactionsp = ncompactions;
  
  return (irtcod);
} /* c_ekktria */

Generated by  Doxygen 1.6.0   Back to index